Skip to content

[PWGEM] Fix potential mismatches between data and MC table size #17411

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

Merged
mhemmer-cern merged 1 commit into from
Aug 11, 2026
+162 −132
Merged

[PWGEM] Fix potential mismatches between data and MC table size #17411

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
221 changes: 128 additions & 93 deletions PWGEM/PhotonMeson/TableProducer/associateMCinfoPhoton.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
// Copyright 2019-2020 CERN and copyright holders of ALICE O2.

Check failure on line 1 in PWGEM/PhotonMeson/TableProducer/associateMCinfoPhoton.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[name/workflow-file] 

Name of a workflow file must match the name of the main struct in it (without the PWG prefix). (Class implementation files should be in "Core" directories.)

Check failure on line 1 in PWGEM/PhotonMeson/TableProducer/associateMCinfoPhoton.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[name/o2-task] 

Specify task name only when it cannot be derived from the struct name. Only append to the default name.
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
// All rights not expressly granted are reserved.
//
Expand Down Expand Up @@ -73,6 +73,8 @@
kElectron = 0x8,
};

static constexpr uint16_t kNoMcParticleMask = (1 << 14); // isNoise bit: no genuine MC particle for this leg/track

Produces<o2::aod::EMMCEvents> mcevents;
Produces<o2::aod::EMMCEventLabels> mceventlabels;
Produces<o2::aod::EMMCParticles> emmcparticles;
Expand All @@ -82,14 +84,15 @@

Produces<o2::aod::BinnedGenPts> binnedGenPt;

Configurable<float> max_eta_gen_secondary{"max_eta_gen_secondary", 0.9, "max eta to store generated information"};

Check failure on line 87 in PWGEM/PhotonMeson/TableProducer/associateMCinfoPhoton.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[name/configurable] 

Use lowerCamelCase for names of configurables and use the same name for the struct member as for the JSON string. (Declare the type and names on the same line.)
Configurable<float> margin_z_gen{"margin_z_gen", 15.f, "margin for Z of true photon conversion point to store generated information"};

Check failure on line 88 in PWGEM/PhotonMeson/TableProducer/associateMCinfoPhoton.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[name/configurable] 

Use lowerCamelCase for names of configurables and use the same name for the struct member as for the JSON string. (Declare the type and names on the same line.)
Configurable<float> max_rxy_gen{"max_rxy_gen", 100, "max rxy to store generated information"};

Check failure on line 89 in PWGEM/PhotonMeson/TableProducer/associateMCinfoPhoton.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[name/configurable] 

Use lowerCamelCase for names of configurables and use the same name for the struct member as for the JSON string. (Declare the type and names on the same line.)
Configurable<bool> requireGammaGammaDecay{"requireGammaGammaDecay", false, "require gamma gamma decay for generated pi0 and eta meson"};
Configurable<float> cfg_max_eta_photon{"cfg_max_eta_gamma", 0.8, "max eta for photons at PV"};

Check failure on line 91 in PWGEM/PhotonMeson/TableProducer/associateMCinfoPhoton.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[name/configurable] 

Use lowerCamelCase for names of configurables and use the same name for the struct member as for the JSON string. (Declare the type and names on the same line.)
Configurable<float> maxPt{"maxPt", 20.f, "max pT for BinnedGenPts table"};
Configurable<float> maxY{"maxY", 0.9f, "max |rapidity| for BinnedGenPts table"};
Configurable<bool> doStoreAllDaughters{"doStoreAllDaughters", false, "flag to enable storing of all photon, pi0, eta, eta' and omega daughters. This will increase the dervied data size!"};
Configurable<bool> useIsSelected{"useIsSelected", true, "flag to enable the usage of isSelected from our EM event selection. Disabling this will increase the derived data size!"};

HistogramRegistry registry{"EMMCEvent"};

Expand Down Expand Up @@ -164,15 +167,34 @@
}
}

template <soa::is_iterator TCollision>
bool passesCollisionGate(TCollision const& collisionIter)
{
if (!collisionIter.has_mcCollision()) {
return false;
}
if (useIsSelected.value && !collisionIter.isSelected()) {
return false;
}
return true;
}

template <o2::soa::is_iterator TMCParticle, o2::soa::is_iterator TMCCollision>
void selectMothersToStore(int motherId, int64_t mcParticleSize, TMCParticle& motherParticle, TMCParticle& daughterIter, TMCCollision const& mcCollisionIter, std::unordered_map<uint64_t, int>& fNewLabels, std::map<uint64_t, int>& fNewLabelsReversed, std::unordered_map<uint64_t, int>& fEventIdx, std::unordered_map<uint64_t, int>& fEventLabels, Counter& fCounter)
{
while (motherId > -1) {
if (motherId < mcParticleSize) { // protect against bad mother indices. why is this needed?
if (motherId < mcParticleSize) {
motherParticle.setCursor(motherId);
int eventIdx = fEventLabels.find(mcCollisionIter.globalIndex())->second;

// if the MC truth particle corresponding to this reconstructed track which is not already written, add it to the skimmed MC stack
auto evIt = fEventLabels.find(mcCollisionIter.globalIndex());
if (evIt == fEventLabels.end()) {
// Can't attribute this ancestor to a known event; stop walking rather than risk UB.
LOG(warning) << "selectMothersToStore: mcCollision " << mcCollisionIter.globalIndex()
<< " not found in fEventLabels, aborting mother-chain walk.";
break;
}
int eventIdx = evIt->second;

if (!fNewLabels.contains(motherParticle.globalIndex())) {
fNewLabels[motherParticle.globalIndex()] = fCounter.particles;
fNewLabelsReversed[fCounter.particles] = motherParticle.globalIndex();
Expand Down Expand Up @@ -232,12 +254,7 @@
for (const auto& collision : collisions) {
registry.fill(HIST("hEventCounter"), 1);

// TODO: investigate the collisions without corresponding mcCollision
if (!collision.has_mcCollision()) {
continue;
}

if (!collision.isSelected()) {
if (!passesCollisionGate(collision)) {
continue;
}

Expand Down Expand Up @@ -306,11 +323,11 @@
mceventlabels(fEventLabels.find(mcCollisionIter.globalIndex())->second, collision.mcMask());

for (const auto& mcParticle : groupedMcParticles) { // store necessary information for denominator of efficiency
if (mcParticle.pt() < 1e-3 || std::fabs(mcParticle.vz()) > 250 || std::sqrt(std::pow(mcParticle.vx(), 2) + std::pow(mcParticle.vy(), 2)) > max_rxy_gen) {

Check failure on line 326 in PWGEM/PhotonMeson/TableProducer/associateMCinfoPhoton.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.
continue;
}
int pdg = mcParticle.pdgCode();
if (std::abs(pdg) > 1e+9) {

Check failure on line 330 in PWGEM/PhotonMeson/TableProducer/associateMCinfoPhoton.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.
continue;
}

Expand Down Expand Up @@ -359,122 +376,135 @@
auto o2TrackIter = o2tracks.begin();
for (const auto& v0 : v0photons) {
collisionIter.setCursor(v0.collisionId());
if (!collisionIter.has_mcCollision()) {
continue;
}
mcCollisionIter.setCursor(collisionIter.mcCollisionId());

ele.setCursor(v0.negTrackId());
pos.setCursor(v0.posTrackId());
// check if the collision is valid AND the legs have valid MC particles
bool validMc = false;
if (passesCollisionGate(collisionIter)) {
mcCollisionIter.setCursor(collisionIter.mcCollisionId());

o2TrackEle.setCursor(ele.trackId());
o2TrackPos.setCursor(pos.trackId());
ele.setCursor(v0.negTrackId());
pos.setCursor(v0.posTrackId());
o2TrackEle.setCursor(ele.trackId());
o2TrackPos.setCursor(pos.trackId());

if (!o2TrackEle.has_mcParticle() || !o2TrackPos.has_mcParticle()) {
continue; // If no MC particle is found, skip the v0
validMc = o2TrackEle.has_mcParticle() && o2TrackPos.has_mcParticle();
}

for (const auto& leg : {pos, ele}) { // be carefull of order {pos, ele}!
o2TrackIter.setCursor(leg.trackId());
mcParticleIter.setCursor(o2TrackIter.mcParticleId());
// LOGF(info, "mcParticleIter.globalIndex() = %d, mcParticleIter.index() = %d", mcParticleIter.globalIndex(), mcParticleIter.index()); // these are exactly the same.

// if the MC truth particle corresponding to this reconstructed track which is not already written, add it to the skimmed MC stack
auto [iter, isNew] = fNewLabels.try_emplace(mcParticleIter.globalIndex(), fCounter.particles);
if (isNew) {
fNewLabelsReversed[fCounter.particles] = mcParticleIter.globalIndex();
fEventIdx[mcParticleIter.globalIndex()] = fEventLabels.find(mcCollisionIter.globalIndex())->second;
fCounter.particles++;
}
v0legmclabels(iter->second, o2TrackIter.mcMask());
for (const auto& leg : {pos, ele}) { // careful: order {pos, ele}, must stay fixed either way
if (validMc) {
o2TrackIter.setCursor(leg.trackId());
mcParticleIter.setCursor(o2TrackIter.mcParticleId());

auto [iter, isNew] = fNewLabels.try_emplace(mcParticleIter.globalIndex(), fCounter.particles);
if (isNew) {
fNewLabelsReversed[fCounter.particles] = mcParticleIter.globalIndex();
auto evIt = fEventLabels.find(mcCollisionIter.globalIndex());
if (evIt != fEventLabels.end()) {
fEventIdx[mcParticleIter.globalIndex()] = evIt->second;
}
fCounter.particles++;
}
v0legmclabels(iter->second, o2TrackIter.mcMask());

// Next, store mother-chain of this reconstructed track.
int motherid = -999; // first mother index
if (mcParticleIter.has_mothers()) {
motherid = mcParticleIter.mothersIds()[0]; // first mother index
int motherid = -999;
if (mcParticleIter.has_mothers()) {
motherid = mcParticleIter.mothersIds()[0];
}
selectMothersToStore(motherid, mcParticles.size(), motherParticle, daughterIter, mcCollisionIter,
fNewLabels, fNewLabelsReversed, fEventIdx, fEventLabels, fCounter);
} else {
v0legmclabels(-1, kNoMcParticleMask); // no MC match for this leg, but the row must still exist
}
selectMothersToStore(motherid, mcParticles.size(), motherParticle, daughterIter, mcCollisionIter, fNewLabels, fNewLabelsReversed, fEventIdx, fEventLabels, fCounter);
} // end of leg loop
} // end of loop over legs {pos, ele}
} // end of v0 loop
}

if constexpr (static_cast<bool>(system & kElectron)) {
auto o2TrackIter = o2tracks.begin();
// auto emprimaryelectrons_coll = emprimaryelectrons.sliceBy(perCollisionEl, collision.globalIndex());
for (const auto& emprimaryelectron : emprimaryelectrons) {
collisionIter.setCursor(emprimaryelectron.collisionId());
if (!collisionIter.has_mcCollision()) {
continue;
}
mcCollisionIter.setCursor(collisionIter.mcCollisionId());

o2TrackIter.setCursor(emprimaryelectron.trackId());
if (!o2TrackIter.has_mcParticle()) {
continue; // If no MC particle is found, skip the dilepton
}
mcParticleIter.setCursor(o2TrackIter.mcParticleId());
int32_t mcLabel = -1;
uint16_t mcMask = kNoMcParticleMask;

// if the MC truth particle corresponding to this reconstructed track which is not already written, add it to the skimmed MC stack
auto [iter, isNew] = fNewLabels.try_emplace(mcParticleIter.globalIndex(), fCounter.particles);
if (isNew) {
fNewLabelsReversed[fCounter.particles] = mcParticleIter.globalIndex();
fEventIdx[mcParticleIter.globalIndex()] = fEventLabels.find(mcCollisionIter.globalIndex())->second;
fCounter.particles++;
}
emprimaryelectronmclabels(iter->second, o2TrackIter.mcMask());
if (passesCollisionGate(collisionIter)) {
mcCollisionIter.setCursor(collisionIter.mcCollisionId());
o2TrackIter.setCursor(emprimaryelectron.trackId());

if (o2TrackIter.has_mcParticle()) {
mcParticleIter.setCursor(o2TrackIter.mcParticleId());

// if the MC truth particle corresponding to this reconstructed track which is not already written, add it to the skimmed MC stack
auto [iter, isNew] = fNewLabels.try_emplace(mcParticleIter.globalIndex(), fCounter.particles);
if (isNew) {
fNewLabelsReversed[fCounter.particles] = mcParticleIter.globalIndex();
auto evIt = fEventLabels.find(mcCollisionIter.globalIndex());
if (evIt != fEventLabels.end()) {
fEventIdx[mcParticleIter.globalIndex()] = evIt->second;
}
fCounter.particles++;
}
mcLabel = iter->second;
mcMask = o2TrackIter.mcMask();

// Next, store mother-chain of this reconstructed track.
int motherid = -999; // first mother index
if (mcParticleIter.has_mothers()) {
motherid = mcParticleIter.mothersIds()[0]; // first mother index
int motherid = -999;
if (mcParticleIter.has_mothers()) {
motherid = mcParticleIter.mothersIds()[0];
}
selectMothersToStore(motherid, mcParticles.size(), motherParticle, daughterIter, mcCollisionIter,
fNewLabels, fNewLabelsReversed, fEventIdx, fEventLabels, fCounter);
} // if (o2TrackIter.has_mcParticle())
}
selectMothersToStore(motherid, mcParticles.size(), motherParticle, daughterIter, mcCollisionIter, fNewLabels, fNewLabelsReversed, fEventIdx, fEventLabels, fCounter);
emprimaryelectronmclabels(mcLabel, mcMask); // always exactly one row per emprimaryelectron
} // end of em primary electron loop
}

if constexpr (static_cast<bool>(system & kEMC)) { // for emc photons
// auto ememcclusters_coll = emcclusters.sliceBy(perCollisionEMC, collision.globalIndex());
for (const auto& emccluster : emcclusters) {
collisionIter.setCursor(emccluster.collisionId());
if (!collisionIter.has_mcCollision()) {
continue;
}
mcCollisionIter.setCursor(collisionIter.mcCollisionId());

// TODO: test
if (emccluster.mcParticleIds().size() <= 0) {
continue;
}
std::vector<int32_t> vEmcMcParticleIds;
std::vector<float> vAmplitudes;

vEmcMcParticleIds.reserve(emccluster.mcParticleIds().size());
vAmplitudes.reserve(emccluster.mcParticleIds().size());

for (size_t iCont = 0; iCont < emccluster.mcParticleIds().size(); iCont++) {
mcPhoton.setCursor(emccluster.mcParticleIds()[iCont]);
if (passesCollisionGate(collisionIter)) {
mcCollisionIter.setCursor(collisionIter.mcCollisionId());

// if the MC truth particle corresponding to this reconstructed track which is not already written, add it to the skimmed MC stack
auto [iter, isNew] = fNewLabels.try_emplace(mcPhoton.globalIndex(), fCounter.particles);
if (isNew) {
fNewLabelsReversed[fCounter.particles] = mcPhoton.globalIndex();
fEventIdx[mcPhoton.globalIndex()] = fEventLabels.find(mcCollisionIter.globalIndex())->second;
fCounter.particles++;
if (emccluster.mcParticleIds().size() != emccluster.amplitude().size()) {
LOG(warning) << "Cluster " << emccluster.globalIndex()
<< ": mcParticleIds/amplitude size mismatch ("
<< emccluster.mcParticleIds().size() << " vs " << emccluster.amplitude().size() << ")";
}
vEmcMcParticleIds.emplace_back(iter->second);
vAmplitudes.emplace_back(emccluster.amplitude()[iCont]);
// ememcclustermclabels(fNewLabels.find(mcPhoton.index())->second);

// Next, store mother-chain of this reconstructed track.
int motherid = -999; // first mother index
if (mcPhoton.has_mothers()) {
motherid = mcPhoton.mothersIds()[0]; // first mother index
size_t nCont = std::min(emccluster.mcParticleIds().size(), emccluster.amplitude().size());

vEmcMcParticleIds.reserve(nCont);
vAmplitudes.reserve(nCont);

for (size_t iCont = 0; iCont < nCont; iCont++) {
mcPhoton.setCursor(emccluster.mcParticleIds()[iCont]);
// if the MC truth particle corresponding to this reconstructed track which is not already written, add it to the skimmed MC stack
auto [iter, isNew] = fNewLabels.try_emplace(mcPhoton.globalIndex(), fCounter.particles);
if (isNew) {
fNewLabelsReversed[fCounter.particles] = mcPhoton.globalIndex();
auto evIt = fEventLabels.find(mcCollisionIter.globalIndex());
if (evIt != fEventLabels.end()) { // keep your earlier guard here too
fEventIdx[mcPhoton.globalIndex()] = evIt->second;
}
fCounter.particles++;
}
vEmcMcParticleIds.emplace_back(iter->second);
vAmplitudes.emplace_back(emccluster.amplitude()[iCont]);

int motherid = -999;
if (mcPhoton.has_mothers()) {
motherid = mcPhoton.mothersIds()[0];
}
selectMothersToStore(motherid, mcParticles.size(), motherParticle, daughterIter, mcCollisionIter,
fNewLabels, fNewLabelsReversed, fEventIdx, fEventLabels, fCounter);
}
selectMothersToStore(motherid, mcParticles.size(), motherParticle, daughterIter, mcCollisionIter, fNewLabels, fNewLabelsReversed, fEventIdx, fEventLabels, fCounter);
} // end of loop over mc particles of the current emc cluster
// whether or not the collision was okay, exactly one row is always written here so later the tables are ensured to have same size!
ememcclustermclabels(vEmcMcParticleIds, vAmplitudes);

} // end of em emc cluster loop
}
}

// Loop over the label map, create the mother/daughter relationships if these exist and write the skimmed MC stack
Expand Down Expand Up @@ -517,7 +547,12 @@
}
}

emmcparticles(fEventIdx.find(oldLabel)->second, mcParticleIter.pdgCode(), mcParticleIter.flags(), mcParticleIter.statusCode(),
auto evIt = fEventIdx.find(oldLabel);
if (evIt == fEventIdx.end()) {
LOG(warning) << "emmcparticles: no event index for particle label " << oldLabel << ", skipping.";
continue;
}
emmcparticles(evIt->second, mcParticleIter.pdgCode(), mcParticleIter.flags(), mcParticleIter.statusCode(),
mothers, daughters,
mcParticleIter.px(), mcParticleIter.py(), mcParticleIter.pz(), mcParticleIter.e(),
mcParticleIter.vx(), mcParticleIter.vy(), mcParticleIter.vz());
Expand All @@ -535,7 +570,7 @@
template <o2::soa::is_table Daughters>
inline bool areTwoPhotonDaughtersAccepted(const Daughters& lDaughters, float maxEta)
{
if (lDaughters.size() != 2) {

Check failure on line 573 in PWGEM/PhotonMeson/TableProducer/associateMCinfoPhoton.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 false;
}

Expand Down Expand Up @@ -635,5 +670,5 @@
WorkflowSpec defineDataProcessing(ConfigContext const& context)
{
return WorkflowSpec{
adaptAnalysisTask<AssociateMCInfoPhoton>(context, TaskName{"associate-mc-info-photon"})};

Check failure on line 673 in PWGEM/PhotonMeson/TableProducer/associateMCinfoPhoton.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[name/o2-task] 

Device names associate-mc-info-photon and associate-m-c-info-photon generated from the specified task name associate-mc-info-photon and from the struct name AssociateMCInfoPhoton, respectively, differ in hyphenation. Consider fixing capitalisation of the struct name to AssociateMcInfoPhoton and removing TaskName.
}
Loading
Loading