From 598928ad66e22bd9377a7556c2937f20e799cf94 Mon Sep 17 00:00:00 2001 From: Marvin Hemmer Date: Wed, 4 Feb 2026 15:36:55 +0100 Subject: [PATCH 1/2] [PWGEM] PhotonMeson: change more cutlibs to use concepts - Use is_table and is_iterator concepts in DalitzEE and V0Photon Cut libs. - Add PCM to flow task --- PWGEM/PhotonMeson/Core/DalitzEECut.h | 15 ++- PWGEM/PhotonMeson/Core/V0PhotonCut.h | 16 +++ PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx | 147 ++++++++++++++++++++- 3 files changed, 166 insertions(+), 12 deletions(-) diff --git a/PWGEM/PhotonMeson/Core/DalitzEECut.h b/PWGEM/PhotonMeson/Core/DalitzEECut.h index c238b5bead3..9eb6b1a0801 100644 --- a/PWGEM/PhotonMeson/Core/DalitzEECut.h +++ b/PWGEM/PhotonMeson/Core/DalitzEECut.h @@ -20,6 +20,7 @@ #include "PWGEM/Dilepton/Utils/PairUtilities.h" #include +#include #include // IWYU pragma: keep #include @@ -72,7 +73,7 @@ class DalitzEECut : public TNamed kTPConly = 1, }; - template + template bool IsSelected(TTrack1 const& t1, TTrack2 const& t2, float bz) const { if (!IsSelectedTrack(t1) || !IsSelectedTrack(t2)) { @@ -86,7 +87,7 @@ class DalitzEECut : public TNamed return true; } - template + template bool IsSelectedPair(TTrack1 const& t1, TTrack2 const& t2, const float bz) const { ROOT::Math::PtEtaPhiMVector v1(t1.pt(), t1.eta(), t1.phi(), o2::constants::physics::MassElectron); @@ -108,7 +109,7 @@ class DalitzEECut : public TNamed return true; } - template + template bool IsSelectedTrack(TTrack const& track, TCollision const& = 0) const { if (!track.hasITS()) { @@ -194,7 +195,7 @@ class DalitzEECut : public TNamed return true; } - template + template bool PassPID(T const& track) const { switch (mPIDScheme) { @@ -212,7 +213,7 @@ class DalitzEECut : public TNamed } } - template + template bool PassTPConly(T const& track) const { bool is_el_included_TPC = mMinTPCNsigmaEl < track.tpcNSigmaEl() && track.tpcNSigmaEl() < mMaxTPCNsigmaEl; @@ -220,7 +221,7 @@ class DalitzEECut : public TNamed return is_el_included_TPC && is_pi_excluded_TPC; } - template + template bool PassTOFif(T const& track) const { bool is_el_included_TPC = mMinTPCNsigmaEl < track.tpcNSigmaEl() && track.tpcNSigmaEl() < mMaxTPCNsigmaEl; @@ -229,7 +230,7 @@ class DalitzEECut : public TNamed return is_el_included_TPC && is_pi_excluded_TPC && is_el_included_TOF; } - template + template bool IsSelectedTrack(T const& track, const DalitzEECuts& cut) const { switch (cut) { diff --git a/PWGEM/PhotonMeson/Core/V0PhotonCut.h b/PWGEM/PhotonMeson/Core/V0PhotonCut.h index 1a1567c454c..1bbcf5c4103 100644 --- a/PWGEM/PhotonMeson/Core/V0PhotonCut.h +++ b/PWGEM/PhotonMeson/Core/V0PhotonCut.h @@ -16,6 +16,7 @@ #ifndef PWGEM_PHOTONMESON_CORE_V0PHOTONCUT_H_ #define PWGEM_PHOTONMESON_CORE_V0PHOTONCUT_H_ +#include "PWGEM/PhotonMeson/Core/EMBitFlags.h" #include "PWGEM/PhotonMeson/Core/EmMlResponsePCM.h" #include "PWGEM/PhotonMeson/Core/V0PhotonCandidate.h" #include "PWGEM/PhotonMeson/Utils/TrackSelection.h" @@ -193,6 +194,21 @@ class V0PhotonCut : public TNamed kNCuts }; + /// \brief check if given v0 photon survives all cuts + /// \param flags EMBitFlags where results will be stored + /// \param v0s v0 photon table to check + template + void AreSelectedRunning(EMBitFlags& flags, TV0 const& v0s) const + { + size_t iV0 = 0; + for (const auto& v0 : v0s) { + if (!IsSelected(v0)) { + flags.set(iV0); + } + ++iV0; + } + } + template bool IsSelected(TV0 const& v0) const { diff --git a/PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx b/PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx index 0b4f156f27b..d61205a6c12 100644 --- a/PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx +++ b/PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx @@ -1240,6 +1240,9 @@ struct TaskPi0FlowEMC { EMBitFlags emcFlags(clusters.size()); fEMCCut.AreSelectedRunning(emcFlags, clusters, matchedPrims, matchedSeconds); + EMBitFlags v0flags(photons.size()); + fV0PhotonCut.AreSelectedRunning(v0flags, photons); + if (cfgDoReverseScaling.value) { energyCorrectionFactor = 1.0505f; } @@ -1267,7 +1270,7 @@ struct TaskPi0FlowEMC { } } for (const auto& [g1, g2] : combinations(CombinationsFullIndexPolicy(photonsEMCPerCollision, photonsPCMPerCollision))) { - if (!(emcFlags.test(g1.globalIndex())) || !(fV0PhotonCut.IsSelected(g2))) { + if (!(emcFlags.test(g1.globalIndex())) || !(v0flags.test(g2.globalIndex()))) { continue; } @@ -1304,7 +1307,7 @@ struct TaskPi0FlowEMC { registry.fill(HIST("hMesonCuts"), 4); continue; } - if (mesonConfig.cfgEnableQA) { + if (mesonConfig.cfgEnableQA.value) { registry.fill(HIST("mesonQA/hInvMassPt"), vMeson.M(), vMeson.Pt()); registry.fill(HIST("mesonQA/hTanThetaPhi"), vMeson.M(), getAngleDegree(std::atan(dTheta / dPhi))); registry.fill(HIST("mesonQA/hAlphaPt"), (v1.E() - v2.E()) / (v1.E() + v2.E()), vMeson.Pt()); @@ -1332,6 +1335,10 @@ struct TaskPi0FlowEMC { float energyCorrectionFactor = 1.f; EMBitFlags emcFlags(clusters.size()); fEMCCut.AreSelectedRunning(emcFlags, clusters, matchedPrims, matchedSeconds); + + EMBitFlags v0flags(pcmPhotons.size()); + fV0PhotonCut.AreSelectedRunning(v0flags, pcmPhotons); + for (const auto& [c1, photonEMC, c2, photonPCM] : pairPCMEMC) { if (!(fEMEventCut.IsSelected(c1)) || !(fEMEventCut.IsSelected(c2))) { // general event selection @@ -1352,7 +1359,7 @@ struct TaskPi0FlowEMC { } registry.fill(HIST("h3DMixingCount"), c1.posZ(), getCentrality(c1), c1.ep2ft0m()); for (const auto& [g1, g2] : combinations(CombinationsFullIndexPolicy(photonEMC, photonPCM))) { - if (!(emcFlags.test(g1.globalIndex())) || !(fV0PhotonCut.IsSelected(g2))) { + if (!(emcFlags.test(g1.globalIndex())) || !(v0flags.test(g2.globalIndex()))) { continue; } // Cut edge clusters away, similar to rotation method to ensure same acceptance is used @@ -1383,7 +1390,7 @@ struct TaskPi0FlowEMC { registry.fill(HIST("hMesonCutsMixed"), 4); continue; } - if (mesonConfig.cfgEnableQA) { + if (mesonConfig.cfgEnableQA.value) { registry.fill(HIST("mesonQA/hInvMassPtMixed"), vMeson.M(), vMeson.Pt()); registry.fill(HIST("mesonQA/hTanThetaPhiMixed"), vMeson.M(), getAngleDegree(std::atan(dTheta / dPhi))); registry.fill(HIST("mesonQA/hAlphaPtMixed"), (v1.E() - v2.E()) / (v1.E() + v2.E()), vMeson.Pt()); @@ -1402,6 +1409,9 @@ struct TaskPi0FlowEMC { // Pi0 from EMCal void processM02(CollsWithQvecs const& collisions, EMCalPhotons const& clusters, MinMTracks const& matchedPrims, MinMSTracks const& matchedSeconds) { + EMBitFlags emcFlags(clusters.size()); + fEMCCut.AreSelectedRunning(emcFlags, clusters, matchedPrims, matchedSeconds); + for (const auto& collision : collisions) { o2::aod::pwgem::photonmeson::utils::eventhistogram::fillEventInfo<0>(®istry, collision); if (!(fEMEventCut.IsSelected(collision))) { @@ -1439,7 +1449,7 @@ struct TaskPi0FlowEMC { } auto matchedPrimsPerCluster = matchedPrims.sliceBy(perEMCClusterMT, photon.globalIndex()); auto matchedSecondsPerCluster = matchedSeconds.sliceBy(perEMCClusterMS, photon.globalIndex()); - if (!(fEMCCut.IsSelected(photon, matchedPrimsPerCluster, matchedSecondsPerCluster))) { + if (!(emcFlags.test(photon.globalIndex()))) { continue; } if (cfgDistanceToEdge.value && (checkEtaPhi1D(photon.eta(), RecoDecay::constrainAngle(photon.phi())) >= cfgEMCalMapLevelSameEvent.value)) { @@ -1472,6 +1482,133 @@ struct TaskPi0FlowEMC { } // processM02 PROCESS_SWITCH(TaskPi0FlowEMC, processM02, "Process single EMCal clusters as function of M02", false); + // Pi0 from EMCal + void processPCM(CollsWithQvecs const& collisions, PCMPhotons const& photons, aod::V0Legs const&) + { + EMBitFlags v0flags(photons.size()); + fV0PhotonCut.AreSelectedRunning(v0flags, photons); + for (const auto& collision : collisions) { + + if (!isFullEventSelected(collision)) { + continue; + } + runNow = collision.runNumber(); + if (runNow != runBefore) { + initCCDB(collision); + runBefore = runNow; + } + + auto photonsPerCollision = photons.sliceBy(perCollisionPCM, collision.globalIndex()); + for (const auto& [g1, g2] : combinations(CombinationsStrictlyUpperIndexPolicy(photonsPerCollision, photonsPerCollision))) { + if (!(v0flags.test(g1.globalIndex())) || !(v0flags.test(g2.globalIndex()))) { + continue; + } + if (correctionConfig.cfgEnableNonLin.value) { + energyCorrectionFactor = fEMCalCorrectionFactor->Eval(g1.e() > MinEnergy ? g1.e() : MinEnergy); + } + if (correctionConfig.cfgEnableNonLin.value) { + energyCorrectionFactor = fEMCalCorrectionFactor->Eval(g2.e() > MinEnergy ? g2.e() : MinEnergy); + } + + ROOT::Math::PtEtaPhiMVector v1(energyCorrectionFactor * g1.pt(), g1.eta(), g1.phi(), 0.); + ROOT::Math::PtEtaPhiMVector v2(energyCorrectionFactor * g2.pt(), g2.eta(), g2.phi(), 0.); + ROOT::Math::PtEtaPhiMVector vMeson = v1 + v2; + + float openingAngle = std::acos(v1.Vect().Dot(v2.Vect()) / (v1.P() * v2.P())); + + registry.fill(HIST("hMesonCuts"), 1); + if (openingAngle <= mesonConfig.minOpenAngle) { + registry.fill(HIST("hMesonCuts"), 2); + continue; + } + if (thnConfigAxisInvMass.value[1] > vMeson.M() || thnConfigAxisInvMass.value.back() < vMeson.M()) { + registry.fill(HIST("hMesonCuts"), 3); + continue; + } + if (thnConfigAxisPt.value[1] > vMeson.Pt() || thnConfigAxisPt.value.back() < vMeson.Pt()) { + registry.fill(HIST("hMesonCuts"), 4); + continue; + } + if (mesonConfig.cfgEnableQA.value) { + registry.fill(HIST("mesonQA/hInvMassPt"), vMeson.M(), vMeson.Pt()); + registry.fill(HIST("mesonQA/hAlphaPt"), (v1.E() - v2.E()) / (v1.E() + v2.E()), vMeson.Pt()); + } + registry.fill(HIST("hMesonCuts"), 6); + runFlowAnalysis<0>(collision, vMeson); + } + } // end of loop over collisions + } + PROCESS_SWITCH(TaskPi0FlowEMC, processPCM, "Process PCM Pi0 candidates", false); + + // PCM-EMCal mixed event + void processPCMMixed(FilteredCollsWithQvecs const& collisions, PCMPhotons const& pcmPhotons, aod::V0Legs const&) + { + + using BinningType = ColumnBinningPolicy>; + BinningType binningMixedEvent{{mixingConfig.cfgVtxBins, mixingConfig.cfgCentBins, mixingConfig.cfgEPBins}, true}; + + auto pcmPhotonTuple = std::make_tuple(pcmPhotons); + SameKindPair pair{binningMixedEvent, mixingConfig.cfgMixingDepth, -1, collisions, pcmPhotonTuple, &cache}; // indicates that 5 events should be mixed and under/overflow (-1) to be ignored + + float energyCorrectionFactor = 1.f; + + EMBitFlags v0flags(pcmPhotons.size()); + fV0PhotonCut.AreSelectedRunning(v0flags, pcmPhotons); + + for (const auto& [c1, photon1, c2, photon2] : pair) { + if (!(fEMEventCut.IsSelected(c1)) || !(fEMEventCut.IsSelected(c2))) { + // general event selection + continue; + } + if (!(eventcuts.cfgFT0COccupancyMin <= c1.ft0cOccupancyInTimeRange() && c1.ft0cOccupancyInTimeRange() < eventcuts.cfgFT0COccupancyMax) || !(eventcuts.cfgFT0COccupancyMin <= c2.ft0cOccupancyInTimeRange() && c2.ft0cOccupancyInTimeRange() < eventcuts.cfgFT0COccupancyMax)) { + // occupancy selection + continue; + } + if (getCentrality(c1) < eventcuts.cfgMinCent || getCentrality(c1) > eventcuts.cfgMaxCent || getCentrality(c2) < eventcuts.cfgMinCent || getCentrality(c2) > eventcuts.cfgMaxCent) { + // event selection + continue; + } + runNow = c1.runNumber(); + if (runNow != runBefore) { + initCCDB(c1); + runBefore = runNow; + } + registry.fill(HIST("h3DMixingCount"), c1.posZ(), getCentrality(c1), c1.ep2ft0m()); + for (const auto& [g1, g2] : combinations(CombinationsFullIndexPolicy(photon1, photon2))) { + if (!(v0flags.test(g1.globalIndex())) || !(v0flags.test(g2.globalIndex()))) { + continue; + } + energyCorrectionFactor = fEMCalCorrectionFactor->Eval(g1.e() > MinEnergy ? g1.e() : MinEnergy); + ROOT::Math::PtEtaPhiMVector v1(energyCorrectionFactor * g1.pt(), g1.eta(), g1.phi(), 0.); + ROOT::Math::PtEtaPhiMVector v2(g2.pt(), g2.eta(), g2.phi(), 0.); + ROOT::Math::PtEtaPhiMVector vMeson = v1 + v2; + + float openingAngle = std::acos(v1.Vect().Dot(v2.Vect()) / (v1.P() * v2.P())); + + registry.fill(HIST("hMesonCutsMixed"), 1); + if (openingAngle <= mesonConfig.minOpenAngle) { + registry.fill(HIST("hMesonCutsMixed"), 2); + continue; + } + if (thnConfigAxisInvMass.value[1] > vMeson.M() || thnConfigAxisInvMass.value.back() < vMeson.M()) { + registry.fill(HIST("hMesonCutsMixed"), 3); + continue; + } + if (thnConfigAxisPt.value[1] > vMeson.Pt() || thnConfigAxisPt.value.back() < vMeson.Pt()) { + registry.fill(HIST("hMesonCutsMixed"), 4); + continue; + } + if (mesonConfig.cfgEnableQA.value) { + registry.fill(HIST("mesonQA/hInvMassPtMixed"), vMeson.M(), vMeson.Pt()); + registry.fill(HIST("mesonQA/hAlphaPtMixed"), (v1.E() - v2.E()) / (v1.E() + v2.E()), vMeson.Pt()); + } + registry.fill(HIST("hMesonCutsMixed"), 6); + runFlowAnalysis<2>(c1, vMeson); + } + } + } + PROCESS_SWITCH(TaskPi0FlowEMC, processPCMMixed, "Process neutral meson flow using PCM-EMC mixed event", false); + }; // End struct TaskPi0FlowEMC WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) From 896777ff007047718ca07b8d8a0c67fe186133ee Mon Sep 17 00:00:00 2001 From: Marvin Hemmer Date: Wed, 4 Feb 2026 15:44:51 +0100 Subject: [PATCH 2/2] [PWGEM] PhotonMeson V0PhotonCuts.h - Fix missing std prefix and using `o2::constants::math::Deg2Rad` instead of ROOT function - Remove `using namespace` in header --- PWGEM/PhotonMeson/Core/V0PhotonCut.h | 41 ++++++++++++++-------------- 1 file changed, 20 insertions(+), 21 deletions(-) diff --git a/PWGEM/PhotonMeson/Core/V0PhotonCut.h b/PWGEM/PhotonMeson/Core/V0PhotonCut.h index 1bbcf5c4103..45fae0defbd 100644 --- a/PWGEM/PhotonMeson/Core/V0PhotonCut.h +++ b/PWGEM/PhotonMeson/Core/V0PhotonCut.h @@ -22,6 +22,7 @@ #include "PWGEM/PhotonMeson/Utils/TrackSelection.h" #include +#include #include #include @@ -41,8 +42,6 @@ #include #include -using namespace o2::pwgem::photonmeson; - namespace o2::analysis { @@ -263,7 +262,7 @@ class V0PhotonCut : public TNamed auto pos = v0.template posTrack_as(); auto ele = v0.template negTrack_as(); - for (auto& track : {pos, ele}) { + for (const auto& track : {pos, ele}) { if (!IsSelectedTrack(track, V0PhotonCuts::kTrackPtRange)) { return false; } @@ -279,10 +278,10 @@ class V0PhotonCut : public TNamed if (!track.hasITS() && !track.hasTPC()) { // track has to be ITSonly or TPConly or ITS-TPC return false; } - if (mDisableITSonly && isITSonlyTrack(track)) { + if (mDisableITSonly && o2::pwgem::photonmeson::isITSonlyTrack(track)) { return false; } - if (mDisableTPConly && isTPConlyTrack(track)) { + if (mDisableTPConly && o2::pwgem::photonmeson::isTPConlyTrack(track)) { return false; } @@ -290,13 +289,13 @@ class V0PhotonCut : public TNamed auto hits_ib = std::count_if(its_ib_Requirement.second.begin(), its_ib_Requirement.second.end(), [&](auto&& requiredLayer) { return track.itsClusterMap() & (1 << requiredLayer); }); auto hits_ob = std::count_if(its_ob_Requirement.second.begin(), its_ob_Requirement.second.end(), [&](auto&& requiredLayer) { return track.itsClusterMap() & (1 << requiredLayer); }); bool its_ob_only = (hits_ib <= its_ib_Requirement.first) && (hits_ob >= its_ob_Requirement.first); - if (isITSonlyTrack(track) && !its_ob_only) { // ITSonly tracks should not have any ITSib hits. + if (o2::pwgem::photonmeson::isITSonlyTrack(track) && !its_ob_only) { // ITSonly tracks should not have any ITSib hits. return false; } auto hits_ob_itstpc = std::count_if(its_ob_Requirement_ITSTPC.second.begin(), its_ob_Requirement_ITSTPC.second.end(), [&](auto&& requiredLayer) { return track.itsClusterMap() & (1 << requiredLayer); }); bool its_ob_only_itstpc = (hits_ib <= its_ib_Requirement.first) && (hits_ob_itstpc >= its_ob_Requirement_ITSTPC.first); - if (isITSTPCTrack(track) && !its_ob_only_itstpc) { // ITSTPC tracks should not have any ITSib hits. + if (o2::pwgem::photonmeson::isITSTPCTrack(track) && !its_ob_only_itstpc) { // ITSTPC tracks should not have any ITSib hits. return false; } } @@ -416,7 +415,7 @@ class V0PhotonCut : public TNamed return v0.mGamma() < mMaxQt * 2.f; case V0PhotonCuts::kAP: - return pow(v0.alpha() / mMaxAlpha, 2) + pow(v0.qtarm() / mMaxQt, 2) < 1.0; + return std::pow(v0.alpha() / mMaxAlpha, 2) + std::pow(v0.qtarm() / mMaxQt, 2) < 1.0; case V0PhotonCuts::kPsiPair: return true; @@ -454,24 +453,24 @@ class V0PhotonCut : public TNamed float y = v0.vy(); // cm, measured secondary vertex of gamma->ee float z = v0.vz(); // cm, measured secondary vertex of gamma->ee - float rxy = sqrt(x * x + y * y); + float rxy = std::sqrt(x * x + y * y); if (rxy < 7.0 || 14.0 < rxy) { return false; } - // r = 0.192 * z + 8.88 (cm) expected wire position in RZ plane.TMath::Tan(10.86 * TMath::DegToRad()) = 0.192 + // r = 0.192 * z + 8.88 (cm) expected wire position in RZ plane.TMath::Tan(10.86 * o2::constants::math::Deg2Rad) = 0.192 if (rxy > 0.192 * z + 14.0) { // upper limit return false; } - float dxy = std::fabs(1.0 * y - x * std::tan(-8.52 * TMath::DegToRad())) / sqrt(pow(1.0, 2) + pow(std::tan(-8.52 * TMath::DegToRad()), 2)); + float dxy = std::fabs(1.0 * y - x * std::tan(-8.52f * o2::constants::math::Deg2Rad)) / std::sqrt(std::pow(1.0, 2) + std::pow(std::tan(-8.52f * o2::constants::math::Deg2Rad), 2)); return !(dxy > margin_xy); } case V0PhotonCuts::kOnWwireOB: { - const float margin_xy = 1.0; // cm - const float rxy_exp = 30.8; // cm - const float x_exp = rxy_exp * std::cos(-1.3 * TMath::DegToRad()); // cm, expected position x of W wire - const float y_exp = rxy_exp * std::sin(-1.3 * TMath::DegToRad()); // cm, expected position y of W wire + const float margin_xy = 1.0; // cm + const float rxy_exp = 30.8; // cm + const float x_exp = rxy_exp * std::cos(-1.3f * o2::constants::math::Deg2Rad); // cm, expected position x of W wire + const float y_exp = rxy_exp * std::sin(-1.3f * o2::constants::math::Deg2Rad); // cm, expected position y of W wire // const float z_min = -47.0; // cm // const float z_max = +47.0; // cm float x = v0.vx(); // cm, measured secondary vertex of gamma->ee @@ -536,26 +535,26 @@ class V0PhotonCut : public TNamed return mMinChi2PerClusterITS < track.itsChi2NCl() && track.itsChi2NCl() < mMaxChi2PerClusterITS; case V0PhotonCuts::kITSClusterSize: { - if (!isITSonlyTrack(track)) { + if (!o2::pwgem::photonmeson::isITSonlyTrack(track)) { return true; } return mMinMeanClusterSizeITS < track.meanClusterSizeITSob() * std::cos(std::atan(track.tgl())) && track.meanClusterSizeITSob() * std::cos(std::atan(track.tgl())) < mMaxMeanClusterSizeITS; } case V0PhotonCuts::kRequireITSTPC: - return isITSTPCTrack(track); + return o2::pwgem::photonmeson::isITSTPCTrack(track); case V0PhotonCuts::kRequireITSonly: - return isITSonlyTrack(track); + return o2::pwgem::photonmeson::isITSonlyTrack(track); case V0PhotonCuts::kRequireTPConly: - return isTPConlyTrack(track); + return o2::pwgem::photonmeson::isTPConlyTrack(track); case V0PhotonCuts::kRequireTPCTRD: - return isTPCTRDTrack(track); + return o2::pwgem::photonmeson::isTPCTRDTrack(track); case V0PhotonCuts::kRequireTPCTOF: - return isTPCTOFTrack(track); + return o2::pwgem::photonmeson::isTPCTOFTrack(track); default: return false;