diff --git a/PWGJE/Tasks/gammaJetTreeProducer.cxx b/PWGJE/Tasks/gammaJetTreeProducer.cxx index de043149abe..c117287038c 100644 --- a/PWGJE/Tasks/gammaJetTreeProducer.cxx +++ b/PWGJE/Tasks/gammaJetTreeProducer.cxx @@ -55,6 +55,7 @@ #include #include #include +#include #include #include @@ -120,9 +121,11 @@ struct GammaJetTreeProducer { std::vector trackEta; std::vector trackPhi; std::vector trackPt; + std::vector trackSourceIndex; std::vector mcParticleEta; std::vector mcParticlePhi; std::vector mcParticlePt; + std::vector mcParticleSourceIndex; TKDTree* trackTree = nullptr; TKDTree* mcParticleTree = nullptr; @@ -265,22 +268,36 @@ struct GammaJetTreeProducer { curMcParticleTreeColIndex = collision.globalIndex(); } - trackEta.clear(); - trackPhi.clear(); - trackPt.clear(); - mcParticleEta.clear(); - mcParticlePhi.clear(); - mcParticlePt.clear(); - // 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>) { + trackEta.clear(); + trackPhi.clear(); + trackPt.clear(); + trackSourceIndex.clear(); + const float maxMatchingDistance = std::max(static_cast(isoR), static_cast(perpConeJetR)); + const float additionalMargin = 0.05f; for (const auto& track : objects) { if (!isTrackSelected(track)) { continue; } + const float phi = RecoDecay::constrainAngle(track.phi(), 0.0); + const int originalIndex = static_cast(trackPt.size()); trackEta.push_back(track.eta()); - trackPhi.push_back(track.phi()); + trackPhi.push_back(phi); trackPt.push_back(track.pt()); + trackSourceIndex.push_back(originalIndex); + if (phi <= (maxMatchingDistance + additionalMargin)) { + trackEta.push_back(track.eta()); + trackPhi.push_back(phi + o2::constants::math::TwoPI); + 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); + trackPt.push_back(track.pt()); + trackSourceIndex.push_back(originalIndex); + } } if (trackEta.size() > 0) { delete trackTree; @@ -288,10 +305,20 @@ struct GammaJetTreeProducer { trackTree->SetData(0, trackEta.data()); trackTree->SetData(1, trackPhi.data()); trackTree->Build(); + } else { + // Keep tree state consistent with the current event content. + delete trackTree; + trackTree = nullptr; } } // 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>) { + mcParticleEta.clear(); + mcParticlePhi.clear(); + mcParticlePt.clear(); + mcParticleSourceIndex.clear(); + const float maxMatchingDistance = std::max(static_cast(isoR), static_cast(perpConeJetR)); + const float additionalMargin = 0.05f; for (const auto& particle : objects) { if (!particle.isPhysicalPrimary()) { continue; @@ -302,9 +329,24 @@ struct GammaJetTreeProducer { if (particle.pt() < trackMinPt) { continue; } + const float phi = RecoDecay::constrainAngle(particle.phi(), 0.0); + const int originalIndex = static_cast(mcParticlePt.size()); mcParticleEta.push_back(particle.eta()); - mcParticlePhi.push_back(particle.phi()); + mcParticlePhi.push_back(phi); mcParticlePt.push_back(particle.pt()); + mcParticleSourceIndex.push_back(originalIndex); + if (phi <= (maxMatchingDistance + additionalMargin)) { + mcParticleEta.push_back(particle.eta()); + mcParticlePhi.push_back(phi + o2::constants::math::TwoPI); + 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); + mcParticlePt.push_back(particle.pt()); + mcParticleSourceIndex.push_back(originalIndex); + } } if (mcParticleEta.size() > 0) { delete mcParticleTree; @@ -312,6 +354,10 @@ struct GammaJetTreeProducer { mcParticleTree->SetData(0, mcParticleEta.data()); mcParticleTree->SetData(1, mcParticlePhi.data()); mcParticleTree->Build(); + } else { + // Keep tree state consistent with the current event content. + delete mcParticleTree; + mcParticleTree = nullptr; } } } @@ -350,7 +396,7 @@ struct GammaJetTreeProducer { { mHistograms.fill(HIST("eventQA"), 0); - if (collision.posZ() > mVertexCut) { + if (std::abs(collision.posZ()) > mVertexCut) { return false; } mHistograms.fill(HIST("eventQA"), 1); @@ -385,6 +431,33 @@ struct GammaJetTreeProducer { return std::abs(pdg->GetParticle(particle.pdgCode())->Charge()) >= 1.; } + /// \brief Sums pt values once per unique source index from KD-tree query indices + /// \param indices KD-tree result indices + /// \param sourceIndexMap Maps KD-tree entry index to original object index + /// \param ptValues pt values of original objects + /// \param uniqueSourceIndices Set used to deduplicate source indices + /// \return Sum of unique pt values for the given indices + double sumUniquePtFromIndices(const std::vector& indices, + const std::vector& sourceIndexMap, + const std::vector& ptValues, + std::unordered_set& uniqueSourceIndices) + { + double ptSum = 0; + for (const auto& index : indices) { + if (index < 0 || static_cast(index) >= sourceIndexMap.size()) { + continue; + } + const int sourceIndex = sourceIndexMap[index]; + if (sourceIndex < 0 || static_cast(sourceIndex) >= ptValues.size()) { + continue; + } + if (uniqueSourceIndices.insert(sourceIndex).second) { + ptSum += ptValues[sourceIndex]; + } + } + return ptSum; + } + /// \brief Calculates the charged particle isolation in a cone of given size using a pre-built kd tree /// \param particle The particle to calculate the isolation for /// \param radius The cone radius @@ -394,15 +467,14 @@ struct GammaJetTreeProducer { 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()}; + float point[2] = {particle.eta(), RecoDecay::constrainAngle(particle.phi(), 0.0)}; std::vector indices; + std::unordered_set uniqueSourceIndices; if (!mcGenIso) { if (trackTree) { trackTree->FindInRange(point, radius, indices); - for (const auto& index : indices) { - iso += trackPt[index]; - } + iso += sumUniquePtFromIndices(indices, trackSourceIndex, trackPt, uniqueSourceIndices); } else { LOG(error) << "Track tree not found"; return 0; @@ -410,9 +482,7 @@ struct GammaJetTreeProducer { } else { if (mcParticleTree) { mcParticleTree->FindInRange(point, radius, indices); - for (const auto& index : indices) { - iso += mcParticlePt[index]; - } + iso += sumUniquePtFromIndices(indices, mcParticleSourceIndex, mcParticlePt, uniqueSourceIndices); } else { LOG(error) << "MC particle tree not found"; return 0; @@ -435,14 +505,16 @@ struct GammaJetTreeProducer { double cPhi = TVector2::Phi_0_2pi(object.phi()); // rotate cone left by 90 degrees - float cPhiLeft = cPhi - o2::constants::math::PIHalf; - float cPhiRight = cPhi + o2::constants::math::PIHalf; + 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::vector indicesLeft; std::vector indicesRight; + std::unordered_set uniqueSourceIndicesLeft; + std::unordered_set uniqueSourceIndicesRight; if (!mcGenIso) { if (trackTree) { @@ -453,12 +525,8 @@ struct GammaJetTreeProducer { return 0; } - for (const auto& index : indicesLeft) { - ptSumLeft += trackPt[index]; - } - for (const auto& index : indicesRight) { - ptSumRight += trackPt[index]; - } + ptSumLeft += sumUniquePtFromIndices(indicesLeft, trackSourceIndex, trackPt, uniqueSourceIndicesLeft); + ptSumRight += sumUniquePtFromIndices(indicesRight, trackSourceIndex, trackPt, uniqueSourceIndicesRight); } else { if (mcParticleTree) { mcParticleTree->FindInRange(pointLeft, radius, indicesLeft); @@ -467,12 +535,8 @@ struct GammaJetTreeProducer { LOG(error) << "MC particle tree not found"; return 0; } - for (const auto& index : indicesLeft) { - ptSumLeft += mcParticlePt[index]; - } - for (const auto& index : indicesRight) { - ptSumRight += mcParticlePt[index]; - } + ptSumLeft += sumUniquePtFromIndices(indicesLeft, mcParticleSourceIndex, mcParticlePt, uniqueSourceIndicesLeft); + ptSumRight += sumUniquePtFromIndices(indicesRight, mcParticleSourceIndex, mcParticlePt, uniqueSourceIndicesRight); } float rho = (ptSumLeft + ptSumRight) / (o2::constants::math::TwoPI * radius * radius); @@ -494,6 +558,46 @@ struct GammaJetTreeProducer { } } + /// \brief Prints the mother chain of a given MC particle. + /// \param particle The particle to print the mother chain for + template + void printMotherChain(const T& particle) const + { + T current = particle; + int depth = 0; + bool continueTracing = true; + while (continueTracing) { + LOG(info) << "Level " << depth + << " | PDG: " << current.pdgCode() + << " | E: " << current.energy() + << " | pt: " << current.pt() + << " | eta: " << current.eta() + << " | phi: " << current.phi() + << " | genStatusCode: " << current.getGenStatusCode() + << " | globalIndex: " << current.globalIndex(); + auto mothers = current.template mothers_as(); + if (mothers.size() == 0) { + break; + } + // 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()) { + selectedMother = 0; + } + if (selectedMother == -1) { + LOG(info) << "Branching in mother chain or unclear mothers, stopping trace."; + break; + } + current = mothers[selectedMother]; + ++depth; + } + } + /// \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 @@ -504,10 +608,36 @@ struct GammaJetTreeProducer { 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) { - iUp = mothers[0].globalIndex(); - currentParticle = mothers[0]; + // sometimes the particle will have two identical mothers, not sure why, but we need to cover this case + // if there are two mothers, but they are not identical, we will stop + bool twoMothersIdentical = false; + int selectedMother = -1; + if (mothers.size() == 1) { + selectedMother = 0; + } else if (mothers.size() == 2) { + twoMothersIdentical = (mothers[0].globalIndex() == mothers[1].globalIndex() && + mothers[0].pdgCode() == mothers[1].pdgCode() && + mothers[0].getGenStatusCode() == mothers[1].getGenStatusCode()); + if (twoMothersIdentical) { + selectedMother = 0; + } + } + while (iUp > 0 && selectedMother >= 0 && mothers[selectedMother].globalIndex() > 0 && mothers[selectedMother].pdgCode() == pdgCode) { + iUp = mothers[selectedMother].globalIndex(); + currentParticle = mothers[selectedMother]; mothers = currentParticle.template mothers_as(); + selectedMother = -1; + twoMothersIdentical = false; + if (mothers.size() == 1) { + selectedMother = 0; + } else if (mothers.size() == 2) { + twoMothersIdentical = (mothers[0].globalIndex() == mothers[1].globalIndex() && + mothers[0].pdgCode() == mothers[1].pdgCode() && + mothers[0].getGenStatusCode() == mothers[1].getGenStatusCode()); + if (twoMothersIdentical) { + selectedMother = 0; + } + } } return currentParticle; } @@ -889,7 +1019,7 @@ struct GammaJetTreeProducer { // check that mother has exactly two daughters which are e+ and e- 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"; } @@ -932,7 +1062,7 @@ struct GammaJetTreeProducer { int nRecCollisions = 0; mHistograms.fill(HIST("mcCollisionsWithRecCollisions"), 0); for (auto const& collision : collisions) { - if (collision.posZ() > mVertexCut) { + if (std::abs(collision.posZ()) > mVertexCut) { continue; } if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits, true, true, rctLabel)) { @@ -1203,7 +1333,7 @@ struct GammaJetTreeProducer { if (jet.pt() < jetPtMin) { return; } - for (auto& jetConstituent : jet.template tracks_as()) { + for (const auto& jetConstituent : jet.template tracks_as()) { fastjetutilities::fillTracks(jetConstituent, jetConstituents, jetConstituent.globalIndex()); } @@ -1305,11 +1435,21 @@ struct GammaJetTreeProducer { if (particle.pdgCode() != PDG_t::kPi0 && particle.pdgCode() != PDG_t::kGamma) { continue; } + // check the origin of the particle uint16_t origin = getMCParticleOrigin(particle); double mcIsolation = ch_iso_in_cone(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 + // leaving this in here + // if (particle.pdgCode() == PDG_t::kGamma) { + // LOG(info) << "Particle PDG: " << particle.pdgCode() << " isPhysicalPrimary: " << particle.isPhysicalPrimary() << " getGenStatusCode: " << particle.getGenStatusCode(); + // LOG(info) << "Origin: " << "kPromptPhoton (" << (origin & (1 << static_cast(gjanalysis::ParticleOrigin::kPromptPhoton))) << "), kDirectPromptPhoton (" << (origin & (1 << static_cast(gjanalysis::ParticleOrigin::kDirectPromptPhoton))) << "), kFragmentationPhoton (" << (origin & (1 << static_cast(gjanalysis::ParticleOrigin::kFragmentationPhoton))) << "), kDecayPhoton (" << (origin & (1 << static_cast(gjanalysis::ParticleOrigin::kDecayPhoton))) << "), kDecayPhotonPi0 (" << (origin & (1 << static_cast(gjanalysis::ParticleOrigin::kDecayPhotonPi0))) << "), kDecayPhotonEta (" << (origin & (1 << static_cast(gjanalysis::ParticleOrigin::kDecayPhotonEta))) << "), kDecayPhotonOther (" << (origin & (1 << static_cast(gjanalysis::ParticleOrigin::kDecayPhotonOther))) << "), kPi0 (" << (origin & (1 << static_cast(gjanalysis::ParticleOrigin::kPi0))) << ")"; + // LOG(info) << "Mother chain: "; + // printMotherChain(particle); + // } + // fill mc gen trigger particle histograms mHistograms.fill(HIST("mcGenTrigger_E"), particle.energy()); mHistograms.fill(HIST("mcGenTrigger_Eta"), particle.eta());