From bf1f2cf3e25e1b56d43646e9af6c6d65fb706e7a Mon Sep 17 00:00:00 2001 From: Florian Jonas Date: Thu, 9 Jul 2026 10:26:04 +0200 Subject: [PATCH 1/3] [PWGJE] Enhance GammaJetAnalysisTree with additional MC information --- PWGJE/DataModel/GammaJetAnalysisTree.h | 21 +++- PWGJE/Tasks/gammaJetTreeProducer.cxx | 138 ++++++++++++++++++++----- 2 files changed, 133 insertions(+), 26 deletions(-) diff --git a/PWGJE/DataModel/GammaJetAnalysisTree.h b/PWGJE/DataModel/GammaJetAnalysisTree.h index 1f527ff8137..d734a8057a5 100644 --- a/PWGJE/DataModel/GammaJetAnalysisTree.h +++ b/PWGJE/DataModel/GammaJetAnalysisTree.h @@ -106,8 +106,16 @@ namespace gjgammamcinfo { DECLARE_SOA_COLUMN(Origin, origin, uint16_t); DECLARE_SOA_COLUMN(LeadingEnergyFraction, leadingEnergyFraction, float); // fraction of energy from the leading MC particle +// further information about the particle that produced this cluster, one could also linnk to mcgen particle table +// but for now this is easier as the MCgen table is filtered and I do not want to deal with ensuring the filered table exists and +// does not filter out a mother +DECLARE_SOA_COLUMN(Eta, eta, float); +DECLARE_SOA_COLUMN(Phi, phi, float); +DECLARE_SOA_COLUMN(Energy, energy, float); +DECLARE_SOA_COLUMN(Pt, pt, float); +DECLARE_SOA_COLUMN(MCIso, mcIso, float); } // namespace gjgammamcinfo -DECLARE_SOA_TABLE(GjGammaMCInfos, "AOD", "GJGAMMAMCINFO", gjgamma::GjEventId, gjgammamcinfo::Origin, gjgammamcinfo::LeadingEnergyFraction) +DECLARE_SOA_TABLE(GjGammaMCInfos, "AOD", "GJGAMMAMCINFO", gjgamma::GjEventId, gjgammamcinfo::Origin, gjgammamcinfo::LeadingEnergyFraction, gjgammamcinfo::Eta, gjgammamcinfo::Phi, gjgammamcinfo::Energy, gjgammamcinfo::Pt, gjgammamcinfo::MCIso) // Generator level particle information from the MC collision that was matched to the reconstructed collision namespace gjmcparticle @@ -178,6 +186,17 @@ DECLARE_SOA_COLUMN(PerpConeRho, perpconerho, float); } // namespace gjmcjet DECLARE_SOA_TABLE(GjMCJets, "AOD", "GJMCJET", gjgamma::GjEventId, gjmcjet::Pt, gjmcjet::Eta, gjmcjet::Phi, gjmcjet::Radius, gjmcjet::Energy, gjmcjet::Mass, gjmcjet::Area, gjmcjet::PerpConeRho) +// MC gen level substructure +namespace gjmcjetsubstructure +{ +DECLARE_SOA_COLUMN(EnergyMother, energyMother, std::vector); //! energy of mother subjet at each splitting +DECLARE_SOA_COLUMN(PtLeading, ptLeading, std::vector); //! pt of leading subjet at each splitting +DECLARE_SOA_COLUMN(PtSubLeading, ptSubLeading, std::vector); //! pt of subleading subjet at each splitting +DECLARE_SOA_COLUMN(Theta, theta, std::vector); //! opening angle theta at each splitting +} // namespace gjmcjetsubstructure +DECLARE_SOA_TABLE(GjMCJetSubstructures, "AOD", "GJMCJETSUBSTR", gjgamma::GjEventId, gjmcjetsubstructure::EnergyMother, gjmcjetsubstructure::PtLeading, gjmcjetsubstructure::PtSubLeading, gjmcjetsubstructure::Theta) + + } // namespace o2::aod #endif // PWGJE_DATAMODEL_GAMMAJETANALYSISTREE_H_ diff --git a/PWGJE/Tasks/gammaJetTreeProducer.cxx b/PWGJE/Tasks/gammaJetTreeProducer.cxx index 65b97d96fa9..64b68692383 100644 --- a/PWGJE/Tasks/gammaJetTreeProducer.cxx +++ b/PWGJE/Tasks/gammaJetTreeProducer.cxx @@ -25,6 +25,7 @@ #include "Common/Core/RecoDecay.h" #include +#include #include #include #include @@ -69,6 +70,7 @@ using namespace o2; using namespace o2::aod; using namespace o2::framework; using namespace o2::framework::expressions; +using namespace o2::constants::physics; using emcClusters = o2::soa::Join; using emcMCClusters = o2::soa::Join; @@ -85,6 +87,7 @@ struct GammaJetTreeProducer { Produces chJetMCInfosTable; // detector level charged jets MC information Produces mcJetsTable; // gen level jets Produces jetSubstructuresTable; // jet substructure observables + Produces mcJetSubstructuresTable; // MC gen-level jet substructure observables HistogramRegistry mHistograms{"GammaJetTreeProducerHisto"}; @@ -95,7 +98,7 @@ struct GammaJetTreeProducer { // --------------- // event cuts - Configurable mVertexCut{"vertexCut", 10.0, "apply z-vertex cut with value in cm"}; + Configurable vertexCut{"vertexCut", 10.0, "apply z-vertex cut with value in cm"}; Configurable eventSelections{"eventSelections", "sel8", "choose event selection"}; Configurable triggerMasks{"triggerMasks", "", "possible JE Trigger masks: fJetChLowPt,fJetChHighPt,fTrackLowPt,fTrackHighPt,fJetD0ChLowPt,fJetD0ChHighPt,fJetLcChLowPt,fJetLcChHighPt,fEMCALReadout,fJetFullHighPt,fJetFullLowPt,fJetNeutralHighPt,fJetNeutralLowPt,fGammaVeryHighPtEMCAL,fGammaVeryHighPtDCAL,fGammaHighPtEMCAL,fGammaHighPtDCAL,fGammaLowPtEMCAL,fGammaLowPtDCAL,fGammaVeryLowPtEMCAL,fGammaVeryLowPtDCAL"}; Configurable trackSelections{"trackSelections", "globalTracks", "set track selections"}; @@ -139,6 +142,9 @@ struct GammaJetTreeProducer { std::vector ptSubLeadingVec; std::vector thetaVec; + // other constants + const float kInvalidValue = -99.0f; + // keeping track of the current collision index // to determine if we need to rebuild the kdTree // this makes loadorder more robust @@ -397,7 +403,7 @@ struct GammaJetTreeProducer { { mHistograms.fill(HIST("eventQA"), 0); - if (std::abs(collision.posZ()) > mVertexCut) { + if (std::abs(collision.posZ()) > vertexCut) { return false; } mHistograms.fill(HIST("eventQA"), 1); @@ -566,8 +572,7 @@ struct GammaJetTreeProducer { { T current = particle; int depth = 0; - bool continueTracing = true; - while (continueTracing) { + while (depth < kMaxRecursionDepth) { LOG(info) << "Level " << depth << " | PDG: " << current.pdgCode() << " | E: " << current.energy() @@ -719,7 +724,7 @@ struct GammaJetTreeProducer { // check if it has mothers that are etas const auto& mothers = particle.template mothers_as(); for (const auto& mother : mothers) { - if (mother.pdgCode() == 221) { + if (mother.pdgCode() == Pdg::kEta) { return true; } } @@ -736,7 +741,7 @@ struct GammaJetTreeProducer { // check if you find a pi0 mother or a eta mother const auto& mothers = particle.template mothers_as(); for (const auto& mother : mothers) { - if (mother.pdgCode() == PDG_t::kPi0 || mother.pdgCode() == 221) { + if (mother.pdgCode() == PDG_t::kPi0 || mother.pdgCode() == Pdg::kEta) { return false; } } @@ -927,15 +932,16 @@ struct GammaJetTreeProducer { /// \brief Gets the origin bitmap for a cluster /// \param cluster The cluster to check /// \param mcParticles The MC particles collection - /// \return A bitmap indicating the cluster's origin + /// \return A a pari containing the bitmap of the origin of the particle, as well as the MC index of the associated MC gen particle template - uint16_t getClusterOrigin(const T& cluster, U const& mcParticles) + std::pair getClusterOrigin(const T& cluster, U const& mcParticles) { uint16_t origin = 0; + int returnMCIndex = -1; auto inducerIDs = cluster.mcParticlesIds(); if (inducerIDs.size() == 0) { SETBIT(origin, static_cast(gjanalysis::ClusterOrigin::kUnknown)); - return origin; + return std::make_pair(origin, returnMCIndex); } // loop over all inducers and print their energy @@ -960,7 +966,7 @@ struct GammaJetTreeProducer { LOG(debug) << "Leading particle primary ID: " << leadingParticlePrimaryID; if (leadingParticlePrimaryID == -1) { SETBIT(origin, static_cast(gjanalysis::ClusterOrigin::kUnknown)); - return origin; + return std::make_pair(origin, returnMCIndex); } const auto& leadingParticlePrimary = mcParticles.iteratorAt(leadingParticlePrimaryID); LOG(debug) << "Leading particle primary PDG: " << leadingParticlePrimary.pdgCode(); @@ -993,17 +999,25 @@ struct GammaJetTreeProducer { LOG(debug) << "Leading particle primary is a decay photon from eta"; } + // for all the cases so far, we can use as MC index the index of the leading particle primary + returnMCIndex = leadingParticlePrimaryID; + // Do checks if a cluster is a merged pi0 decay // we classify a cluster as merged pi0 if the leading and subleading contribution to a cluster come from two photons that are part of a pi0 decay if (isMergedFromPDGDecay(cluster, mcParticles, PDG_t::kPi0)) { SETBIT(origin, static_cast(gjanalysis::ClusterOrigin::kMergedPi0)); LOG(debug) << "Cluster is a merged pi0"; + // if this is a merged pion decay, it should return the MC index of the pi0 mother (it is ensured that this is properly returned otherwise we would not be in this if statement) + returnMCIndex = getIndexMotherChain(mcParticles.iteratorAt(inducerIDs[0]), mcParticles, PDG_t::kPi0); + } - if (isMergedFromPDGDecay(cluster, mcParticles, 221)) { + if (isMergedFromPDGDecay(cluster, mcParticles, Pdg::kEta)) { SETBIT(origin, static_cast(gjanalysis::ClusterOrigin::kMergedEta)); LOG(debug) << "Cluster is a merged eta"; + returnMCIndex = getIndexMotherChain(mcParticles.iteratorAt(inducerIDs[0]), mcParticles, Pdg::kEta); } + // check if photon conversion // check that leading contribution is an electron or positron LOG(debug) << "Checking if cluster is a converted photon"; @@ -1023,6 +1037,8 @@ struct GammaJetTreeProducer { if ((daughters.iteratorAt(0).pdgCode() == PDG_t::kElectron && daughters.iteratorAt(1).pdgCode() == PDG_t::kPositron) || (daughters.iteratorAt(0).pdgCode() == PDG_t::kPositron && daughters.iteratorAt(1).pdgCode() == PDG_t::kElectron)) { SETBIT(origin, static_cast(gjanalysis::ClusterOrigin::kConvertedPhoton)); LOG(debug) << "Cluster is a converted photon"; + // if we found a converted photon, it should use for the the MC index the one of the photon + returnMCIndex = mother.globalIndex(); } } } @@ -1031,7 +1047,7 @@ struct GammaJetTreeProducer { } // display bit origin LOG(debug) << "Origin bits: " << std::bitset<16>(origin); - return origin; + return std::make_pair(origin, returnMCIndex); } // --------------------- @@ -1046,7 +1062,6 @@ struct GammaJetTreeProducer { mcCollisionsMultiRecCollisions.clear(); mcJetIndexMapping.clear(); } - PROCESS_SWITCH(GammaJetTreeProducer, processClearMaps, "process function that clears all the maps in each dataframe", true); // WARNING: This function always has to run second in the processing chain /// \brief Processes MC event matching QA @@ -1063,7 +1078,7 @@ struct GammaJetTreeProducer { int nRecCollisions = 0; mHistograms.fill(HIST("mcCollisionsWithRecCollisions"), 0); for (auto const& collision : collisions) { - if (std::abs(collision.posZ()) > mVertexCut) { + if (std::abs(collision.posZ()) > vertexCut) { continue; } if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits, true, true, rctLabel)) { @@ -1099,7 +1114,6 @@ struct GammaJetTreeProducer { mHistograms.fill(HIST("numberRecCollisionsVsPhotonPt"), nRecCollisions, particle.pt()); } } - PROCESS_SWITCH(GammaJetTreeProducer, processMCCollisionsMatching, "Process MC event matching QA", false); /// \brief Processes data events in data fill event table /// \param collision The collision to process @@ -1113,7 +1127,6 @@ struct GammaJetTreeProducer { eventsTable(collision.multFT0M(), collision.centFT0M(), collision.rho(), collision.eventSel(), collision.trackOccupancyInTimeRange(), collision.alias_raw()); collisionMapping[collision.globalIndex()] = eventsTable.lastIndex(); } - PROCESS_SWITCH(GammaJetTreeProducer, processEventData, "Process event data", true); using MCCol = o2::soa::Join; @@ -1154,7 +1167,6 @@ struct GammaJetTreeProducer { } mcEventsTable(eventsTable.lastIndex(), mcCollision.weight(), mcCollision.rho(), isMultipleAssigned); } - PROCESS_SWITCH(GammaJetTreeProducer, processEventMC, "Process MC event MC", false); // --------------------- // Processing functions can be safely added below this line @@ -1224,8 +1236,6 @@ struct GammaJetTreeProducer { mHistograms.fill(HIST("trackPt"), track.pt()); } } - PROCESS_SWITCH(GammaJetTreeProducer, processClusters, "Process EMCal clusters", true); - /// \brief Processes MC cluster information (rec level) /// \param collision The collision to process /// \param mcClusters The MC clusters to process @@ -1240,7 +1250,7 @@ struct GammaJetTreeProducer { // TODO: add weights for (const auto& mcCluster : mcClusters) { mHistograms.fill(HIST("clusterMC_E_All"), mcCluster.energy()); - uint16_t origin = getClusterOrigin(mcCluster, mcParticles); + auto [origin, mcIndex] = getClusterOrigin(mcCluster, mcParticles); float leadingEnergyFraction = mcCluster.amplitudeA()[0] / mcCluster.energy(); // Fill MC origin QA histograms if (origin & (1 << static_cast(gjanalysis::ClusterOrigin::kPhoton))) { @@ -1283,11 +1293,26 @@ struct GammaJetTreeProducer { mHistograms.fill(HIST("clusterMC_E_ConvertedPhoton"), mcCluster.energy()); mHistograms.fill(HIST("clusterMC_m02_ConvertedPhoton"), mcCluster.m02()); } + + // determine properties fo particle that produced the given cluster + float mcIsolation = -1; + float mcEnergy = kInvalidValue; + float mcPt = kInvalidValue; + float mcEta = kInvalidValue; + float mcPhi = kInvalidValue; + if (mcIndex != -1) { + const auto& mcParticle = mcParticles.iteratorAt(mcIndex); + mcIsolation = ch_iso_in_cone(mcParticle, isoR, true); + mcEnergy = mcParticle.energy(); + mcPt = mcParticle.pt(); + mcEta = mcParticle.eta(); + mcPhi = mcParticle.phi(); + } + // fill table - gammaMCInfosTable(storedColIndex, origin, leadingEnergyFraction); + gammaMCInfosTable(storedColIndex, origin, leadingEnergyFraction, mcEta, mcPhi, mcEnergy, mcPt, mcIsolation); } } - PROCESS_SWITCH(GammaJetTreeProducer, processClustersMCInfo, "Process MC cluster information", false); /// \brief Fills the charged jet table with jet information and calculates jet properties /// \param storedColIndex The stored collision index @@ -1369,6 +1394,58 @@ struct GammaJetTreeProducer { jetSubstructuresTable(storedColIndex, energyMotherVec, ptLeadingVec, ptSubLeadingVec, thetaVec); } + /// \brief Fills the MC gen-level substructure table for each splitting in the jet + /// \param storedColIndex The stored collision index + /// \param jet The MC particle-level jet to process + /// \param mcparticles The MC particles collection + template + void fillMCSubstructureTable(int32_t storedColIndex, T const& jet, U const& /*mcparticles*/) + { + // adjust settings according to the jet radius + jetReclusterer.jetR = jet.r() / 100.0; + // clear vectors + energyMotherVec.clear(); + ptLeadingVec.clear(); + ptSubLeadingVec.clear(); + thetaVec.clear(); + jetReclustered.clear(); + jetConstituents.clear(); + + if (jet.pt() < jetPtMin) { + return; + } + for (const auto& jetConstituent : jet.template tracks_as()) { + fastjetutilities::fillTracks(jetConstituent, jetConstituents, jetConstituent.globalIndex()); + } + + fastjet::ClusterSequenceArea clusterSeq(jetReclusterer.findJets(jetConstituents, jetReclustered)); + jetReclustered = sorted_by_pt(jetReclustered); + + if (jetReclustered.size() == 0) { + return; + } + + fastjet::PseudoJet daughterSubJet = jetReclustered[0]; + fastjet::PseudoJet parentSubJet1; + fastjet::PseudoJet parentSubJet2; + + while (daughterSubJet.has_parents(parentSubJet1, parentSubJet2)) { + if (parentSubJet1.perp() < parentSubJet2.perp()) { + std::swap(parentSubJet1, parentSubJet2); + } + + energyMotherVec.push_back(daughterSubJet.e()); + ptLeadingVec.push_back(parentSubJet1.pt()); + ptSubLeadingVec.push_back(parentSubJet2.pt()); + thetaVec.push_back(parentSubJet1.delta_R(parentSubJet2)); + + // Continue with harder parent + daughterSubJet = parentSubJet1; + } + + mcJetSubstructuresTable(storedColIndex, energyMotherVec, ptLeadingVec, ptSubLeadingVec, thetaVec); + } + Filter jetCuts = aod::jet::pt > jetPtMin; /// \brief Processes charged jets and fills jet table /// \param collision The collision to process @@ -1395,7 +1472,6 @@ struct GammaJetTreeProducer { } } } - PROCESS_SWITCH(GammaJetTreeProducer, processChargedJetsData, "Process charged jets", true); Preslice ParticlesPerMCCollisions = aod::jmcparticle::mcCollisionId; /// \brief Processes MC particles and fills MC particle table @@ -1482,7 +1558,6 @@ struct GammaJetTreeProducer { } } } - PROCESS_SWITCH(GammaJetTreeProducer, processMCParticles, "Process MC particles", false); // NOTE: The KD tree is now built lazily in each function that needs it, so execution order is no longer critical Preslice PJetsPerMCCollisions = aod::jmcparticle::mcCollisionId; @@ -1515,9 +1590,12 @@ struct GammaJetTreeProducer { mcJetIndexMapping[pjet.globalIndex()] = localIndex; localIndex++; mHistograms.fill(HIST("mcpJetPt"), pjet.pt()); + + if (calculateJetSubstructure) { + fillMCSubstructureTable(storedColIndex, pjet, mcgenparticles); + } } } - PROCESS_SWITCH(GammaJetTreeProducer, processChargedJetsMCP, "Process MC particle level jets", false); // NOTE: It is important that this function runs after the processChargedJetsMCP function (where the mc jet index mapping is built) using JetMCPTable = soa::Filtered>; @@ -1564,6 +1642,16 @@ struct GammaJetTreeProducer { chJetMCInfosTable(storedColIndex, iLocalIndexGeo, iLocalIndexPt); } } + // definition of all process switches, ordering is important! + PROCESS_SWITCH(GammaJetTreeProducer, processClearMaps, "process function that clears all the maps in each dataframe", true); // always run first + PROCESS_SWITCH(GammaJetTreeProducer, processMCCollisionsMatching, "Process MC event matching QA", false); + PROCESS_SWITCH(GammaJetTreeProducer, processEventData, "Process event data", true); + PROCESS_SWITCH(GammaJetTreeProducer, processEventMC, "Process MC event MC", false); + PROCESS_SWITCH(GammaJetTreeProducer, processClusters, "Process EMCal clusters", true); + PROCESS_SWITCH(GammaJetTreeProducer, processMCParticles, "Process MC particles", false); // has to run before processClustersMCInfo + PROCESS_SWITCH(GammaJetTreeProducer, processClustersMCInfo, "Process MC cluster information", false); + PROCESS_SWITCH(GammaJetTreeProducer, processChargedJetsData, "Process charged jets", true); + PROCESS_SWITCH(GammaJetTreeProducer, processChargedJetsMCP, "Process MC particle level jets", false); PROCESS_SWITCH(GammaJetTreeProducer, processChargedJetsMCD, "Process MC detector level jets", false); }; From d671b282a06c473b43bf32f9ea28cf2a3665e3da Mon Sep 17 00:00:00 2001 From: Florian Jonas Date: Thu, 9 Jul 2026 13:36:53 +0200 Subject: [PATCH 2/3] [PWGJE] changes to satisfy linter --- PWGJE/DataModel/GammaJetAnalysisTree.h | 8 +- PWGJE/Tasks/gammaJetTreeProducer.cxx | 361 ++++++++++++------------- 2 files changed, 181 insertions(+), 188 deletions(-) diff --git a/PWGJE/DataModel/GammaJetAnalysisTree.h b/PWGJE/DataModel/GammaJetAnalysisTree.h index d734a8057a5..5fb142469a2 100644 --- a/PWGJE/DataModel/GammaJetAnalysisTree.h +++ b/PWGJE/DataModel/GammaJetAnalysisTree.h @@ -9,9 +9,8 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -/// +/// \file GammaJetAnalysisTree.h /// \brief Table definitions for gamma-jet analyses -/// /// \author Florian Jonas #ifndef PWGJE_DATAMODEL_GAMMAJETANALYSISTREE_H_ @@ -106,8 +105,8 @@ namespace gjgammamcinfo { DECLARE_SOA_COLUMN(Origin, origin, uint16_t); DECLARE_SOA_COLUMN(LeadingEnergyFraction, leadingEnergyFraction, float); // fraction of energy from the leading MC particle -// further information about the particle that produced this cluster, one could also linnk to mcgen particle table -// but for now this is easier as the MCgen table is filtered and I do not want to deal with ensuring the filered table exists and +// further information about the particle that produced this cluster, one could also linnk to mcgen particle table +// but for now this is easier as the MCgen table is filtered and I do not want to deal with ensuring the filered table exists and // does not filter out a mother DECLARE_SOA_COLUMN(Eta, eta, float); DECLARE_SOA_COLUMN(Phi, phi, float); @@ -196,7 +195,6 @@ DECLARE_SOA_COLUMN(Theta, theta, std::vector); //! opening } // namespace gjmcjetsubstructure DECLARE_SOA_TABLE(GjMCJetSubstructures, "AOD", "GJMCJETSUBSTR", gjgamma::GjEventId, gjmcjetsubstructure::EnergyMother, gjmcjetsubstructure::PtLeading, gjmcjetsubstructure::PtSubLeading, gjmcjetsubstructure::Theta) - } // namespace o2::aod #endif // PWGJE_DATAMODEL_GAMMAJETANALYSISTREE_H_ diff --git a/PWGJE/Tasks/gammaJetTreeProducer.cxx b/PWGJE/Tasks/gammaJetTreeProducer.cxx index 64b68692383..9529f2e4684 100644 --- a/PWGJE/Tasks/gammaJetTreeProducer.cxx +++ b/PWGJE/Tasks/gammaJetTreeProducer.cxx @@ -51,6 +51,7 @@ #include #include +#include #include #include #include @@ -71,27 +72,27 @@ using namespace o2::aod; using namespace o2::framework; using namespace o2::framework::expressions; using namespace o2::constants::physics; -using emcClusters = o2::soa::Join; -using emcMCClusters = o2::soa::Join; +using EmcClusters = o2::soa::Join; +using EmcMCClusters = o2::soa::Join; struct GammaJetTreeProducer { // analysis tree // charged jets // photon candidates - Produces chargedJetsTable; // detector level jets - Produces eventsTable; // rec events - Produces gammasTable; // detector level clusters - Produces mcEventsTable; // mc collisions information - Produces mcParticlesTable; // gen level particles (photons and pi0) - Produces gammaMCInfosTable; // detector level clusters MC information - Produces chJetMCInfosTable; // detector level charged jets MC information - Produces mcJetsTable; // gen level jets - Produces jetSubstructuresTable; // jet substructure observables + Produces chargedJetsTable; // detector level jets + Produces eventsTable; // rec events + Produces gammasTable; // detector level clusters + Produces mcEventsTable; // mc collisions information + Produces mcParticlesTable; // gen level particles (photons and pi0) + Produces gammaMCInfosTable; // detector level clusters MC information + Produces chJetMCInfosTable; // detector level charged jets MC information + Produces mcJetsTable; // gen level jets + Produces jetSubstructuresTable; // jet substructure observables Produces mcJetSubstructuresTable; // MC gen-level jet substructure observables HistogramRegistry mHistograms{"GammaJetTreeProducerHisto"}; - Service pdg; + Service pdg{}; // --------------- // Configureables @@ -114,7 +115,9 @@ struct GammaJetTreeProducer { int mRunNumber = 0; std::vector eventSelectionBits; int trackSelection = -1; - const int kMaxRecursionDepth = 100; + static constexpr int MaxRecursionDepth = 100; + static constexpr int PromptMaxStatus = 90; + static constexpr int MinFSRStatus = 40; std::unordered_map collisionMapping; std::unordered_map mcJetIndexMapping; // maps the global index to the index in the mc jets table (per event). This is because later we want to later construct all trees on a per event level, and we need to know at what position in the table per event this is stored @@ -143,7 +146,7 @@ struct GammaJetTreeProducer { std::vector thetaVec; // other constants - const float kInvalidValue = -99.0f; + static constexpr float InvalidValue = -99.0f; // keeping track of the current collision index // to determine if we need to rebuild the kdTree @@ -153,8 +156,8 @@ struct GammaJetTreeProducer { void init(InitContext const&) { - using o2HistType = HistType; - using o2Axis = AxisSpec; + using O2HistType = HistType; + using O2Axis = AxisSpec; eventSelectionBits = jetderiveddatautilities::initialiseEventSelectionBits(static_cast(eventSelections)); triggerMaskBits = jetderiveddatautilities::initialiseTriggerMaskBits(triggerMasks); @@ -168,86 +171,86 @@ struct GammaJetTreeProducer { // create histograms LOG(info) << "Creating histograms"; - const o2Axis ptAxis{100, 0, 200, "p_{T} (GeV/c)"}; - const o2Axis ptRecAxis{100, 0, 200, "p_{T}^{rec} (GeV/c)"}; - const o2Axis ptGenAxis{100, 0, 200, "p_{T}^{gen} (GeV/c)"}; - const o2Axis energyAxis{100, 0, 100, "E (GeV)"}; - const o2Axis m02Axis{100, 0, 3, "m02"}; - const o2Axis etaAxis{100, -1, 1, "#eta"}; - const o2Axis phiAxis{100, 0, o2::constants::math::TwoPI, "#phi"}; - const o2Axis dRAxis{100, 0, 1, "dR"}; - const o2Axis occupancyAxis{300, 0, 30000, "occupancy"}; - const o2Axis nCollisionsAxis{10, -0.5, 9.5, "nCollisions"}; - mHistograms.add("clusterE", "Energy of cluster", o2HistType::kTH1F, {energyAxis}); - mHistograms.add("trackPt", "pT of track", o2HistType::kTH1F, {ptAxis}); - mHistograms.add("chjetPt", "pT of charged jet", o2HistType::kTH1F, {ptAxis}); - mHistograms.add("chjetPtEtaPhi", "pT of charged jet", o2HistType::kTHnSparseF, {ptAxis, etaAxis, phiAxis}); - mHistograms.add("chjetpt_vs_constpt", "pT of charged jet vs pT of constituents", o2HistType::kTH2F, {ptRecAxis, ptGenAxis}); + const O2Axis ptAxis{100, 0, 200, "p_{T} (GeV/c)"}; + const O2Axis ptRecAxis{100, 0, 200, "p_{T}^{rec} (GeV/c)"}; + const O2Axis ptGenAxis{100, 0, 200, "p_{T}^{gen} (GeV/c)"}; + const O2Axis energyAxis{100, 0, 100, "E (GeV)"}; + const O2Axis m02Axis{100, 0, 3, "m02"}; + const O2Axis etaAxis{100, -1, 1, "#eta"}; + const O2Axis phiAxis{100, 0, o2::constants::math::TwoPI, "#phi"}; + const O2Axis dRAxis{100, 0, 1, "dR"}; + const O2Axis occupancyAxis{300, 0, 30000, "occupancy"}; + const O2Axis nCollisionsAxis{10, -0.5, 9.5, "nCollisions"}; + mHistograms.add("clusterE", "Energy of cluster", O2HistType::kTH1F, {energyAxis}); + mHistograms.add("trackPt", "pT of track", O2HistType::kTH1F, {ptAxis}); + mHistograms.add("chjetPt", "pT of charged jet", O2HistType::kTH1F, {ptAxis}); + mHistograms.add("chjetPtEtaPhi", "pT of charged jet", O2HistType::kTHnSparseF, {ptAxis, etaAxis, phiAxis}); + mHistograms.add("chjetpt_vs_constpt", "pT of charged jet vs pT of constituents", O2HistType::kTH2F, {ptRecAxis, ptGenAxis}); // track QA THnSparse - mHistograms.add("trackPtEtaPhi", "Track QA", o2HistType::kTHnSparseF, {ptAxis, etaAxis, phiAxis}); - mHistograms.add("trackPtEtaOccupancy", "Track QA vs occupancy", o2HistType::kTHnSparseF, {ptAxis, etaAxis, occupancyAxis}); + mHistograms.add("trackPtEtaPhi", "Track QA", O2HistType::kTHnSparseF, {ptAxis, etaAxis, phiAxis}); + mHistograms.add("trackPtEtaOccupancy", "Track QA vs occupancy", O2HistType::kTHnSparseF, {ptAxis, etaAxis, occupancyAxis}); // QA for MC collisions to rec collision matching // number of reconstructed and matched collisions for each MC collision vs mc gen photon energy - mHistograms.add("numberRecCollisionsVsPhotonPt", "Number of rec collisions vs photon energy", o2HistType::kTH2F, {nCollisionsAxis, energyAxis}); + mHistograms.add("numberRecCollisionsVsPhotonPt", "Number of rec collisions vs photon energy", O2HistType::kTH2F, {nCollisionsAxis, energyAxis}); // Cluster MC histograms - mHistograms.add("clusterMC_E_All", "Cluster energy for photons", o2HistType::kTH1F, {energyAxis}); - mHistograms.add("clusterMC_E_Photon", "Cluster energy for photons", o2HistType::kTH1F, {energyAxis}); - mHistograms.add("clusterMC_E_PromptPhoton", "Cluster energy for prompt photons", o2HistType::kTH1F, {energyAxis}); - mHistograms.add("clusterMC_E_DirectPromptPhoton", "Cluster energy for direct prompt photons", o2HistType::kTH1F, {energyAxis}); - mHistograms.add("clusterMC_E_FragmentationPhoton", "Cluster energy for fragmentation photons", o2HistType::kTH1F, {energyAxis}); - mHistograms.add("clusterMC_E_DecayPhoton", "Cluster energy for decay photons", o2HistType::kTH1F, {energyAxis}); - mHistograms.add("clusterMC_E_DecayPhotonPi0", "Cluster energy for decay photons from pi0", o2HistType::kTH1F, {energyAxis}); - mHistograms.add("clusterMC_E_DecayPhotonEta", "Cluster energy for decay photons from eta", o2HistType::kTH1F, {energyAxis}); - mHistograms.add("clusterMC_E_MergedPi0", "Cluster energy for merged pi0s", o2HistType::kTH1F, {energyAxis}); - mHistograms.add("clusterMC_E_MergedEta", "Cluster energy for merged etas", o2HistType::kTH1F, {energyAxis}); - mHistograms.add("clusterMC_E_ConvertedPhoton", "Cluster energy for converted photons", o2HistType::kTH1F, {energyAxis}); - mHistograms.add("clusterMC_m02_Photon", "M02 for photons", o2HistType::kTH1F, {m02Axis}); - mHistograms.add("clusterMC_m02_PromptPhoton", "M02 for prompt photons", o2HistType::kTH1F, {m02Axis}); - mHistograms.add("clusterMC_m02_DirectPromptPhoton", "M02 for direct prompt photons", o2HistType::kTH1F, {m02Axis}); - mHistograms.add("clusterMC_m02_FragmentationPhoton", "M02 for fragmentation photons", o2HistType::kTH1F, {m02Axis}); - mHistograms.add("clusterMC_m02_DecayPhoton", "M02 for decay photons", o2HistType::kTH1F, {m02Axis}); - mHistograms.add("clusterMC_m02_DecayPhotonPi0", "M02 for decay photons from pi0", o2HistType::kTH1F, {m02Axis}); - mHistograms.add("clusterMC_m02_DecayPhotonEta", "M02 for decay photons from eta", o2HistType::kTH1F, {m02Axis}); - mHistograms.add("clusterMC_m02_MergedPi0", "M02 for merged pi0s", o2HistType::kTH1F, {m02Axis}); - mHistograms.add("clusterMC_m02_MergedEta", "M02 for merged etas", o2HistType::kTH1F, {m02Axis}); - mHistograms.add("clusterMC_m02_ConvertedPhoton", "M02 for converted photons", o2HistType::kTH1F, {m02Axis}); + mHistograms.add("clusterMC_E_All", "Cluster energy for photons", O2HistType::kTH1F, {energyAxis}); + mHistograms.add("clusterMC_E_Photon", "Cluster energy for photons", O2HistType::kTH1F, {energyAxis}); + mHistograms.add("clusterMC_E_PromptPhoton", "Cluster energy for prompt photons", O2HistType::kTH1F, {energyAxis}); + mHistograms.add("clusterMC_E_DirectPromptPhoton", "Cluster energy for direct prompt photons", O2HistType::kTH1F, {energyAxis}); + mHistograms.add("clusterMC_E_FragmentationPhoton", "Cluster energy for fragmentation photons", O2HistType::kTH1F, {energyAxis}); + mHistograms.add("clusterMC_E_DecayPhoton", "Cluster energy for decay photons", O2HistType::kTH1F, {energyAxis}); + mHistograms.add("clusterMC_E_DecayPhotonPi0", "Cluster energy for decay photons from pi0", O2HistType::kTH1F, {energyAxis}); + mHistograms.add("clusterMC_E_DecayPhotonEta", "Cluster energy for decay photons from eta", O2HistType::kTH1F, {energyAxis}); + mHistograms.add("clusterMC_E_MergedPi0", "Cluster energy for merged pi0s", O2HistType::kTH1F, {energyAxis}); + mHistograms.add("clusterMC_E_MergedEta", "Cluster energy for merged etas", O2HistType::kTH1F, {energyAxis}); + mHistograms.add("clusterMC_E_ConvertedPhoton", "Cluster energy for converted photons", O2HistType::kTH1F, {energyAxis}); + mHistograms.add("clusterMC_m02_Photon", "M02 for photons", O2HistType::kTH1F, {m02Axis}); + mHistograms.add("clusterMC_m02_PromptPhoton", "M02 for prompt photons", O2HistType::kTH1F, {m02Axis}); + mHistograms.add("clusterMC_m02_DirectPromptPhoton", "M02 for direct prompt photons", O2HistType::kTH1F, {m02Axis}); + mHistograms.add("clusterMC_m02_FragmentationPhoton", "M02 for fragmentation photons", O2HistType::kTH1F, {m02Axis}); + mHistograms.add("clusterMC_m02_DecayPhoton", "M02 for decay photons", O2HistType::kTH1F, {m02Axis}); + mHistograms.add("clusterMC_m02_DecayPhotonPi0", "M02 for decay photons from pi0", O2HistType::kTH1F, {m02Axis}); + mHistograms.add("clusterMC_m02_DecayPhotonEta", "M02 for decay photons from eta", O2HistType::kTH1F, {m02Axis}); + mHistograms.add("clusterMC_m02_MergedPi0", "M02 for merged pi0s", O2HistType::kTH1F, {m02Axis}); + mHistograms.add("clusterMC_m02_MergedEta", "M02 for merged etas", O2HistType::kTH1F, {m02Axis}); + mHistograms.add("clusterMC_m02_ConvertedPhoton", "M02 for converted photons", O2HistType::kTH1F, {m02Axis}); // MC Gen trigger particle histograms - mHistograms.add("mcGenTrigger_Eta", "eta of mc gen trigger particle", o2HistType::kTH1F, {etaAxis}); - mHistograms.add("mcGenTrigger_Phi", "phi of mc gen trigger particle", o2HistType::kTH1F, {phiAxis}); - mHistograms.add("mcGenTrigger_Pt", "pT of mc gen trigger particle", o2HistType::kTH1F, {ptAxis}); - mHistograms.add("mcGenTrigger_E", "E of mc gen trigger particle", o2HistType::kTH1F, {energyAxis}); - mHistograms.add("mcGenTrigger_E_PromptPhoton", "E of mc gen trigger prompt photon", o2HistType::kTH1F, {energyAxis}); - mHistograms.add("mcGenTrigger_E_DirectPromptPhoton", "E of mc gen trigger direct prompt photon", o2HistType::kTH1F, {energyAxis}); - mHistograms.add("mcGenTrigger_E_FragmentationPhoton", "E of mc gen trigger fragmentation photon", o2HistType::kTH1F, {energyAxis}); - mHistograms.add("mcGenTrigger_E_DecayPhoton", "E of mc gen trigger decay photon", o2HistType::kTH1F, {energyAxis}); - mHistograms.add("mcGenTrigger_E_DecayPhotonPi0", "E of mc gen trigger decay photon from pi0", o2HistType::kTH1F, {energyAxis}); - mHistograms.add("mcGenTrigger_E_DecayPhotonEta", "E of mc gen trigger decay photon from eta", o2HistType::kTH1F, {energyAxis}); - mHistograms.add("mcGenTrigger_E_DecayPhotonOther", "E of mc gen trigger decay photon from other", o2HistType::kTH1F, {energyAxis}); - mHistograms.add("mcGenTrigger_E_Pi0", "E of mc gen trigger pi0", o2HistType::kTH1F, {energyAxis}); + mHistograms.add("mcGenTrigger_Eta", "eta of mc gen trigger particle", O2HistType::kTH1F, {etaAxis}); + mHistograms.add("mcGenTrigger_Phi", "phi of mc gen trigger particle", O2HistType::kTH1F, {phiAxis}); + mHistograms.add("mcGenTrigger_Pt", "pT of mc gen trigger particle", O2HistType::kTH1F, {ptAxis}); + mHistograms.add("mcGenTrigger_E", "E of mc gen trigger particle", O2HistType::kTH1F, {energyAxis}); + mHistograms.add("mcGenTrigger_E_PromptPhoton", "E of mc gen trigger prompt photon", O2HistType::kTH1F, {energyAxis}); + mHistograms.add("mcGenTrigger_E_DirectPromptPhoton", "E of mc gen trigger direct prompt photon", O2HistType::kTH1F, {energyAxis}); + mHistograms.add("mcGenTrigger_E_FragmentationPhoton", "E of mc gen trigger fragmentation photon", O2HistType::kTH1F, {energyAxis}); + mHistograms.add("mcGenTrigger_E_DecayPhoton", "E of mc gen trigger decay photon", O2HistType::kTH1F, {energyAxis}); + mHistograms.add("mcGenTrigger_E_DecayPhotonPi0", "E of mc gen trigger decay photon from pi0", O2HistType::kTH1F, {energyAxis}); + mHistograms.add("mcGenTrigger_E_DecayPhotonEta", "E of mc gen trigger decay photon from eta", O2HistType::kTH1F, {energyAxis}); + mHistograms.add("mcGenTrigger_E_DecayPhotonOther", "E of mc gen trigger decay photon from other", O2HistType::kTH1F, {energyAxis}); + mHistograms.add("mcGenTrigger_E_Pi0", "E of mc gen trigger pi0", O2HistType::kTH1F, {energyAxis}); // MC Particle level jet histograms - mHistograms.add("mcpJetPt", "pT of mc particle level jet", o2HistType::kTH1F, {ptAxis}); + mHistograms.add("mcpJetPt", "pT of mc particle level jet", O2HistType::kTH1F, {ptAxis}); // MC Detector level jet matching jet histograms - mHistograms.add("mcdJetPtVsTrueJetPtMatchingGeo", "pT rec (x-axis) of detector level jets vs pT true (y-axis) of mc particle level jet (geo matching)", o2HistType::kTH2F, {ptRecAxis, ptGenAxis}); - mHistograms.add("mcdJetPtVsTrueJetPtMatchingPt", "pT rec (x-axis) of detector level jets vs pT true (y-axis) of mc particle level jet (pt matching)", o2HistType::kTH2F, {ptRecAxis, ptGenAxis}); + mHistograms.add("mcdJetPtVsTrueJetPtMatchingGeo", "pT rec (x-axis) of detector level jets vs pT true (y-axis) of mc particle level jet (geo matching)", O2HistType::kTH2F, {ptRecAxis, ptGenAxis}); + mHistograms.add("mcdJetPtVsTrueJetPtMatchingPt", "pT rec (x-axis) of detector level jets vs pT true (y-axis) of mc particle level jet (pt matching)", O2HistType::kTH2F, {ptRecAxis, ptGenAxis}); // Event QA histogram const int nEventBins = 8; - const TString eventLabels[nEventBins] = {"All", "AfterVertexCut", "AfterCollisionSelection", "AfterTriggerSelection", "AfterEMCALSelection", "AfterClusterESelection", "Has MC collision", "is not MB Gap"}; - mHistograms.add("eventQA", "Event QA", o2HistType::kTH1F, {{nEventBins, -0.5, 7.5}}); + const std::array eventLabels = {"All", "AfterVertexCut", "AfterCollisionSelection", "AfterTriggerSelection", "AfterEMCALSelection", "AfterClusterESelection", "Has MC collision", "is not MB Gap"}; + mHistograms.add("eventQA", "Event QA", O2HistType::kTH1F, {{nEventBins, -0.5, 7.5}}); for (int iBin = 0; iBin < nEventBins; iBin++) { mHistograms.get(HIST("eventQA"))->GetXaxis()->SetBinLabel(iBin + 1, eventLabels[iBin]); } // MC collisions QA histograms) const int nRecCollisionBins = 4; - const TString recCollisionLabels[nRecCollisionBins] = {"All", "1 Rec collision", "More than 1 rec collisions", "No rec collisions"}; - mHistograms.add("mcCollisionsWithRecCollisions", "MC collisions with rec collisions", o2HistType::kTH1F, {{nRecCollisionBins, -0.5, 3.5}}); + const std::array recCollisionLabels = {"All", "1 Rec collision", "More than 1 rec collisions", "No rec collisions"}; + mHistograms.add("mcCollisionsWithRecCollisions", "MC collisions with rec collisions", O2HistType::kTH1F, {{nRecCollisionBins, -0.5, 3.5}}); for (int iBin = 0; iBin < nRecCollisionBins; iBin++) { mHistograms.get(HIST("mcCollisionsWithRecCollisions"))->GetXaxis()->SetBinLabel(iBin + 1, recCollisionLabels[iBin]); } @@ -295,13 +298,13 @@ struct GammaJetTreeProducer { trackSourceIndex.push_back(originalIndex); if (phi <= (maxMatchingDistance + additionalMargin)) { trackEta.push_back(track.eta()); - trackPhi.push_back(phi + o2::constants::math::TwoPI); + trackPhi.push_back(phi + o2::constants::math::TwoPI); // o2-linter: disable=two-pi-add-subtract (periodic KD-tree ghost copy) trackPt.push_back(track.pt()); trackSourceIndex.push_back(originalIndex); } if (phi >= (o2::constants::math::TwoPI - (maxMatchingDistance + additionalMargin))) { trackEta.push_back(track.eta()); - trackPhi.push_back(phi - o2::constants::math::TwoPI); + trackPhi.push_back(phi - o2::constants::math::TwoPI); // o2-linter: disable=two-pi-add-subtract (periodic KD-tree ghost copy) trackPt.push_back(track.pt()); trackSourceIndex.push_back(originalIndex); } @@ -344,13 +347,13 @@ struct GammaJetTreeProducer { mcParticleSourceIndex.push_back(originalIndex); if (phi <= (maxMatchingDistance + additionalMargin)) { mcParticleEta.push_back(particle.eta()); - mcParticlePhi.push_back(phi + o2::constants::math::TwoPI); + mcParticlePhi.push_back(phi + o2::constants::math::TwoPI); // o2-linter: disable=two-pi-add-subtract (periodic KD-tree ghost copy) mcParticlePt.push_back(particle.pt()); mcParticleSourceIndex.push_back(originalIndex); } if (phi >= (o2::constants::math::TwoPI - (maxMatchingDistance + additionalMargin))) { mcParticleEta.push_back(particle.eta()); - mcParticlePhi.push_back(phi - o2::constants::math::TwoPI); + mcParticlePhi.push_back(phi - o2::constants::math::TwoPI); // o2-linter: disable=two-pi-add-subtract (periodic KD-tree ghost copy) mcParticlePt.push_back(particle.pt()); mcParticleSourceIndex.push_back(originalIndex); } @@ -471,16 +474,16 @@ struct GammaJetTreeProducer { /// \param mcGenIso Whether to use the mc gen particle tree (if false, use the track tree) /// \return The charged particle isolation template - double ch_iso_in_cone(const T& particle, float radius = 0.4, bool mcGenIso = false) + double chIsoInCone(const T& particle, float radius = 0.4, bool mcGenIso = false) { double iso = 0; - float point[2] = {particle.eta(), RecoDecay::constrainAngle(particle.phi(), 0.0)}; + std::array point = {particle.eta(), RecoDecay::constrainAngle(particle.phi(), 0.0)}; std::vector indices; std::unordered_set uniqueSourceIndices; if (!mcGenIso) { if (trackTree) { - trackTree->FindInRange(point, radius, indices); + trackTree->FindInRange(point.data(), radius, indices); iso += sumUniquePtFromIndices(indices, trackSourceIndex, trackPt, uniqueSourceIndices); } else { LOG(error) << "Track tree not found"; @@ -488,7 +491,7 @@ struct GammaJetTreeProducer { } } else { if (mcParticleTree) { - mcParticleTree->FindInRange(point, radius, indices); + mcParticleTree->FindInRange(point.data(), radius, indices); iso += sumUniquePtFromIndices(indices, mcParticleSourceIndex, mcParticlePt, uniqueSourceIndices); } else { LOG(error) << "MC particle tree not found"; @@ -504,7 +507,7 @@ struct GammaJetTreeProducer { /// \param radius The cone radius for density calculation /// \return The average charged particle density in the perpendicular cones template - double ch_perp_cone_rho(const T& object, float radius = 0.4, bool mcGenIso = false) + double chPerpConeRho(const T& object, float radius = 0.4, bool mcGenIso = false) { double ptSumLeft = 0; double ptSumRight = 0; @@ -515,8 +518,8 @@ struct GammaJetTreeProducer { float cPhiLeft = RecoDecay::constrainAngle(cPhi - o2::constants::math::PIHalf, 0.0); float cPhiRight = RecoDecay::constrainAngle(cPhi + o2::constants::math::PIHalf, 0.0); - float pointLeft[2] = {object.eta(), cPhiLeft}; - float pointRight[2] = {object.eta(), cPhiRight}; + std::array pointLeft = {object.eta(), cPhiLeft}; + std::array pointRight = {object.eta(), cPhiRight}; std::vector indicesLeft; std::vector indicesRight; @@ -525,8 +528,8 @@ struct GammaJetTreeProducer { if (!mcGenIso) { if (trackTree) { - trackTree->FindInRange(pointLeft, radius, indicesLeft); - trackTree->FindInRange(pointRight, radius, indicesRight); + trackTree->FindInRange(pointLeft.data(), radius, indicesLeft); + trackTree->FindInRange(pointRight.data(), radius, indicesRight); } else { LOG(error) << "Track tree not found"; return 0; @@ -536,8 +539,8 @@ struct GammaJetTreeProducer { ptSumRight += sumUniquePtFromIndices(indicesRight, trackSourceIndex, trackPt, uniqueSourceIndicesRight); } else { if (mcParticleTree) { - mcParticleTree->FindInRange(pointLeft, radius, indicesLeft); - mcParticleTree->FindInRange(pointRight, radius, indicesRight); + mcParticleTree->FindInRange(pointLeft.data(), radius, indicesLeft); + mcParticleTree->FindInRange(pointRight.data(), radius, indicesRight); } else { LOG(error) << "MC particle tree not found"; return 0; @@ -572,7 +575,7 @@ struct GammaJetTreeProducer { { T current = particle; int depth = 0; - while (depth < kMaxRecursionDepth) { + while (depth < MaxRecursionDepth) { LOG(info) << "Level " << depth << " | PDG: " << current.pdgCode() << " | E: " << current.energy() @@ -587,12 +590,10 @@ struct GammaJetTreeProducer { } // Handle case with potentially duplicate mothers int selectedMother = -1; - if (mothers.size() == 1) { - selectedMother = 0; - } else if (mothers.size() == 2 && - mothers[0].globalIndex() == mothers[1].globalIndex() && - mothers[0].pdgCode() == mothers[1].pdgCode() && - mothers[0].getGenStatusCode() == mothers[1].getGenStatusCode()) { + if (mothers.size() == 1 || (mothers.size() == 2 && // o2-linter: disable=magic-number (it is just counting number of mothers) + mothers[0].globalIndex() == mothers[1].globalIndex() && + mothers[0].pdgCode() == mothers[1].pdgCode() && + mothers[0].getGenStatusCode() == mothers[1].getGenStatusCode())) { selectedMother = 0; } if (selectedMother == -1) { @@ -620,7 +621,7 @@ struct GammaJetTreeProducer { int selectedMother = -1; if (mothers.size() == 1) { selectedMother = 0; - } else if (mothers.size() == 2) { + } else if (mothers.size() == 2) { // o2-linter: disable=magic-number (it is just counting number of mothers) twoMothersIdentical = (mothers[0].globalIndex() == mothers[1].globalIndex() && mothers[0].pdgCode() == mothers[1].pdgCode() && mothers[0].getGenStatusCode() == mothers[1].getGenStatusCode()); @@ -636,7 +637,7 @@ struct GammaJetTreeProducer { twoMothersIdentical = false; if (mothers.size() == 1) { selectedMother = 0; - } else if (mothers.size() == 2) { + } else if (mothers.size() == 2) { // o2-linter: disable=magic-number (it is just counting number of mothers) twoMothersIdentical = (mothers[0].globalIndex() == mothers[1].globalIndex() && mothers[0].pdgCode() == mothers[1].pdgCode() && mothers[0].getGenStatusCode() == mothers[1].getGenStatusCode()); @@ -653,10 +654,7 @@ struct GammaJetTreeProducer { /// \return true if particle is a prompt photon, false otherwise bool isPromptPhoton(const auto& particle) { - if (particle.pdgCode() == PDG_t::kGamma && particle.isPhysicalPrimary() && std::abs(particle.getGenStatusCode()) < 90) { - return true; - } - return false; + return static_cast(particle.pdgCode() == PDG_t::kGamma && particle.isPhysicalPrimary() && std::abs(particle.getGenStatusCode()) < PromptMaxStatus); } /// \brief Checks if a particle is a direct prompt photon /// \param particle The particle to check @@ -664,10 +662,10 @@ struct GammaJetTreeProducer { bool isDirectPromptPhoton(const auto& particle) { // check if particle isa prompt photon - if (particle.pdgCode() == PDG_t::kGamma && particle.isPhysicalPrimary() && std::abs(particle.getGenStatusCode()) < 90) { + if (particle.pdgCode() == PDG_t::kGamma && particle.isPhysicalPrimary() && std::abs(particle.getGenStatusCode()) < PromptMaxStatus) { // find the top carbon copy auto topCopy = iTopCopy(particle); - if (topCopy.pdgCode() == PDG_t::kGamma && std::abs(topCopy.getGenStatusCode()) < 40) { // < 40 is particle directly produced in hard scattering + if (topCopy.pdgCode() == PDG_t::kGamma && std::abs(topCopy.getGenStatusCode()) < MinFSRStatus) { // < 40 is particle directly produced in hard scattering return true; } } @@ -678,10 +676,10 @@ struct GammaJetTreeProducer { /// \return true if particle is a fragmentation photon, false otherwise bool isFragmentationPhoton(const auto& particle) { - if (particle.pdgCode() == PDG_t::kGamma && particle.isPhysicalPrimary() && std::abs(particle.getGenStatusCode()) < 90) { + if (particle.pdgCode() == PDG_t::kGamma && particle.isPhysicalPrimary() && std::abs(particle.getGenStatusCode()) < PromptMaxStatus) { // find the top carbon copy auto topCopy = iTopCopy(particle); - if (topCopy.pdgCode() == PDG_t::kGamma && std::abs(topCopy.getGenStatusCode()) >= 40) { // frag photon + if (topCopy.pdgCode() == PDG_t::kGamma && std::abs(topCopy.getGenStatusCode()) >= MinFSRStatus) { // frag photon return true; } } @@ -692,10 +690,7 @@ struct GammaJetTreeProducer { /// \return true if particle is a decay photon, false otherwise bool isDecayPhoton(const auto& particle) { - if (particle.pdgCode() == PDG_t::kGamma && particle.isPhysicalPrimary() && std::abs(particle.getGenStatusCode()) >= 90) { - return true; - } - return false; + return static_cast(particle.pdgCode() == PDG_t::kGamma && particle.isPhysicalPrimary() && std::abs(particle.getGenStatusCode()) >= PromptMaxStatus); } /// \brief Checks if a particle is a decay photon from pi0 /// \param particle The particle to check @@ -703,7 +698,7 @@ struct GammaJetTreeProducer { template bool isDecayPhotonPi0(const T& particle) { - if (particle.pdgCode() == PDG_t::kGamma && particle.isPhysicalPrimary() && std::abs(particle.getGenStatusCode()) >= 90) { + if (particle.pdgCode() == PDG_t::kGamma && particle.isPhysicalPrimary() && std::abs(particle.getGenStatusCode()) >= PromptMaxStatus) { // check if it has mothers that are pi0s const auto& mothers = particle.template mothers_as(); for (const auto& mother : mothers) { @@ -720,7 +715,7 @@ struct GammaJetTreeProducer { template bool isDecayPhotonEta(const T& particle) { - if (particle.pdgCode() == PDG_t::kGamma && particle.isPhysicalPrimary() && std::abs(particle.getGenStatusCode()) >= 90) { + if (particle.pdgCode() == PDG_t::kGamma && particle.isPhysicalPrimary() && std::abs(particle.getGenStatusCode()) >= PromptMaxStatus) { // check if it has mothers that are etas const auto& mothers = particle.template mothers_as(); for (const auto& mother : mothers) { @@ -737,7 +732,7 @@ struct GammaJetTreeProducer { template bool isDecayPhotonOther(const T& particle) { - if (particle.pdgCode() == PDG_t::kGamma && particle.isPhysicalPrimary() && std::abs(particle.getGenStatusCode()) >= 90) { + if (particle.pdgCode() == PDG_t::kGamma && particle.isPhysicalPrimary() && std::abs(particle.getGenStatusCode()) >= PromptMaxStatus) { // check if you find a pi0 mother or a eta mother const auto& mothers = particle.template mothers_as(); for (const auto& mother : mothers) { @@ -754,10 +749,7 @@ struct GammaJetTreeProducer { /// \return true if particle is a pi0, false otherwise bool isPi0(const auto& particle) { - if (particle.pdgCode() == PDG_t::kPi0) { - return true; - } - return false; + return static_cast(particle.pdgCode() == PDG_t::kPi0); } /// \brief Gets the bitmap for a MC particle that indicated what type of particle it is @@ -802,16 +794,15 @@ struct GammaJetTreeProducer { int getIndexMotherChain(const T& particle, aod::JMcParticles const& mcParticles, int pdgCode, int depth = 0) { // Limit recursion depth to avoid infinite loops - if (depth > kMaxRecursionDepth) { // 100 generations should be more than enough + if (depth > MaxRecursionDepth) { // 100 generations should be more than enough return -1; } const auto& mothers = particle.template mothers_as(); for (const auto& mother : mothers) { if (mother.pdgCode() == pdgCode) { return mother.globalIndex(); - } else { - return getIndexMotherChain(mother, mcParticles, pdgCode, depth + 1); } + return getIndexMotherChain(mother, mcParticles, pdgCode, depth + 1); } return -1; } @@ -823,7 +814,7 @@ struct GammaJetTreeProducer { void getDaughtersInChain(const T& particle, std::vector& daughters, int depth = 0) { // Limit recursion depth to avoid infinite loops - if (depth > kMaxRecursionDepth) { // 100 generations should be more than enough + if (depth > MaxRecursionDepth) { // 100 generations should be more than enough return; } @@ -836,7 +827,6 @@ struct GammaJetTreeProducer { daughters.push_back(daughter.globalIndex()); getDaughtersInChain(daughter, daughters, depth + 1); } - return; } /// \brief Finds the first physical primary particle in the decay chain (upwards) /// \param particle The particle to start from @@ -845,7 +835,7 @@ struct GammaJetTreeProducer { int findPhysicalPrimaryInChain(const T& particle, int depth = 0) { // Limit recursion depth to avoid infinite loops - if (depth > kMaxRecursionDepth) { // 100 generations should be more than enough + if (depth > MaxRecursionDepth) { // 100 generations should be more than enough return -1; } @@ -855,13 +845,15 @@ struct GammaJetTreeProducer { } // check if the particle has mothers - if (!particle.has_mothers()) + if (!particle.has_mothers()) { return -1; + } // now get mothers const auto mothers = particle.template mothers_as(); - if (mothers.size() == 0) + if (mothers.size() == 0) { return -1; + } // get first mother for (const auto& mother : mothers) { @@ -884,7 +876,7 @@ struct GammaJetTreeProducer { bool isMergedFromPDGDecay(const T& cluster, U const& mcParticles, int pdgCode) { auto inducerIDs = cluster.mcParticlesIds(); - if (inducerIDs.size() < 2) { // it can not me "merged" if it has less than 2 inducers + if (inducerIDs.size() < 2) { // o2-linter: disable=magic-number (it is just counting number of inducers) return false; } @@ -896,7 +888,7 @@ struct GammaJetTreeProducer { // get daughters of pi0 mother auto daughtersMother = mother.template daughters_as(); // check if there are two daughters that are both photons - if (daughtersMother.size() == 2) { + if (daughtersMother.size() == 2) { // o2-linter: disable=magic-number (it is just counting number of daughters) const auto& daughter1 = daughtersMother.iteratorAt(0); const auto& daughter2 = daughtersMother.iteratorAt(1); if (daughter1.pdgCode() == PDG_t::kGamma && daughter2.pdgCode() == PDG_t::kGamma) { @@ -1009,7 +1001,6 @@ struct GammaJetTreeProducer { LOG(debug) << "Cluster is a merged pi0"; // if this is a merged pion decay, it should return the MC index of the pi0 mother (it is ensured that this is properly returned otherwise we would not be in this if statement) returnMCIndex = getIndexMotherChain(mcParticles.iteratorAt(inducerIDs[0]), mcParticles, PDG_t::kPi0); - } if (isMergedFromPDGDecay(cluster, mcParticles, Pdg::kEta)) { SETBIT(origin, static_cast(gjanalysis::ClusterOrigin::kMergedEta)); @@ -1017,7 +1008,6 @@ struct GammaJetTreeProducer { returnMCIndex = getIndexMotherChain(mcParticles.iteratorAt(inducerIDs[0]), mcParticles, Pdg::kEta); } - // check if photon conversion // check that leading contribution is an electron or positron LOG(debug) << "Checking if cluster is a converted photon"; @@ -1032,7 +1022,7 @@ struct GammaJetTreeProducer { LOG(debug) << "Got the mother with PDG 22 and daughters"; const auto& daughters = mother.template daughters_as(); // check that mother has exactly two daughters which are e+ and e- - if (daughters.size() == 2) { + if (daughters.size() == 2) { // o2-linter: disable=magic-number (it is just counting number of daughters) LOG(debug) << "Got the daughters"; if ((daughters.iteratorAt(0).pdgCode() == PDG_t::kElectron && daughters.iteratorAt(1).pdgCode() == PDG_t::kPositron) || (daughters.iteratorAt(0).pdgCode() == PDG_t::kPositron && daughters.iteratorAt(1).pdgCode() == PDG_t::kElectron)) { SETBIT(origin, static_cast(gjanalysis::ClusterOrigin::kConvertedPhoton)); @@ -1118,7 +1108,7 @@ struct GammaJetTreeProducer { /// \brief Processes data events in data fill event table /// \param collision The collision to process /// \param clusters The EMCAL clusters in the event - void processEventData(soa::Join::iterator const& collision, emcClusters const& clusters) + void processEventData(soa::Join::iterator const& collision, EmcClusters const& clusters) { if (!isEventAccepted(collision, clusters)) { return; @@ -1134,7 +1124,7 @@ struct GammaJetTreeProducer { /// \param collision The collision to process /// \param clusters The EMCAL clusters in the event /// \param mcCollisions The MC collisions collection - void processEventMC(soa::Join::iterator const& collision, emcClusters const& clusters, MCCol const&) + void processEventMC(soa::Join::iterator const& collision, EmcClusters const& clusters, MCCol const&) { if (!isEventAccepted(collision, clusters)) { return; @@ -1175,19 +1165,20 @@ struct GammaJetTreeProducer { // define cluster filter. It selects only those clusters which are of the type // sadly passing of the string at runtime is not possible for technical region so cluster definition is // an integer instead - PresliceUnsorted EMCTrackPerTrack = aod::jemctrack::trackId; + PresliceUnsorted emcTrackPerTrack = aod::jemctrack::trackId; // Process clusters /// \brief Processes clusters and fills cluster table /// \param collision The collision to process /// \param clusters The EMCAL clusters to process /// \param tracks The tracks collection /// \param emctracks The EMCAL tracks collection from track matching - void processClusters(soa::Join::iterator const& collision, emcClusters const& clusters, aod::JetTracks const& tracks, aod::JEMCTracks const& emctracks) + void processClusters(soa::Join::iterator const& collision, EmcClusters const& clusters, aod::JetTracks const& tracks, aod::JEMCTracks const& emctracks) { // event selection int32_t storedColIndex = getStoredColIndex(collision); - if (storedColIndex == -1) + if (storedColIndex == -1) { return; + } // loop over tracks one time for QA runTrackQA(collision, tracks); @@ -1201,8 +1192,8 @@ struct GammaJetTreeProducer { // fill histograms mHistograms.fill(HIST("clusterE"), cluster.energy()); - double isoraw = ch_iso_in_cone(cluster, isoR, false); - double perpconerho = ch_perp_cone_rho(cluster, isoR, false); + double isoraw = chIsoInCone(cluster, isoR, false); + double perpconerho = chPerpConeRho(cluster, isoR, false); // find closest matched track double dEta = 0; @@ -1216,17 +1207,16 @@ struct GammaJetTreeProducer { if (!isTrackSelected(track)) { continue; } - auto emcTracksPerTrack = emctracks.sliceBy(EMCTrackPerTrack, track.globalIndex()); + auto emcTracksPerTrack = emctracks.sliceBy(emcTrackPerTrack, track.globalIndex()); auto emcTrack = emcTracksPerTrack.iteratorAt(0); // find closest track that still has E/p < trackMatchingEoverP if (cluster.energy() / track.p() > trackMatchingEoverP) { continue; - } else { - dEta = cluster.eta() - emcTrack.etaEmcal(); - dPhi = RecoDecay::constrainAngle(RecoDecay::constrainAngle(emcTrack.phiEmcal(), -o2::constants::math::PI) - RecoDecay::constrainAngle(cluster.phi(), -o2::constants::math::PI), -o2::constants::math::PI); - p = track.p(); - break; } + dEta = cluster.eta() - emcTrack.etaEmcal(); + dPhi = RecoDecay::constrainAngle(RecoDecay::constrainAngle(emcTrack.phiEmcal(), -o2::constants::math::PI) - RecoDecay::constrainAngle(cluster.phi(), -o2::constants::math::PI), -o2::constants::math::PI); + p = track.p(); + break; } gammasTable(storedColIndex, cluster.energy(), cluster.definition(), cluster.eta(), cluster.phi(), cluster.m02(), cluster.m20(), cluster.nCells(), cluster.time(), cluster.isExotic(), cluster.distanceToBadChannel(), cluster.nlm(), isoraw, perpconerho, dPhi, dEta, p); } @@ -1240,12 +1230,13 @@ struct GammaJetTreeProducer { /// \param collision The collision to process /// \param mcClusters The MC clusters to process /// \param mcParticles The MC particles collection - void processClustersMCInfo(soa::Join::iterator const& collision, emcMCClusters const& mcClusters, aod::JMcParticles const& mcParticles) + void processClustersMCInfo(soa::Join::iterator const& collision, EmcMCClusters const& mcClusters, aod::JMcParticles const& mcParticles) { // event selection int32_t storedColIndex = getStoredColIndex(collision); - if (storedColIndex == -1) + if (storedColIndex == -1) { return; + } // loop over mcClusters // TODO: add weights for (const auto& mcCluster : mcClusters) { @@ -1253,56 +1244,56 @@ struct GammaJetTreeProducer { auto [origin, mcIndex] = getClusterOrigin(mcCluster, mcParticles); float leadingEnergyFraction = mcCluster.amplitudeA()[0] / mcCluster.energy(); // Fill MC origin QA histograms - if (origin & (1 << static_cast(gjanalysis::ClusterOrigin::kPhoton))) { + if ((origin & (1 << static_cast(gjanalysis::ClusterOrigin::kPhoton))) != 0) { mHistograms.fill(HIST("clusterMC_E_Photon"), mcCluster.energy()); mHistograms.fill(HIST("clusterMC_m02_Photon"), mcCluster.m02()); } - if (origin & (1 << static_cast(gjanalysis::ClusterOrigin::kPromptPhoton))) { + if ((origin & (1 << static_cast(gjanalysis::ClusterOrigin::kPromptPhoton))) != 0) { mHistograms.fill(HIST("clusterMC_E_PromptPhoton"), mcCluster.energy()); mHistograms.fill(HIST("clusterMC_m02_PromptPhoton"), mcCluster.m02()); } - if (origin & (1 << static_cast(gjanalysis::ClusterOrigin::kDirectPromptPhoton))) { + if ((origin & (1 << static_cast(gjanalysis::ClusterOrigin::kDirectPromptPhoton))) != 0) { mHistograms.fill(HIST("clusterMC_E_DirectPromptPhoton"), mcCluster.energy()); mHistograms.fill(HIST("clusterMC_m02_DirectPromptPhoton"), mcCluster.m02()); } - if (origin & (1 << static_cast(gjanalysis::ClusterOrigin::kFragmentationPhoton))) { + if ((origin & (1 << static_cast(gjanalysis::ClusterOrigin::kFragmentationPhoton))) != 0) { mHistograms.fill(HIST("clusterMC_E_FragmentationPhoton"), mcCluster.energy()); mHistograms.fill(HIST("clusterMC_m02_FragmentationPhoton"), mcCluster.m02()); } - if (origin & (1 << static_cast(gjanalysis::ClusterOrigin::kDecayPhoton))) { + if ((origin & (1 << static_cast(gjanalysis::ClusterOrigin::kDecayPhoton))) != 0) { mHistograms.fill(HIST("clusterMC_E_DecayPhoton"), mcCluster.energy()); mHistograms.fill(HIST("clusterMC_m02_DecayPhoton"), mcCluster.m02()); } - if (origin & (1 << static_cast(gjanalysis::ClusterOrigin::kDecayPhotonPi0))) { + if ((origin & (1 << static_cast(gjanalysis::ClusterOrigin::kDecayPhotonPi0))) != 0) { mHistograms.fill(HIST("clusterMC_E_DecayPhotonPi0"), mcCluster.energy()); mHistograms.fill(HIST("clusterMC_m02_DecayPhotonPi0"), mcCluster.m02()); } - if (origin & (1 << static_cast(gjanalysis::ClusterOrigin::kDecayPhotonEta))) { + if ((origin & (1 << static_cast(gjanalysis::ClusterOrigin::kDecayPhotonEta))) != 0) { mHistograms.fill(HIST("clusterMC_E_DecayPhotonEta"), mcCluster.energy()); mHistograms.fill(HIST("clusterMC_m02_DecayPhotonEta"), mcCluster.m02()); } - if (origin & (1 << static_cast(gjanalysis::ClusterOrigin::kMergedPi0))) { + if ((origin & (1 << static_cast(gjanalysis::ClusterOrigin::kMergedPi0))) != 0) { mHistograms.fill(HIST("clusterMC_E_MergedPi0"), mcCluster.energy()); mHistograms.fill(HIST("clusterMC_m02_MergedPi0"), mcCluster.m02()); } - if (origin & (1 << static_cast(gjanalysis::ClusterOrigin::kMergedEta))) { + if ((origin & (1 << static_cast(gjanalysis::ClusterOrigin::kMergedEta))) != 0) { mHistograms.fill(HIST("clusterMC_E_MergedEta"), mcCluster.energy()); mHistograms.fill(HIST("clusterMC_m02_MergedEta"), mcCluster.m02()); } - if (origin & (1 << static_cast(gjanalysis::ClusterOrigin::kConvertedPhoton))) { + if ((origin & (1 << static_cast(gjanalysis::ClusterOrigin::kConvertedPhoton))) != 0) { mHistograms.fill(HIST("clusterMC_E_ConvertedPhoton"), mcCluster.energy()); mHistograms.fill(HIST("clusterMC_m02_ConvertedPhoton"), mcCluster.m02()); } // determine properties fo particle that produced the given cluster float mcIsolation = -1; - float mcEnergy = kInvalidValue; - float mcPt = kInvalidValue; - float mcEta = kInvalidValue; - float mcPhi = kInvalidValue; + float mcEnergy = InvalidValue; + float mcPt = InvalidValue; + float mcEta = InvalidValue; + float mcPhi = InvalidValue; if (mcIndex != -1) { const auto& mcParticle = mcParticles.iteratorAt(mcIndex); - mcIsolation = ch_iso_in_cone(mcParticle, isoR, true); + mcIsolation = chIsoInCone(mcParticle, isoR, true); mcEnergy = mcParticle.energy(); mcPt = mcParticle.pt(); mcEta = mcParticle.eta(); @@ -1333,7 +1324,7 @@ struct GammaJetTreeProducer { leadingTrackPt = constituent.pt(); } } - double perpconerho = ch_perp_cone_rho(jet, perpConeJetR, false); + double perpconerho = chPerpConeRho(jet, perpConeJetR, false); chargedJetsTable(storedColIndex, jet.pt(), jet.eta(), jet.phi(), jet.r(), jet.energy(), jet.mass(), jet.area(), leadingTrackPt, perpconerho, nconst); mHistograms.fill(HIST("chjetPtEtaPhi"), jet.pt(), jet.eta(), jet.phi()); mHistograms.fill(HIST("chjetPt"), jet.pt()); @@ -1455,8 +1446,9 @@ struct GammaJetTreeProducer { { // event selection int32_t storedColIndex = getStoredColIndex(collision); - if (storedColIndex == -1) + if (storedColIndex == -1) { return; + } // build kd tree for tracks (needed for perpendicular cone rho calculation) buildKdTree(collision, tracks); @@ -1473,7 +1465,7 @@ struct GammaJetTreeProducer { } } - Preslice ParticlesPerMCCollisions = aod::jmcparticle::mcCollisionId; + Preslice particlesPerMcCollisions = aod::jmcparticle::mcCollisionId; /// \brief Processes MC particles and fills MC particle table /// \param collision The collision to process /// \param mcgenparticles The MC particles to process @@ -1481,15 +1473,16 @@ struct GammaJetTreeProducer { { // event selection int32_t storedColIndex = getStoredColIndex(collision); - if (storedColIndex == -1) + if (storedColIndex == -1) { return; + } if (!collision.has_mcCollision()) { return; } // only storing MC particles if we found a reconstructed collision - auto particlesPerMcCollision = mcgenparticles.sliceBy(ParticlesPerMCCollisions, collision.mcCollisionId()); + auto particlesPerMcCollision = mcgenparticles.sliceBy(particlesPerMcCollisions, collision.mcCollisionId()); // build kd tree for mc particles buildKdTree(collision, particlesPerMcCollision); @@ -1504,7 +1497,7 @@ struct GammaJetTreeProducer { // Particles produced in the collision including products of strong and // electromagnetic decay and excluding feed-down from weak decays of strange // particles. - if (!(particle.isPhysicalPrimary() || particle.pdgCode() == PDG_t::kPi0)) { + if (!particle.isPhysicalPrimary() && particle.pdgCode() != PDG_t::kPi0) { continue; } @@ -1515,7 +1508,7 @@ struct GammaJetTreeProducer { // check the origin of the particle uint16_t origin = getMCParticleOrigin(particle); - double mcIsolation = ch_iso_in_cone(particle, isoR, true); + double mcIsolation = chIsoInCone(particle, isoR, true); mcParticlesTable(storedColIndex, particle.energy(), particle.eta(), particle.phi(), particle.pt(), particle.pdgCode(), mcIsolation, origin); // DEBUGGING for photons. If it is a photon, print the origin and then print the chain @@ -1532,35 +1525,35 @@ struct GammaJetTreeProducer { mHistograms.fill(HIST("mcGenTrigger_Eta"), particle.eta()); mHistograms.fill(HIST("mcGenTrigger_Phi"), particle.phi()); mHistograms.fill(HIST("mcGenTrigger_Pt"), particle.pt()); - if (origin & (1 << static_cast(gjanalysis::ParticleOrigin::kPromptPhoton))) { + if ((origin & (1 << static_cast(gjanalysis::ParticleOrigin::kPromptPhoton))) != 0) { mHistograms.fill(HIST("mcGenTrigger_E_PromptPhoton"), particle.energy()); } - if (origin & (1 << static_cast(gjanalysis::ParticleOrigin::kDirectPromptPhoton))) { + if ((origin & (1 << static_cast(gjanalysis::ParticleOrigin::kDirectPromptPhoton))) != 0) { mHistograms.fill(HIST("mcGenTrigger_E_DirectPromptPhoton"), particle.energy()); } - if (origin & (1 << static_cast(gjanalysis::ParticleOrigin::kFragmentationPhoton))) { + if ((origin & (1 << static_cast(gjanalysis::ParticleOrigin::kFragmentationPhoton))) != 0) { mHistograms.fill(HIST("mcGenTrigger_E_FragmentationPhoton"), particle.energy()); } - if (origin & (1 << static_cast(gjanalysis::ParticleOrigin::kDecayPhoton))) { + if ((origin & (1 << static_cast(gjanalysis::ParticleOrigin::kDecayPhoton))) != 0) { mHistograms.fill(HIST("mcGenTrigger_E_DecayPhoton"), particle.energy()); } - if (origin & (1 << static_cast(gjanalysis::ParticleOrigin::kDecayPhotonPi0))) { + if ((origin & (1 << static_cast(gjanalysis::ParticleOrigin::kDecayPhotonPi0))) != 0) { mHistograms.fill(HIST("mcGenTrigger_E_DecayPhotonPi0"), particle.energy()); } - if (origin & (1 << static_cast(gjanalysis::ParticleOrigin::kDecayPhotonEta))) { + if ((origin & (1 << static_cast(gjanalysis::ParticleOrigin::kDecayPhotonEta))) != 0) { mHistograms.fill(HIST("mcGenTrigger_E_DecayPhotonEta"), particle.energy()); } - if (origin & (1 << static_cast(gjanalysis::ParticleOrigin::kDecayPhotonOther))) { + if ((origin & (1 << static_cast(gjanalysis::ParticleOrigin::kDecayPhotonOther))) != 0) { mHistograms.fill(HIST("mcGenTrigger_E_DecayPhotonOther"), particle.energy()); } - if (origin & (1 << static_cast(gjanalysis::ParticleOrigin::kPi0))) { + if ((origin & (1 << static_cast(gjanalysis::ParticleOrigin::kPi0))) != 0) { mHistograms.fill(HIST("mcGenTrigger_E_Pi0"), particle.energy()); } } } // NOTE: The KD tree is now built lazily in each function that needs it, so execution order is no longer critical - Preslice PJetsPerMCCollisions = aod::jmcparticle::mcCollisionId; + Preslice pJetsPerMcCollisions = aod::jmcparticle::mcCollisionId; /// \brief Processes MC particle level charged jets and fills MC jet table /// \param collision The collision to process /// \param chargedJets The MC particle level charged jets to process @@ -1570,22 +1563,23 @@ struct GammaJetTreeProducer { { // event selection int32_t storedColIndex = getStoredColIndex(collision); - if (storedColIndex == -1) + if (storedColIndex == -1) { return; + } // loop over charged jets if (!collision.has_mcCollision()) { return; } // build kd tree for mc particles (needed for perpendicular cone rho calculation) - auto particlesPerMcCollision = mcgenparticles.sliceBy(ParticlesPerMCCollisions, collision.mcCollisionId()); + auto particlesPerMcCollision = mcgenparticles.sliceBy(particlesPerMcCollisions, collision.mcCollisionId()); buildKdTree(collision, particlesPerMcCollision); int localIndex = 0; - auto pjetsPerMcCollision = chargedJets.sliceBy(PJetsPerMCCollisions, collision.mcCollisionId()); + auto pjetsPerMcCollision = chargedJets.sliceBy(pJetsPerMcCollisions, collision.mcCollisionId()); for (const auto& pjet : pjetsPerMcCollision) { // fill MC particle level jet table - float perpconerho = ch_perp_cone_rho(pjet, perpConeJetR, true); + float perpconerho = chPerpConeRho(pjet, perpConeJetR, true); mcJetsTable(storedColIndex, pjet.pt(), pjet.eta(), pjet.phi(), pjet.r(), pjet.energy(), pjet.mass(), pjet.area(), perpconerho); mcJetIndexMapping[pjet.globalIndex()] = localIndex; localIndex++; @@ -1609,8 +1603,9 @@ struct GammaJetTreeProducer { { // event selection int32_t storedColIndex = getStoredColIndex(collision); - if (storedColIndex == -1) + if (storedColIndex == -1) { return; + } // build kd tree for tracks (needed for perpendicular cone rho calculation) buildKdTree(collision, tracks); @@ -1655,9 +1650,9 @@ struct GammaJetTreeProducer { PROCESS_SWITCH(GammaJetTreeProducer, processChargedJetsMCD, "Process MC detector level jets", false); }; -WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +WorkflowSpec defineDataProcessing(ConfigContext const& context) { WorkflowSpec workflow{ - adaptAnalysisTask(cfgc, TaskName{"gamma-jet-tree-producer"})}; + adaptAnalysisTask(context)}; return workflow; } From e82cafc5e9bd71acaf1d12d7ab2d61886404bd58 Mon Sep 17 00:00:00 2001 From: Florian Jonas Date: Thu, 9 Jul 2026 14:19:48 +0200 Subject: [PATCH 3/3] [PWGJE] introducing TESTBIT and .empty() --- PWGJE/Tasks/gammaJetTreeProducer.cxx | 46 ++++++++++++++-------------- 1 file changed, 23 insertions(+), 23 deletions(-) diff --git a/PWGJE/Tasks/gammaJetTreeProducer.cxx b/PWGJE/Tasks/gammaJetTreeProducer.cxx index 9529f2e4684..732eb1f9106 100644 --- a/PWGJE/Tasks/gammaJetTreeProducer.cxx +++ b/PWGJE/Tasks/gammaJetTreeProducer.cxx @@ -309,7 +309,7 @@ struct GammaJetTreeProducer { trackSourceIndex.push_back(originalIndex); } } - if (trackEta.size() > 0) { + if (!trackEta.empty()) { delete trackTree; trackTree = new TKDTree(trackEta.size(), 2, 1); trackTree->SetData(0, trackEta.data()); @@ -358,7 +358,7 @@ struct GammaJetTreeProducer { mcParticleSourceIndex.push_back(originalIndex); } } - if (mcParticleEta.size() > 0) { + if (!mcParticleEta.empty()) { delete mcParticleTree; mcParticleTree = new TKDTree(mcParticleEta.size(), 2, 1); mcParticleTree->SetData(0, mcParticleEta.data()); @@ -924,7 +924,7 @@ struct GammaJetTreeProducer { /// \brief Gets the origin bitmap for a cluster /// \param cluster The cluster to check /// \param mcParticles The MC particles collection - /// \return A a pari containing the bitmap of the origin of the particle, as well as the MC index of the associated MC gen particle + /// \return A pair containing the bitmap of the origin of the particle, as well as the MC index of the associated MC gen particle template std::pair getClusterOrigin(const T& cluster, U const& mcParticles) { @@ -1244,43 +1244,43 @@ struct GammaJetTreeProducer { auto [origin, mcIndex] = getClusterOrigin(mcCluster, mcParticles); float leadingEnergyFraction = mcCluster.amplitudeA()[0] / mcCluster.energy(); // Fill MC origin QA histograms - if ((origin & (1 << static_cast(gjanalysis::ClusterOrigin::kPhoton))) != 0) { + if (TESTBIT(origin, static_cast(gjanalysis::ClusterOrigin::kPhoton))) { mHistograms.fill(HIST("clusterMC_E_Photon"), mcCluster.energy()); mHistograms.fill(HIST("clusterMC_m02_Photon"), mcCluster.m02()); } - if ((origin & (1 << static_cast(gjanalysis::ClusterOrigin::kPromptPhoton))) != 0) { + if (TESTBIT(origin, static_cast(gjanalysis::ClusterOrigin::kPromptPhoton))) { mHistograms.fill(HIST("clusterMC_E_PromptPhoton"), mcCluster.energy()); mHistograms.fill(HIST("clusterMC_m02_PromptPhoton"), mcCluster.m02()); } - if ((origin & (1 << static_cast(gjanalysis::ClusterOrigin::kDirectPromptPhoton))) != 0) { + if (TESTBIT(origin, static_cast(gjanalysis::ClusterOrigin::kDirectPromptPhoton))) { mHistograms.fill(HIST("clusterMC_E_DirectPromptPhoton"), mcCluster.energy()); mHistograms.fill(HIST("clusterMC_m02_DirectPromptPhoton"), mcCluster.m02()); } - if ((origin & (1 << static_cast(gjanalysis::ClusterOrigin::kFragmentationPhoton))) != 0) { + if (TESTBIT(origin, static_cast(gjanalysis::ClusterOrigin::kFragmentationPhoton))) { mHistograms.fill(HIST("clusterMC_E_FragmentationPhoton"), mcCluster.energy()); mHistograms.fill(HIST("clusterMC_m02_FragmentationPhoton"), mcCluster.m02()); } - if ((origin & (1 << static_cast(gjanalysis::ClusterOrigin::kDecayPhoton))) != 0) { + if (TESTBIT(origin, static_cast(gjanalysis::ClusterOrigin::kDecayPhoton))) { mHistograms.fill(HIST("clusterMC_E_DecayPhoton"), mcCluster.energy()); mHistograms.fill(HIST("clusterMC_m02_DecayPhoton"), mcCluster.m02()); } - if ((origin & (1 << static_cast(gjanalysis::ClusterOrigin::kDecayPhotonPi0))) != 0) { + if (TESTBIT(origin, static_cast(gjanalysis::ClusterOrigin::kDecayPhotonPi0))) { mHistograms.fill(HIST("clusterMC_E_DecayPhotonPi0"), mcCluster.energy()); mHistograms.fill(HIST("clusterMC_m02_DecayPhotonPi0"), mcCluster.m02()); } - if ((origin & (1 << static_cast(gjanalysis::ClusterOrigin::kDecayPhotonEta))) != 0) { + if (TESTBIT(origin, static_cast(gjanalysis::ClusterOrigin::kDecayPhotonEta))) { mHistograms.fill(HIST("clusterMC_E_DecayPhotonEta"), mcCluster.energy()); mHistograms.fill(HIST("clusterMC_m02_DecayPhotonEta"), mcCluster.m02()); } - if ((origin & (1 << static_cast(gjanalysis::ClusterOrigin::kMergedPi0))) != 0) { + if (TESTBIT(origin, static_cast(gjanalysis::ClusterOrigin::kMergedPi0))) { mHistograms.fill(HIST("clusterMC_E_MergedPi0"), mcCluster.energy()); mHistograms.fill(HIST("clusterMC_m02_MergedPi0"), mcCluster.m02()); } - if ((origin & (1 << static_cast(gjanalysis::ClusterOrigin::kMergedEta))) != 0) { + if (TESTBIT(origin, static_cast(gjanalysis::ClusterOrigin::kMergedEta))) { mHistograms.fill(HIST("clusterMC_E_MergedEta"), mcCluster.energy()); mHistograms.fill(HIST("clusterMC_m02_MergedEta"), mcCluster.m02()); } - if ((origin & (1 << static_cast(gjanalysis::ClusterOrigin::kConvertedPhoton))) != 0) { + if (TESTBIT(origin, static_cast(gjanalysis::ClusterOrigin::kConvertedPhoton))) { mHistograms.fill(HIST("clusterMC_E_ConvertedPhoton"), mcCluster.energy()); mHistograms.fill(HIST("clusterMC_m02_ConvertedPhoton"), mcCluster.m02()); } @@ -1357,7 +1357,7 @@ struct GammaJetTreeProducer { fastjet::ClusterSequenceArea clusterSeq(jetReclusterer.findJets(jetConstituents, jetReclustered)); jetReclustered = sorted_by_pt(jetReclustered); - if (jetReclustered.size() == 0) { + if (jetReclustered.empty()) { return; } @@ -1412,7 +1412,7 @@ struct GammaJetTreeProducer { fastjet::ClusterSequenceArea clusterSeq(jetReclusterer.findJets(jetConstituents, jetReclustered)); jetReclustered = sorted_by_pt(jetReclustered); - if (jetReclustered.size() == 0) { + if (jetReclustered.empty()) { return; } @@ -1525,28 +1525,28 @@ struct GammaJetTreeProducer { mHistograms.fill(HIST("mcGenTrigger_Eta"), particle.eta()); mHistograms.fill(HIST("mcGenTrigger_Phi"), particle.phi()); mHistograms.fill(HIST("mcGenTrigger_Pt"), particle.pt()); - if ((origin & (1 << static_cast(gjanalysis::ParticleOrigin::kPromptPhoton))) != 0) { + if (TESTBIT(origin, static_cast(gjanalysis::ParticleOrigin::kPromptPhoton))) { mHistograms.fill(HIST("mcGenTrigger_E_PromptPhoton"), particle.energy()); } - if ((origin & (1 << static_cast(gjanalysis::ParticleOrigin::kDirectPromptPhoton))) != 0) { + if (TESTBIT(origin, static_cast(gjanalysis::ParticleOrigin::kDirectPromptPhoton))) { mHistograms.fill(HIST("mcGenTrigger_E_DirectPromptPhoton"), particle.energy()); } - if ((origin & (1 << static_cast(gjanalysis::ParticleOrigin::kFragmentationPhoton))) != 0) { + if (TESTBIT(origin, static_cast(gjanalysis::ParticleOrigin::kFragmentationPhoton))) { mHistograms.fill(HIST("mcGenTrigger_E_FragmentationPhoton"), particle.energy()); } - if ((origin & (1 << static_cast(gjanalysis::ParticleOrigin::kDecayPhoton))) != 0) { + if (TESTBIT(origin, static_cast(gjanalysis::ParticleOrigin::kDecayPhoton))) { mHistograms.fill(HIST("mcGenTrigger_E_DecayPhoton"), particle.energy()); } - if ((origin & (1 << static_cast(gjanalysis::ParticleOrigin::kDecayPhotonPi0))) != 0) { + if (TESTBIT(origin, static_cast(gjanalysis::ParticleOrigin::kDecayPhotonPi0))) { mHistograms.fill(HIST("mcGenTrigger_E_DecayPhotonPi0"), particle.energy()); } - if ((origin & (1 << static_cast(gjanalysis::ParticleOrigin::kDecayPhotonEta))) != 0) { + if (TESTBIT(origin, static_cast(gjanalysis::ParticleOrigin::kDecayPhotonEta))) { mHistograms.fill(HIST("mcGenTrigger_E_DecayPhotonEta"), particle.energy()); } - if ((origin & (1 << static_cast(gjanalysis::ParticleOrigin::kDecayPhotonOther))) != 0) { + if (TESTBIT(origin, static_cast(gjanalysis::ParticleOrigin::kDecayPhotonOther))) { mHistograms.fill(HIST("mcGenTrigger_E_DecayPhotonOther"), particle.energy()); } - if ((origin & (1 << static_cast(gjanalysis::ParticleOrigin::kPi0))) != 0) { + if (TESTBIT(origin, static_cast(gjanalysis::ParticleOrigin::kPi0))) { mHistograms.fill(HIST("mcGenTrigger_E_Pi0"), particle.energy()); } }