From f37b741aba8122a4e4e7dc09ddaa8917f539f074 Mon Sep 17 00:00:00 2001 From: Marvin Hemmer Date: Mon, 10 Aug 2026 10:06:35 +0200 Subject: [PATCH] [PWGEM,Photon] Adding EMCal ML response plus initCCDB bugfix - Add EMCal ML response which can be used to identify conversion cluster pairs - Fix a bug in the initCCDB function in multiple files where the member variable that stores the run number was not updated triggering the update of the CCDB for every collision including logging - Update `TruthClass` in `emcalPhotonMcTask.cxx` - Update `MCUtilities.h` the get Origin functions to use `TMCProcess` enums from ROOT - Add `GetMesonInChain` function which uses `FindMotherInChain` but instead of returning the ID of the daughter of the meson it returns the meson ID - Fix clang-tidy warnings and errors --- .../PhotonMeson/Core/EMCConversionCandidate.h | 63 +++ PWGEM/PhotonMeson/Core/EMCPhotonCut.h | 4 + .../Core/EmMlResponseEMCConversion.h | 102 ++++ PWGEM/PhotonMeson/Core/Pi0EtaToGammaGamma.h | 5 +- PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h | 5 +- PWGEM/PhotonMeson/DataModel/gammaTables.h | 4 +- .../TableProducer/photonconversionbuilder.cxx | 4 +- PWGEM/PhotonMeson/Tasks/MaterialBudgetMC.cxx | 2 +- PWGEM/PhotonMeson/Tasks/SinglePhotonMC.cxx | 2 - PWGEM/PhotonMeson/Tasks/TagAndProbe.cxx | 2 +- PWGEM/PhotonMeson/Tasks/emcalPhotonMcTask.cxx | 474 +++++++++++++----- PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx | 23 +- PWGEM/PhotonMeson/Utils/MCUtilities.h | 162 ++++-- PWGEM/PhotonMeson/Utils/ParticleOrigin.h | 46 ++ 14 files changed, 711 insertions(+), 187 deletions(-) create mode 100644 PWGEM/PhotonMeson/Core/EMCConversionCandidate.h create mode 100644 PWGEM/PhotonMeson/Core/EmMlResponseEMCConversion.h create mode 100644 PWGEM/PhotonMeson/Utils/ParticleOrigin.h diff --git a/PWGEM/PhotonMeson/Core/EMCConversionCandidate.h b/PWGEM/PhotonMeson/Core/EMCConversionCandidate.h new file mode 100644 index 00000000000..07eac2a030a --- /dev/null +++ b/PWGEM/PhotonMeson/Core/EMCConversionCandidate.h @@ -0,0 +1,63 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file EMCConversionCandidate.h +/// \brief Header file that defines EMCConversionCandidate, a struct to be used with the EMCConversion ML model +/// \author Marvin Hemmer + +#ifndef PWGEM_PHOTONMESON_CORE_EMCCONVERSIONCANDIDATE_H_ +#define PWGEM_PHOTONMESON_CORE_EMCCONVERSIONCANDIDATE_H_ + +#include + +namespace o2::analysis::em +{ + +// Requires T to expose the specific named getters getInputFeatures() needs, +// each returning something convertible to float. +template +concept IsEmcConversionCandidate = requires(T const& c) { + { c.minv() } -> std::convertible_to; + { c.deltaEta() } -> std::convertible_to; + { c.deltaR() } -> std::convertible_to; + { c.phiv() } -> std::convertible_to; + { c.rConv() } -> std::convertible_to; + { c.totE() } -> std::convertible_to; + { c.e2() } -> std::convertible_to; + { c.e1() } -> std::convertible_to; + { c.deltaPhi() } -> std::convertible_to; +}; + +struct EMCConversionCandidate { + float mMinv; + float mDeltaEta; + float mDeltaR; + float mPhiv; + float mRConv; + float mTotE; + float mE2; + float mE1; + float mDeltaPhi; + + [[nodiscard]] float minv() const { return mMinv; } + [[nodiscard]] float deltaEta() const { return mDeltaEta; } + [[nodiscard]] float deltaR() const { return mDeltaR; } + [[nodiscard]] float phiv() const { return mPhiv; } + [[nodiscard]] float rConv() const { return mRConv; } + [[nodiscard]] float totE() const { return mTotE; } + [[nodiscard]] float e2() const { return mE2; } + [[nodiscard]] float e1() const { return mE1; } + [[nodiscard]] float deltaPhi() const { return mDeltaPhi; } +}; + +} // namespace o2::analysis::em + +#endif // PWGEM_PHOTONMESON_CORE_EMCCONVERSIONCANDIDATE_H_ diff --git a/PWGEM/PhotonMeson/Core/EMCPhotonCut.h b/PWGEM/PhotonMeson/Core/EMCPhotonCut.h index 5567ca37c5f..1e93046b573 100644 --- a/PWGEM/PhotonMeson/Core/EMCPhotonCut.h +++ b/PWGEM/PhotonMeson/Core/EMCPhotonCut.h @@ -443,6 +443,10 @@ class EMCPhotonCut /// \return true if cluster survives cut else false bool IsSelectedEMCalRunning(const EMCPhotonCuts& cut, o2::soa::is_iterator auto const& cluster, IsFullTrackIterator auto& matchedTrackIter, int64_t const nMatchedTracks, o2::framework::HistogramRegistry* fRegistry = nullptr) const { + if (nMatchedTracks == 0) { + // there are not tracks to match with, so its true + return true; + } switch (cut) { case EMCPhotonCuts::kTM: return checkTrackMatching(cluster, matchedTrackIter, nMatchedTracks, true, [this](float pt) { return GetTrackMatchingEta(pt); }, [this](float pt) { return GetTrackMatchingPhi(pt); }, fRegistry, TrackType::kPrimary); diff --git a/PWGEM/PhotonMeson/Core/EmMlResponseEMCConversion.h b/PWGEM/PhotonMeson/Core/EmMlResponseEMCConversion.h new file mode 100644 index 00000000000..c60e3343264 --- /dev/null +++ b/PWGEM/PhotonMeson/Core/EmMlResponseEMCConversion.h @@ -0,0 +1,102 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file EmMlResponseEMCConversion.h +/// \brief Class to compute the ML response for EMC conversion selections +/// \author Marvin Hemmer + +#ifndef PWGEM_PHOTONMESON_CORE_EMMLRESPONSEEMCCONVERSION_H_ +#define PWGEM_PHOTONMESON_CORE_EMMLRESPONSEEMCCONVERSION_H_ + +#include "PWGEM/PhotonMeson/Core/EMCConversionCandidate.h" + +#include "Tools/ML/MlResponse.h" + +#include +#include + +// Fill the map of available input features +// the key is the feature's name (std::string) +// the value is the corresponding value in EnumInputFeatures +#define FILL_MAP_EMC_CONV(FEATURE) \ + { \ + #FEATURE, static_cast(InputFeaturesEMCConversion::FEATURE)} + +// Check if the index of mCachedIndices (index associated to a FEATURE) +// matches the entry in EnumInputFeatures associated to this FEATURE +// if so, the inputFeatures vector is filled with the FEATURE's value +// by calling the corresponding GETTER from OBJECT +// NOLINTNEXTLINE(cppcoreguidelines-macro-usage) +#define CHECK_AND_FILL_VEC_EMC_CONV(GETTER) \ + case static_cast(InputFeaturesEMCConversion::GETTER): { \ + inputFeatures.emplace_back(candidate.GETTER()); \ + break; \ + } + +namespace o2::analysis::em::emcconv +{ + +// input feature used in the ml model +enum class InputFeaturesEMCConversion : uint8_t { + minv, + deltaEta, + deltaR, + phiv, + rConv, + totE, + e2, + e1, + deltaPhi +}; + +template +class EmMlResponseEMCConversion : public MlResponse +{ + public: + EmMlResponseEMCConversion() = default; + virtual ~EmMlResponseEMCConversion() = default; + + template + std::vector getInputFeatures(TCandidate const& candidate) + { + std::vector inputFeatures; + for (const auto& idx : MlResponse::mCachedIndices) { + switch (idx) { + CHECK_AND_FILL_VEC_EMC_CONV(minv) + CHECK_AND_FILL_VEC_EMC_CONV(deltaEta) + CHECK_AND_FILL_VEC_EMC_CONV(deltaR) + CHECK_AND_FILL_VEC_EMC_CONV(phiv) + CHECK_AND_FILL_VEC_EMC_CONV(rConv) + CHECK_AND_FILL_VEC_EMC_CONV(totE) + CHECK_AND_FILL_VEC_EMC_CONV(e2) + CHECK_AND_FILL_VEC_EMC_CONV(e1) + CHECK_AND_FILL_VEC_EMC_CONV(deltaPhi) + } + } + return inputFeatures; + } + + protected: + void setAvailableInputFeatures() + { + MlResponse::mAvailableInputFeatures = { + FILL_MAP_EMC_CONV(minv), FILL_MAP_EMC_CONV(deltaEta), FILL_MAP_EMC_CONV(deltaR), + FILL_MAP_EMC_CONV(phiv), FILL_MAP_EMC_CONV(rConv), FILL_MAP_EMC_CONV(totE), + FILL_MAP_EMC_CONV(e2), FILL_MAP_EMC_CONV(e1), FILL_MAP_EMC_CONV(deltaPhi)}; + } +}; + +} // namespace o2::analysis::em::emcconv + +#undef FILL_MAP_EMC_CONV +#undef CHECK_AND_FILL_VEC_EMC_CONV + +#endif // PWGEM_PHOTONMESON_CORE_EMMLRESPONSEEMCCONVERSION_H_ diff --git a/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGamma.h b/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGamma.h index e1db2a28d7b..45e4f7c26d7 100644 --- a/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGamma.h +++ b/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGamma.h @@ -458,10 +458,8 @@ struct Pi0EtaToGammaGamma { if (mRunNumber == collision.runNumber()) { return; } + mRunNumber = collision.runNumber(); - if (mRunNumber == collision.runNumber()) { - return; - } // In case override, don't proceed, please - no CCDB access required if (d_bz_input > -990) { // o2-linter: disable=magic-number (override value) d_bz = d_bz_input; @@ -470,7 +468,6 @@ struct Pi0EtaToGammaGamma { grpmag.setL3Current(30000.f / (d_bz / 5.0f)); // o2-linter: disable=magic-number (override value) } o2::base::Propagator::initFieldFromGRP(&grpmag); - mRunNumber = collision.runNumber(); return; } diff --git a/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h b/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h index ea43bd54900..79c7b4308ae 100644 --- a/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h +++ b/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h @@ -313,10 +313,8 @@ struct Pi0EtaToGammaGammaMC { if (mRunNumber == collision.runNumber()) { return; } + mRunNumber = collision.runNumber(); - if (mRunNumber == collision.runNumber()) { - return; - } // In case override, don't proceed, please - no CCDB access required if (d_bz_input > -990) { // o2-linter: disable=magic-number (override value) d_bz = d_bz_input; @@ -325,7 +323,6 @@ struct Pi0EtaToGammaGammaMC { grpmag.setL3Current(30000.f / (d_bz / 5.0f)); // o2-linter: disable=magic-number (override value) } o2::base::Propagator::initFieldFromGRP(&grpmag); - mRunNumber = collision.runNumber(); return; } diff --git a/PWGEM/PhotonMeson/DataModel/gammaTables.h b/PWGEM/PhotonMeson/DataModel/gammaTables.h index 33dca54b611..1a5ae7e4a84 100644 --- a/PWGEM/PhotonMeson/DataModel/gammaTables.h +++ b/PWGEM/PhotonMeson/DataModel/gammaTables.h @@ -647,10 +647,10 @@ DECLARE_SOA_INDEX_COLUMN(EmEmcCluster, emEmcCluster); //! } // namespace trackmatching DECLARE_SOA_TABLE(EmEmcMTracks, "AOD", "EMEMCMTRACK", //! - trackmatching::EmEmcClusterId, emctm::DeltaPhi, emctm::DeltaEta, emctm::TrackP, emctm::TrackPt); + o2::soa::Index<>, trackmatching::EmEmcClusterId, emctm::DeltaPhi, emctm::DeltaEta, emctm::TrackP, emctm::TrackPt); DECLARE_SOA_TABLE(EmEmcMSTracks, "AOD", "EMEMCMSTRACK", //! - trackmatching::EmEmcClusterId, emctm::DeltaPhi, emctm::DeltaEta, emctm::TrackP, emctm::TrackPt); + o2::soa::Index<>, trackmatching::EmEmcClusterId, emctm::DeltaPhi, emctm::DeltaEta, emctm::TrackP, emctm::TrackPt); DECLARE_SOA_TABLE(EMCEMEventIds_000, "AOD", "EMCEMEVENTID", emccluster::EMEventId); // To be joined with SkimEMCClusters table at analysis level. DECLARE_SOA_TABLE_VERSIONED(EMCEMEventIds_001, "AOD", "EMCEMEVENTID", 1, emccluster::PMEventId); // To be joined with SkimEMCClusters table at analysis level. diff --git a/PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.cxx b/PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.cxx index b8d3229c6df..01429a688b1 100644 --- a/PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.cxx +++ b/PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.cxx @@ -468,6 +468,7 @@ struct PhotonConversionBuilder { if (mRunNumber == bc.runNumber()) { return; } + mRunNumber = bc.runNumber(); // In case override, don't proceed, please - no CCDB access required if (d_bz_input > -990) { // o2-linter: disable=magic-number (override value) d_bz = d_bz_input; @@ -476,7 +477,7 @@ struct PhotonConversionBuilder { grpmag.setL3Current(30000.f / (d_bz / 5.0f)); // o2-linter: disable=magic-number (override value) } o2::base::Propagator::initFieldFromGRP(&grpmag); - mRunNumber = bc.runNumber(); + return; } @@ -485,7 +486,6 @@ struct PhotonConversionBuilder { // Fetch magnetic field from ccdb for current collision d_bz = bc.grpMagField().getNominalL3Field(); LOG(info) << "Retrieved GRP for timestamp " << run3grp_timestamp << " with magnetic field of " << d_bz << " kZG"; - mRunNumber = bc.runNumber(); if (useMatCorrType == 2) { // o2-linter: disable=magic-number (material budget correction) // setMatLUT only after magfield has been initalized (setMatLUT has implicit and problematic init field call if not) diff --git a/PWGEM/PhotonMeson/Tasks/MaterialBudgetMC.cxx b/PWGEM/PhotonMeson/Tasks/MaterialBudgetMC.cxx index 30bf1e1d469..4db16a9e817 100644 --- a/PWGEM/PhotonMeson/Tasks/MaterialBudgetMC.cxx +++ b/PWGEM/PhotonMeson/Tasks/MaterialBudgetMC.cxx @@ -131,7 +131,7 @@ struct MaterialBudgetMC { auto* list_pair_subsys_photoncut = dynamic_cast(list_pair_subsys->FindObject(photon_cut_name.data())); for (const auto& cut3 : cuts3) { - std::string pair_cut_name = cut3.getName(); + std::string const& pair_cut_name = cut3.getName(); o2::aod::pwgem::photon::histogram::AddHistClass(list_pair_subsys_photoncut, pair_cut_name.data()); auto* list_pair_subsys_paircut = dynamic_cast(list_pair_subsys_photoncut->FindObject(pair_cut_name.data())); o2::aod::pwgem::photon::histogram::DefineHistograms(list_pair_subsys_paircut, "material_budget_study", "Pair"); diff --git a/PWGEM/PhotonMeson/Tasks/SinglePhotonMC.cxx b/PWGEM/PhotonMeson/Tasks/SinglePhotonMC.cxx index 40918d54757..d5619fd4edb 100644 --- a/PWGEM/PhotonMeson/Tasks/SinglePhotonMC.cxx +++ b/PWGEM/PhotonMeson/Tasks/SinglePhotonMC.cxx @@ -17,9 +17,7 @@ #include "PWGEM/Dilepton/Utils/MCUtilities.h" #include "PWGEM/PhotonMeson/Core/CutsLibrary.h" -#include "PWGEM/PhotonMeson/Core/EMCPhotonCut.h" #include "PWGEM/PhotonMeson/Core/HistogramsLibrary.h" -#include "PWGEM/PhotonMeson/Core/PHOSPhotonCut.h" #include "PWGEM/PhotonMeson/Core/V0PhotonCut.h" #include "PWGEM/PhotonMeson/DataModel/EventTables.h" #include "PWGEM/PhotonMeson/DataModel/gammaTables.h" diff --git a/PWGEM/PhotonMeson/Tasks/TagAndProbe.cxx b/PWGEM/PhotonMeson/Tasks/TagAndProbe.cxx index 8de6e48e99e..6bb206de9dc 100644 --- a/PWGEM/PhotonMeson/Tasks/TagAndProbe.cxx +++ b/PWGEM/PhotonMeson/Tasks/TagAndProbe.cxx @@ -140,7 +140,7 @@ struct TagAndProbe { THashList* list_pair_subsys_photoncut = o2::aod::pwgem::photon::histogram::AddHistClass(list_pair, photon_cut_name.data()); for (auto& cut3 : paircuts) { - std::string pair_cut_name = cut3.getName(); + std::string const& pair_cut_name = cut3.getName(); o2::aod::pwgem::photon::histogram::AddHistClass(list_pair_subsys_photoncut, pair_cut_name.data()); auto* list_pair_subsys_paircut = dynamic_cast(list_pair_subsys_photoncut->FindObject(pair_cut_name.data())); o2::aod::pwgem::photon::histogram::DefineHistograms(list_pair_subsys_paircut, "tag_and_probe", pairname.data()); diff --git a/PWGEM/PhotonMeson/Tasks/emcalPhotonMcTask.cxx b/PWGEM/PhotonMeson/Tasks/emcalPhotonMcTask.cxx index 27f17c4d352..164d8a67b06 100644 --- a/PWGEM/PhotonMeson/Tasks/emcalPhotonMcTask.cxx +++ b/PWGEM/PhotonMeson/Tasks/emcalPhotonMcTask.cxx @@ -14,18 +14,22 @@ /// \author M. Hemmer, marvin.hemmer@cern.ch #include "PWGEM/PhotonMeson/Core/EMBitFlags.h" +#include "PWGEM/PhotonMeson/Core/EMCConversionCandidate.h" #include "PWGEM/PhotonMeson/Core/EMCPhotonCut.h" #include "PWGEM/PhotonMeson/Core/EMPhotonEventCut.h" +#include "PWGEM/PhotonMeson/Core/EmMlResponseEMCConversion.h" #include "PWGEM/PhotonMeson/DataModel/ConversionMl.h" #include "PWGEM/PhotonMeson/DataModel/EventTables.h" #include "PWGEM/PhotonMeson/DataModel/GammaTablesRedux.h" #include "PWGEM/PhotonMeson/DataModel/gammaTables.h" #include "PWGEM/PhotonMeson/Utils/EventHistograms.h" #include "PWGEM/PhotonMeson/Utils/MCUtilities.h" +#include "PWGEM/PhotonMeson/Utils/ParticleOrigin.h" #include "Common/Core/RecoDecay.h" +#include "Tools/ML/MlResponse.h" -#include +#include #include #include #include @@ -33,6 +37,7 @@ #include #include #include +#include #include #include #include @@ -49,6 +54,7 @@ #include #include #include +#include #include #include @@ -61,6 +67,7 @@ #include #include #include +#include #include #include @@ -72,10 +79,10 @@ using namespace o2::soa; using namespace o2::aod::pwgem::photon; using namespace o2::constants::physics; using namespace o2::aod::pwgem::photonmeson::utils::mcutil; +using namespace o2::analysis::em; 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 - +constexpr float MinAmpFraction = 0.6f; // Minimum fraction of particle energy that it needs to deposit in cluster to count as significant enum CentralityEstimator { None = 0, CFT0A, @@ -85,29 +92,19 @@ enum CentralityEstimator { }; enum class TruthClass { - Conversion = 0, // 0: true e+/e- pair from the same conversion - - 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, // 4: photon + electron cluster, same Pi0 - PhotonElectronDiffPi0, // 5: photon + electron cluster, different Pi0s - PhotonElectronOnePi0, // 6: photon + electron cluster, only one from a Pi0 - BSPhotonElectron, // 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 - - SplitPhotonCluster, // 11: one photon producing two clusters - SplitLeptonCluster, // 12: one lepton producing two clusters - PhotonBSPhotonPair, // 13: photon + photon from Bremsstrahlung - ElectronBSPhotonPair, // 14: one cluster from Bremsstrahlung and one electron cluster except case BSPhotonElectron - BSPhotonPair, // 15: both photons from Bremsstrahlung - - Background, // 16: else / uncorrelated - + Conversion, // two electron legs, same conversion vertex + GammaGammaSamePi0, // both clusters are DIRECT daughter photons of the same generator-level meson + GammaGammaAnnihilation, // two photons from the same e+e- annihilation vertex + BSPhotonElectron, // a bremsstrahlung photon paired with the specific lepton that radiated it + PhotonComptonElectronPair, // a Compton-scattered photon paired with its own recoil electron + ElectronPairSamePi0, // two lepton clusters, directly from same meson + CrossConvertedSiblings, // two lepton clusters, same meson cross-converted siblings + DalitzDecaySiblings, // two lepton clusters, same meson from Dalitz + SplitPhotonCluster, // one physical photon shower reconstructed as two clusters + SplitLeptonCluster, // one physical lepton shower reconstructed as two clusters + IndirGammaGammaSamePi0, // both clusters are indirect daughter photons of the same generator-level meson + SameMesonIndirect, // both clusters trace to the SAME generator-level meson, but at least one path passes through extra generations (e.g. a further conversion/BS/scatter) before reaching the cluster -- NOT a clean two-body relationship + Background, // no common meson ancestor at all -- genuinely uncorrelated combinatorics NClasses }; @@ -127,52 +124,73 @@ enum class TagDecision { }; struct ClusterMcInfo { - bool isLepton = false; - bool isPhoton = false; - bool isFromConv = false; - bool isMergedConv = false; - bool isFromPi0 = false; - bool isFromBremsstrahlung = false; - int convMotherId = -1; - int photonId = -1; - float purity = 0; + bool isLepton = false; // is cluster from a lepton + bool isPhoton = false; // is cluster from a photon + bool isMergedConv = false; // is cluster from a merged conversion + LeptonOrigin leptonOrigin = LeptonOrigin::Other; // origin of lepton in case isLepton == true + PhotonOrigin photonOrigin = PhotonOrigin::Other; // origin of photon in case isPhoton == true + int photonMotherId = -1; // for leptons: the photon this lepton traces to + int mesonId = -1; // for Decay photons: the pi0/eta/omega/etaprime id | for leptons from decay photons: the pi0/eta/omega/etaprime id + int hardPartonId = -1; // for Direct photons: the quark/gluon id + int decayPhotonId = -1; // unified: the photon (self, if isPhoton; or photonMotherId, if isLepton) whose immediate mother is a meson. -1 if not applicable. + float purity = 0.f; // fraction of energy the main particle contributed to the cluster + float radius = 0.f; // radius in xy from where the main contributor to the cluster originated }; -template -ClusterMcInfo classifyCluster(const TGroup& g, TIter& mcCluster, TIter& mcClusterLooper, TIter& mcClusterLooper2, McParticles const& mcParticles) +template +ClusterMcInfo classifyCluster(const TCluster& g, TIter& mcCluster, TIter& mcClusterLooper, TIter& mcClusterLooper2, McParticles const& mcParticles) { + static const std::array kMesonPdgs{PDG_t::kPi0, Pdg::kEta, Pdg::kOmega, Pdg::kEtaPrime}; ClusterMcInfo info; mcCluster.setCursor(g.emmcparticleIds()[0]); - info.isFromBremsstrahlung = isFromBremsstrahlung(mcCluster, mcClusterLooper); // particle has to be a photon and it has to have a e+ or e- as mother! + info.radius = std::hypot(mcCluster.vx(), mcCluster.vy()); float leadingAmplitude = g.amplitude()[0]; info.purity = leadingAmplitude; + info.mesonId = o2::aod::pwgem::photonmeson::utils::mcutil::GetMesonInChain(mcCluster, mcParticles, kMesonPdgs); if (std::abs(mcCluster.pdgCode()) == PDG_t::kElectron) { info.isLepton = true; - info.convMotherId = getMotherIndexFromChain(mcCluster, mcClusterLooper, PDG_t::kGamma); - info.isFromConv = info.convMotherId >= 0; - - 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 && mcClusterLooper.pdgCode() == -1 * mcCluster.pdgCode()) { - int32_t otherConvMotherId = getMotherIndexFromChain(mcClusterLooper, mcClusterLooper2, PDG_t::kGamma); - if (otherConvMotherId == info.convMotherId) { - if (g.amplitude()[i] >= leadingAmplitude - MaxAmpDiff && g.amplitude()[i] > MinAmpThreshold) { - info.isMergedConv = true; - } - break; - } - } + info.leptonOrigin = getLeptonOriginType(mcCluster, mcClusterLooper2, kMesonPdgs); + mcClusterLooper.setCursor(g.emmcparticleIds()[0]); + if (info.leptonOrigin == LeptonOrigin::Conversion) { + if (!mcCluster.has_mothers()) [[unlikely]] { + // conersion with no mother does not make any sense + info.leptonOrigin = LeptonOrigin::Other; + } else { + info.photonMotherId = mcCluster.mothersIds()[0]; + // since this is a real conversion where both daughters exist, check if the other conversion leg entered this cluster as well + for (size_t i = 1; i < g.emmcparticleIds().size(); ++i) { + mcClusterLooper.setCursor(g.emmcparticleIds()[i]); + if (std::abs(mcClusterLooper.pdgCode()) == PDG_t::kElectron && mcClusterLooper.pdgCode() == -1 * mcCluster.pdgCode()) { + mcClusterLooper2.setCursor(mcClusterLooper.globalIndex()); + int32_t otherConvMotherId = getMotherIndexFromChain(mcClusterLooper2, PDG_t::kGamma); + if (otherConvMotherId == info.photonMotherId) { + float energyFraction = (g.amplitude()[i] * g.e()) / mcClusterLooper.e(); + if (energyFraction >= MinAmpFraction && g.amplitude()[i] > MinAmpThreshold) { + info.isMergedConv = true; + } + // we found the sibling leg no need to search further + break; + } // if (otherConvMotherId == info.photonMotherId) + } // if (std::abs(mcClusterLooper.pdgCode()) == PDG_t::kElectron && mcClusterLooper.pdgCode() == -1 * mcCluster.pdgCode()) + } // end of loop over other cluster contributions + } // particle has mothers + } // if(info.leptonOrigin == LeptonOrigin::Conversion) + if (info.photonMotherId >= 0) { + mcClusterLooper.setCursor(info.photonMotherId); + info.photonOrigin = getPhotonOriginType(mcClusterLooper, mcClusterLooper2, kMesonPdgs, info.hardPartonId); + if (info.photonOrigin == PhotonOrigin::Decay) { + info.decayPhotonId = info.photonMotherId; + info.mesonId = o2::aod::pwgem::photonmeson::utils::mcutil::GetMesonInChain(mcClusterLooper, mcParticles, kMesonPdgs); } } } if (std::abs(mcCluster.pdgCode()) == PDG_t::kGamma) { info.isPhoton = true; + info.photonOrigin = getPhotonOriginType(mcCluster, mcClusterLooper, kMesonPdgs, info.hardPartonId); + if (info.photonOrigin == PhotonOrigin::Decay) { + info.mesonId = o2::aod::pwgem::photonmeson::utils::mcutil::GetMesonInChain(mcCluster, mcParticles, kMesonPdgs); + } } - - 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; } @@ -181,17 +199,36 @@ struct EmcalPhotonMcTask { static constexpr float PhiVUndefined = -999.f; static constexpr float Epsilon = 1.e-6f; + static constexpr std::array, 1> defaultCutsMl{{{0.0, 0.25}}}; + static constexpr std::array(TruthClass::NClasses)> kTruthClassNames = { - "Conversion", "PhotonPairSamePi0", "PhotonPairDiffPi0", "PhotonPairOnePi0", - "PhotonElectronSamePi0", "PhotonElectronDiffPi0", "PhotonElectronOnePi0", "BSPhotonElectron", - "ElectronPairSamePi0", "ElectronPairDiffPi0", "ElectronPairOnePi0", - "SplitPhotonCluster", "SplitLeptonCluster", "Background"}; + "Conversion", "GammaGammaSamePi0", "GammaGammaAnnihilation", "BSPhotonElectron", + "PhotonComptonElectronPair", "ElectronPairSamePi0", "CrossConvertedSiblings", "DalitzDecaySiblings", + "SplitPhotonCluster", "SplitLeptonCluster", "IndirGammaGammaSamePi0", "SameMesonIndirect", + "Background"}; Produces convTagCandidates; Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; Configurable writeTable{"writeTable", true, "write table for ML."}; - Configurable> classPrescale{"classPrescale", {1, 1, 700, 25, 1, 350, 15, 1, 1, 35, 2, 1, 1, 1, 1, 1, 1000}, "prescale factor per TruthClass, indexed 0..10 matching the enum order"}; + Configurable> classPrescale{"classPrescale", + { + 1, // Conversion + 1, // GammaGammaSamePi0 + 1, // GammaGammaAnnihilation + 1, // BSPhotonElectron + 1, // PhotonComptonElectronPair + 1, // ElectronPairSamePi0 + 1, // CrossConvertedSiblings + 1, // DalitzDecaySiblings + 1, // SplitPhotonCluster + 1, // SplitLeptonCluster + 1, // IndirGammaGammaSamePi0 + 1, // SameMesonIndirect + 5000, // Background + }, + "prescale factor per TruthClass, indexed 0..12 matching the enum order"}; + Configurable bkgPrescaleSeed{"bkgPrescaleSeed", 42, "seed for the background-prescale RNG"}; // configurable axis @@ -259,6 +296,25 @@ struct EmcalPhotonMcTask { Configurable cfgEnableQA{"cfgEnableQA", false, "flag to turn QA plots on/off"}; } mesonConfig; + struct : ConfigurableGroup { + std::string prefix = "mlConfig"; + Configurable useMlTagging{"useMlTagging", false, "use ML score instead of box cut for conversion tagging"}; + Configurable loadModelsFromCCDB{"loadModelsFromCCDB", false, "load ML model from CCDB"}; + Configurable> mlInputFeatures{ + "mlInputFeatures", + {"minv", "deltaEta", "deltaR", "phiv", "rConv", "totE", "e2", "e1", "deltaPhi"}, + "input feature names -- content and order must match the Python training FEATURES list"}; + Configurable mlModelPathLocal{"mlModelPathLocal", "/data/mhemmer/O2ML/code/conversion_tagging_bdt_conversion_splits_brems.onnx", "local ONNX model path"}; + Configurable> modelPathsCCDB{"modelPathsCCDB", std::vector{"Users/m/mhemmer/EM/ML/"}, "Paths of models on CCDB"}; + Configurable> onnxFileNames{"onnxFileNames", std::vector{"conversion_tagging_bdt_conversion_splits_brems.onnx"}, "ONNX file names for each pT bin (if not from CCDB full path)"}; + Configurable mlThreshold{"mlThreshold", 0.5f, "positive-class score threshold for tagging"}; + Configurable> cutsMl{"cutsMl", {defaultCutsMl[0].data(), 1, 2, {"pT bin 0"}, { + "score photon pairs", + "score conversion pairs", + }}, + "ML selections per pT bin"}; + } mlConfig; + SliceCache cache; using EMCalPhotons = soa::Join; @@ -276,13 +332,15 @@ struct EmcalPhotonMcTask { int8_t bTruthLabel{}; - Service ccdb{}; + o2::ccdb::CcdbApi ccdbApi; int mRunNumber{0}; float dBz{0.f}; std::mt19937 mRandGen; std::uniform_int_distribution mPrescaleDist; + o2::analysis::em::emcconv::EmMlResponseEMCConversion mMlResponse; + void defineEMEventCut() { fEMEventCut = EMPhotonEventCut("fEMEventCut", "fEMEventCut"); @@ -325,11 +383,6 @@ struct EmcalPhotonMcTask { mRunNumber = 0; dBz = 0; - ccdb->setURL(ccdbUrl); - ccdb->setCaching(true); - ccdb->setLocalObjectValidityChecking(); - ccdb->setFatalWhenNull(false); - defineEMEventCut(); defineEMCCut(); fEMCCut.addQAHistograms(®istry); @@ -339,8 +392,8 @@ struct EmcalPhotonMcTask { const AxisSpec thnAxisPtRec{thnConfigAxisPt, "#it{p}_{T} (GeV/#it{c})"}; const AxisSpec thnAxisInvMass{thnConfigAxisInvMass, "#it{M}_{#gamma#gamma} (GeV/#it{c}^{2})"}; - const AxisSpec thnAxisrConvRec{100, 0, 500, "#it{R}_{rec}"}; - const AxisSpec thnAxisrConvGen{100, 0, 500, "#it{R}_{gen}"}; + const AxisSpec thnAxisrConvRec{1000, 0, 500, "#it{R}_{rec}"}; + const AxisSpec thnAxisrConvGen{1000, 0, 500, "#it{R}_{gen}"}; const AxisSpec thnAxisDeltaEta{thnConfigAxisDeltaEta, "#Delta#it{eta}"}; const AxisSpec thnAxisDeltaPhi{thnConfigAxisDeltaPhi, "#Delta#it{#varphi} (rad)"}; @@ -348,6 +401,8 @@ struct EmcalPhotonMcTask { const AxisSpec thnAxisTagging{static_cast(TagDecision::NTags), -0.5, static_cast(TagDecision::NTags) - 0.5, ""}; const AxisSpec thnAxisClasses{static_cast(TruthClass::NClasses), -0.5, static_cast(TruthClass::NClasses) - 0.5, ""}; + const AxisSpec thnAxisM02{250, 0., 2.5, "#it{M}_{02}"}; + AxisSpec thnAxisCentOrMult{1, 0., 1., "Centrality/Multiplicity"}; // placeholder, overwritten in init if (useCent.value) { // PbPb: use centrality @@ -361,21 +416,17 @@ struct EmcalPhotonMcTask { // set bin labels once at init, so histogram is human-readable without decoding the enum hTruthLabel->GetXaxis()->SetBinLabel(static_cast(TruthClass::Conversion) + 1, "Conversion"); - hTruthLabel->GetXaxis()->SetBinLabel(static_cast(TruthClass::PhotonPairSamePi0) + 1, "PhotonPairSamePi0"); - hTruthLabel->GetXaxis()->SetBinLabel(static_cast(TruthClass::PhotonPairDiffPi0) + 1, "PhotonPairDiffPi0"); - hTruthLabel->GetXaxis()->SetBinLabel(static_cast(TruthClass::PhotonPairOnePi0) + 1, "PhotonPairOnePi0"); - 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::GammaGammaSamePi0) + 1, "GammaGammaSamePi0"); + hTruthLabel->GetXaxis()->SetBinLabel(static_cast(TruthClass::GammaGammaAnnihilation) + 1, "GammaGammaAnnihilation"); hTruthLabel->GetXaxis()->SetBinLabel(static_cast(TruthClass::BSPhotonElectron) + 1, "BSPhotonElectron"); + hTruthLabel->GetXaxis()->SetBinLabel(static_cast(TruthClass::PhotonComptonElectronPair) + 1, "PhotonComptonElectronPair"); 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::CrossConvertedSiblings) + 1, "CrossConvertedSiblings"); + hTruthLabel->GetXaxis()->SetBinLabel(static_cast(TruthClass::DalitzDecaySiblings) + 1, "DalitzDecaySiblings"); hTruthLabel->GetXaxis()->SetBinLabel(static_cast(TruthClass::SplitPhotonCluster) + 1, "SplitPhotonCluster"); hTruthLabel->GetXaxis()->SetBinLabel(static_cast(TruthClass::SplitLeptonCluster) + 1, "SplitLeptonCluster"); - hTruthLabel->GetXaxis()->SetBinLabel(static_cast(TruthClass::PhotonBSPhotonPair) + 1, "PhotonBSPhotonPair"); - hTruthLabel->GetXaxis()->SetBinLabel(static_cast(TruthClass::ElectronBSPhotonPair) + 1, "ElectronBSPhotonPair"); - hTruthLabel->GetXaxis()->SetBinLabel(static_cast(TruthClass::BSPhotonPair) + 1, "BSPhotonPair"); + hTruthLabel->GetXaxis()->SetBinLabel(static_cast(TruthClass::IndirGammaGammaSamePi0) + 1, "IndirGammaGammaSamePi0"); + hTruthLabel->GetXaxis()->SetBinLabel(static_cast(TruthClass::SameMesonIndirect) + 1, "SameMesonIndirect"); 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}}); @@ -392,6 +443,43 @@ struct EmcalPhotonMcTask { hConfusionMatrixConversionTagging->GetYaxis()->SetBinLabel(4, "background"); hConfusionMatrixConversionTagging->GetYaxis()->SetBinLabel(5, "#gamma"); + registry.add("hBSRadius", "Radius of BS photons;;counts", HistType::kTH1D, {thnAxisrConvGen}); + registry.add("hLeptonRadius", "Radius of leptons;;counts", HistType::kTH1D, {thnAxisrConvGen}); + registry.add("hConvLeptonRadius", "Radius of leptons from conversions;;counts", HistType::kTH1D, {thnAxisrConvGen}); + + registry.add("Photon/M02", "M02 distribution;;counts", HistType::kTH1D, {thnAxisM02}); + auto hPhotonProcess = registry.add("Photon/hProcess", "Production process type", HistType::kTH1D, {{kMaxMCProcess, -0.5, kMaxMCProcess - 0.5}}); + for (int i = 0; i < kMaxMCProcess; ++i) { + hPhotonProcess->GetXaxis()->SetBinLabel(i + 1, TMCProcessName[i]); + } + registry.addClone("Photon/", "Electron/"); + registry.addClone("Photon/", "Positron/"); + registry.addClone("Photon/", "BSPhoton/"); + registry.addClone("Photon/", "MergedConv/"); + registry.addClone("Photon/", "ConvElectron/"); + registry.addClone("Photon/", "ConvPositron/"); + registry.addClone("Photon/", "Other/"); + registry.addClone("Photon/", "Lepton/"); + + auto hClusterType = registry.add("hClusterType", "Truth label distribution;;Counts", HistType::kTH2D, {{8, -0.5, 7.5}, thnAxisPtRec}); + hClusterType->GetXaxis()->SetBinLabel(1, "Photon"); + hClusterType->GetXaxis()->SetBinLabel(2, "Electron"); + hClusterType->GetXaxis()->SetBinLabel(3, "Positron"); + hClusterType->GetXaxis()->SetBinLabel(4, "BSPhoton"); + hClusterType->GetXaxis()->SetBinLabel(5, "MergedConv"); + hClusterType->GetXaxis()->SetBinLabel(6, "Conv electron"); + hClusterType->GetXaxis()->SetBinLabel(7, "Conv positron"); + hClusterType->GetXaxis()->SetBinLabel(8, "Other"); + + auto hBSLeptonFate = registry.add("hBSLeptonFate", "Fate of the Bremsstrahlungsphotons mother lepton;;Counts", HistType::kTH1D, {{3, -0.5, 2.5}}); + hBSLeptonFate->GetXaxis()->SetBinLabel(1, "dominant"); + hBSLeptonFate->GetXaxis()->SetBinLabel(2, "non dominant"); + hBSLeptonFate->GetXaxis()->SetBinLabel(3, "absent"); + + if (mlConfig.useMlTagging.value) { + registry.add("hMlScore", "BDT score;;Counts", HistType::kTH1D, {{100, -10, 10}}); + } + mRandGen.seed(bkgPrescaleSeed.value); if (classPrescale.value.size() != kTruthClassNames.size()) { @@ -403,6 +491,7 @@ struct EmcalPhotonMcTask { for (size_t i = 0; i < kTruthClassNames.size(); ++i) { LOG(info) << " [" << i << "] " << kTruthClassNames[i] << " -> prescale = " << classPrescale.value[i]; } + }; // end init template @@ -413,10 +502,34 @@ struct EmcalPhotonMcTask { } mRunNumber = collision.runNumber(); - auto run3grp_timestamp = collision.timestamp(); + auto timestamp = collision.timestamp(); // Fetch magnetic field from ccdb for current collision dBz = collision.grpMagField().getNominalL3Field(); - LOG(info) << "Retrieved GRP for timestamp " << run3grp_timestamp << " with magnetic field of " << dBz << " kZG"; + LOG(info) << "Retrieved GRP for timestamp " << timestamp << " with magnetic field of " << dBz << " kZG"; + + if (mlConfig.useMlTagging.value) { + // single bin, full range -- see earlier discussion: pT/energy-binned + // thresholds are a straightforward future extension of this same + // machinery if ever needed, not used right now + std::vector binsLimits = {0., 1000.}; + std::vector cutDir = { + static_cast(o2::cuts_ml::CutDirection::CutNot), // class 0 (negative-class prob) -- no cut + static_cast(o2::cuts_ml::CutDirection::CutSmaller), // class 1 (positive-class prob) -- reject if score < threshold + }; + + mMlResponse.configure(binsLimits, mlConfig.cutsMl, cutDir, /*nClasses=*/2); + mMlResponse.cacheInputFeaturesIndices(mlConfig.mlInputFeatures.value); + if (mlConfig.loadModelsFromCCDB) { + ccdbApi.init(ccdbUrl); + mMlResponse.setModelPathsCCDB(mlConfig.onnxFileNames, ccdbApi, mlConfig.modelPathsCCDB.value, timestamp); + } else { + mMlResponse.setModelPathsLocal({mlConfig.mlModelPathLocal.value}); + } + mMlResponse.init(); + + LOG(info) << "ML conversion tagging enabled -- model: " << mlConfig.mlModelPathLocal.value + << ", threshold: " << mlConfig.mlThreshold.value; + } } template @@ -499,6 +612,7 @@ struct EmcalPhotonMcTask { EMBitFlags emcFlagsFromTrueMesonSameGamma(clusters.size()); EMBitFlags emcFlagsFromTrueConversion(clusters.size()); EMBitFlags emcFlagsTagging(clusters.size()); + EMBitFlags emcFlagsMlTagging(clusters.size()); EMBitFlags emcFlags(clusters.size()); if (clusters.size() > 0) { fEMCCut.AreSelectedRunning(emcFlags, clusters, matchedPrims, matchedSeconds, ®istry); @@ -513,7 +627,9 @@ struct EmcalPhotonMcTask { for (const auto& collision : collisions) { initCCDB(collision); - isFullEventSelected(collision, true); + if (!isFullEventSelected(collision, true)) { + continue; + } float centOrMult = getCentralityOrMultiplicity(collision); @@ -565,18 +681,30 @@ struct EmcalPhotonMcTask { emcFlagsTagging.set(g2.globalIndex()); } + if (mlConfig.useMlTagging.value) { + o2::analysis::em::EMCConversionCandidate candidate{ + .mMinv = static_cast(vMeson.M()), .mDeltaEta = deltaEta, .mDeltaR = std::hypot(deltaEta, deltaPhi), .mPhiv = phiV, .mRConv = rConv, .mTotE = (g2.e() + g1.e()), .mE2 = g2.e(), .mE1 = g1.e(), .mDeltaPhi = deltaPhi}; + std::vector mlInput = mMlResponse.getInputFeatures(candidate); + std::vector mlOutput; + bool isTagged = mMlResponse.isSelectedMl(mlInput, 0.f, mlOutput); + if (isTagged) { + emcFlagsMlTagging.set(g1.globalIndex()); + emcFlagsMlTagging.set(g2.globalIndex()); + } + registry.fill(HIST("hMlScore"), mlOutput[1]); // positive-class score, always, tagged or not + } + // set MC particle cursors to the largest cluster contributor mcCluster1.setCursor(g1.emmcparticleIds()[0]); mcCluster2.setCursor(g2.emmcparticleIds()[0]); - bool areFromSamePi0 = false; + auto c1 = classifyCluster(g1, mcCluster1, mcClusterLooper, mcClusterLooper2, mcParticles); + auto c2 = classifyCluster(g2, mcCluster2, mcClusterLooper, mcClusterLooper2, mcParticles); + + const bool areFromSamePi0 = c1.mesonId >= 0 && c1.mesonId == c2.mesonId; bool areConversionLegs = false; bool areSplitPhotonCluster = false; bool areSplitLeptonCluster = false; - bool areBSPhotonElectron = false; - - auto c1 = classifyCluster(g1, mcCluster1, mcClusterLooper, mcClusterLooper2, mcParticles); - auto c2 = classifyCluster(g2, mcCluster2, mcClusterLooper, mcClusterLooper2, mcParticles); // 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 @@ -590,34 +718,29 @@ struct EmcalPhotonMcTask { } } - const bool isAnyBSPhoton = c1.isFromBremsstrahlung || c2.isFromBremsstrahlung; - const bool areBSPhotons = c1.isFromBremsstrahlung && c2.isFromBremsstrahlung; - // if they are not a split cluster check for proper conversion pair - if (!isSameDominantParticle && c1.isFromConv && c2.isFromConv && c1.convMotherId == c2.convMotherId) { + if (!isSameDominantParticle && c1.leptonOrigin == LeptonOrigin::Conversion && c2.leptonOrigin == LeptonOrigin::Conversion && c1.photonMotherId == c2.photonMotherId) { emcFlagsFromTrueConversion.set(g1.globalIndex()); emcFlagsFromTrueConversion.set(g2.globalIndex()); areConversionLegs = true; } // 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); + if (!isSameDominantParticle && c1.decayPhotonId >= 0 && c2.decayPhotonId >= 0) { + mcPhoton1.setCursor(c1.decayPhotonId); + mcPhoton2.setCursor(c2.decayPhotonId); mcMother.setCursor(mcPhoton1.mothersIds()[0]); if (mcMother.producedByGenerator()) { - if (c1.photonId == c2.photonId) { + if (c1.mesonId == c2.mesonId) { // 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); + const bool photonIsBS = (c1.isPhoton && c1.photonOrigin == PhotonOrigin::Bremsstrahlung) || (c2.isPhoton && c2.photonOrigin == PhotonOrigin::Bremsstrahlung); if (photonIsBS && ((c1.isLepton && c2.isPhoton) || (c2.isLepton && c1.isPhoton))) { - areBSPhotonElectron = true; + // nothing } 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()); emcFlagsFromTrueMeson.set(g2.globalIndex()); } @@ -625,46 +748,45 @@ struct EmcalPhotonMcTask { } bTruthLabel = static_cast(TruthClass::Background); + if (!mcCluster1.has_mothers() || !mcCluster2.has_mothers()) { + registry.fill(HIST("hTruthLabel"), bTruthLabel, vMeson.Pt()); + + // final tree values plus filling + const int prescale = classPrescale.value[static_cast(bTruthLabel)]; + const bool keepThisRow = (prescale <= 1) || (std::uniform_int_distribution(0, prescale - 1)(mRandGen) == 0); + if (writeTable.value && keepThisRow) { + convTagCandidates(collision.globalIndex(), vMeson.M(), harmonicET, deltaEta, deltaPhi, phiV, g1.e(), g2.e(), g1.m02(), g2.m02(), g1.time(), g2.time(), g1.nCells(), g2.nCells(), c1.purity, c2.purity, bTruthLabel, centOrMult); + } + continue; + } 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 (areBSPhotonElectron) { + } else if (c1.photonOrigin == PhotonOrigin::Annihilation && c2.photonOrigin == PhotonOrigin::Annihilation && mcCluster1.mothersIds()[0] == mcCluster2.mothersIds()[0]) { + bTruthLabel = static_cast(TruthClass::GammaGammaAnnihilation); + } else if ((c1.photonOrigin == PhotonOrigin::Bremsstrahlung && mcCluster1.mothersIds()[0] == mcCluster2.globalIndex()) || (c2.photonOrigin == PhotonOrigin::Bremsstrahlung && mcCluster2.mothersIds()[0] == mcCluster1.globalIndex())) { bTruthLabel = static_cast(TruthClass::BSPhotonElectron); - } else if (areBSPhotons && (c1.isFromPi0 || c2.isFromPi0)) { - bTruthLabel = static_cast(TruthClass::BSPhotonPair); - } else if (isAnyBSPhoton && (c1.isFromPi0 || c2.isFromPi0) && ((c1.isPhoton && c2.isLepton) || (c2.isPhoton && c1.isLepton))) { - bTruthLabel = static_cast(TruthClass::ElectronBSPhotonPair); - } else if (isAnyBSPhoton && (c1.isFromPi0 || c2.isFromPi0) && (c1.isPhoton && c2.isPhoton)) { - bTruthLabel = static_cast(TruthClass::PhotonBSPhotonPair); - } 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) { + } else if ((c1.leptonOrigin == LeptonOrigin::Compton && mcCluster1.mothersIds()[0] == mcCluster2.globalIndex()) || (c2.leptonOrigin == LeptonOrigin::Compton && mcCluster2.mothersIds()[0] == mcCluster1.globalIndex())) { + bTruthLabel = static_cast(TruthClass::PhotonComptonElectronPair); + } else if (c1.leptonOrigin == LeptonOrigin::DirectMesonDecay && c2.leptonOrigin == LeptonOrigin::DirectMesonDecay && mcCluster1.mothersIds()[0] == mcCluster2.mothersIds()[0]) { // both cluster are leptons that come from the same meson decay + mcMother.setCursor(mcCluster1.mothersIds()[0]); + if (mcMother.daughtersIds().size() == 2) { bTruthLabel = static_cast(TruthClass::ElectronPairSamePi0); + } else if (mcMother.daughtersIds().size() == 3) { + bTruthLabel = static_cast(TruthClass::DalitzDecaySiblings); } - } 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 (c1.leptonOrigin == LeptonOrigin::Conversion && c2.leptonOrigin == LeptonOrigin::Conversion && mcCluster1.mothersIds()[0] != mcCluster2.mothersIds()[0] && areFromSamePi0) { // both cluster are leptons that come from different conversions that come from the same meson + bTruthLabel = static_cast(TruthClass::CrossConvertedSiblings); + } else if (c1.photonOrigin == PhotonOrigin::Decay && c2.photonOrigin == PhotonOrigin::Decay && areFromSamePi0) { // both clusters are photons from decay from the same meson + bTruthLabel = static_cast(TruthClass::GammaGammaSamePi0); + } else if (c1.isPhoton && c2.isPhoton && areFromSamePi0) { + bTruthLabel = static_cast(TruthClass::IndirGammaGammaSamePi0); + } else if (areFromSamePi0) { // both cluster do not fit into one of the categories above, but they share a common meson ancestry + bTruthLabel = static_cast(TruthClass::SameMesonIndirect); } - registry.fill(HIST("hTruthLabel"), bTruthLabel, vMeson.Pt()); // final tree values plus filling @@ -674,11 +796,58 @@ struct EmcalPhotonMcTask { convTagCandidates(collision.globalIndex(), vMeson.M(), harmonicET, deltaEta, deltaPhi, phiV, g1.e(), g2.e(), g1.m02(), g2.m02(), g1.time(), g2.time(), g1.nCells(), g2.nCells(), c1.purity, c2.purity, bTruthLabel, centOrMult); } } // pair loop + + // key: MC particle global index -> list of (cluster global index, contributor rank) + std::unordered_map>> particleToClusterContributions; + for (const auto& cluster : photonsEMCPerCollision) { + if (!emcFlags.test(cluster.globalIndex())) { + continue; + } + const auto& ids = cluster.emmcparticleIds(); + for (size_t i = 0; i < ids.size(); ++i) { + particleToClusterContributions[ids[i]].emplace_back(cluster.globalIndex(), i); + } + } // cluster loop + + for (const auto& cluster : photonsEMCPerCollision) { + if (!emcFlags.test(cluster.globalIndex())) { + continue; + } + auto c1 = classifyCluster(cluster, mcCluster1, mcClusterLooper, mcClusterLooper2, mcParticles); + + // NEW: bremsstrahlung sibling-fate check, mirrors the conversion one above + if (c1.photonOrigin == PhotonOrigin::Bremsstrahlung) { + // mcCluster1 is currently sitting on the bremsstrahlung photon itself + // (classifyCluster leaves it there for the photon branch) -- its + // immediate mother is the radiating lepton we want to look up. + if (mcCluster1.has_mothers()) { + const int radiatingLeptonId = mcCluster1.mothersIds()[0]; + + auto it = particleToClusterContributions.find(radiatingLeptonId); + if (it == particleToClusterContributions.end()) { + registry.fill(HIST("hBSLeptonFate"), 2); // radiating lepton absent from any cluster + } else { + bool isDominantSomewhere = false; + for (const auto& [clusterId, rank] : it->second) { + if (clusterId == cluster.globalIndex()) { + continue; // skip itself (shouldn't normally match, but same safety as before) + } + if (rank == 0) { + isDominantSomewhere = true; + break; + } + } + registry.fill(HIST("hBSLeptonFate"), isDominantSomewhere ? 0 : 1); // 0=dominant elsewhere, 1=leakage-only + } + } + } + } // cluster loop } // collision loop std::vector photonSeen(mcParticles.size(), false); // this decay photon has >=1 resolved cluster std::vector photonTagged(mcParticles.size(), false); // >=1 of those clusters got conversion-tagged auto collision = collisions.begin(); + for (const auto& cluster : clusters) { if (!(emcFlags.test(cluster.globalIndex()))) { continue; @@ -690,6 +859,43 @@ struct EmcalPhotonMcTask { if (cluster.pmeventId() > collision.globalIndex()) { collision.setCursor(cluster.pmeventId()); } + if (!isFullEventSelected(collision, false)) { + continue; + } + + auto clusterMcInfo = classifyCluster(cluster, mcCluster1, mcClusterLooper, mcClusterLooper2, mcParticles); + if (clusterMcInfo.photonOrigin == PhotonOrigin::Bremsstrahlung) { + registry.fill(HIST("hBSRadius"), clusterMcInfo.radius); + registry.fill(HIST("BSPhoton/M02"), cluster.m02()); + registry.fill(HIST("hClusterType"), 3, cluster.e()); + } else if (clusterMcInfo.isMergedConv) { + registry.fill(HIST("MergedConv/M02"), cluster.m02()); + registry.fill(HIST("hClusterType"), 4, cluster.e()); + } else if (clusterMcInfo.leptonOrigin == LeptonOrigin::Conversion) { + registry.fill(HIST("hConvLeptonRadius"), clusterMcInfo.radius); + if (mcCluster1.pdgCode() == PDG_t::kElectron) { + registry.fill(HIST("ConvElectron/M02"), cluster.m02()); + registry.fill(HIST("hClusterType"), 5, cluster.e()); + } else if (mcCluster1.pdgCode() == PDG_t::kPositron) { + registry.fill(HIST("ConvPositron/M02"), cluster.m02()); + registry.fill(HIST("hClusterType"), 6, cluster.e()); + } + } else if (clusterMcInfo.isLepton) { + if (mcCluster1.pdgCode() == PDG_t::kElectron) { + registry.fill(HIST("Electron/M02"), cluster.m02()); + registry.fill(HIST("hClusterType"), 1, cluster.e()); + } else if (mcCluster1.pdgCode() == PDG_t::kPositron) { + registry.fill(HIST("Positron/M02"), cluster.m02()); + registry.fill(HIST("hClusterType"), 2, cluster.e()); + } + registry.fill(HIST("hLeptonRadius"), clusterMcInfo.radius); + } else if (clusterMcInfo.isPhoton) { + registry.fill(HIST("hClusterType"), 0, cluster.e()); + registry.fill(HIST("Photon/M02"), cluster.m02()); + } else { + registry.fill(HIST("hClusterType"), 7, cluster.e()); + registry.fill(HIST("Other/M02"), cluster.m02()); + } mcCluster1.setCursor(cluster.emmcparticleIds()[0]); int photonid1 = o2::aod::pwgem::photonmeson::utils::mcutil::FindMotherInChain(mcCluster1, mcParticles, std::vector{PDG_t::kPi0, Pdg::kEta, Pdg::kOmega, Pdg::kEtaPrime}); @@ -709,6 +915,10 @@ struct EmcalPhotonMcTask { if (mcCluster1.pdgCode() == PDG_t::kGamma) { registry.fill(HIST("EMCal/ConfusionMatrixConversionTagging"), emcFlagsTagging.test(cluster.globalIndex()) ? 0 : 1, static_cast(ClusterTruthClass::Photon)); + registry.fill(HIST("Photon/hProcess"), mcCluster1.getProcess()); + } + if (std::abs(mcCluster1.pdgCode()) == PDG_t::kElectron) { + registry.fill(HIST("Lepton/hProcess"), mcCluster1.getProcess()); } if (!emcFlagsFromTrueConversion.test(cluster.globalIndex())) { registry.fill(HIST("EMCal/ConfusionMatrixConversionTagging"), emcFlagsTagging.test(cluster.globalIndex()) ? 0 : 1, static_cast(ClusterTruthClass::Conversion)); @@ -756,7 +966,7 @@ struct EmcalPhotonMcTask { if (lost) { registry.fill(HIST("EMCal/hPi0BothResolvedLost"), 1.0); // "lost" bin } - } + } // end of loop over mc particles } PROCESS_SWITCH(EmcalPhotonMcTask, processEmcal, "Process for pcm and emcal photons", true); diff --git a/PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx b/PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx index db7281218ba..146f38007f2 100644 --- a/PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx +++ b/PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx @@ -173,7 +173,9 @@ 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"}; + Configurable useEMCal{"useEMCal", false, "flag to use EMCal clusters"}; + Configurable useDCal{"useDCal", false, "flag to use DCal clusters"}; + Configurable useCrosspairs{"useCrosspairs", true, "flag to allow pairing of EMCal with DCal clusters. If this is set, useEMCal and useDCal are ignored!"}; } emccuts; V0PhotonCut fV0PhotonCut; @@ -253,8 +255,8 @@ struct TaskPi0FlowEMC { int runNow = 0; int runBefore = -1; - static constexpr float MaxPhiEMCal = 3.5f; - static constexpr uint16_t MaxPhiEMCalUint = static_cast(INT16_MAX); // Maximum value currently useable for partitions for some weird reason, but luckily enough + static constexpr float MaxPhiEMCal = 3.9f; // exatly the middle between EMCal and DCal + static constexpr uint16_t MaxPhiEMCalUint = static_cast(39000u); // exatly the middle between EMCal and DCal but as uint16_t that is used for storing phi values in derived data // 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); @@ -1093,7 +1095,6 @@ struct TaskPi0FlowEMC { fEMCCut.AreSelectedRunning(flags, clusters, matchedPrims, matchedSeconds, ®istry); for (const auto& collision : collisions) { - if (!isFullEventSelected(collision, true)) { continue; } @@ -1121,10 +1122,13 @@ struct TaskPi0FlowEMC { registry.fill(HIST("clusterQA/hClusterEtaPhiAfter"), photon.phi(), photon.eta()); // after cuts } } - if (emccuts.separateEMCalDCal.value) { + if (emccuts.useEMCal.value && !emccuts.useCrosspairs.value) { runPairingLoop(collision, emcalPhotonsPerCollision, emcalPhotonsPerCollision, flags, flags); + } + if (emccuts.useDCal.value && !emccuts.useCrosspairs.value) { runPairingLoop(collision, dcalPhotonsPerCollision, dcalPhotonsPerCollision, flags, flags); - } else { + } + if (emccuts.useCrosspairs.value) { runPairingLoop(collision, photonsPerCollision, photonsPerCollision, flags, flags); } if (rotationConfig.cfgDoRotation.value) { @@ -1181,8 +1185,11 @@ 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 + if (emccuts.useEMCal.value && !emccuts.useCrosspairs.value && (!isEMCalRegion(g1.phi()) || !isEMCalRegion(g2.phi()))) { + continue; + } + if (emccuts.useDCal.value && !emccuts.useCrosspairs.value && (isEMCalRegion(g1.phi()) || isEMCalRegion(g2.phi()))) { + continue; } // Cut edge clusters away, similar to rotation method to ensure same acceptance is used diff --git a/PWGEM/PhotonMeson/Utils/MCUtilities.h b/PWGEM/PhotonMeson/Utils/MCUtilities.h index 5dd76b4ec8b..4a31a33e51b 100644 --- a/PWGEM/PhotonMeson/Utils/MCUtilities.h +++ b/PWGEM/PhotonMeson/Utils/MCUtilities.h @@ -16,9 +16,12 @@ #ifndef PWGEM_PHOTONMESON_UTILS_MCUTILITIES_H_ #define PWGEM_PHOTONMESON_UTILS_MCUTILITIES_H_ +#include "PWGEM/PhotonMeson/Utils/ParticleOrigin.h" + #include #include +#include #include #include @@ -31,6 +34,9 @@ //_______________________________________________________________________ namespace o2::aod::pwgem::photonmeson::utils::mcutil { + +constexpr float kVertexEps = 1e-4f; // cm + template bool IsPhysicalPrimary(TTrack const& mctrack) { @@ -89,6 +95,7 @@ int IsXFromY(T const& mctrack, TMCs const& mcTracks, const int pdgX, const int p } return -1; } + //_______________________________________________________________________ // Go up the decay chain of a mcparticle looking for a mother with the given pdg codes, if found return this mothers daughter // E.g. Find the gamma that was created in a pi0 or eta decay @@ -106,6 +113,51 @@ int FindMotherInChain(T const& mcparticle, TMCs const& mcparticles, TTargetPDGs } return FindMotherInChain(mother, mcparticles, motherpdgs, Depth - 1); } + +//_______________________________________________________________________ +/// \brief Go up the decay chain of a mcparticle looking for a mother with the given pdg codes, +/// and return that MOTHER's own index (unlike FindMotherInChain, which returns its daughter). +/// Two different particles that share the same meson ancestor will resolve to the same value here. +/// \param mcparticle iterator of McParticles +/// \param mcparticles table of McParticles +/// \param motherpdgs ranges of mother PDG values to compare against +/// \param Depth how many links should this go up +template +int GetMesonInChain(T const& mcparticle, TMCs const& mcparticles, TTargetPDGs const& motherpdgs, const int Depth = 15) +{ + int decayChildIdx = FindMotherInChain(mcparticle, mcparticles, motherpdgs, Depth); + if (decayChildIdx < 0) { + return -1; + } + auto decayChild = mcparticles.iteratorAt(decayChildIdx); + return decayChild.mothersIds()[0]; // the meson itself, not its daughter +} + +//_______________________________________________________________________ +/// \brief Go up the decay chain of a mcparticle looking for a mother with the given pdg codes, if found return this mothers daughter +/// E.g. Find the gamma that was created in a pi0 or eta decay +/// \param mcIter iterator of mcparticle -- WILL BE MODIFIED/CONSUMED by this function +/// \param motherPdgs target mother PDG values +/// \param depth how many steps in the chain this check should go maximum before failing +template +int findMotherInChain(T& mcIter, TTargetPDGs const& motherPdgs, const int depth = 50) +{ + int currentIndex = mcIter.globalIndex(); // the node whose immediate mother we're about to test + + for (int d = 0; d < depth; ++d) { + if (!mcIter.has_mothers()) { + return -1; + } + const int motherId = mcIter.mothersIds()[0]; + mcIter.setCursor(motherId); + if (std::find(motherPdgs.begin(), motherPdgs.end(), mcIter.pdgCode()) != motherPdgs.end()) { + return currentIndex; // mother matches -- return the node directly below it + } + currentIndex = motherId; // no match -- this mother becomes "current" for the next step up + } + return -1; +} + //_______________________________________________________________________ template int IsEleFromPC(T const& mctrack, TMCs const& mcTracks) @@ -321,7 +373,7 @@ bool isMotherPDG(const T& mcparticle, T& mcparticleWorking, const int motherPDG, /// \param mcCursor iterator of mcparticle /// \param iter shared iterator used to walk to the mother template -bool isFromBremsstrahlung(TIter const& mcCursor, TIter& iter) +bool isFromBremsstrahlung(TIter const& mcCursor, TIter& iter, int& motherId) { if (!mcCursor.has_mothers()) { return false; @@ -329,7 +381,10 @@ bool isFromBremsstrahlung(TIter const& mcCursor, TIter& iter) if (mcCursor.pdgCode() != PDG_t::kGamma) { return false; // only a photon can itself be a bremsstrahlung emission } - const int motherId = mcCursor.mothersIds()[0]; + if (mcCursor.mothersIds().size() != 1) { + return false; // mother can be only a single lepton, otherwise it might be e+e- annihilation or something else + } + motherId = mcCursor.mothersIds()[0]; iter.setCursor(motherId); return std::abs(iter.pdgCode()) == PDG_t::kElectron; } @@ -337,48 +392,93 @@ bool isFromBremsstrahlung(TIter const& mcCursor, TIter& iter) //_______________________________________________________________________ /// \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 -/// \param mcparticle iterator of mxparticle, WILL BE CHANGED by this function! +/// \param mcParticle iterator of mxparticle, WILL BE CHANGED by this function! /// \param motherPDG target mother PDG value /// \param depth how many steps in the chain this check should go maximum before failing template -int32_t getMotherIndexFromChain(T& mcparticle, const int motherPDG, const int depth = 10) // o2-linter: disable=pdg/explicit-code (false positive) +int32_t getMotherIndexFromChain(T& mcParticle, const int motherPDG, const int depth = 10) { - if (!mcparticle.has_mothers() || depth < 1) { - return -1; - } - - int32_t motherid = mcparticle.mothersIds()[0]; - mcparticle.setCursor(motherid); - if (mcparticle.pdgCode() == motherPDG) { - return motherid; // The mother has the required pdg code, so return its daughters global mc particle code. + for (int d = 0; d < depth; ++d) { + if (!mcParticle.has_mothers()) { + return -1; + } + const int32_t motherid = mcParticle.mothersIds()[0]; + mcParticle.setCursor(motherid); + if (mcParticle.pdgCode() == motherPDG) { + return motherid; + } } - return getMotherIndexFromChain(mcparticle, motherPDG, depth - 1); + return -1; } //_______________________________________________________________________ -/// \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 the photon's id, if primary electron return -1 -/// \param mcparticle iterator of mcparticle, NOT modified by this function -/// \param mcparticleWorking a second iterator of the SAME table, used as scratch space to walk up the chain -- caller must supply this so the function doesn't construct its own -/// \param motherPDG target mother PDG value -/// \param depth how many steps in the chain this check should go maximum before failing -template -int32_t getMotherIndexFromChain(const T& mcparticle, T& mcparticleWorking, const int motherPDG, const int depth = 10) // o2-linter: disable=pdg/explicit-code (false positive) +/// \brief Obtains given photon mcpartlices origin type +/// \param mcPhoton mcparticle iterator of photon +/// \param iter mcparticle iterator used to walk to the mother +/// \param motherDPGs list of pdg values of mothers that would be excepted for PhotonOrigin::Decay +/// \return given photon mcpartlices origin type +template +o2::analysis::em::PhotonOrigin getPhotonOriginType(TIter const& mcPhoton, TIter& mcIter, TTargetPDGs const& motherPdgs, int& motherId) { - if (!mcparticle.has_mothers() || depth < 1) { - return -1; + switch (mcPhoton.getProcess()) { + case TMCProcess::kPBrem: + return o2::analysis::em::PhotonOrigin::Bremsstrahlung; + case TMCProcess::kPAnnihilation: + return o2::analysis::em::PhotonOrigin::Annihilation; + case TMCProcess::kPHadronic: + return o2::analysis::em::PhotonOrigin::Hadronic; + default: + break; } - - int32_t motherid = mcparticle.mothersIds()[0]; - mcparticleWorking.setCursor(motherid); - if (mcparticleWorking.pdgCode() == motherPDG) { - return motherid; + if (!mcPhoton.has_mothers()) { + return o2::analysis::em::PhotonOrigin::Other; + } + mcIter.setCursor(mcPhoton.globalIndex()); + motherId = findMotherInChain(mcIter, motherPdgs); + if (motherId >= 0) { + return o2::analysis::em::PhotonOrigin::Decay; + } + motherId = mcPhoton.mothersIds()[0]; + mcIter.setCursor(motherId); + if ((std::abs(mcIter.pdgCode()) >= PDG_t::kDown && std::abs(mcIter.pdgCode()) <= PDG_t::kTop) || std::abs(mcIter.pdgCode()) == PDG_t::kGluon) { + return o2::analysis::em::PhotonOrigin::Direct; } - return getMotherIndexFromChain(mcparticleWorking, mcparticleWorking, motherPDG, depth - 1); + return o2::analysis::em::PhotonOrigin::Other; } //_______________________________________________________________________ +/// \brief Obtains given lepton mcpartlices origin type +/// \param mcLepton mcparticle iterator of lepton +/// \param iter mcparticle iterator used to walk to the mother +/// \param motherDPGs list of pdg values of mothers that would be excepted for PhotonOrigin::Decay +/// \return given lepton mcpartlices origin type +template +o2::analysis::em::LeptonOrigin getLeptonOriginType(TIter const& mcLepton, TIter& mcIter, TTargetPDGs const& motherPdgs) +{ + switch (mcLepton.getProcess()) { + case TMCProcess::kPPair: + return o2::analysis::em::LeptonOrigin::Conversion; + case TMCProcess::kPCompton: + return o2::analysis::em::LeptonOrigin::Compton; + case TMCProcess::kPPhotoelectric: + return o2::analysis::em::LeptonOrigin::PhotoElectric; + case TMCProcess::kPDeltaRay: + return o2::analysis::em::LeptonOrigin::DeltaRay; + case TMCProcess::kPDecay: { + if (!mcLepton.has_mothers()) { + return o2::analysis::em::LeptonOrigin::Other; + } + mcIter.setCursor(mcLepton.mothersIds()[0]); + if (std::find(motherPdgs.begin(), motherPdgs.end(), mcIter.pdgCode()) != motherPdgs.end()) { + return o2::analysis::em::LeptonOrigin::DirectMesonDecay; + } + return o2::analysis::em::LeptonOrigin::Other; // decay, but not from your target meson list + } + default: + return o2::analysis::em::LeptonOrigin::Other; + } +} + } // namespace o2::aod::pwgem::photonmeson::utils::mcutil -//_______________________________________________________________________ -//_______________________________________________________________________ + #endif // PWGEM_PHOTONMESON_UTILS_MCUTILITIES_H_ diff --git a/PWGEM/PhotonMeson/Utils/ParticleOrigin.h b/PWGEM/PhotonMeson/Utils/ParticleOrigin.h new file mode 100644 index 00000000000..86e0c3b1470 --- /dev/null +++ b/PWGEM/PhotonMeson/Utils/ParticleOrigin.h @@ -0,0 +1,46 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file ParticleOrigin.h +/// \brief commonly used enums for particle origins. +/// \author marvin.hemmer@cern.ch + +#ifndef PWGEM_PHOTONMESON_UTILS_PARTICLEORIGIN_H_ +#define PWGEM_PHOTONMESON_UTILS_PARTICLEORIGIN_H_ + +#include + +namespace o2::analysis::em +{ + +// Classifies the production history of a lepton +enum class LeptonOrigin : uint8_t { + Conversion = 0, + Compton = 1, + PhotoElectric = 2, + DeltaRay = 3, + DirectMesonDecay = 4, + Other = 5 +}; + +// Classifies the production history of a photon +enum class PhotonOrigin : uint8_t { + Direct = 0, // direct photons from quarks or gluons from inital scattering + Decay = 1, // decay photons from + Bremsstrahlung = 2, // photon from bremsstrahlung + Annihilation = 3, // e+e- -> gammagamma + Hadronic = 4, // hadronic interactions like charged pions with material + Other = 5 // anything else +}; + +} // namespace o2::analysis::em + +#endif // PWGEM_PHOTONMESON_UTILS_PARTICLEORIGIN_H_