diff --git a/PWGJE/DataModel/GammaJetAnalysisTree.h b/PWGJE/DataModel/GammaJetAnalysisTree.h index 5b756178cc9..02a9602ec6e 100644 --- a/PWGJE/DataModel/GammaJetAnalysisTree.h +++ b/PWGJE/DataModel/GammaJetAnalysisTree.h @@ -17,43 +17,45 @@ #ifndef PWGJE_DATAMODEL_GAMMAJETANALYSISTREE_H_ #define PWGJE_DATAMODEL_GAMMAJETANALYSISTREE_H_ -#include "Framework/AnalysisDataModel.h" -#include "PWGJE/DataModel/EMCALClusters.h" #include "PWGJE/Core/JetDerivedDataUtilities.h" +#include "PWGJE/DataModel/EMCALClusters.h" #include "PWGJE/DataModel/Jet.h" -namespace o2::aod::gjanalysis{ - enum class ClusterOrigin{ - kUnknown = 0, - kPhoton, // dominant amount of energy from the cluster is from a photon - kPromptPhoton, - kDirectPromptPhoton, - kFragmentationPhoton, - kDecayPhoton, // the particle that produced the cluster is a decay product - kDecayPhotonPi0, // the cluster was produced by a pi0 decay - kDecayPhotonEta, // the cluster was produced by a eta decay - kMergedPi0, // the cluster was produced by a merged pi0, i.e. two photons contribute to the cluster that both come from pi0 decay - kMergedEta, // the cluster was produced by a merged eta, i.e. two photons contribute to the cluster that both come from eta decay - kConvertedPhoton, // the cluster was produced by a converted photon, i.e. a photon that converted to an electron-positron pair and one of the electrons was detected in the cluster - }; - enum class ParticleOrigin{ - kUnknown = 0, - kPromptPhoton, - kDirectPromptPhoton, - kFragmentationPhoton, - kDecayPhoton, - kDecayPhotonPi0, - kDecayPhotonEta, - kDecayPhotonOther, - kPi0 - }; -} +#include "Framework/AnalysisDataModel.h" + +namespace o2::aod::gjanalysis +{ +enum class ClusterOrigin { + kUnknown = 0, + kPhoton, // dominant amount of energy from the cluster is from a photon + kPromptPhoton, + kDirectPromptPhoton, + kFragmentationPhoton, + kDecayPhoton, // the particle that produced the cluster is a decay product + kDecayPhotonPi0, // the cluster was produced by a pi0 decay + kDecayPhotonEta, // the cluster was produced by a eta decay + kMergedPi0, // the cluster was produced by a merged pi0, i.e. two photons contribute to the cluster that both come from pi0 decay + kMergedEta, // the cluster was produced by a merged eta, i.e. two photons contribute to the cluster that both come from eta decay + kConvertedPhoton, // the cluster was produced by a converted photon, i.e. a photon that converted to an electron-positron pair and one of the electrons was detected in the cluster +}; +enum class ParticleOrigin { + kUnknown = 0, + kPromptPhoton, + kDirectPromptPhoton, + kFragmentationPhoton, + kDecayPhoton, + kDecayPhotonPi0, + kDecayPhotonEta, + kDecayPhotonOther, + kPi0 +}; +} // namespace o2::aod::gjanalysis namespace o2::aod { // Collision level information namespace gjevent -{ //! event index +{ //! event index DECLARE_SOA_COLUMN(Multiplicity, multiplicity, float); DECLARE_SOA_COLUMN(Centrality, centrality, float); DECLARE_SOA_COLUMN(Rho, rho, float); @@ -70,7 +72,7 @@ namespace gjmcevent { DECLARE_SOA_INDEX_COLUMN(GjEvent, gjevent); DECLARE_SOA_COLUMN(Weight, weight, double); -DECLARE_SOA_COLUMN(Rho, rho, float); // gen level rho +DECLARE_SOA_COLUMN(Rho, rho, float); // gen level rho DECLARE_SOA_COLUMN(IsMultipleAssigned, isMultipleAssigned, bool); // if the corresponding MC collision matched to this rec collision was also matched to other rec collisions (allows to skip those on analysis level ) } // namespace gjmcevent DECLARE_SOA_TABLE(GjMCEvents, "AOD", "GJMCEVENT", gjmcevent::GjEventId, gjmcevent::Weight, gjmcevent::Rho, gjmcevent::IsMultipleAssigned) @@ -106,7 +108,6 @@ DECLARE_SOA_COLUMN(LeadingEnergyFraction, leadingEnergyFraction, float); // frac } // namespace gjgammamcinfo DECLARE_SOA_TABLE(GjGammaMCInfos, "AOD", "GJGAMMAMCINFO", gjgamma::GjEventId, gjgammamcinfo::Origin, gjgammamcinfo::LeadingEnergyFraction) - // Generator level particle information from the MC collision that was matched to the reconstructed collision namespace gjmcparticle { @@ -115,10 +116,10 @@ DECLARE_SOA_COLUMN(Energy, energy, float); DECLARE_SOA_COLUMN(Eta, eta, float); DECLARE_SOA_COLUMN(Phi, phi, float); DECLARE_SOA_COLUMN(Pt, pt, float); -DECLARE_SOA_COLUMN(PdgCode, pdgCode, ushort); // TODO also add smoe origin of particle? maybe only save original pi0 and eta and photon (not decay photons) +DECLARE_SOA_COLUMN(PdgCode, pdgCode, ushort); // TODO also add smoe origin of particle? maybe only save original pi0 and eta and photon (not decay photons) DECLARE_SOA_COLUMN(MCIsolation, mcIsolation, float); // isolation in cone on mc gen level -DECLARE_SOA_COLUMN(Origin, origin, uint16_t); // origin of particle -} +DECLARE_SOA_COLUMN(Origin, origin, uint16_t); // origin of particle +} // namespace gjmcparticle DECLARE_SOA_TABLE(GjMCParticles, "AOD", "GJMCPARTICLE", gjmcparticle::GjEventId, gjmcparticle::Energy, gjmcparticle::Eta, gjmcparticle::Phi, gjmcparticle::Pt, gjmcparticle::PdgCode, gjmcparticle::MCIsolation, gjmcparticle::Origin) // Reconstructed charged jet information @@ -161,4 +162,4 @@ DECLARE_SOA_TABLE(GjMCJets, "AOD", "GJMCJET", gjgamma::GjEventId, gjmcjet::Pt, g } // namespace o2::aod -#endif // PWGJE_DATAMODEL_GAMMAJETANALYSISTREE_H_ \ No newline at end of file +#endif // PWGJE_DATAMODEL_GAMMAJETANALYSISTREE_H_ diff --git a/PWGJE/Tasks/gammaJetTreeProducer.cxx b/PWGJE/Tasks/gammaJetTreeProducer.cxx index 6bb997dfe2d..8d828ad66b1 100644 --- a/PWGJE/Tasks/gammaJetTreeProducer.cxx +++ b/PWGJE/Tasks/gammaJetTreeProducer.cxx @@ -15,41 +15,40 @@ /// \since 02.08.2024 // C++ system headers first +#include + #include #include #include -#include // Framework and other headers after -#include "Framework/ASoA.h" -#include "Framework/AnalysisDataModel.h" -#include "Framework/AnalysisTask.h" -#include "Framework/HistogramRegistry.h" -#include "Framework/O2DatabasePDGPlugin.h" +#include "PWGJE/Core/FastJetUtilities.h" +#include "PWGJE/Core/JetDerivedDataUtilities.h" +#include "PWGJE/Core/JetUtilities.h" +#include "PWGJE/DataModel/EMCALClusters.h" +#include "PWGJE/DataModel/GammaJetAnalysisTree.h" +#include "PWGJE/DataModel/Jet.h" #include "Common/Core/RecoDecay.h" #include "Common/Core/TrackSelection.h" #include "Common/Core/TrackSelectionDefaults.h" #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/TrackSelectionTables.h" +#include "EventFiltering/filterTables.h" -#include "PWGJE/Core/FastJetUtilities.h" -#include "PWGJE/Core/JetDerivedDataUtilities.h" -#include "PWGJE/Core/JetUtilities.h" -#include "PWGJE/DataModel/Jet.h" -#include "PWGJE/DataModel/GammaJetAnalysisTree.h" - -#include "EMCALBase/Geometry.h" -#include "EMCALCalib/BadChannelMap.h" -#include "PWGJE/DataModel/EMCALClusters.h" +#include "CommonDataFormat/InteractionRecord.h" +#include "DataFormatsEMCAL/AnalysisCluster.h" #include "DataFormatsEMCAL/Cell.h" #include "DataFormatsEMCAL/Constants.h" -#include "DataFormatsEMCAL/AnalysisCluster.h" -#include "TVector2.h" - -#include "CommonDataFormat/InteractionRecord.h" +#include "EMCALBase/Geometry.h" +#include "EMCALCalib/BadChannelMap.h" +#include "Framework/ASoA.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/AnalysisTask.h" +#include "Framework/HistogramRegistry.h" +#include "Framework/O2DatabasePDGPlugin.h" -#include "EventFiltering/filterTables.h" +#include "TVector2.h" // \struct GammaJetTreeProducer /// \brief Task to produce a tree for gamma-jet analysis, including photons (and information of isolation) and charged and full jets @@ -69,14 +68,14 @@ 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 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 mcJetsTable; // gen level jets HistogramRegistry mHistograms{"GammaJetTreeProducerHisto"}; @@ -152,7 +151,7 @@ struct GammaJetTreeProducer { 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 + // 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}); @@ -200,26 +199,21 @@ struct GammaJetTreeProducer { 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"}; + 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}}); - for (int iBin = 0; iBin < nEventBins; iBin++){ + 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"}; + 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}}); - for (int iBin = 0; iBin < nRecCollisionBins; iBin++){ + for (int iBin = 0; iBin < nRecCollisionBins; iBin++) { mHistograms.get(HIST("mcCollisionsWithRecCollisions"))->GetXaxis()->SetBinLabel(iBin + 1, recCollisionLabels[iBin]); } - - - - } // --------------------- @@ -228,7 +222,9 @@ struct GammaJetTreeProducer { /// \brief Builds the kd tree for the tracks or mc particles (used for fast isolation calculation) /// \param objects The objects to build the kd tree for (tracks or mc particles) - template void buildKdTree(const T& objects){ + template + void buildKdTree(const T& objects) + { trackEta.clear(); trackPhi.clear(); trackPt.clear(); @@ -238,15 +234,15 @@ struct GammaJetTreeProducer { // if the track type is aod::JetTracks, we need to build the kd tree for the tracks if constexpr (std::is_same_v, aod::JetTracks>) { - for(auto track : objects){ - if(!isTrackSelected(track)){ + for (auto track : objects) { + if (!isTrackSelected(track)) { continue; } trackEta.push_back(track.eta()); trackPhi.push_back(track.phi()); trackPt.push_back(track.pt()); } - if(trackEta.size() > 0){ + if (trackEta.size() > 0) { delete trackTree; trackTree = new TKDTree(trackEta.size(), 2, 1); trackTree->SetData(0, trackEta.data()); @@ -256,11 +252,11 @@ struct GammaJetTreeProducer { } // if the track type is aod::JetParticles, we need to build the kd tree for the mc particles if constexpr (std::is_same_v, aod::JetParticles>) { - for(auto particle : objects){ - if(!particle.isPhysicalPrimary()){ + for (auto particle : objects) { + if (!particle.isPhysicalPrimary()) { continue; } - if(!isCharged(particle)){ + if (!isCharged(particle)) { continue; } if (particle.pt() < trackMinPt) { @@ -270,7 +266,7 @@ struct GammaJetTreeProducer { mcParticlePhi.push_back(particle.phi()); mcParticlePt.push_back(particle.pt()); } - if(mcParticleEta.size() > 0){ + if (mcParticleEta.size() > 0) { delete mcParticleTree; mcParticleTree = new TKDTree(mcParticleEta.size(), 2, 1); mcParticleTree->SetData(0, mcParticleEta.data()); @@ -336,7 +332,6 @@ struct GammaJetTreeProducer { if (cluster.energy() > minClusterETrigger) { mHistograms.fill(HIST("eventQA"), 5); return true; - } } return false; @@ -345,7 +340,8 @@ struct GammaJetTreeProducer { /// \brief Checks if a particle is charged /// \param particle The MC particle to check /// \return true if particle has non-zero charge, false otherwise - bool isCharged(const auto& particle){ + bool isCharged(const auto& particle) + { return std::abs(pdg->GetParticle(particle.pdgCode())->Charge()) >= 1.; } @@ -354,12 +350,14 @@ struct GammaJetTreeProducer { /// \param radius The cone radius /// \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){ + template + double ch_iso_in_cone(const T& particle, float radius = 0.4, bool mcGenIso = false) + { double iso = 0; float point[2] = {particle.eta(), particle.phi()}; std::vector indices; - if(!mcGenIso){ + if (!mcGenIso) { if (trackTree) { trackTree->FindInRange(point, radius, indices); for (const auto& index : indices) { @@ -382,13 +380,14 @@ struct GammaJetTreeProducer { } return iso; } - + /// \brief Calculates the charged particle density in perpendicular cones /// \param object The reference object (cluster or jet) /// \param tracks The tracks to check /// \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) + template + double ch_perp_cone_rho(const T& object, float radius = 0.4, bool mcGenIso = false) { double ptSumLeft = 0; double ptSumRight = 0; @@ -405,7 +404,7 @@ struct GammaJetTreeProducer { std::vector indicesLeft; std::vector indicesRight; - if(!mcGenIso){ + if (!mcGenIso) { if (trackTree) { trackTree->FindInRange(pointLeft, radius, indicesLeft); trackTree->FindInRange(pointRight, radius, indicesRight); @@ -455,65 +454,69 @@ struct GammaJetTreeProducer { } } - /// \brief Finds the top-most copy of a particle in the decay chain (following carbon copies) /// \param particle The particle to start from /// \return The top-most copy of the particle - template T iTopCopy(const T& particle) const { + template + T iTopCopy(const T& particle) const + { int iUp = particle.globalIndex(); T currentParticle = particle; int pdgCode = particle.pdgCode(); auto mothers = particle.template mothers_as(); - while ( iUp > 0 && mothers.size() == 1 && mothers[0].globalIndex() > 0 && mothers[0].pdgCode() == pdgCode){ + while (iUp > 0 && mothers.size() == 1 && mothers[0].globalIndex() > 0 && mothers[0].pdgCode() == pdgCode) { iUp = mothers[0].globalIndex(); currentParticle = mothers[0]; mothers = currentParticle.template mothers_as(); } return currentParticle; } - /// \brief Checks if a particle is a prompt photon /// \param particle The MC particle to check /// \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; + bool isPromptPhoton(const auto& particle) + { + if (particle.pdgCode() == PDG_t::kGamma && particle.isPhysicalPrimary() && std::abs(particle.getGenStatusCode()) < 90) { + return true; } return false; } /// \brief Checks if a particle is a direct prompt photon /// \param particle The particle to check /// \return true if particle is a direct prompt photon, false otherwise - bool isDirectPromptPhoton(const auto& particle){ + bool isDirectPromptPhoton(const auto& particle) + { // check if particle isa prompt photon - if(particle.pdgCode() == PDG_t::kGamma && particle.isPhysicalPrimary() && std::abs(particle.getGenStatusCode()) < 90){ - // 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 - return true; - } + if (particle.pdgCode() == PDG_t::kGamma && particle.isPhysicalPrimary() && std::abs(particle.getGenStatusCode()) < 90) { + // 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 + return true; + } } return false; } /// \brief Checks if a particle is a fragmentation photon /// \param particle The particle to check /// \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){ - // find the top carbon copy - auto topCopy = iTopCopy(particle); - if(topCopy.pdgCode() == PDG_t::kGamma && std::abs(topCopy.getGenStatusCode()) >= 40){ // frag photon - return true; - } + bool isFragmentationPhoton(const auto& particle) + { + if (particle.pdgCode() == PDG_t::kGamma && particle.isPhysicalPrimary() && std::abs(particle.getGenStatusCode()) < 90) { + // find the top carbon copy + auto topCopy = iTopCopy(particle); + if (topCopy.pdgCode() == PDG_t::kGamma && std::abs(topCopy.getGenStatusCode()) >= 40) { // frag photon + return true; + } } return false; } /// \brief Checks if a particle is a decay photon /// \param particle The particle to check /// \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){ + bool isDecayPhoton(const auto& particle) + { + if (particle.pdgCode() == PDG_t::kGamma && particle.isPhysicalPrimary() && std::abs(particle.getGenStatusCode()) >= 90) { return true; } return false; @@ -521,12 +524,14 @@ struct GammaJetTreeProducer { /// \brief Checks if a particle is a decay photon from pi0 /// \param particle The particle to check /// \return true if particle is a decay photon from pi0, false otherwise - template bool isDecayPhotonPi0(const T& particle){ - if(particle.pdgCode() == PDG_t::kGamma && particle.isPhysicalPrimary() && std::abs(particle.getGenStatusCode()) >= 90){ + template + bool isDecayPhotonPi0(const T& particle) + { + if (particle.pdgCode() == PDG_t::kGamma && particle.isPhysicalPrimary() && std::abs(particle.getGenStatusCode()) >= 90) { // check if it has mothers that are pi0s - const auto &mothers = particle.template mothers_as(); - for(auto mother : mothers){ - if(mother.pdgCode() == PDG_t::kPi0){ + const auto& mothers = particle.template mothers_as(); + for (auto mother : mothers) { + if (mother.pdgCode() == PDG_t::kPi0) { return true; } } @@ -536,12 +541,14 @@ struct GammaJetTreeProducer { /// \brief Checks if a particle is a decay photon from eta /// \param particle The particle to check /// \return true if particle is a decay photon from eta, false otherwise - template bool isDecayPhotonEta(const T& particle){ - if(particle.pdgCode() == PDG_t::kGamma && particle.isPhysicalPrimary() && std::abs(particle.getGenStatusCode()) >= 90){ + template + bool isDecayPhotonEta(const T& particle) + { + if (particle.pdgCode() == PDG_t::kGamma && particle.isPhysicalPrimary() && std::abs(particle.getGenStatusCode()) >= 90) { // check if it has mothers that are etas - const auto &mothers = particle.template mothers_as(); - for(auto mother : mothers){ - if(mother.pdgCode() == 221){ + const auto& mothers = particle.template mothers_as(); + for (auto mother : mothers) { + if (mother.pdgCode() == 221) { return true; } } @@ -551,12 +558,14 @@ struct GammaJetTreeProducer { /// \brief Checks if a particle is a decay photon from other sources /// \param particle The particle to check /// \return true if particle is a decay photon from other sources, false otherwise - template bool isDecayPhotonOther(const T& particle){ - if(particle.pdgCode() == PDG_t::kGamma && particle.isPhysicalPrimary() && std::abs(particle.getGenStatusCode()) >= 90){ + template + bool isDecayPhotonOther(const T& particle) + { + if (particle.pdgCode() == PDG_t::kGamma && particle.isPhysicalPrimary() && std::abs(particle.getGenStatusCode()) >= 90) { // check if you find a pi0 mother or a eta mother - const auto &mothers = particle.template mothers_as(); - for(auto mother : mothers){ - if(mother.pdgCode() == PDG_t::kPi0 || mother.pdgCode() == 221){ + const auto& mothers = particle.template mothers_as(); + for (auto mother : mothers) { + if (mother.pdgCode() == PDG_t::kPi0 || mother.pdgCode() == 221) { return false; } } @@ -567,8 +576,9 @@ struct GammaJetTreeProducer { /// \brief Checks if a particle is a pi0 /// \param particle The particle to check /// \return true if particle is a pi0, false otherwise - bool isPi0(const auto& particle){ - if(particle.pdgCode() == PDG_t::kPi0){ + bool isPi0(const auto& particle) + { + if (particle.pdgCode() == PDG_t::kPi0) { return true; } return false; @@ -577,30 +587,31 @@ struct GammaJetTreeProducer { /// \brief Gets the bitmap for a MC particle that indicated what type of particle it is /// \param particle The particle to check /// \return A bitmap indicating the particle's origin - uint16_t getMCParticleOrigin(const auto& particle){ + uint16_t getMCParticleOrigin(const auto& particle) + { uint16_t origin = 0; - if(isPromptPhoton(particle)){ + if (isPromptPhoton(particle)) { SETBIT(origin, static_cast(gjanalysis::ParticleOrigin::kPromptPhoton)); } - if(isDirectPromptPhoton(particle)){ + if (isDirectPromptPhoton(particle)) { SETBIT(origin, static_cast(gjanalysis::ParticleOrigin::kDirectPromptPhoton)); } - if(isFragmentationPhoton(particle)){ + if (isFragmentationPhoton(particle)) { SETBIT(origin, static_cast(gjanalysis::ParticleOrigin::kFragmentationPhoton)); } - if(isDecayPhoton(particle)){ + if (isDecayPhoton(particle)) { SETBIT(origin, static_cast(gjanalysis::ParticleOrigin::kDecayPhoton)); } - if(isDecayPhotonPi0(particle)){ + if (isDecayPhotonPi0(particle)) { SETBIT(origin, static_cast(gjanalysis::ParticleOrigin::kDecayPhotonPi0)); } - if(isDecayPhotonEta(particle)){ + if (isDecayPhotonEta(particle)) { SETBIT(origin, static_cast(gjanalysis::ParticleOrigin::kDecayPhotonEta)); } - if(isDecayPhotonOther(particle)){ + if (isDecayPhotonOther(particle)) { SETBIT(origin, static_cast(gjanalysis::ParticleOrigin::kDecayPhotonOther)); } - if(isPi0(particle)){ + if (isPi0(particle)) { SETBIT(origin, static_cast(gjanalysis::ParticleOrigin::kPi0)); } return origin; @@ -611,14 +622,16 @@ struct GammaJetTreeProducer { /// \param mcParticles The MC particles collection /// \param pdgCode The PDG code to search for /// \return The index of the mother particle, or -1 if not found - template int getIndexMotherChain(const T& particle, aod::JMcParticles const& mcParticles, int pdgCode, int depth = 0){ + template + 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 return -1; } - const auto &mothers = particle.template mothers_as(); - for(auto mother : mothers){ - if(mother.pdgCode() == pdgCode){ + const auto& mothers = particle.template mothers_as(); + for (auto mother : mothers) { + if (mother.pdgCode() == pdgCode) { return mother.globalIndex(); } else { return getIndexMotherChain(mother, mcParticles, pdgCode, depth + 1); @@ -630,16 +643,18 @@ struct GammaJetTreeProducer { /// \brief Gets all daughter particle IDs in the decay chain /// \param particle The particle to start from /// \return Vector of daughter particle IDs - template std::vector getDaughtersInChain(const T& particle, int depth = 0){ + template + std::vector getDaughtersInChain(const T& particle, int depth = 0) + { std::vector daughters; // Limit recursion depth to avoid infinite loops if (depth > kMaxRecursionDepth) { // 100 generations should be more than enough return daughters; } T currentParticle = particle; - while(currentParticle.has_daughters()){ - const auto &daughtersIDs = currentParticle.template daughters_as(); - for(auto daughter : daughtersIDs){ + while (currentParticle.has_daughters()) { + const auto& daughtersIDs = currentParticle.template daughters_as(); + for (auto daughter : daughtersIDs) { daughters.push_back(daughter.globalIndex()); } currentParticle = daughtersIDs.iteratorAt(0); @@ -653,30 +668,32 @@ struct GammaJetTreeProducer { /// \brief Finds the first physical primary particle in the decay chain (upwards) /// \param particle The particle to start from /// \return The index of the first physical primary particle, or -1 if not found - template int findPhysicalPrimaryInChain(const T& particle, int depth = 0){ + template + 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 return -1; } // first check if current particle is physical primary - if(particle.isPhysicalPrimary()){ + if (particle.isPhysicalPrimary()) { return particle.globalIndex(); } // 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(auto mother : mothers){ + for (auto mother : mothers) { int primaryIndex = findPhysicalPrimaryInChain(mother, depth + 1); - if(primaryIndex >= 0) { + if (primaryIndex >= 0) { return primaryIndex; } break; // only check first mother @@ -690,24 +707,26 @@ struct GammaJetTreeProducer { /// \param mcParticles The MC particles collection /// \param pdgCode The PDG code to check for /// \return true if cluster is merged from the specified decay, false otherwise - template bool isMergedFromPDGDecay(const T& cluster, U const& mcParticles, int pdgCode){ + template + 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) { // it can not me "merged" if it has less than 2 inducers return false; } bool isMerged = false; int motherIndex = getIndexMotherChain(mcParticles.iteratorAt(inducerIDs[0]), mcParticles, pdgCode); - if(motherIndex != -1){ - const auto &mother = mcParticles.iteratorAt(motherIndex); + if (motherIndex != -1) { + const auto& mother = mcParticles.iteratorAt(motherIndex); // 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){ - const auto &daughter1 = daughtersMother.iteratorAt(0); - const auto &daughter2 = daughtersMother.iteratorAt(1); - if(daughter1.pdgCode() == PDG_t::kGamma && daughter2.pdgCode() == PDG_t::kGamma){ + if (daughtersMother.size() == 2) { + const auto& daughter1 = daughtersMother.iteratorAt(0); + const auto& daughter2 = daughtersMother.iteratorAt(1); + if (daughter1.pdgCode() == PDG_t::kGamma && daughter2.pdgCode() == PDG_t::kGamma) { // get the full stack of particles that these daughters create auto fullDecayChain1 = getDaughtersInChain(daughter1); auto fullDecayChain2 = getDaughtersInChain(daughter2); @@ -715,17 +734,17 @@ struct GammaJetTreeProducer { bool photon2Found = false; // check if any of the particles in the fullDecayChain are leading or subleading in the cluster - for(auto particleID : fullDecayChain1){ - if(particleID == inducerIDs[0] || particleID == inducerIDs[1]){ + for (auto particleID : fullDecayChain1) { + if (particleID == inducerIDs[0] || particleID == inducerIDs[1]) { photon1Found = true; } } - for(auto particleID : fullDecayChain2){ - if(particleID == inducerIDs[0] || particleID == inducerIDs[1]){ + for (auto particleID : fullDecayChain2) { + if (particleID == inducerIDs[0] || particleID == inducerIDs[1]) { photon2Found = true; } } - if(photon1Found && photon2Found){ + if (photon1Found && photon2Found) { isMerged = true; } } @@ -734,16 +753,17 @@ struct GammaJetTreeProducer { return isMerged; } - // determine cluster origin /// \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 - template uint16_t getClusterOrigin(const T& cluster, U const& mcParticles){ + template + uint16_t getClusterOrigin(const T& cluster, U const& mcParticles) + { uint16_t origin = 0; auto inducerIDs = cluster.mcParticlesIds(); - if(inducerIDs.size() == 0){ + if (inducerIDs.size() == 0) { SETBIT(origin, static_cast(gjanalysis::ClusterOrigin::kUnknown)); return origin; } @@ -752,10 +772,10 @@ struct GammaJetTreeProducer { LOG(debug) << "Cluster with energy: " << cluster.energy() << " and nInducers: " << inducerIDs.size(); LOG(debug) << "Number of stored amplitudes: " << cluster.amplitudeA().size(); int aCounter = 0; - for(auto inducerID : inducerIDs){ - const auto &inducer = mcParticles.iteratorAt(inducerID); + for (auto inducerID : inducerIDs) { + const auto& inducer = mcParticles.iteratorAt(inducerID); int motherPDG = -1; - if(inducer.has_mothers()){ + if (inducer.has_mothers()) { motherPDG = inducer.template mothers_as()[0].pdgCode(); } LOG(debug) << "Inducer energy: " << inducer.energy() << " amplitude: " << cluster.amplitudeA()[aCounter] << " and PDG: " << inducer.pdgCode() << " isPhysicalPrimary: " << inducer.isPhysicalPrimary() << " motherPDG: " << motherPDG; @@ -763,53 +783,53 @@ struct GammaJetTreeProducer { } // check if leading energy contribution is from a photon - const auto &leadingParticle = mcParticles.iteratorAt(inducerIDs[0]); + const auto& leadingParticle = mcParticles.iteratorAt(inducerIDs[0]); LOG(debug) << "Leading particle: PDG" << leadingParticle.pdgCode(); - //leading particle primary ID + // leading particle primary ID int leadingParticlePrimaryID = findPhysicalPrimaryInChain(leadingParticle); LOG(debug) << "Leading particle primary ID: " << leadingParticlePrimaryID; - if(leadingParticlePrimaryID == -1){ + if (leadingParticlePrimaryID == -1) { SETBIT(origin, static_cast(gjanalysis::ClusterOrigin::kUnknown)); return origin; } - const auto &leadingParticlePrimary = mcParticles.iteratorAt(leadingParticlePrimaryID); + const auto& leadingParticlePrimary = mcParticles.iteratorAt(leadingParticlePrimaryID); LOG(debug) << "Leading particle primary PDG: " << leadingParticlePrimary.pdgCode(); - if(leadingParticlePrimary.pdgCode() == PDG_t::kGamma){ + if (leadingParticlePrimary.pdgCode() == PDG_t::kGamma) { LOG(debug) << "Leading particle primary is a photon"; SETBIT(origin, static_cast(gjanalysis::ClusterOrigin::kPhoton)); } - if(isPromptPhoton(leadingParticlePrimary)){ + if (isPromptPhoton(leadingParticlePrimary)) { SETBIT(origin, static_cast(gjanalysis::ClusterOrigin::kPromptPhoton)); LOG(debug) << "Leading particle primary is a prompt photon"; } - if(isDirectPromptPhoton(leadingParticlePrimary)){ + if (isDirectPromptPhoton(leadingParticlePrimary)) { SETBIT(origin, static_cast(gjanalysis::ClusterOrigin::kDirectPromptPhoton)); LOG(debug) << "Leading particle primary is a direct prompt photon"; } - if(isFragmentationPhoton(leadingParticlePrimary)){ + if (isFragmentationPhoton(leadingParticlePrimary)) { SETBIT(origin, static_cast(gjanalysis::ClusterOrigin::kFragmentationPhoton)); LOG(debug) << "Leading particle primary is a fragmentation photon"; } - if(isDecayPhoton(leadingParticlePrimary)){ + if (isDecayPhoton(leadingParticlePrimary)) { SETBIT(origin, static_cast(gjanalysis::ClusterOrigin::kDecayPhoton)); LOG(debug) << "Leading particle primary is a decay photon"; } - if(isDecayPhotonPi0(leadingParticlePrimary)){ + if (isDecayPhotonPi0(leadingParticlePrimary)) { SETBIT(origin, static_cast(gjanalysis::ClusterOrigin::kDecayPhotonPi0)); LOG(debug) << "Leading particle primary is a decay photon from pi0"; } - if(isDecayPhotonEta(leadingParticlePrimary)){ + if (isDecayPhotonEta(leadingParticlePrimary)) { SETBIT(origin, static_cast(gjanalysis::ClusterOrigin::kDecayPhotonEta)); LOG(debug) << "Leading particle primary is a decay photon from eta"; } // 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)){ + if (isMergedFromPDGDecay(cluster, mcParticles, PDG_t::kPi0)) { SETBIT(origin, static_cast(gjanalysis::ClusterOrigin::kMergedPi0)); LOG(debug) << "Cluster is a merged pi0"; - } - if(isMergedFromPDGDecay(cluster, mcParticles, 221)){ + } + if (isMergedFromPDGDecay(cluster, mcParticles, 221)) { SETBIT(origin, static_cast(gjanalysis::ClusterOrigin::kMergedEta)); LOG(debug) << "Cluster is a merged eta"; } @@ -817,20 +837,20 @@ struct GammaJetTreeProducer { // check if photon conversion // check that leading contribution is an electron or positron LOG(debug) << "Checking if cluster is a converted photon"; - if(std::abs(leadingParticle.pdgCode()) == PDG_t::kElectron){ + if (std::abs(leadingParticle.pdgCode()) == PDG_t::kElectron) { // make sure this electron is not a physicsl primary and has mothers - if(!leadingParticle.isPhysicalPrimary() && leadingParticle.has_mothers()){ + if (!leadingParticle.isPhysicalPrimary() && leadingParticle.has_mothers()) { const auto mothers = leadingParticle.template mothers_as(); - if(mothers.size() > 0) { + if (mothers.size() > 0) { LOG(debug) << "Got the mother"; - const auto &mother = mothers[0]; - if(mother.pdgCode() == PDG_t::kGamma && mother.has_daughters()){ + const auto& mother = mothers[0]; + if (mother.pdgCode() == PDG_t::kGamma && mother.has_daughters()) { 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) { 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)){ + 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"; } @@ -844,7 +864,6 @@ struct GammaJetTreeProducer { return origin; } - // --------------------- // Processing functions // --------------------- @@ -864,15 +883,16 @@ struct GammaJetTreeProducer { /// \param mcCollision The MC collision to process /// \param collisions The rec collisions collection /// \param mcgenparticles The MC particles collection - void processMCCollisionsMatching(aod::JetMcCollision const& mcCollision, soa::SmallGroups const& collisions, aod::JetParticles const& mcgenparticles){ - if(mcCollision.weight() == 0){ + void processMCCollisionsMatching(aod::JetMcCollision const& mcCollision, soa::SmallGroups const& collisions, aod::JetParticles const& mcgenparticles) + { + if (mcCollision.weight() == 0) { return; } // determine number of rec collisions int nRecCollisions = 0; mHistograms.fill(HIST("mcCollisionsWithRecCollisions"), 0); - for(auto const& collision : collisions){ + for (auto const& collision : collisions) { if (collision.posZ() > mVertexCut) { continue; } @@ -887,33 +907,33 @@ struct GammaJetTreeProducer { } nRecCollisions++; } - if(nRecCollisions == 0){ + if (nRecCollisions == 0) { mHistograms.fill(HIST("mcCollisionsWithRecCollisions"), 3); } - if(nRecCollisions == 1){ + if (nRecCollisions == 1) { mHistograms.fill(HIST("mcCollisionsWithRecCollisions"), 1); } - if(nRecCollisions > 1){ + if (nRecCollisions > 1) { mHistograms.fill(HIST("mcCollisionsWithRecCollisions"), 2); } - if(nRecCollisions > 1){ + if (nRecCollisions > 1) { mcCollisionsMultiRecCollisions.push_back(mcCollision.globalIndex()); } // loop over mcgenparticles - for(auto const& particle : mcgenparticles){ - if(!particle.isPhysicalPrimary()){ + for (auto const& particle : mcgenparticles) { + if (!particle.isPhysicalPrimary()) { continue; - } - if(particle.pdgCode() != 22){ + } + if (particle.pdgCode() != 22) { continue; - } + } 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 + /// \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) @@ -927,32 +947,33 @@ struct GammaJetTreeProducer { } PROCESS_SWITCH(GammaJetTreeProducer, processEventData, "Process event data", true); - using MCCol = o2::soa::Join; - + using MCCol = o2::soa::Join; + /// \brief Processes MC events and fills rec and MC event tables (disable processEventData) /// \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; } // check that this event has a MC collision - if(!collision.has_mcCollision()){ + if (!collision.has_mcCollision()) { return; } mHistograms.fill(HIST("eventQA"), 6); // check if this event is not MB gap event - if(collision.subGeneratorId() == jetderiveddatautilities::JCollisionSubGeneratorId::mbGap){ + if (collision.subGeneratorId() == jetderiveddatautilities::JCollisionSubGeneratorId::mbGap) { return; } mHistograms.fill(HIST("eventQA"), 7); // fill rec collision table eventsTable(collision.multiplicity(), collision.centrality(), collision.rho(), collision.eventSel(), collision.trackOccupancyInTimeRange(), collision.alias_raw()); - + // fill collision mapping collisionMapping[collision.globalIndex()] = eventsTable.lastIndex(); @@ -960,14 +981,13 @@ struct GammaJetTreeProducer { bool isMultipleAssigned = false; // check if we are dealing with a rec collision matched to a MC collision that was matched to multiple rec collisions - if(std::find(mcCollisionsMultiRecCollisions.begin(), mcCollisionsMultiRecCollisions.end(), mcCollision.globalIndex()) != mcCollisionsMultiRecCollisions.end()){ + if (std::find(mcCollisionsMultiRecCollisions.begin(), mcCollisionsMultiRecCollisions.end(), mcCollision.globalIndex()) != mcCollisionsMultiRecCollisions.end()) { isMultipleAssigned = true; } 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 // --------------------- @@ -1042,14 +1062,15 @@ 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) return; // loop over mcClusters // TODO: add weights - for(auto mcCluster : mcClusters){ + for (auto mcCluster : mcClusters) { mHistograms.fill(HIST("clusterMC_E_All"), mcCluster.energy()); uint16_t origin = getClusterOrigin(mcCluster, mcParticles); float leadingEnergyFraction = mcCluster.amplitudeA()[0] / mcCluster.energy(); @@ -1058,39 +1079,39 @@ struct GammaJetTreeProducer { 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))) { 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))) { 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))) { 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))) { 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))) { 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))) { 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))) { 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))) { 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))) { mHistograms.fill(HIST("clusterMC_E_ConvertedPhoton"), mcCluster.energy()); mHistograms.fill(HIST("clusterMC_m02_ConvertedPhoton"), mcCluster.m02()); } @@ -1104,16 +1125,18 @@ struct GammaJetTreeProducer { /// \param storedColIndex The stored collision index /// \param jet The jet to process /// \param tracks The tracks collection - template void fillChargedJetTable(int32_t storedColIndex, T const& jet, U const& /*tracks*/){ - if(jet.pt() < jetPtMin){ + template + void fillChargedJetTable(int32_t storedColIndex, T const& jet, U const& /*tracks*/) + { + if (jet.pt() < jetPtMin) { return; - } - ushort nconst = 0; - float leadingTrackPt = 0; - for(auto& constituent : jet.template tracks_as()){ + } + ushort nconst = 0; + float leadingTrackPt = 0; + for (auto& constituent : jet.template tracks_as()) { mHistograms.fill(HIST("chjetpt_vs_constpt"), jet.pt(), constituent.pt()); nconst++; - if(constituent.pt() > leadingTrackPt){ + if (constituent.pt() > leadingTrackPt) { leadingTrackPt = constituent.pt(); } } @@ -1145,13 +1168,14 @@ struct GammaJetTreeProducer { /// \brief Processes MC particles and fills MC particle table /// \param collision The collision to process /// \param mcgenparticles The MC particles to process - void processMCParticles(soa::Join::iterator const& collision, aod::JetParticles const& mcgenparticles, MCCol const&){ + void processMCParticles(soa::Join::iterator const& collision, aod::JetParticles const& mcgenparticles, MCCol const&) + { // event selection int32_t storedColIndex = getStoredColIndex(collision); if (storedColIndex == -1) return; - if(!collision.has_mcCollision()){ + if (!collision.has_mcCollision()) { return; } @@ -1162,21 +1186,21 @@ struct GammaJetTreeProducer { buildKdTree(particlesPerMcCollision); // Now we want to store every pi0 and every prompt photon that we find on generator level - for(auto particle : particlesPerMcCollision){ + for (auto particle : particlesPerMcCollision) { // only store particles above a given threshold - if(particle.pt() < minMCGenPt){ + if (particle.pt() < minMCGenPt) { continue; } // Test if a particle is a physical primary according to the following definition: // 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; } // only store photons and pi0s in mcgen stack - if(particle.pdgCode() != PDG_t::kPi0 && particle.pdgCode() != PDG_t::kGamma){ + if (particle.pdgCode() != PDG_t::kPi0 && particle.pdgCode() != PDG_t::kGamma) { continue; } // check the origin of the particle @@ -1189,32 +1213,31 @@ 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))) { mHistograms.fill(HIST("mcGenTrigger_E_PromptPhoton"), particle.energy()); } - if(origin & (1 << static_cast(gjanalysis::ParticleOrigin::kDirectPromptPhoton))){ + if (origin & (1 << static_cast(gjanalysis::ParticleOrigin::kDirectPromptPhoton))) { mHistograms.fill(HIST("mcGenTrigger_E_DirectPromptPhoton"), particle.energy()); } - if(origin & (1 << static_cast(gjanalysis::ParticleOrigin::kFragmentationPhoton))){ + if (origin & (1 << static_cast(gjanalysis::ParticleOrigin::kFragmentationPhoton))) { mHistograms.fill(HIST("mcGenTrigger_E_FragmentationPhoton"), particle.energy()); } - if(origin & (1 << static_cast(gjanalysis::ParticleOrigin::kDecayPhoton))){ + if (origin & (1 << static_cast(gjanalysis::ParticleOrigin::kDecayPhoton))) { mHistograms.fill(HIST("mcGenTrigger_E_DecayPhoton"), particle.energy()); } - if(origin & (1 << static_cast(gjanalysis::ParticleOrigin::kDecayPhotonPi0))){ + if (origin & (1 << static_cast(gjanalysis::ParticleOrigin::kDecayPhotonPi0))) { mHistograms.fill(HIST("mcGenTrigger_E_DecayPhotonPi0"), particle.energy()); } - if(origin & (1 << static_cast(gjanalysis::ParticleOrigin::kDecayPhotonEta))){ + if (origin & (1 << static_cast(gjanalysis::ParticleOrigin::kDecayPhotonEta))) { mHistograms.fill(HIST("mcGenTrigger_E_DecayPhotonEta"), particle.energy()); } - if(origin & (1 << static_cast(gjanalysis::ParticleOrigin::kDecayPhotonOther))){ + if (origin & (1 << static_cast(gjanalysis::ParticleOrigin::kDecayPhotonOther))) { mHistograms.fill(HIST("mcGenTrigger_E_DecayPhotonOther"), particle.energy()); } - if(origin & (1 << static_cast(gjanalysis::ParticleOrigin::kPi0))){ + if (origin & (1 << static_cast(gjanalysis::ParticleOrigin::kPi0))) { mHistograms.fill(HIST("mcGenTrigger_E_Pi0"), particle.energy()); } } - } PROCESS_SWITCH(GammaJetTreeProducer, processMCParticles, "Process MC particles", false); @@ -1224,18 +1247,19 @@ struct GammaJetTreeProducer { /// \param collision The collision to process /// \param chargedJets The MC particle level charged jets to process /// \param mcCollisions The MC collisions collection - void processChargedJetsMCP(soa::Join::iterator const& collision, soa::Filtered> const& chargedJets, MCCol const&){ + void processChargedJetsMCP(soa::Join::iterator const& collision, soa::Filtered> const& chargedJets, MCCol const&) + { // event selection int32_t storedColIndex = getStoredColIndex(collision); if (storedColIndex == -1) return; // loop over charged jets - if(!collision.has_mcCollision()){ + if (!collision.has_mcCollision()) { return; } int localIndex = 0; auto pjetsPerMcCollision = chargedJets.sliceBy(PJetsPerMCCollisions, collision.mcCollisionId()); - for(auto pjet : pjetsPerMcCollision){ + for (auto pjet : pjetsPerMcCollision) { // fill MC particle level jet table float perpconerho = ch_perp_cone_rho(pjet, perpConeJetR, true); mcJetsTable(storedColIndex, pjet.pt(), pjet.eta(), pjet.phi(), pjet.r(), pjet.energy(), pjet.mass(), pjet.area(), perpconerho); @@ -1254,7 +1278,7 @@ struct GammaJetTreeProducer { /// \param chargedJets The MC detector level charged jets to process /// \param tracks The tracks collection /// \param pjets The MC particle level jets collection (just loaded to have subscription to the table) - void processChargedJetsMCD(soa::Join::iterator const& collision, soa::Filtered> const& chargedJets, aod::JetTracks const& tracks, JetMCPTable const& /*pjets*/) + void processChargedJetsMCD(soa::Join::iterator const& collision, soa::Filtered> const& chargedJets, aod::JetTracks const& tracks, JetMCPTable const& /*pjets*/) { // event selection int32_t storedColIndex = getStoredColIndex(collision); @@ -1268,12 +1292,12 @@ struct GammaJetTreeProducer { int iLocalIndexGeo = -1; int iLocalIndexPt = -1; // We will always store the information for both in our tree - if(jet.has_matchedJetGeo()){ + if (jet.has_matchedJetGeo()) { const auto& pjet = jet.template matchedJetGeo_first_as(); iLocalIndexGeo = mcJetIndexMapping[pjet.globalIndex()]; mHistograms.fill(HIST("mcdJetPtVsTrueJetPtMatchingGeo"), jet.pt(), pjet.pt()); } - if(jet.has_matchedJetPt()){ + if (jet.has_matchedJetPt()) { const auto& pjet = jet.template matchedJetPt_first_as(); iLocalIndexPt = mcJetIndexMapping[pjet.globalIndex()]; mHistograms.fill(HIST("mcdJetPtVsTrueJetPtMatchingPt"), jet.pt(), pjet.pt()); @@ -1289,4 +1313,5 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) WorkflowSpec workflow{ adaptAnalysisTask(cfgc, TaskName{"gamma-jet-tree-producer"})}; return workflow; -} \ No newline at end of file +} + \ No newline at end of file