From 53eaaa735f95382c460141593de799c323900874 Mon Sep 17 00:00:00 2001 From: Marvin Hemmer Date: Tue, 3 Feb 2026 15:04:13 +0100 Subject: [PATCH 1/2] [PWGEM] PhotonMeson: Add basic EMCal QA to EMCPhotonCut - Add basic cluster QA to EMCPhotonCut class when using `AreSelectedRunning` function. QA is not enabled for `IsSelected` since IsSelected is considered legecy code that should slowly be removed. --- PWGEM/PhotonMeson/Core/EMCPhotonCut.cxx | 6 + PWGEM/PhotonMeson/Core/EMCPhotonCut.h | 151 ++++++++++++++++++++- PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx | 5 +- 3 files changed, 154 insertions(+), 8 deletions(-) diff --git a/PWGEM/PhotonMeson/Core/EMCPhotonCut.cxx b/PWGEM/PhotonMeson/Core/EMCPhotonCut.cxx index 1ccca2a2be7..3d911df19e3 100644 --- a/PWGEM/PhotonMeson/Core/EMCPhotonCut.cxx +++ b/PWGEM/PhotonMeson/Core/EMCPhotonCut.cxx @@ -82,6 +82,12 @@ void EMCPhotonCut::SetUseSecondaryTM(bool flag) LOG(info) << "EM Photon Cluster Cut, using secondary TM cut is set to : " << mUseTM; } +void EMCPhotonCut::SetDoQA(bool flag) +{ + mDoQA = flag; + LOG(info) << "EM Photon Cluster Cut, QA is set to: " << mUseTM; +} + void EMCPhotonCut::print() const { LOG(info) << "EMCal Photon Cut:"; diff --git a/PWGEM/PhotonMeson/Core/EMCPhotonCut.h b/PWGEM/PhotonMeson/Core/EMCPhotonCut.h index 64913f86d3e..d34c1003315 100644 --- a/PWGEM/PhotonMeson/Core/EMCPhotonCut.h +++ b/PWGEM/PhotonMeson/Core/EMCPhotonCut.h @@ -20,6 +20,7 @@ #include "PWGEM/PhotonMeson/DataModel/gammaTables.h" #include +#include #include @@ -92,7 +93,7 @@ class EMCPhotonCut : public TNamed EMCPhotonCut() = default; EMCPhotonCut(const char* name, const char* title) : TNamed(name, title) {} - enum class EMCPhotonCuts : int { + enum class EMCPhotonCuts : std::uint8_t { // cluster cut kDefinition = 0, kEnergy, @@ -105,6 +106,12 @@ class EMCPhotonCut : public TNamed kNCuts }; + enum class TrackType : std::uint8_t { + // cluster cut + kPrimary = 0, + kSecondary, + }; + static const char* mCutNames[static_cast(EMCPhotonCuts::kNCuts)]; constexpr auto getClusterId(o2::soa::is_iterator auto const& t) const @@ -126,7 +133,7 @@ class EMCPhotonCut : public TNamed /// \param GetPhiCut lambda to get the phi cut value /// \param applyEoverP bool to check if E/p should be checked (for secondaries we do not check this!) bool checkTrackMatching(o2::soa::is_iterator auto const& cluster, IsTrackIterator auto& emcmatchedtrack, o2::soa::RowViewSentinel const emcmatchedtrackEnd, - bool applyEoverP, auto GetEtaCut, auto GetPhiCut) const + bool applyEoverP, auto GetEtaCut, auto GetPhiCut, o2::framework::HistogramRegistry* fRegistry = nullptr, TrackType trackType = TrackType::kPrimary) const { // advance to cluster while (emcmatchedtrack != emcmatchedtrackEnd && getClusterId(emcmatchedtrack) < cluster.globalIndex()) { @@ -150,6 +157,22 @@ class EMCPhotonCut : public TNamed (dPhi > GetPhiCut(trackpt)) || (applyEoverP && cluster.e() / trackp >= mMinEoverP); if (!fail) { + if (mDoQA && fRegistry != nullptr) { + switch (trackType) { + case TrackType::kPrimary: + fRegistry->fill(HIST("Cluster/hTrackdEtadPhi"), dEta, dPhi); + fRegistry->fill(HIST("Cluster/hTrackdEtaPt"), dEta, trackpt); + fRegistry->fill(HIST("Cluster/hTrackdPhiPt"), dPhi, trackpt); + break; + case TrackType::kSecondary: + fRegistry->fill(HIST("Cluster/hSecTrackdEtadPhi"), dEta, dPhi); + fRegistry->fill(HIST("Cluster/hSecTrackdEtaPt"), dEta, trackpt); + fRegistry->fill(HIST("Cluster/hSecTrackdPhiPt"), dPhi, trackpt); + break; + default: + break; + } + } return false; // cluster got a track matche to it } ++emcmatchedtrack; @@ -157,22 +180,114 @@ class EMCPhotonCut : public TNamed return true; // all tracks checked, cluster survives } + void fillBeforeClusterHistogram(o2::soa::is_iterator auto const& cluster, o2::framework::HistogramRegistry* fRegistry = nullptr) + { + + if (mDoQA == false || fRegistry == nullptr) { + return; + } + + fRegistry->fill(HIST("Cluster/hClusterQualityCuts"), 0, cluster.e()); + fRegistry->fill(HIST("Cluster/before/hE"), 0, cluster.e()); + fRegistry->fill(HIST("Cluster/before/hPt"), 0, cluster.pt()); + fRegistry->fill(HIST("Cluster/before/hEtaPhi"), 0, cluster.eta(), cluster.phi()); + fRegistry->fill(HIST("Cluster/before/hNCell"), 0, cluster.nCells(), cluster.e()); + fRegistry->fill(HIST("Cluster/before/hM02"), 0, cluster.m02(), cluster.e()); + fRegistry->fill(HIST("Cluster/before/hTime"), 0, cluster.time(), cluster.e()); + } + + void fillAfterClusterHistogram(o2::soa::is_iterator auto const& cluster, o2::framework::HistogramRegistry* fRegistry = nullptr) + { + + if (mDoQA == false || fRegistry == nullptr) { + return; + } + + fRegistry->fill(HIST("Cluster/hClusterQualityCuts"), 0, cluster.e()); + fRegistry->fill(HIST("Cluster/before/hE"), 0, cluster.e()); + fRegistry->fill(HIST("Cluster/before/hPt"), 0, cluster.pt()); + fRegistry->fill(HIST("Cluster/before/hEtaPhi"), 0, cluster.eta(), cluster.phi()); + fRegistry->fill(HIST("Cluster/before/hNCell"), 0, cluster.nCells(), cluster.e()); + fRegistry->fill(HIST("Cluster/before/hM02"), 0, cluster.m02(), cluster.e()); + fRegistry->fill(HIST("Cluster/before/hTime"), 0, cluster.time(), cluster.e()); + } + /// \brief check if given clusters survives all cuts /// \param flags EMBitFlags where results will be stored /// \param cluster cluster table to check /// \param matchedTracks matched primary tracks table /// \param matchedSecondaries matched secondary tracks table - void AreSelectedRunning(EMBitFlags& flags, o2::soa::is_table auto const& clusters, IsTrackContainer auto const& emcmatchedtracks, IsTrackContainer auto const& secondaries) const + /// \param fRegistry HistogramRegistry pointer of the main task + void AreSelectedRunning(EMBitFlags& flags, o2::soa::is_table auto const& clusters, IsTrackContainer auto const& emcmatchedtracks, IsTrackContainer auto const& secondaries, o2::framework::HistogramRegistry* fRegistry = nullptr) { auto emcmatchedtrackIter = emcmatchedtracks.begin(); auto emcmatchedtrackEnd = emcmatchedtracks.end(); auto secondaryIter = secondaries.begin(); auto secondaryEnd = secondaries.end(); size_t iCluster = 0; + + // Check if we want to do QA and provided proper histogram registry + if (mDoQA && fRegistry != nullptr) { + const o2::framework::AxisSpec thAxisClusterEnergy{500, 0, 50, "#it{E}_{cls} (GeV)"}; + const o2::framework::AxisSpec thAxisMomentum{250, 0., 25., "#it{p}_{T} (GeV/#it{c})"}; + const o2::framework::AxisSpec thAxisDEta{200, -0.1, 0.1, "#Delta#eta"}; + const o2::framework::AxisSpec thAxisDPhi{200, -0.1, 0.1, "#Delta#varphi (rad)"}; + const o2::framework::AxisSpec thAxisEnergy{500, 0., 50., "#it{E} (GeV)"}; + const o2::framework::AxisSpec thAxisEta{320, -0.8, 0.8, "#eta"}; + const o2::framework::AxisSpec thAxisPhi{500, 0, o2::constants::math::TwoPI, "#varphi (rad)"}; + const o2::framework::AxisSpec thAxisNCell{51, -0.5, 50.5, "#it{N}_{cell}"}; + const o2::framework::AxisSpec thAxisM02{200, 0, 2.0, "#it{M}_{02}"}; + const o2::framework::AxisSpec thAxisTime{300, -150, +150, "#it{t}_{cls} (ns)"}; + const o2::framework::AxisSpec thAxisEoverP{400, 0, 10., "#it{E}_{cls}/#it{p}_{track} (#it{c})"}; + + fRegistry->add("Cluster/before/hE", "E_{cluster};#it{E}_{cluster} (GeV);#it{N}_{cluster}", o2::framework::kTH1F, {thAxisClusterEnergy}, true); + fRegistry->add("Cluster/before/hPt", "Transverse momenta of clusters;#it{p}_{T} (GeV/c);#it{N}_{cluster}", o2::framework::kTH1F, {thAxisClusterEnergy}, true); + fRegistry->add("Cluster/before/hNgamma", "Number of #gamma candidates per collision;#it{N}_{#gamma} per collision;#it{N}_{collisions}", o2::framework::kTH1F, {{1001, -0.5f, 1000.5f}}, true); + fRegistry->add("Cluster/before/hEtaPhi", "#eta vs #varphi;#eta;#varphi (rad.)", o2::framework::kTH2F, {thAxisEta, thAxisPhi}, true); + fRegistry->add("Cluster/before/hNCell", "#it{N}_{cells};N_{cells} (GeV);#it{E}_{cluster} (GeV)", o2::framework::kTH2F, {thAxisNCell, thAxisClusterEnergy}, true); + fRegistry->add("Cluster/before/hM02", "Long ellipse axis;#it{M}_{02} (cm);#it{E}_{cluster} (GeV)", o2::framework::kTH2F, {thAxisM02, thAxisClusterEnergy}, true); + fRegistry->add("Cluster/before/hTime", "Cluster time;#it{t}_{cls} (ns);#it{E}_{cluster} (GeV)", o2::framework::kTH2F, {thAxisTime, thAxisClusterEnergy}, true); + + fRegistry->addClone("Cluster/before/", "Cluster/after/"); + + auto hClusterQualityCuts = fRegistry->add("Cluster/hClusterQualityCuts", "Energy at which clusters are removed by a given cut;;#it{E} (GeV)", o2::framework::kTH2F, {{static_cast(EMCPhotonCut::EMCPhotonCuts::kNCuts) + 2, -0.5, static_cast(EMCPhotonCut::EMCPhotonCuts::kNCuts) + 1.5}, thAxisClusterEnergy}, true); + hClusterQualityCuts->GetXaxis()->SetBinLabel(1, "In"); + hClusterQualityCuts->GetXaxis()->SetBinLabel(2, "Definition"); + hClusterQualityCuts->GetXaxis()->SetBinLabel(3, "Energy"); + hClusterQualityCuts->GetXaxis()->SetBinLabel(4, "NCell"); + hClusterQualityCuts->GetXaxis()->SetBinLabel(5, "M02"); + hClusterQualityCuts->GetXaxis()->SetBinLabel(6, "Timing"); + hClusterQualityCuts->GetXaxis()->SetBinLabel(7, "TM"); + hClusterQualityCuts->GetXaxis()->SetBinLabel(8, "Sec. TM"); + hClusterQualityCuts->GetXaxis()->SetBinLabel(9, "Exotic"); + hClusterQualityCuts->GetXaxis()->SetBinLabel(10, "Out"); + + fRegistry->add("Cluster/hTrackdEtadPhi", "d#eta vs. d#varphi of matched tracks;d#eta;d#varphi (rad.)", o2::framework::kTH2F, {thAxisDEta, thAxisDPhi}, true); + fRegistry->add("Cluster/hTrackdEtaPt", "d#eta vs. track pT of matched tracks;d#eta;d#varphi (rad.)", o2::framework::kTH2F, {thAxisDEta, thAxisMomentum}, true); + fRegistry->add("Cluster/hTrackdPhiPt", "d#varphi vs. track pT of matched tracks;d#eta;d#varphi (rad.)", o2::framework::kTH2F, {thAxisDPhi, thAxisMomentum}, true); + fRegistry->add("Cluster/hSecTrackdEtadPhi", "d#eta vs. d#varphi of matched secondary tracks;d#eta;d#varphi (rad.)", o2::framework::kTH2F, {thAxisDEta, thAxisDPhi}, true); + fRegistry->add("Cluster/hSecTrackdEtaPt", "d#eta vs. track pT of matched secondary tracks;d#eta;d#varphi (rad.)", o2::framework::kTH2F, {thAxisDEta, thAxisMomentum}, true); + fRegistry->add("Cluster/hSecTrackdPhiPt", "d#varphi vs. track pT of matched secondary tracks;d#eta;d#varphi (rad.)", o2::framework::kTH2F, {thAxisDPhi, thAxisMomentum}, true); + } + + nTotClusterPerColl = 0; + currentCollID = clusters.iteratorAt(0).emeventId(); for (const auto& cluster : clusters) { + if (currentCollID == cluster.emeventId()) { + ++nTotClusterPerColl; + } else { + fRegistry->fill(HIST("Cluster/before/hNgamma"), nTotClusterPerColl); + nTotClusterPerColl = 0; + } + if (mDoQA == true || fRegistry != nullptr) { + fillBeforeClusterHistogram(cluster, fRegistry); + } if (!IsSelectedRunning(cluster, emcmatchedtrackIter, emcmatchedtrackEnd, secondaryIter, secondaryEnd)) { flags.set(iCluster); + } else if (mDoQA == true || fRegistry != nullptr) { + fillAfterClusterHistogram(cluster, fRegistry); } + currentCollID = cluster.emeventId(); ++iCluster; } } @@ -184,32 +299,47 @@ class EMCPhotonCut : public TNamed /// \param secondaryIter current iterator of matched secondary tracks /// \param secondaryEnd end iterator of matched secondary tracks /// \return true if cluster survives all cuts else false - bool IsSelectedRunning(o2::soa::is_iterator auto const& cluster, IsTrackIterator auto& emcmatchedtrackIter, o2::soa::RowViewSentinel const emcmatchedtrackEnd, IsTrackIterator auto& secondaryIter, o2::soa::RowViewSentinel const secondaryEnd) const + bool IsSelectedRunning(o2::soa::is_iterator auto const& cluster, IsTrackIterator auto& emcmatchedtrackIter, o2::soa::RowViewSentinel const emcmatchedtrackEnd, IsTrackIterator auto& secondaryIter, o2::soa::RowViewSentinel const secondaryEnd, o2::framework::HistogramRegistry* fRegistry = nullptr) { if (!IsSelectedEMCalRunning(EMCPhotonCuts::kDefinition, cluster)) { + fRegistry->fill(HIST("Cluster/hClusterQualityCuts"), static_cast(EMCPhotonCuts::kDefinition) + 1, cluster.e()); return false; } if (!IsSelectedEMCalRunning(EMCPhotonCuts::kEnergy, cluster)) { + fRegistry->fill(HIST("Cluster/hClusterQualityCuts"), static_cast(EMCPhotonCuts::kEnergy) + 1, cluster.e()); return false; } if (!IsSelectedEMCalRunning(EMCPhotonCuts::kNCell, cluster)) { + fRegistry->fill(HIST("Cluster/hClusterQualityCuts"), static_cast(EMCPhotonCuts::kNCell) + 1, cluster.e()); return false; } if (!IsSelectedEMCalRunning(EMCPhotonCuts::kM02, cluster)) { + fRegistry->fill(HIST("Cluster/hClusterQualityCuts"), static_cast(EMCPhotonCuts::kM02) + 1, cluster.e()); return false; } if (!IsSelectedEMCalRunning(EMCPhotonCuts::kTiming, cluster)) { + fRegistry->fill(HIST("Cluster/hClusterQualityCuts"), static_cast(EMCPhotonCuts::kTiming) + 1, cluster.e()); return false; } if (mUseTM && (!IsSelectedEMCalRunning(EMCPhotonCuts::kTM, cluster, emcmatchedtrackIter, emcmatchedtrackEnd))) { + fRegistry->fill(HIST("Cluster/hClusterQualityCuts"), static_cast(EMCPhotonCuts::kTM) + 1, cluster.e()); return false; } if (mUseSecondaryTM && (!IsSelectedEMCalRunning(EMCPhotonCuts::kSecondaryTM, cluster, secondaryIter, secondaryEnd))) { + fRegistry->fill(HIST("Cluster/hClusterQualityCuts"), static_cast(EMCPhotonCuts::kSecondaryTM) + 1, cluster.e()); return false; } if (!IsSelectedEMCalRunning(EMCPhotonCuts::kExotic, cluster)) { + fRegistry->fill(HIST("Cluster/hClusterQualityCuts"), static_cast(EMCPhotonCuts::kExotic) + 1, cluster.e()); return false; } + fRegistry->fill(HIST("Cluster/hClusterQualityCuts"), static_cast(EMCPhotonCuts::kNCuts) + 1, cluster.e()); + if (currentCollID == cluster.emeventId()) { + ++nAccClusterPerColl; + } else { + fRegistry->fill(HIST("Cluster/before/hNgamma"), nAccClusterPerColl); + nAccClusterPerColl = 0; + } return true; } @@ -439,6 +569,10 @@ class EMCPhotonCut : public TNamed /// \param flag flag to use secondary track matching void SetUseSecondaryTM(bool flag = false); + /// \brief Set flag to do QA + /// \param flag flag to do QA + void SetDoQA(bool flag = false); + /// \brief Set parameters for track matching delta eta = a + (pT + b)^c /// \param a a in a + (pT + b)^c /// \param b b in a + (pT + b)^c @@ -517,8 +651,13 @@ class EMCPhotonCut : public TNamed float mMaxTime{25.f}; ///< maximum cluster timing float mMinEoverP{1.75f}; ///< minimum cluster energy over track momentum ratio needed for the pair to be considered matched bool mUseExoticCut{true}; ///< flag to decide if the exotic cluster cut is to be checked or not - bool mUseTM{true}; ///< flag to decide if track matching cut is to be checek or not - bool mUseSecondaryTM{false}; ///< flag to decide if seconary track matching cut is to be checek or not + bool mUseTM{true}; ///< flag to decide if track matching cut is to be check or not + bool mUseSecondaryTM{false}; ///< flag to decide if seconary track matching cut is to be check or not + bool mDoQA{false}; ///< flag to decide if QA should be done or not + + uint nTotClusterPerColl{0}; ///< running number of all cluster per collision used for QA + uint nAccClusterPerColl{0}; ///< running number of accepted cluster per collision used for QA + int currentCollID{-1}; ///< running collision ID of clusters used for QA TrackMatchingParams mTrackMatchingEtaParams = {-1.f, 0.f, 0.f}; TrackMatchingParams mTrackMatchingPhiParams = {-1.f, 0.f, 0.f}; diff --git a/PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx b/PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx index 0b4f156f27b..9d2a3b7295e 100644 --- a/PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx +++ b/PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx @@ -348,6 +348,7 @@ struct TaskPi0FlowEMC { fEMCCut.SetClusterizer(emccuts.clusterDefinition); fEMCCut.SetUseTM(emccuts.cfgEMCUseTM.value); // disables or enables TM fEMCCut.SetUseSecondaryTM(emccuts.emcUseSecondaryTM.value); // disables or enables secondary TM + fEMCCut.SetDoQA(emccuts.cfgEnableQA.value); } void definePCMCut() @@ -1095,7 +1096,7 @@ struct TaskPi0FlowEMC { void processEMCal(CollsWithQvecs const& collisions, EMCalPhotons const& clusters, MinMTracks const& matchedPrims, MinMSTracks const& matchedSeconds) { EMBitFlags flags(clusters.size()); - fEMCCut.AreSelectedRunning(flags, clusters, matchedPrims, matchedSeconds); + fEMCCut.AreSelectedRunning(flags, clusters, matchedPrims, matchedSeconds, ®istry); if (cfgDoReverseScaling.value) { energyCorrectionFactor = 1.0505f; @@ -1238,7 +1239,7 @@ struct TaskPi0FlowEMC { void processEMCalPCMC(CollsWithQvecs const& collisions, EMCalPhotons const& clusters, PCMPhotons const& photons, aod::V0Legs const&, MinMTracks const& matchedPrims, MinMSTracks const& matchedSeconds) { EMBitFlags emcFlags(clusters.size()); - fEMCCut.AreSelectedRunning(emcFlags, clusters, matchedPrims, matchedSeconds); + fEMCCut.AreSelectedRunning(emcFlags, clusters, matchedPrims, matchedSeconds, ®istry); if (cfgDoReverseScaling.value) { energyCorrectionFactor = 1.0505f; From 243024cd7ff1d0ed4e7afabc0b54bb442fe98223 Mon Sep 17 00:00:00 2001 From: Marvin Hemmer Date: Wed, 4 Feb 2026 09:46:44 +0100 Subject: [PATCH 2/2] fix missing argument in calibETaskEmc.cxx --- PWGEM/PhotonMeson/Tasks/calibTaskEmc.cxx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/PWGEM/PhotonMeson/Tasks/calibTaskEmc.cxx b/PWGEM/PhotonMeson/Tasks/calibTaskEmc.cxx index 91469523448..13e4cf66313 100644 --- a/PWGEM/PhotonMeson/Tasks/calibTaskEmc.cxx +++ b/PWGEM/PhotonMeson/Tasks/calibTaskEmc.cxx @@ -702,7 +702,7 @@ struct CalibTaskEmc { int nColl = 1; EMBitFlags emcFlags(clusters.size()); - fEMCCut.AreSelectedRunning(emcFlags, clusters, matchedPrims, matchedSeconds); + fEMCCut.AreSelectedRunning(emcFlags, clusters, matchedPrims, matchedSeconds, ®istry); for (const auto& collision : collisions) { auto photonsPerCollision = clusters.sliceBy(perCollisionEMC, collision.globalIndex()); @@ -812,7 +812,7 @@ struct CalibTaskEmc { int nColl = 1; EMBitFlags emcFlags(clusters.size()); - fEMCCut.AreSelectedRunning(emcFlags, clusters, matchedPrims, matchedSeconds); + fEMCCut.AreSelectedRunning(emcFlags, clusters, matchedPrims, matchedSeconds, ®istry); for (const auto& collision : collisions) { o2::aod::pwgem::photonmeson::utils::eventhistogram::fillEventInfo<0>(®istry, collision); @@ -917,7 +917,7 @@ struct CalibTaskEmc { { float energyCorrectionFactor = 1.f; EMBitFlags emcFlags(clusters.size()); - fEMCCut.AreSelectedRunning(emcFlags, clusters, matchedPrims, matchedSeconds); + fEMCCut.AreSelectedRunning(emcFlags, clusters, matchedPrims, matchedSeconds, ®istry); SameKindPair pair{binningOnPositions, mixingConfig.cfgMixingDepth, -1, &cache}; // indicates that 5 events should be mixed and under/overflow (-1) to be ignored @@ -998,7 +998,7 @@ struct CalibTaskEmc { { float energyCorrectionFactor = 1.f; EMBitFlags emcFlags(clusters.size()); - fEMCCut.AreSelectedRunning(emcFlags, clusters, matchedPrims, matchedSeconds); + fEMCCut.AreSelectedRunning(emcFlags, clusters, matchedPrims, matchedSeconds, ®istry); for (const auto& [c1, photonEMC, c2, photonPCM] : pairPCMEMC) { if (!(fEMEventCut.IsSelected(c1)) || !(fEMEventCut.IsSelected(c2))) {