Skip to content

[PWGJE] isolation fixes & MC prompt photon information fixes #16308

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
fjonasALICE wants to merge 3 commits into
base: master
Choose a base branch
from
+177 −37
Open

[PWGJE] isolation fixes & MC prompt photon information fixes #16308

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
40085e2
[PWGJE] isolation fixes & MC prompt photon information fixes
fjonasALICE May 19, 2026
5bf6a78
Please consider the following formatting changes
alibuild May 19, 2026
6550a08
Merge pull request #18 from alibuild/alibot-cleanup-16308
fjonasALICE May 20, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
214 changes: 177 additions & 37 deletions PWGJE/Tasks/gammaJetTreeProducer.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -55,6 +55,7 @@
#include <string>
#include <type_traits>
#include <unordered_map>
#include <unordered_set>
#include <utility>
#include <vector>

Expand Down Expand Up @@ -120,9 +121,11 @@
std::vector<float> trackEta;
std::vector<float> trackPhi;
std::vector<float> trackPt;
std::vector<int> trackSourceIndex;
std::vector<float> mcParticleEta;
std::vector<float> mcParticlePhi;
std::vector<float> mcParticlePt;
std::vector<int> mcParticleSourceIndex;
TKDTree<int, float>* trackTree = nullptr;
TKDTree<int, float>* mcParticleTree = nullptr;

Expand Down Expand Up @@ -265,33 +268,57 @@
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<typename std::decay_t<T>, aod::JetTracks>) {
trackEta.clear();
trackPhi.clear();
trackPt.clear();
trackSourceIndex.clear();
const float maxMatchingDistance = std::max(static_cast<float>(isoR), static_cast<float>(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<int>(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);

Check failure on line 291 in PWGJE/Tasks/gammaJetTreeProducer.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[two-pi-add-subtract] 

Use RecoDecay::constrainAngle to restrict angle to a given range.
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);

Check failure on line 297 in PWGJE/Tasks/gammaJetTreeProducer.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[two-pi-add-subtract] 

Use RecoDecay::constrainAngle to restrict angle to a given range.
trackPt.push_back(track.pt());
trackSourceIndex.push_back(originalIndex);
}
}
if (trackEta.size() > 0) {
delete trackTree;
trackTree = new TKDTree<int, float>(trackEta.size(), 2, 1);
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<typename std::decay_t<T>, aod::JetParticles>) {
mcParticleEta.clear();
mcParticlePhi.clear();
mcParticlePt.clear();
mcParticleSourceIndex.clear();
const float maxMatchingDistance = std::max(static_cast<float>(isoR), static_cast<float>(perpConeJetR));
const float additionalMargin = 0.05f;
for (const auto& particle : objects) {
if (!particle.isPhysicalPrimary()) {
continue;
Expand All @@ -302,16 +329,35 @@
if (particle.pt() < trackMinPt) {
continue;
}
const float phi = RecoDecay::constrainAngle(particle.phi(), 0.0);
const int originalIndex = static_cast<int>(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);

Check failure on line 340 in PWGJE/Tasks/gammaJetTreeProducer.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[two-pi-add-subtract] 

Use RecoDecay::constrainAngle to restrict angle to a given range.
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);

Check failure on line 346 in PWGJE/Tasks/gammaJetTreeProducer.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[two-pi-add-subtract] 

Use RecoDecay::constrainAngle to restrict angle to a given range.
mcParticlePt.push_back(particle.pt());
mcParticleSourceIndex.push_back(originalIndex);
}
}
if (mcParticleEta.size() > 0) {
delete mcParticleTree;
mcParticleTree = new TKDTree<int, float>(mcParticleEta.size(), 2, 1);
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;
}
}
}
Expand Down Expand Up @@ -350,7 +396,7 @@
{
mHistograms.fill(HIST("eventQA"), 0);

if (collision.posZ() > mVertexCut) {
if (std::abs(collision.posZ()) > mVertexCut) {
return false;
}
mHistograms.fill(HIST("eventQA"), 1);
Expand Down Expand Up @@ -385,6 +431,33 @@
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<int>& indices,
const std::vector<int>& sourceIndexMap,
const std::vector<float>& ptValues,
std::unordered_set<int>& uniqueSourceIndices)
{
double ptSum = 0;
for (const auto& index : indices) {
if (index < 0 || static_cast<size_t>(index) >= sourceIndexMap.size()) {
continue;
}
const int sourceIndex = sourceIndexMap[index];
if (sourceIndex < 0 || static_cast<size_t>(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
Expand All @@ -394,25 +467,22 @@
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<int> indices;
std::unordered_set<int> 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;
}
} 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;
Expand All @@ -435,14 +505,16 @@
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<int> indicesLeft;
std::vector<int> indicesRight;
std::unordered_set<int> uniqueSourceIndicesLeft;
std::unordered_set<int> uniqueSourceIndicesRight;

if (!mcGenIso) {
if (trackTree) {
Expand All @@ -453,12 +525,8 @@
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);
Expand All @@ -467,12 +535,8 @@
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);
Expand All @@ -494,6 +558,46 @@
}
}

/// \brief Prints the mother chain of a given MC particle.
/// \param particle The particle to print the mother chain for
template <typename T>
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<aod::JMcParticles>();
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 &&

Check failure on line 586 in PWGJE/Tasks/gammaJetTreeProducer.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
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
Expand All @@ -504,10 +608,36 @@
T currentParticle = particle;
int pdgCode = particle.pdgCode();
auto mothers = particle.template mothers_as<aod::JMcParticles>();
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) {

Check failure on line 617 in PWGJE/Tasks/gammaJetTreeProducer.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
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<aod::JMcParticles>();
selectedMother = -1;
twoMothersIdentical = false;
if (mothers.size() == 1) {
selectedMother = 0;
} else if (mothers.size() == 2) {

Check failure on line 633 in PWGJE/Tasks/gammaJetTreeProducer.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
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;
}
Expand All @@ -517,7 +647,7 @@
/// \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) {

Check failure on line 650 in PWGJE/Tasks/gammaJetTreeProducer.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
return true;
}
return false;
Expand Down Expand Up @@ -588,7 +718,7 @@
// check if it has mothers that are etas
const auto& mothers = particle.template mothers_as<aod::JMcParticles>();
for (const auto& mother : mothers) {
if (mother.pdgCode() == 221) {

Check failure on line 721 in PWGJE/Tasks/gammaJetTreeProducer.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[pdg/explicit-code] 

Avoid hard-coded PDG codes. Use named values from PDG_t or o2::constants::physics::Pdg instead.
return true;
}
}
Expand All @@ -605,7 +735,7 @@
// check if you find a pi0 mother or a eta mother
const auto& mothers = particle.template mothers_as<aod::JMcParticles>();
for (const auto& mother : mothers) {
if (mother.pdgCode() == PDG_t::kPi0 || mother.pdgCode() == 221) {

Check failure on line 738 in PWGJE/Tasks/gammaJetTreeProducer.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[pdg/explicit-code] 

Avoid hard-coded PDG codes. Use named values from PDG_t or o2::constants::physics::Pdg instead.
return false;
}
}
Expand Down Expand Up @@ -889,7 +1019,7 @@
// 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<uint16_t>(gjanalysis::ClusterOrigin::kConvertedPhoton));
LOG(debug) << "Cluster is a converted photon";
}
Expand Down Expand Up @@ -932,7 +1062,7 @@
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)) {
Expand Down Expand Up @@ -1203,7 +1333,7 @@
if (jet.pt() < jetPtMin) {
return;
}
for (auto& jetConstituent : jet.template tracks_as<U>()) {
for (const auto& jetConstituent : jet.template tracks_as<U>()) {
fastjetutilities::fillTracks(jetConstituent, jetConstituents, jetConstituent.globalIndex());
}

Expand Down Expand Up @@ -1305,11 +1435,21 @@
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<uint16_t>(gjanalysis::ParticleOrigin::kPromptPhoton))) << "), kDirectPromptPhoton (" << (origin & (1 << static_cast<uint16_t>(gjanalysis::ParticleOrigin::kDirectPromptPhoton))) << "), kFragmentationPhoton (" << (origin & (1 << static_cast<uint16_t>(gjanalysis::ParticleOrigin::kFragmentationPhoton))) << "), kDecayPhoton (" << (origin & (1 << static_cast<uint16_t>(gjanalysis::ParticleOrigin::kDecayPhoton))) << "), kDecayPhotonPi0 (" << (origin & (1 << static_cast<uint16_t>(gjanalysis::ParticleOrigin::kDecayPhotonPi0))) << "), kDecayPhotonEta (" << (origin & (1 << static_cast<uint16_t>(gjanalysis::ParticleOrigin::kDecayPhotonEta))) << "), kDecayPhotonOther (" << (origin & (1 << static_cast<uint16_t>(gjanalysis::ParticleOrigin::kDecayPhotonOther))) << "), kPi0 (" << (origin & (1 << static_cast<uint16_t>(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());
Expand Down
Loading