From 7ba00fa36890b8ef0afafc82626158c3c25ec4df Mon Sep 17 00:00:00 2001 From: Marvin Hemmer Date: Wed, 29 Jul 2026 11:10:48 +0200 Subject: [PATCH 1/2] [PWGEM] Add option to separate EMCal and DCal cluster pairing for flow task [PWGEM] PM: Add amplitude of McParticles to derived data - Add column of vectors to EmEmcMcParticle table which stores the amplitude fraction that each EmMcParticle contributed to the total cluster energy - Add three new TruthClasses in `emcalPhotonMcTask.cxx`: - PhotonElectronBS - SplitPhotonCluster - SplitLeptonCluster - Add new function `isFromBremsstrahlung` to `MCUtilities.h` --- .../PhotonMeson/DataModel/GammaTablesRedux.h | 2 +- PWGEM/PhotonMeson/DataModel/gammaTables.h | 25 +-- .../TableProducer/associateMCinfoPhoton.cxx | 36 ++-- .../TableProducer/skimmerGammaCalo.cxx | 13 +- PWGEM/PhotonMeson/Tasks/emcalPhotonMcTask.cxx | 164 ++++++++++++------ PWGEM/PhotonMeson/Tasks/taskFlowReso.cxx | 10 +- PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx | 32 +++- PWGEM/PhotonMeson/Tasks/testTaskEmc.cxx | 4 +- PWGEM/PhotonMeson/Utils/MCUtilities.h | 18 ++ 9 files changed, 208 insertions(+), 96 deletions(-) diff --git a/PWGEM/PhotonMeson/DataModel/GammaTablesRedux.h b/PWGEM/PhotonMeson/DataModel/GammaTablesRedux.h index d666265af2a..d9d2485206d 100644 --- a/PWGEM/PhotonMeson/DataModel/GammaTablesRedux.h +++ b/PWGEM/PhotonMeson/DataModel/GammaTablesRedux.h @@ -47,7 +47,7 @@ enum Observable { }; // Values in tables are stored in downscaled format to save disk space -const std::array downscalingFactors{ +constexpr std::array downscalingFactors{ 1E0, // Cluster definition 1E3, // Cluster energy 1E4, // Cluster eta diff --git a/PWGEM/PhotonMeson/DataModel/gammaTables.h b/PWGEM/PhotonMeson/DataModel/gammaTables.h index 9ef4b45727b..ca46c23329f 100644 --- a/PWGEM/PhotonMeson/DataModel/gammaTables.h +++ b/PWGEM/PhotonMeson/DataModel/gammaTables.h @@ -102,23 +102,31 @@ using V0LegMCLabel = V0LegMCLabels::iterator; // 3. EMEMCClusterMCLabels: Vector of EM MC particle ID of largest contributor to cluster namespace emcclustermclabel { -DECLARE_SOA_ARRAY_INDEX_COLUMN(McParticle, emmcparticle) //! +DECLARE_SOA_ARRAY_INDEX_COLUMN(McParticle, mcParticle) //! Array of indice of McParticles +DECLARE_SOA_COLUMN(Amplitude, amplitude, std::vector); //! vector of aplitudes stored in cluster. Ordering is identical to the ordering of EMMCParticleId } // namespace emcclustermclabel // NOTE: MC labels. This table has one vector of global mc particle ids for each reconstructed emc cluster (joinable with emccluster table) DECLARE_SOA_TABLE(EMCClusterMCLabels, "AOD", "EMCClsMCLABEL", //! emcclustermclabel::McParticleIds); -using EMCClusterMCLabel = EMCClusterMCLabels::iterator; + +DECLARE_SOA_TABLE_VERSIONED(EMCClusterMCLabels_001, "AOD", "EMCClsMCLABEL", 1, //! + emcclustermclabel::McParticleIds, emcclustermclabel::Amplitude); +using EMCClusterMCLabel = EMCClusterMCLabels_001::iterator; namespace ememcclustermclabel { -DECLARE_SOA_ARRAY_INDEX_COLUMN(EMMCParticle, emmcparticle); //! +DECLARE_SOA_ARRAY_INDEX_COLUMN(EMMCParticle, emmcparticle); //! Array of indice of EmMcParticles +DECLARE_SOA_COLUMN(Amplitude, amplitude, std::vector); //! vector of aplitudes stored in cluster. Ordering is identical to the ordering of EMMCParticleId } // namespace ememcclustermclabel // NOTE: MC labels. This table has a vector of entries for each reconstructed emc cluster (joinable with emccluster table) DECLARE_SOA_TABLE(EMEMCClusterMCLabels, "AOD", "EMEMCClsMCLABEL", //! ememcclustermclabel::EMMCParticleIds); -using EMEMCClusterMCLabel = EMEMCClusterMCLabels::iterator; + +DECLARE_SOA_TABLE_VERSIONED(EMEMCClusterMCLabels_001, "AOD", "EMEMCClsMCLABEL", 1, //! + ememcclustermclabel::EMMCParticleIds, ememcclustermclabel::Amplitude); +using EMEMCClusterMCLabel = EMEMCClusterMCLabels_001::iterator; namespace v0leg { @@ -146,9 +154,8 @@ DECLARE_SOA_DYNAMIC_COLUMN(MeanClusterSizeITS, meanClusterSizeITS, [](uint32_t i } if (nl > 0) { return static_cast(total_cluster_size) / static_cast(nl); - } else { - return 0; } + return 0; }); DECLARE_SOA_DYNAMIC_COLUMN(MeanClusterSizeITSib, meanClusterSizeITSib, [](uint32_t itsClusterSizes) -> float { int total_cluster_size = 0, nl = 0; @@ -161,9 +168,8 @@ DECLARE_SOA_DYNAMIC_COLUMN(MeanClusterSizeITSib, meanClusterSizeITSib, [](uint32 } if (nl > 0) { return static_cast(total_cluster_size) / static_cast(nl); - } else { - return 0; } + return 0; }); DECLARE_SOA_DYNAMIC_COLUMN(MeanClusterSizeITSob, meanClusterSizeITSob, [](uint32_t itsClusterSizes) -> float { int total_cluster_size = 0, nl = 0; @@ -176,9 +182,8 @@ DECLARE_SOA_DYNAMIC_COLUMN(MeanClusterSizeITSob, meanClusterSizeITSob, [](uint32 } if (nl > 0) { return static_cast(total_cluster_size) / static_cast(nl); - } else { - return 0; } + return 0; }); } // namespace v0leg DECLARE_SOA_TABLE(V0Legs_000, "AOD", "V0LEG", //! diff --git a/PWGEM/PhotonMeson/TableProducer/associateMCinfoPhoton.cxx b/PWGEM/PhotonMeson/TableProducer/associateMCinfoPhoton.cxx index 0034d0dc3b6..30e3cb4cd3c 100644 --- a/PWGEM/PhotonMeson/TableProducer/associateMCinfoPhoton.cxx +++ b/PWGEM/PhotonMeson/TableProducer/associateMCinfoPhoton.cxx @@ -10,8 +10,8 @@ // or submit itself to any jurisdiction. /// \file associateMCinfoPhoton.cxx -/// \brief This code produces reduced events for photon analyses -/// \author Daiki Sekihata (daiki.sekihata@cern.ch) +/// \brief This code produces EmMc tables were the McParticleIds get reshuffled due to not storing all McParticles +/// \author Daiki Sekihata (daiki.sekihata@cern.ch), Marvin Hemmer (marvin.hemmer@cern.ch), Nicolas Strangmann (nicolas.strangmann@cern.ch) #include "PWGEM/PhotonMeson/DataModel/GammaTablesRedux.h" #include "PWGEM/PhotonMeson/DataModel/gammaTables.h" @@ -57,7 +57,7 @@ using namespace o2::constants::physics; using MyCollisionsMC = soa::Join; using TracksMC = soa::Join; using FwdTracksMC = soa::Join; -using MyEMCClusters = soa::Join; +using MyEMCClusters = soa::Join; struct Counter { int particles{0}; @@ -77,7 +77,7 @@ struct AssociateMCInfoPhoton { Produces emmcparticles; Produces v0legmclabels; Produces emprimaryelectronmclabels; - Produces ememcclustermclabels; + Produces ememcclustermclabels; Produces binnedGenPt; @@ -100,6 +100,7 @@ struct AssociateMCInfoPhoton { // !!Don't change pt,eta,y binning. These binnings have to be consistent with binned data at analysis.!! std::vector ptbins; + ptbins.reserve(72); for (int i = 0; i < 2; i++) { // o2-linter: disable=magic-number (just numbers for binning) ptbins.emplace_back(0.05 * (i - 0) + 0.0); // from 0 to 0.05 GeV/c, every 0.05 GeV/c } @@ -116,7 +117,7 @@ struct AssociateMCInfoPhoton { const AxisSpec axisRapidity{{0.0, +0.8, +0.9}, "rapidity |y|"}; static constexpr uint NParticleNames = 9; - static constexpr std::string_view ParticleNames[NParticleNames] = { + static constexpr std::array ParticleNames = { "Gamma", "Pi0", "Eta", "Omega", "Phi", "ChargedPion", "ChargedKaon", "K0S", "Lambda"}; @@ -267,8 +268,9 @@ struct AssociateMCInfoPhoton { genGamma[binNumber]++; break; case PDG_t::kPi0: - if (requireGammaGammaDecay && !isGammaGammaDecay(mcParticle, mcParticles)) + if (requireGammaGammaDecay && !isGammaGammaDecay(mcParticle, mcParticles)) { continue; + } registry.fill(HIST("Generated/h2PtY_Pi0"), mcParticle.pt(), std::fabs(mcParticle.y())); genPi0[binNumber]++; if (isMesonAccepted) { @@ -276,8 +278,9 @@ struct AssociateMCInfoPhoton { } break; case Pdg::kEta: - if (requireGammaGammaDecay && !isGammaGammaDecay(mcParticle, mcParticles)) + if (requireGammaGammaDecay && !isGammaGammaDecay(mcParticle, mcParticles)) { continue; + } registry.fill(HIST("Generated/h2PtY_Eta"), mcParticle.pt(), std::fabs(mcParticle.y())); genEta[binNumber]++; if (isMesonAccepted) { @@ -291,7 +294,7 @@ struct AssociateMCInfoPhoton { } // end of mc track loop // make an entry for this MC event only if it was not already added to the table - if (!(fEventLabels.find(mcCollisionIter.globalIndex()) != fEventLabels.end())) { + if (!fEventLabels.contains(mcCollisionIter.globalIndex())) { mcevents(mcCollisionIter.globalIndex(), mcCollisionIter.generatorsID(), mcCollisionIter.posX(), mcCollisionIter.posY(), mcCollisionIter.posZ(), mcCollisionIter.impactParameter(), mcCollisionIter.eventPlaneAngle()); fEventLabels[mcCollisionIter.globalIndex()] = fCounter.events; fCounter.events++; @@ -438,15 +441,17 @@ struct AssociateMCInfoPhoton { mcCollisionIter.setCursor(collisionIter.mcCollisionId()); // TODO: test - if (emccluster.emmcparticleIds().size() <= 0) { + if (emccluster.mcParticleIds().size() <= 0) { continue; } std::vector vEmcMcParticleIds; + std::vector vAmplitudes; - vEmcMcParticleIds.reserve(emccluster.emmcparticleIds().size()); + vEmcMcParticleIds.reserve(emccluster.mcParticleIds().size()); + vAmplitudes.reserve(emccluster.mcParticleIds().size()); - for (const auto& emcParticleId : emccluster.emmcparticleIds()) { - mcPhoton.setCursor(emcParticleId); + for (size_t iCont = 0; iCont < emccluster.mcParticleIds().size(); iCont++) { + mcPhoton.setCursor(emccluster.mcParticleIds()[iCont]); // if the MC truth particle corresponding to this reconstructed track which is not already written, add it to the skimmed MC stack auto [iter, isNew] = fNewLabels.try_emplace(mcPhoton.globalIndex(), fCounter.particles); @@ -456,6 +461,7 @@ struct AssociateMCInfoPhoton { fCounter.particles++; } vEmcMcParticleIds.emplace_back(iter->second); + vAmplitudes.emplace_back(emccluster.amplitude()[iCont]); // ememcclustermclabels(fNewLabels.find(mcPhoton.index())->second); // Next, store mother-chain of this reconstructed track. @@ -465,7 +471,7 @@ struct AssociateMCInfoPhoton { } selectMothersToStore(motherid, mcParticles.size(), motherParticle, daughterIter, mcCollisionIter, fNewLabels, fNewLabelsReversed, fEventIdx, fEventLabels, fCounter); } // end of loop over mc particles of the current emc cluster - ememcclustermclabels(vEmcMcParticleIds); + ememcclustermclabels(vEmcMcParticleIds, vAmplitudes); } // end of em emc cluster loop } @@ -625,8 +631,8 @@ struct AssociateMCInfoPhoton { PROCESS_SWITCH(AssociateMCInfoPhoton, processDummy, "processDummy", true); }; -WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +WorkflowSpec defineDataProcessing(ConfigContext const& context) { return WorkflowSpec{ - adaptAnalysisTask(cfgc, TaskName{"associate-mc-info-photon"})}; + adaptAnalysisTask(context, TaskName{"associate-mc-info-photon"})}; } diff --git a/PWGEM/PhotonMeson/TableProducer/skimmerGammaCalo.cxx b/PWGEM/PhotonMeson/TableProducer/skimmerGammaCalo.cxx index 3c54a08a4bb..872df52551c 100644 --- a/PWGEM/PhotonMeson/TableProducer/skimmerGammaCalo.cxx +++ b/PWGEM/PhotonMeson/TableProducer/skimmerGammaCalo.cxx @@ -53,7 +53,7 @@ struct SkimmerGammaCalo { Preslice psMSperCluster = o2::aod::emcalclustercell::emcalclusterId; Produces tableGammaEMCReco; - Produces tableEMCClusterMCLabels; + Produces tableEMCClusterMCLabels; Produces tableCellEMCReco; Produces tableEmEmcClusters; @@ -331,15 +331,16 @@ struct SkimmerGammaCalo { continue; } std::vector mcLabels; + std::vector amplitudes; + mcLabels.reserve(emccluster.amplitudeA().size()); + amplitudes.reserve(emccluster.amplitudeA().size()); for (size_t iCont = 0; iCont < emccluster.amplitudeA().size(); iCont++) { mcLabels.push_back(emccluster.mcParticleIds()[iCont]); + amplitudes.push_back(emccluster.amplitudeA()[iCont]); } - // LOGF(info, "---- New Cluster ---"); - // for (unsigned long int iCont = 0; iCont < mcLabels.size(); iCont++) { - // LOGF(info, "iCont = %d, mcParticle = %d, amplitudeA = %.5f", iCont, mcLabels.at(iCont), emccluster.amplitudeA()[iCont]); - // } - tableEMCClusterMCLabels(mcLabels); + tableEMCClusterMCLabels(mcLabels, amplitudes); mcLabels.clear(); + amplitudes.clear(); } } PROCESS_SWITCH(SkimmerGammaCalo, processMC, "process MC info", false); // Run this in addition to processRec for MCs to copy the cluster mc labels from the EMCALMCClusters to the skimmed EMCClusterMCLabels table diff --git a/PWGEM/PhotonMeson/Tasks/emcalPhotonMcTask.cxx b/PWGEM/PhotonMeson/Tasks/emcalPhotonMcTask.cxx index b83cc325bc7..f1feeb608a8 100644 --- a/PWGEM/PhotonMeson/Tasks/emcalPhotonMcTask.cxx +++ b/PWGEM/PhotonMeson/Tasks/emcalPhotonMcTask.cxx @@ -73,6 +73,9 @@ using namespace o2::aod::pwgem::photon; using namespace o2::constants::physics; using namespace o2::aod::pwgem::photonmeson::utils::mcutil; +constexpr float MinAmpThreshold = 0.2f; // Minimum cluster amplitude threshold to count as significant +constexpr float MaxAmpDiff = 0.1f; // Maximum cluster amplitude difference to leading cluster contribution to count as significant + enum CentralityEstimator { None = 0, CFT0A, @@ -82,21 +85,25 @@ enum CentralityEstimator { }; enum class TruthClass { - Conversion = 0, // 1: true e+/e- pair from the same conversion + Conversion = 0, // 0: true e+/e- pair from the same conversion - PhotonPairSamePi0, // 2: two photon clusters, same Pi0 - PhotonPairDiffPi0, // 3: two photon clusters, different Pi0s - PhotonPairOnePi0, // 4: two photon clusters, only one from a Pi0 + PhotonPairSamePi0, // 1: two photon clusters, same Pi0 + PhotonPairDiffPi0, // 2: two photon clusters, different Pi0s + PhotonPairOnePi0, // 3: two photon clusters, only one from a Pi0 - PhotonElectronSamePi0, // 5: photon + electron cluster, same Pi0 - PhotonElectronDiffPi0, // 6: photon + electron cluster, different Pi0s - PhotonElectronOnePi0, // 7: photon + electron cluster, only one from a Pi0 + PhotonElectronSamePi0, // 4: photon + electron cluster, same Pi0 + PhotonElectronDiffPi0, // 5: photon + electron cluster, different Pi0s + PhotonElectronOnePi0, // 6: photon + electron cluster, only one from a Pi0 + PhotonElectronBS, // 7: photon + electron cluster, from Bremsstrahlung ElectronPairSamePi0, // 8: e+e cluster pair, same Pi0 (conversion and/or Dalitz) ElectronPairDiffPi0, // 9: e+e cluster pair, different Pi0s ElectronPairOnePi0, // 10: e+e cluster pair, only one from a Pi0 - Background, // 11: else / uncorrelated + SplitPhotonCluster, // 11: one photon producing two clusters + SplitLeptonCluster, // 12: one lepton producing two clusters + + Background, // 13: else / uncorrelated NClasses }; @@ -122,6 +129,7 @@ struct ClusterMcInfo { bool isFromConv = false; bool isMergedConv = false; bool isFromPi0 = false; + bool isFromBremsstrahlung = false; int convMotherId = -1; int photonId = -1; }; @@ -131,6 +139,8 @@ ClusterMcInfo classifyCluster(const TGroup& g, TIter& mcCluster, TIter& mcCluste { ClusterMcInfo info; mcCluster.setCursor(g.emmcparticleIds()[0]); + info.isFromBremsstrahlung = isFromBremsstrahlung(mcCluster, mcClusterLooper); + float leadingAmplitude = g.amplitude()[0]; if (std::abs(mcCluster.pdgCode()) == PDG_t::kElectron) { info.isLepton = true; info.convMotherId = getMotherIndexFromChain(mcCluster, mcClusterLooper, PDG_t::kGamma); @@ -139,12 +149,12 @@ ClusterMcInfo classifyCluster(const TGroup& g, TIter& mcCluster, TIter& mcCluste if (mcCluster.mothersIds().size() > 0 && info.isFromConv) { for (size_t i = 1; i < g.emmcparticleIds().size(); ++i) { mcClusterLooper.setCursor(g.emmcparticleIds()[i]); - if (std::abs(mcClusterLooper.pdgCode()) == PDG_t::kElectron) { + if (std::abs(mcClusterLooper.pdgCode()) == PDG_t::kElectron && mcClusterLooper.pdgCode() == -1 * mcCluster.pdgCode()) { int32_t otherConvMotherId = getMotherIndexFromChain(mcClusterLooper, mcClusterLooper2, PDG_t::kGamma); if (otherConvMotherId == info.convMotherId) { - info.isMergedConv = true; - info.isLepton = false; - info.isPhoton = true; + if (g.amplitude()[i] >= leadingAmplitude - MaxAmpDiff && g.amplitude()[i] > MinAmpThreshold) { + info.isMergedConv = true; + } break; } } @@ -155,14 +165,22 @@ ClusterMcInfo classifyCluster(const TGroup& g, TIter& mcCluster, TIter& mcCluste info.isPhoton = true; } - info.photonId = o2::aod::pwgem::photonmeson::utils::mcutil::FindMotherInChain(mcCluster, mcParticles, std::vector{PDG_t::kPi0, Pdg::kEta}); + info.photonId = o2::aod::pwgem::photonmeson::utils::mcutil::FindMotherInChain(mcCluster, mcParticles, std::vector{PDG_t::kPi0, Pdg::kEta, Pdg::kOmega, Pdg::kEtaPrime}); info.isFromPi0 = info.photonId >= 0; return info; } + struct EmcalPhotonMcTask { static constexpr float EMCALRadius = 440.f; static constexpr float PhiVUndefined = -999.f; + static constexpr float Epsilon = 1.e-6f; + + static constexpr std::array(TruthClass::NClasses)> kTruthClassNames = { + "Conversion", "PhotonPairSamePi0", "PhotonPairDiffPi0", "PhotonPairOnePi0", + "PhotonElectronSamePi0", "PhotonElectronDiffPi0", "PhotonElectronOnePi0", "PhotonElectronBS", + "ElectronPairSamePi0", "ElectronPairDiffPi0", "ElectronPairOnePi0", + "SplitPhotonCluster", "SplitLeptonCluster", "Background"}; Produces convTagCandidates; @@ -171,7 +189,7 @@ struct EmcalPhotonMcTask { Configurable grpmagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"}; Configurable skipGRPOquery{"skipGRPOquery", true, "skip grpo query"}; Configurable writeTable{"writeTable", true, "write table for ML."}; - Configurable> classPrescale{"classPrescale", {1, 1, 700, 25, 1, 350, 15, 1, 35, 2, 1000}, "prescale factor per TruthClass, indexed 0..10 matching the enum order"}; + Configurable> classPrescale{"classPrescale", {1, 1, 700, 25, 1, 350, 15, 1, 1, 35, 2, 1, 1, 1000}, "prescale factor per TruthClass, indexed 0..10 matching the enum order"}; Configurable bkgPrescaleSeed{"bkgPrescaleSeed", 42, "seed for the background-prescale RNG"}; // configurable axis @@ -181,7 +199,7 @@ struct EmcalPhotonMcTask { ConfigurableAxis thnConfigAxisMult{"thnConfigAxisMult", {60, 0., 60000.}, "multiplicity axis for the current event"}; ConfigurableAxis thnConfigAxisDeltaEta{"thnConfigAxisDeltaEta", {100, -1, 1}, "delta eta axis"}; ConfigurableAxis thnConfigAxisDeltaPhi{"thnConfigAxisDeltaPhi", {100, -1, 1}, "delta phi axis"}; - Configurable useCent{"useCent", 0, "flag to enable usage of centrality instead of multiplicity as axis."}; + Configurable useCent{"useCent", false, "flag to enable usage of centrality instead of multiplicity as axis."}; struct : ConfigurableGroup { std::string prefix = "conversiontagging"; @@ -241,14 +259,14 @@ struct EmcalPhotonMcTask { SliceCache cache; - using EMCalPhotons = soa::Join; + using EMCalPhotons = soa::Join; using Colls = soa::Join; using McColls = o2::soa::Join; using McParticles = EMMCParticles; - PresliceOptional perCollisionEMC = o2::aod::emccluster::pmeventId; + PresliceOptional perCollisionEMC = o2::aod::emccluster::pmeventId; PresliceOptional perEMCClusterMT = o2::aod::mintm::minClusterId; PresliceOptional perEMCClusterMS = o2::aod::mintm::minClusterId; @@ -346,9 +364,12 @@ struct EmcalPhotonMcTask { hTruthLabel->GetXaxis()->SetBinLabel(static_cast(TruthClass::PhotonElectronSamePi0) + 1, "PhotonElectronSamePi0"); hTruthLabel->GetXaxis()->SetBinLabel(static_cast(TruthClass::PhotonElectronDiffPi0) + 1, "PhotonElectronDiffPi0"); hTruthLabel->GetXaxis()->SetBinLabel(static_cast(TruthClass::PhotonElectronOnePi0) + 1, "PhotonElectronOnePi0"); + hTruthLabel->GetXaxis()->SetBinLabel(static_cast(TruthClass::PhotonElectronBS) + 1, "PhotonElectronBS"); hTruthLabel->GetXaxis()->SetBinLabel(static_cast(TruthClass::ElectronPairSamePi0) + 1, "ElectronPairSamePi0"); hTruthLabel->GetXaxis()->SetBinLabel(static_cast(TruthClass::ElectronPairDiffPi0) + 1, "ElectronPairDiffPi0"); hTruthLabel->GetXaxis()->SetBinLabel(static_cast(TruthClass::ElectronPairOnePi0) + 1, "ElectronPairOnePi0"); + hTruthLabel->GetXaxis()->SetBinLabel(static_cast(TruthClass::SplitPhotonCluster) + 1, "SplitPhotonCluster"); + hTruthLabel->GetXaxis()->SetBinLabel(static_cast(TruthClass::SplitLeptonCluster) + 1, "SplitLeptonCluster"); hTruthLabel->GetXaxis()->SetBinLabel(static_cast(TruthClass::Background) + 1, "Background"); auto hPi0BothResolvedLost = registry.add("EMCal/hPi0BothResolvedLost", "Confusion matrix for conversion tagging", HistType::kTH1D, {{2, -0.5, 1.5}}); @@ -366,6 +387,16 @@ struct EmcalPhotonMcTask { hConfusionMatrixConversionTagging->GetYaxis()->SetBinLabel(5, "#gamma"); mRandGen.seed(bkgPrescaleSeed.value); + + if (classPrescale.value.size() != kTruthClassNames.size()) { + LOG(fatal) << "classPrescale has " << classPrescale.value.size() << " entries, " + << "but TruthClass has " << kTruthClassNames.size() << " members -- update classPrescale!"; + } + + LOG(info) << "=== classPrescale configuration ==="; + for (size_t i = 0; i < kTruthClassNames.size(); ++i) { + LOG(info) << " [" << i << "] " << kTruthClassNames[i] << " -> prescale = " << classPrescale.value[i]; + } }; // end init template @@ -522,7 +553,7 @@ struct EmcalPhotonMcTask { ROOT::Math::XYZVector nRef = u.Cross(zAxis); // normal to the plane containing u and z float phiV = PhiVUndefined; // sentinel for degenerate geometry - if (nDecay.R() > 1e-6f && nRef.R() > 1e-6f) { + if (nDecay.R() > Epsilon && nRef.R() > Epsilon) { float cosPhiV = static_cast(nDecay.Unit().Dot(nRef.Unit())); cosPhiV = std::clamp(cosPhiV, -1.f, 1.f); phiV = std::acos(cosPhiV); @@ -547,28 +578,50 @@ struct EmcalPhotonMcTask { mcCluster1.setCursor(g1.emmcparticleIds()[0]); mcCluster2.setCursor(g2.emmcparticleIds()[0]); - bool areFromSamePi0 = false; // both clusters are from same pi0 or eta - bool areConversionLegs = false; // both clusters are e+ + e- from one photon + bool areFromSamePi0 = false; + bool areConversionLegs = false; + bool areSplitPhotonCluster = false; + bool areSplitLeptonCluster = false; + bool arePhotonElectronBS = false; auto c1 = classifyCluster(g1, mcCluster1, mcClusterLooper, mcClusterLooper2, mcParticles); auto c2 = classifyCluster(g2, mcCluster2, mcClusterLooper, mcClusterLooper2, mcParticles); - if (c1.isFromConv && c2.isFromConv && c1.convMotherId == c2.convMotherId) { + // split-cluster check MUST run first and take priority over everything else -- + // if both clusters share the same dominant MC particle, this is one physical + // shower reconstructed as two clusters, not a genuine pair of anything. + const bool isSameDominantParticle = (g1.emmcparticleIds()[0] == g2.emmcparticleIds()[0]); + if (isSameDominantParticle) { + if (c1.isPhoton) { + areSplitPhotonCluster = true; + } else if (c1.isLepton) { + areSplitLeptonCluster = true; + } + } + + // if they are not a split cluster check for proper conversion pair + if (!isSameDominantParticle && c1.isFromConv && c2.isFromConv && c1.convMotherId == c2.convMotherId) { emcFlagsFromTrueConversion.set(g1.globalIndex()); emcFlagsFromTrueConversion.set(g2.globalIndex()); areConversionLegs = true; } - if (c1.isFromPi0 && c2.isFromPi0) { - // set the cursors to the photonId + // if they are not a split cluster check for neutral meson connection + if (!isSameDominantParticle && c1.isFromPi0 && c2.isFromPi0) { mcPhoton1.setCursor(c1.photonId); mcPhoton2.setCursor(c2.photonId); mcMother.setCursor(mcPhoton1.mothersIds()[0]); if (mcMother.producedByGenerator()) { if (c1.photonId == c2.photonId) { - areFromSamePi0 = true; - emcFlagsFromTrueMesonSameGamma.set(g1.globalIndex()); - emcFlagsFromTrueMesonSameGamma.set(g2.globalIndex()); + // bremsstrahlung: one side is a photon born from the other side's lepton lineage + const bool photonIsBS = (c1.isPhoton && c1.isFromBremsstrahlung) || (c2.isPhoton && c2.isFromBremsstrahlung); + if (photonIsBS && ((c1.isLepton && c2.isPhoton) || (c2.isLepton && c1.isPhoton))) { + arePhotonElectronBS = true; + } else { + areFromSamePi0 = true; + emcFlagsFromTrueMesonSameGamma.set(g1.globalIndex()); + emcFlagsFromTrueMesonSameGamma.set(g2.globalIndex()); + } } else if (mcPhoton1.mothersIds()[0] == mcPhoton2.mothersIds()[0]) { areFromSamePi0 = true; emcFlagsFromTrueMeson.set(g1.globalIndex()); @@ -576,34 +629,39 @@ struct EmcalPhotonMcTask { } } } + bTruthLabel = static_cast(TruthClass::Background); - if (areConversionLegs) { + if (areSplitPhotonCluster) { + bTruthLabel = static_cast(TruthClass::SplitPhotonCluster); + } else if (areSplitLeptonCluster) { + bTruthLabel = static_cast(TruthClass::SplitLeptonCluster); + } else if (areConversionLegs) { bTruthLabel = static_cast(TruthClass::Conversion); - } else { - if (areFromSamePi0) { - if ((c1.isLepton && c2.isPhoton) || (c2.isLepton && c1.isPhoton)) { - bTruthLabel = static_cast(TruthClass::PhotonElectronSamePi0); - } else if (c1.isPhoton && c2.isPhoton) { - bTruthLabel = static_cast(TruthClass::PhotonPairSamePi0); - } else if (c1.isLepton && c2.isLepton) { - bTruthLabel = static_cast(TruthClass::ElectronPairSamePi0); - } - } else if (c1.isFromPi0 && c2.isFromPi0) { - if ((c1.isLepton && c2.isPhoton) || (c2.isLepton && c1.isPhoton)) { - bTruthLabel = static_cast(TruthClass::PhotonElectronDiffPi0); - } else if (c1.isPhoton && c2.isPhoton) { - bTruthLabel = static_cast(TruthClass::PhotonPairDiffPi0); - } else if (c1.isLepton && c2.isLepton) { - bTruthLabel = static_cast(TruthClass::ElectronPairDiffPi0); - } - } else if ((c1.isFromPi0 && !c2.isFromPi0) || (!c1.isFromPi0 && c2.isFromPi0)) { - if ((c1.isLepton && c2.isPhoton) || (c2.isLepton && c1.isPhoton)) { - bTruthLabel = static_cast(TruthClass::PhotonElectronOnePi0); - } else if (c1.isPhoton && c2.isPhoton) { - bTruthLabel = static_cast(TruthClass::PhotonPairOnePi0); - } else if (c1.isLepton && c2.isLepton) { - bTruthLabel = static_cast(TruthClass::ElectronPairOnePi0); - } + } else if (arePhotonElectronBS) { + bTruthLabel = static_cast(TruthClass::PhotonElectronBS); + } else if (areFromSamePi0) { + if ((c1.isLepton && c2.isPhoton) || (c2.isLepton && c1.isPhoton)) { + bTruthLabel = static_cast(TruthClass::PhotonElectronSamePi0); + } else if (c1.isPhoton && c2.isPhoton) { + bTruthLabel = static_cast(TruthClass::PhotonPairSamePi0); + } else if (c1.isLepton && c2.isLepton) { + bTruthLabel = static_cast(TruthClass::ElectronPairSamePi0); + } + } else if (c1.isFromPi0 && c2.isFromPi0) { + if ((c1.isLepton && c2.isPhoton) || (c2.isLepton && c1.isPhoton)) { + bTruthLabel = static_cast(TruthClass::PhotonElectronDiffPi0); + } else if (c1.isPhoton && c2.isPhoton) { + bTruthLabel = static_cast(TruthClass::PhotonPairDiffPi0); + } else if (c1.isLepton && c2.isLepton) { + bTruthLabel = static_cast(TruthClass::ElectronPairDiffPi0); + } + } else if ((c1.isFromPi0 && !c2.isFromPi0) || (!c1.isFromPi0 && c2.isFromPi0)) { + if ((c1.isLepton && c2.isPhoton) || (c2.isLepton && c1.isPhoton)) { + bTruthLabel = static_cast(TruthClass::PhotonElectronOnePi0); + } else if (c1.isPhoton && c2.isPhoton) { + bTruthLabel = static_cast(TruthClass::PhotonPairOnePi0); + } else if (c1.isLepton && c2.isLepton) { + bTruthLabel = static_cast(TruthClass::ElectronPairOnePi0); } } @@ -634,7 +692,7 @@ struct EmcalPhotonMcTask { } mcCluster1.setCursor(cluster.emmcparticleIds()[0]); - int photonid1 = o2::aod::pwgem::photonmeson::utils::mcutil::FindMotherInChain(mcCluster1, mcParticles, std::vector{PDG_t::kPi0, Pdg::kEta}); + int photonid1 = o2::aod::pwgem::photonmeson::utils::mcutil::FindMotherInChain(mcCluster1, mcParticles, std::vector{PDG_t::kPi0, Pdg::kEta, Pdg::kOmega, Pdg::kEtaPrime}); int motherId = -1; if (photonid1 >= 0) { mcPhoton1.setCursor(photonid1); diff --git a/PWGEM/PhotonMeson/Tasks/taskFlowReso.cxx b/PWGEM/PhotonMeson/Tasks/taskFlowReso.cxx index 0918ba11d8b..65b3ab804df 100644 --- a/PWGEM/PhotonMeson/Tasks/taskFlowReso.cxx +++ b/PWGEM/PhotonMeson/Tasks/taskFlowReso.cxx @@ -51,7 +51,7 @@ using namespace o2::aod::pwgem::photon; enum QvecEstimator { FT0M = 0, - FT0A = 1, + FT0A, FT0C, TPCPos, TPCNeg, @@ -61,7 +61,7 @@ enum QvecEstimator { enum CentralityEstimator { None = 0, - CFT0A = 1, + CFT0A, CFT0C, CFT0M, NCentralityEstimators @@ -383,7 +383,7 @@ struct TaskFlowReso { // no selection on the centrality is applied on purpose to allow for the resolution study in post-processing return; } - if (!(eventcuts.cfgFT0COccupancyMin <= collision.ft0cOccupancyInTimeRange() && collision.ft0cOccupancyInTimeRange() < eventcuts.cfgFT0COccupancyMax)) { + if (!(eventcuts.cfgFT0COccupancyMin <= collision.ft0cOccupancyInTimeRange()) || !(collision.ft0cOccupancyInTimeRange() < eventcuts.cfgFT0COccupancyMax)) { return; } float cent = getCentrality(collision); @@ -531,7 +531,7 @@ struct TaskFlowReso { }; // End struct TaskFlowReso -WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +WorkflowSpec defineDataProcessing(ConfigContext const& context) { - return WorkflowSpec{adaptAnalysisTask(cfgc)}; + return WorkflowSpec{adaptAnalysisTask(context)}; } diff --git a/PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx b/PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx index c8a097a79cf..af98c86de85 100644 --- a/PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx +++ b/PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx @@ -173,6 +173,7 @@ struct TaskPi0FlowEMC { Configurable cfgEMCUseTM{"cfgEMCUseTM", false, "flag to use EMCal track matching cut or not"}; Configurable emcUseSecondaryTM{"emcUseSecondaryTM", false, "flag to use EMCal secondary track matching cut or not"}; Configurable cfgEnableQA{"cfgEnableQA", false, "flag to turn QA plots on/off"}; + Configurable separateEMCalDCal{"separateEMCalDCal", false, "flag to only pair EMCal with EMCal and DCal with DCal clusters"}; } emccuts; V0PhotonCut fV0PhotonCut; @@ -252,6 +253,8 @@ struct TaskPi0FlowEMC { int runNow = 0; int runBefore = -1; + static constexpr float MaxPhiEMCal = 3.5f; + // Filter clusterFilter = aod::skimmedcluster::time >= emccuts.cfgEMCminTime && aod::skimmedcluster::time <= emccuts.cfgEMCmaxTime && aod::skimmedcluster::m02 >= emccuts.cfgEMCminM02 && aod::skimmedcluster::m02 <= emccuts.cfgEMCmaxM02 && aod::skimmedcluster::e >= emccuts.cfgEMCminE; Filter collisionFilter = (nabs(aod::collision::posZ) <= eventcuts.cfgZvtxMax) && (aod::evsel::ft0cOccupancyInTimeRange <= eventcuts.cfgFT0COccupancyMax) && (aod::evsel::ft0cOccupancyInTimeRange >= eventcuts.cfgFT0COccupancyMin); // using FilteredEMCalPhotons = soa::Filtered>; @@ -261,10 +264,13 @@ struct TaskPi0FlowEMC { using CollsWithQvecs = soa::Join; using Colls = soa::Join; + Partition emcalPhotons = aod::mincluster::storedPhi < std::lround(MaxPhiEMCal * emcdownscaling::downscalingFactors[emcdownscaling::kPhi]); + Partition dcalPhotons = aod::mincluster::storedPhi >= std::lround(MaxPhiEMCal * emcdownscaling::downscalingFactors[emcdownscaling::kPhi]); + static constexpr std::size_t NQVecEntries = 6; - PresliceOptional perCollisionEMC = o2::aod::emccluster::pmeventId; - PresliceOptional perCollisionPCM = aod::v0photonkf::pmeventId; + PresliceOptional perCollisionEMC = o2::aod::emccluster::pmeventId; + PresliceOptional perCollisionPCM = aod::v0photonkf::pmeventId; PresliceOptional perEMCClusterMT = o2::aod::mintm::minClusterId; PresliceOptional perEMCClusterMS = o2::aod::mintm::minClusterId; @@ -482,6 +488,13 @@ struct TaskPi0FlowEMC { }; // end init + /// \brief Check whether a photon (by its phi) falls in the EMCal or DCal acceptance + /// \param phi azimuthal angle of the photon + static bool isEMCalRegion(float phi) + { + return phi < MaxPhiEMCal; + } + /// Change radians to degree /// \param angle in radians /// \return angle in degree @@ -1048,7 +1061,7 @@ struct TaskPi0FlowEMC { continue; } if (rotationConfig.cfgDoRotation.value && nColl % rotationConfig.cfgDownsampling == 0) { - rotationBackground(vMeson, v1, v2, photons1, g1.globalIndex(), g2.globalIndex(), collision); + rotationBackground(vMeson, v1, v2, photons1, g1.globalIndex(), g2.globalIndex(), collision); } if (thnConfigAxisInvMass.value[1] > vMeson.M() || thnConfigAxisInvMass.value.back() < vMeson.M()) { registry.fill(HIST("hMesonCuts"), 3); @@ -1094,6 +1107,8 @@ struct TaskPi0FlowEMC { } auto photonsPerCollision = clusters.sliceBy(perCollisionEMC, collision.globalIndex()); + auto emcalPhotonsPerCollision = emcalPhotons.sliceBy(perCollisionEMC, collision.globalIndex()); + auto dcalPhotonsPerCollision = dcalPhotons.sliceBy(perCollisionEMC, collision.globalIndex()); if (emccuts.cfgEnableQA.value) { for (const auto& photon : photonsPerCollision) { @@ -1109,7 +1124,12 @@ struct TaskPi0FlowEMC { registry.fill(HIST("clusterQA/hClusterEtaPhiAfter"), photon.phi(), photon.eta()); // after cuts } } - runPairingLoop(collision, photonsPerCollision, photonsPerCollision, flags, flags); + if (emccuts.separateEMCalDCal.value) { + runPairingLoop(collision, emcalPhotonsPerCollision, emcalPhotonsPerCollision, flags, flags); + runPairingLoop(collision, dcalPhotonsPerCollision, dcalPhotonsPerCollision, flags, flags); + } else { + runPairingLoop(collision, photonsPerCollision, photonsPerCollision, flags, flags); + } if (rotationConfig.cfgDoRotation.value) { if (nColl % rotationConfig.cfgDownsampling == 0) { nColl = 1; // reset counter @@ -1164,6 +1184,10 @@ struct TaskPi0FlowEMC { if (!(flags.test(g1.globalIndex())) || !(flags.test(g2.globalIndex()))) { continue; } + if (emccuts.separateEMCalDCal.value && isEMCalRegion(g1.phi()) != isEMCalRegion(g2.phi())) { + continue; // only pair EMCal-EMCal or DCal-DCal + } + // Cut edge clusters away, similar to rotation method to ensure same acceptance is used if (cfgDistanceToEdge.value > 0) { if (checkEtaPhi1D(g1.eta(), RecoDecay::constrainAngle(g1.phi())) >= cfgEMCalMapLevelBackground.value) { diff --git a/PWGEM/PhotonMeson/Tasks/testTaskEmc.cxx b/PWGEM/PhotonMeson/Tasks/testTaskEmc.cxx index c10b30fc3a6..ad2b5594c15 100644 --- a/PWGEM/PhotonMeson/Tasks/testTaskEmc.cxx +++ b/PWGEM/PhotonMeson/Tasks/testTaskEmc.cxx @@ -172,7 +172,7 @@ struct TestTaskEmc { PROCESS_SWITCH(TestTaskEmc, processEMCalCalib, "Process EMCal calibration same event", true); }; // End struct TestTaskEmc -WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +WorkflowSpec defineDataProcessing(ConfigContext const& context) { - return WorkflowSpec{adaptAnalysisTask(cfgc)}; + return WorkflowSpec{adaptAnalysisTask(context)}; } diff --git a/PWGEM/PhotonMeson/Utils/MCUtilities.h b/PWGEM/PhotonMeson/Utils/MCUtilities.h index 94fceb68445..60557218784 100644 --- a/PWGEM/PhotonMeson/Utils/MCUtilities.h +++ b/PWGEM/PhotonMeson/Utils/MCUtilities.h @@ -315,6 +315,24 @@ bool isMotherPDG(const T& mcparticle, T& mcparticleWorking, const int motherPDG, return isMotherPDG(mcparticleWorking, mcparticleWorking, motherPDG, depth - 1); } +//_______________________________________________________________________ +/// \brief Check if given particle is from Bremsstrahlung +/// \param mcCursor iterator of mcparticle +/// \param iter shared iterator used to walk to the mother +template +bool isFromBremsstrahlung(TIter const& mcCursor, TIter& iter) +{ + if (!mcCursor.has_mothers()) { + return false; + } + if (mcCursor.pdgCode() != PDG_t::kGamma) { + return false; // only a photon can itself be a bremsstrahlung emission + } + const int motherId = mcCursor.mothersIds()[0]; + iter.setCursor(motherId); + return std::abs(iter.pdgCode()) == PDG_t::kElectron; +} + //_______________________________________________________________________ /// \brief Go up the decay chain of a mcparticle looking for a mother with the given pdg codes, if found return id else -1 /// E.g. if electron cluster is coming from a photon return true, if primary electron return false From f09b919486984c3eb1e2a3943378d70e6344b90c Mon Sep 17 00:00:00 2001 From: Marvin Hemmer Date: Fri, 31 Jul 2026 12:57:49 +0200 Subject: [PATCH 2/2] Fix Missing include --- PWGEM/PhotonMeson/TableProducer/associateMCinfoPhoton.cxx | 1 + 1 file changed, 1 insertion(+) diff --git a/PWGEM/PhotonMeson/TableProducer/associateMCinfoPhoton.cxx b/PWGEM/PhotonMeson/TableProducer/associateMCinfoPhoton.cxx index 30e3cb4cd3c..7eb9d556d85 100644 --- a/PWGEM/PhotonMeson/TableProducer/associateMCinfoPhoton.cxx +++ b/PWGEM/PhotonMeson/TableProducer/associateMCinfoPhoton.cxx @@ -38,6 +38,7 @@ #include #include +#include #include #include #include