From 05d94f404526c30a42bfa3dac8a27357cbd256e4 Mon Sep 17 00:00:00 2001 From: mattiaisgro Date: Sat, 14 Jun 2025 21:46:55 +0200 Subject: [PATCH 1/2] Add PDGCode enumeration for MC particle codes The PDGCode enumeration contains common particle codes which can be used to make code clearer and less error prone, instead of hard-coding particle codes. --- .../dataframe/FCCAnalyses/MCParticleCodes.h | 46 +++++++++++++ analyzers/dataframe/src/JetTaggingUtils.cc | 65 +++++++++++++------ analyzers/dataframe/src/MCParticle.cc | 13 ++-- 3 files changed, 99 insertions(+), 25 deletions(-) create mode 100644 analyzers/dataframe/FCCAnalyses/MCParticleCodes.h diff --git a/analyzers/dataframe/FCCAnalyses/MCParticleCodes.h b/analyzers/dataframe/FCCAnalyses/MCParticleCodes.h new file mode 100644 index 00000000000..da0730b257a --- /dev/null +++ b/analyzers/dataframe/FCCAnalyses/MCParticleCodes.h @@ -0,0 +1,46 @@ + +#ifndef MCPARTICLECODES_H +#define MCPARTICLECODES_H + +// This file contains the particle codes used in the FCC Analyses framework. +namespace FCCAnalyses { + + /// Particle codes following the PDG Monte Carlo particle numbering scheme. + enum class PDGCode { + + UNKNOWN = 0, + + // Leptons + ELECTRON = 11, + MUON = 13, + TAU = 15, + E_NEUTRINO = 12, + MU_NEUTRINO = 14, + TAU_NEUTRINO = 16, + + // Partons + QUARK_UP = 1, + QUARK_DOWN = 2, + QUARK_CHARM = 4, + QUARK_STRANGE = 3, + QUARK_TOP = 6, + QUARK_BOTTOM = 5, + GLUON = 21, + + // Bosons + PHOTON = 22, + Z = 23, + W = 24, + HIGGS = 25, + + // Baryons + PION = 211, + KAON = 321, + PROTON = 2212, + NEUTRON = 2112, + }; + +} // namespace FCCAnalyses + + +#endif diff --git a/analyzers/dataframe/src/JetTaggingUtils.cc b/analyzers/dataframe/src/JetTaggingUtils.cc index 7ed529c18a3..c874295d884 100644 --- a/analyzers/dataframe/src/JetTaggingUtils.cc +++ b/analyzers/dataframe/src/JetTaggingUtils.cc @@ -1,4 +1,5 @@ #include "FCCAnalyses/JetTaggingUtils.h" +#include "FCCAnalyses/PDGCodes.h" namespace FCCAnalyses { @@ -11,13 +12,17 @@ get_flavour(ROOT::VecOps::RVec in, int loopcount = 0; for (size_t i = 0; i < MCin.size(); ++i) { + auto &parton = MCin[i]; + // Select partons only (for pythia8 71-79, for pythia6 2): if ((parton.generatorStatus > 80 || parton.generatorStatus < 70) && parton.generatorStatus != 2) continue; - if (std::abs(parton.PDG) > 5 && parton.PDG != 21) + + if (std::abs(parton.PDG) > PDGCode::QUARK_BOTTOM && parton.PDG != PDGCode::GLUON) continue; + ROOT::Math::PxPyPzMVector lv(parton.momentum.x, parton.momentum.y, parton.momentum.z, parton.mass); @@ -39,14 +44,19 @@ get_flavour(ROOT::VecOps::RVec in, Float_t angle = acos(dot / norm); if (angle <= 0.3) { - if (result[j] == 21 or result[j] == 0) { + + if (result[j] == PDGCode::GLUON or result[j] == PDGCode::UNKNOWN) { + // if no match before, or matched to gluon, match to // this particle (favour quarks over gluons) result[j] = std::abs(parton.PDG); - } else if (parton.PDG != 21) { + + } else if (parton.PDG != PDGCode::GLUON) { + // if matched to quark, and this is a quark, favour // heavier flavours result[j] = std::max(result[j], std::abs(parton.PDG)); + } else { // if matched to quark, and this is a gluon, keep // previous result (favour quark) @@ -66,15 +76,20 @@ ROOT::VecOps::RVec get_btag(ROOT::VecOps::RVec in, float efficiency, ROOT::VecOps::RVec result(in.size(), 0); for (size_t j = 0; j < in.size(); ++j) { - if (in.at(j) == 5 && gRandom->Uniform() <= efficiency) + + if (in.at(j) == PDGCode::QUARK_BOTTOM && gRandom->Uniform() <= efficiency) result[j] = 1; - if (in.at(j) == 4 && gRandom->Uniform() <= mistag_c) + + if (in.at(j) == PDGCode::QUARK_CHARM && gRandom->Uniform() <= mistag_c) result[j] = 1; - if (in.at(j) < 4 && gRandom->Uniform() <= mistag_l) + + if (in.at(j) < PDGCode::QUARK_CHARM && gRandom->Uniform() <= mistag_l) result[j] = 1; - if (in.at(j) == 21 && gRandom->Uniform() <= mistag_g) + + if (in.at(j) == PDGCode::GLUON && gRandom->Uniform() <= mistag_g) result[j] = 1; } + return result; } @@ -85,13 +100,17 @@ ROOT::VecOps::RVec get_ctag(ROOT::VecOps::RVec in, float efficiency, ROOT::VecOps::RVec result(in.size(), 0); for (size_t j = 0; j < in.size(); ++j) { - if (in.at(j) == 4 && gRandom->Uniform() <= efficiency) + + if (in.at(j) == PDGCode::QUARK_CHARM && gRandom->Uniform() <= efficiency) result[j] = 1; - if (in.at(j) == 5 && gRandom->Uniform() <= mistag_b) + + if (in.at(j) == PDGCode::QUARK_BOTTOM && gRandom->Uniform() <= mistag_b) result[j] = 1; - if (in.at(j) < 4 && gRandom->Uniform() <= mistag_l) + + if (in.at(j) < PDGCode::QUARK_CHARM && gRandom->Uniform() <= mistag_l) result[j] = 1; - if (in.at(j) == 21 && gRandom->Uniform() <= mistag_g) + + if (in.at(j) == PDGCode::GLUON && gRandom->Uniform() <= mistag_g) result[j] = 1; } return result; @@ -104,13 +123,17 @@ ROOT::VecOps::RVec get_ltag(ROOT::VecOps::RVec in, float efficiency, ROOT::VecOps::RVec result(in.size(), 0); for (size_t j = 0; j < in.size(); ++j) { - if (in.at(j) < 4 && gRandom->Uniform() <= efficiency) + + if (in.at(j) < PDGCode::QUARK_CHARM && gRandom->Uniform() <= efficiency) result[j] = 1; - if (in.at(j) == 5 && gRandom->Uniform() <= mistag_b) + + if (in.at(j) == PDGCode::QUARK_BOTTOM && gRandom->Uniform() <= mistag_b) result[j] = 1; - if (in.at(j) == 4 && gRandom->Uniform() <= mistag_c) + + if (in.at(j) == PDGCode::QUARK_CHARM && gRandom->Uniform() <= mistag_c) result[j] = 1; - if (in.at(j) == 21 && gRandom->Uniform() <= mistag_g) + + if (in.at(j) == PDGCode::GLUON && gRandom->Uniform() <= mistag_g) result[j] = 1; } return result; @@ -123,13 +146,17 @@ ROOT::VecOps::RVec get_gtag(ROOT::VecOps::RVec in, float efficiency, ROOT::VecOps::RVec result(in.size(), 0); for (size_t j = 0; j < in.size(); ++j) { - if (in.at(j) == 21 && gRandom->Uniform() <= efficiency) + + if (in.at(j) == PDGCode::GLUON && gRandom->Uniform() <= efficiency) result[j] = 1; - if (in.at(j) == 5 && gRandom->Uniform() <= mistag_b) + + if (in.at(j) == PDGCode::QUARK_BOTTOM && gRandom->Uniform() <= mistag_b) result[j] = 1; - if (in.at(j) == 4 && gRandom->Uniform() <= mistag_c) + + if (in.at(j) == PDGCode::QUARK_CHARM && gRandom->Uniform() <= mistag_c) result[j] = 1; - if (in.at(j) < 4 && gRandom->Uniform() <= mistag_l) + + if (in.at(j) < PDGCode::QUARK_CHARM && gRandom->Uniform() <= mistag_l) result[j] = 1; } return result; diff --git a/analyzers/dataframe/src/MCParticle.cc b/analyzers/dataframe/src/MCParticle.cc index 009d0d4c255..f5b05c3d2cf 100644 --- a/analyzers/dataframe/src/MCParticle.cc +++ b/analyzers/dataframe/src/MCParticle.cc @@ -1,4 +1,5 @@ #include "FCCAnalyses/MCParticle.h" +#include "FCCAnalyses/PDGCodes.h" #include #include #include @@ -696,7 +697,7 @@ int get_lepton_origin(const edm4hep::MCParticleData &p, // std::cout << std::endl << " enter in MCParticle::get_lepton_origin PDG = " << p.PDG << std::endl; int pdg = std::abs( p.PDG ) ; - if ( pdg != 11 && pdg != 13 && pdg != 15 ) return -1; + if ( pdg != PDGCode::ELECTRON && pdg != PDGCode::MUON && pdg != PDGCode::TAU ) return -1; int result = 0; @@ -706,25 +707,25 @@ int get_lepton_origin(const edm4hep::MCParticleData &p, int pdg_parent = in.at(index).PDG ; // std::cout << " parent has pdg = " << in.at(index).PDG << " status = " << in.at(index).generatorStatus << std::endl; - if ( abs( pdg_parent ) == 23 || abs( pdg_parent ) == 24 ) { + if ( abs( pdg_parent ) == PDGCode::W || abs( pdg_parent ) == PDGCode::Z ) { result = pdg_parent ; //std::cout << " ... Lepton is from W or Z , return code = " << result << std::endl; break; } - if ( abs( pdg_parent ) == 22 ) { + if ( abs( pdg_parent ) == PDGCode::PHOTON ) { result = pdg_parent ; //std::cout << " ... Lepton is from a virtual photon , return code = " << result << std::endl; break; } - if ( abs( pdg_parent ) == 15 ) { + if ( abs( pdg_parent ) == PDGCode::TAU ) { result = pdg_parent ; //std::cout << " ... Lepton is from a tau, return code = " << result << std::endl; break; } - if ( abs( pdg_parent ) == 11 ) { // beam particle ? + if ( abs( pdg_parent ) == PDGCode::ELECTRON ) { // beam particle ? // beam particles should have generatorStatus = 4, // but that is not the case in files produced from Whizard + p6 if ( in.at(index).generatorStatus == 4 || ind.at ( in.at(index).parents_begin ) == 0 ) { @@ -734,7 +735,7 @@ int get_lepton_origin(const edm4hep::MCParticleData &p, } } - if ( pdg == 11 && abs( pdg_parent ) == 13 ) { // mu -> e + if ( pdg == PDGCode::ELECTRON && abs( pdg_parent ) == PDGCode::MUON ) { // mu -> e result = pdg_parent; //std::cout << " ... Electron from a muon decay, return code = " << result << std::endl; break; From 0af032e5fe324242918d477bae11c3b0d9b01c72 Mon Sep 17 00:00:00 2001 From: mattiaisgro Date: Wed, 18 Jun 2025 14:09:23 +0200 Subject: [PATCH 2/2] Change PDGCode from scoped enum to C enum --- .../dataframe/FCCAnalyses/MCParticleCodes.h | 46 +++++++++---------- analyzers/dataframe/src/JetTaggingUtils.cc | 40 ++++++++-------- analyzers/dataframe/src/MCParticle.cc | 14 +++--- 3 files changed, 48 insertions(+), 52 deletions(-) diff --git a/analyzers/dataframe/FCCAnalyses/MCParticleCodes.h b/analyzers/dataframe/FCCAnalyses/MCParticleCodes.h index da0730b257a..6b0a7b4eadb 100644 --- a/analyzers/dataframe/FCCAnalyses/MCParticleCodes.h +++ b/analyzers/dataframe/FCCAnalyses/MCParticleCodes.h @@ -6,38 +6,34 @@ namespace FCCAnalyses { /// Particle codes following the PDG Monte Carlo particle numbering scheme. - enum class PDGCode { + enum PDGCode : int { - UNKNOWN = 0, + PDG_UNKNOWN = 0, // Leptons - ELECTRON = 11, - MUON = 13, - TAU = 15, - E_NEUTRINO = 12, - MU_NEUTRINO = 14, - TAU_NEUTRINO = 16, + PDG_ELECTRON = 11, + PDG_MUON = 13, + PDG_TAU = 15, + PDG_E_NEUTRINO = 12, + PDG_MU_NEUTRINO = 14, + PDG_TAU_NEUTRINO = 16, // Partons - QUARK_UP = 1, - QUARK_DOWN = 2, - QUARK_CHARM = 4, - QUARK_STRANGE = 3, - QUARK_TOP = 6, - QUARK_BOTTOM = 5, - GLUON = 21, + PDG_QUARK_UP = 1, + PDG_QUARK_DOWN = 2, + PDG_QUARK_CHARM = 4, + PDG_QUARK_STRANGE = 3, + PDG_QUARK_TOP = 6, + PDG_QUARK_BOTTOM = 5, + PDG_GLUON = 21, // Bosons - PHOTON = 22, - Z = 23, - W = 24, - HIGGS = 25, - - // Baryons - PION = 211, - KAON = 321, - PROTON = 2212, - NEUTRON = 2112, + PDG_PHOTON = 22, + PDG_Z = 23, + PDG_W = 24, + PDG_HIGGS = 25, + + // ... }; } // namespace FCCAnalyses diff --git a/analyzers/dataframe/src/JetTaggingUtils.cc b/analyzers/dataframe/src/JetTaggingUtils.cc index c874295d884..fc560342e47 100644 --- a/analyzers/dataframe/src/JetTaggingUtils.cc +++ b/analyzers/dataframe/src/JetTaggingUtils.cc @@ -1,5 +1,5 @@ #include "FCCAnalyses/JetTaggingUtils.h" -#include "FCCAnalyses/PDGCodes.h" +#include "FCCAnalyses/MCParticleCodes.h" namespace FCCAnalyses { @@ -20,7 +20,7 @@ get_flavour(ROOT::VecOps::RVec in, parton.generatorStatus != 2) continue; - if (std::abs(parton.PDG) > PDGCode::QUARK_BOTTOM && parton.PDG != PDGCode::GLUON) + if (std::abs(parton.PDG) > PDG_QUARK_BOTTOM && parton.PDG != PDG_GLUON) continue; ROOT::Math::PxPyPzMVector lv(parton.momentum.x, parton.momentum.y, @@ -45,13 +45,13 @@ get_flavour(ROOT::VecOps::RVec in, if (angle <= 0.3) { - if (result[j] == PDGCode::GLUON or result[j] == PDGCode::UNKNOWN) { + if (result[j] == PDG_GLUON or result[j] == PDG_UNKNOWN) { // if no match before, or matched to gluon, match to // this particle (favour quarks over gluons) result[j] = std::abs(parton.PDG); - } else if (parton.PDG != PDGCode::GLUON) { + } else if (parton.PDG != PDG_GLUON) { // if matched to quark, and this is a quark, favour // heavier flavours @@ -77,16 +77,16 @@ ROOT::VecOps::RVec get_btag(ROOT::VecOps::RVec in, float efficiency, for (size_t j = 0; j < in.size(); ++j) { - if (in.at(j) == PDGCode::QUARK_BOTTOM && gRandom->Uniform() <= efficiency) + if (in.at(j) == PDG_QUARK_BOTTOM && gRandom->Uniform() <= efficiency) result[j] = 1; - if (in.at(j) == PDGCode::QUARK_CHARM && gRandom->Uniform() <= mistag_c) + if (in.at(j) == PDG_QUARK_CHARM && gRandom->Uniform() <= mistag_c) result[j] = 1; - if (in.at(j) < PDGCode::QUARK_CHARM && gRandom->Uniform() <= mistag_l) + if (in.at(j) < PDG_QUARK_CHARM && gRandom->Uniform() <= mistag_l) result[j] = 1; - if (in.at(j) == PDGCode::GLUON && gRandom->Uniform() <= mistag_g) + if (in.at(j) == PDG_GLUON && gRandom->Uniform() <= mistag_g) result[j] = 1; } @@ -101,16 +101,16 @@ ROOT::VecOps::RVec get_ctag(ROOT::VecOps::RVec in, float efficiency, for (size_t j = 0; j < in.size(); ++j) { - if (in.at(j) == PDGCode::QUARK_CHARM && gRandom->Uniform() <= efficiency) + if (in.at(j) == PDG_QUARK_CHARM && gRandom->Uniform() <= efficiency) result[j] = 1; - if (in.at(j) == PDGCode::QUARK_BOTTOM && gRandom->Uniform() <= mistag_b) + if (in.at(j) == PDG_QUARK_BOTTOM && gRandom->Uniform() <= mistag_b) result[j] = 1; - if (in.at(j) < PDGCode::QUARK_CHARM && gRandom->Uniform() <= mistag_l) + if (in.at(j) < PDG_QUARK_CHARM && gRandom->Uniform() <= mistag_l) result[j] = 1; - if (in.at(j) == PDGCode::GLUON && gRandom->Uniform() <= mistag_g) + if (in.at(j) == PDG_GLUON && gRandom->Uniform() <= mistag_g) result[j] = 1; } return result; @@ -124,16 +124,16 @@ ROOT::VecOps::RVec get_ltag(ROOT::VecOps::RVec in, float efficiency, for (size_t j = 0; j < in.size(); ++j) { - if (in.at(j) < PDGCode::QUARK_CHARM && gRandom->Uniform() <= efficiency) + if (in.at(j) < PDG_QUARK_CHARM && gRandom->Uniform() <= efficiency) result[j] = 1; - if (in.at(j) == PDGCode::QUARK_BOTTOM && gRandom->Uniform() <= mistag_b) + if (in.at(j) == PDG_QUARK_BOTTOM && gRandom->Uniform() <= mistag_b) result[j] = 1; - if (in.at(j) == PDGCode::QUARK_CHARM && gRandom->Uniform() <= mistag_c) + if (in.at(j) == PDG_QUARK_CHARM && gRandom->Uniform() <= mistag_c) result[j] = 1; - if (in.at(j) == PDGCode::GLUON && gRandom->Uniform() <= mistag_g) + if (in.at(j) == PDG_GLUON && gRandom->Uniform() <= mistag_g) result[j] = 1; } return result; @@ -147,16 +147,16 @@ ROOT::VecOps::RVec get_gtag(ROOT::VecOps::RVec in, float efficiency, for (size_t j = 0; j < in.size(); ++j) { - if (in.at(j) == PDGCode::GLUON && gRandom->Uniform() <= efficiency) + if (in.at(j) == PDG_GLUON && gRandom->Uniform() <= efficiency) result[j] = 1; - if (in.at(j) == PDGCode::QUARK_BOTTOM && gRandom->Uniform() <= mistag_b) + if (in.at(j) == PDG_QUARK_BOTTOM && gRandom->Uniform() <= mistag_b) result[j] = 1; - if (in.at(j) == PDGCode::QUARK_CHARM && gRandom->Uniform() <= mistag_c) + if (in.at(j) == PDG_QUARK_CHARM && gRandom->Uniform() <= mistag_c) result[j] = 1; - if (in.at(j) < PDGCode::QUARK_CHARM && gRandom->Uniform() <= mistag_l) + if (in.at(j) < PDG_QUARK_CHARM && gRandom->Uniform() <= mistag_l) result[j] = 1; } return result; diff --git a/analyzers/dataframe/src/MCParticle.cc b/analyzers/dataframe/src/MCParticle.cc index f5b05c3d2cf..b8f59977608 100644 --- a/analyzers/dataframe/src/MCParticle.cc +++ b/analyzers/dataframe/src/MCParticle.cc @@ -1,5 +1,5 @@ #include "FCCAnalyses/MCParticle.h" -#include "FCCAnalyses/PDGCodes.h" +#include "FCCAnalyses/MCParticleCodes.h" #include #include #include @@ -697,7 +697,7 @@ int get_lepton_origin(const edm4hep::MCParticleData &p, // std::cout << std::endl << " enter in MCParticle::get_lepton_origin PDG = " << p.PDG << std::endl; int pdg = std::abs( p.PDG ) ; - if ( pdg != PDGCode::ELECTRON && pdg != PDGCode::MUON && pdg != PDGCode::TAU ) return -1; + if ( pdg != PDG_ELECTRON && pdg != PDG_MUON && pdg != PDG_TAU ) return -1; int result = 0; @@ -707,25 +707,25 @@ int get_lepton_origin(const edm4hep::MCParticleData &p, int pdg_parent = in.at(index).PDG ; // std::cout << " parent has pdg = " << in.at(index).PDG << " status = " << in.at(index).generatorStatus << std::endl; - if ( abs( pdg_parent ) == PDGCode::W || abs( pdg_parent ) == PDGCode::Z ) { + if ( abs( pdg_parent ) == PDG_W || abs( pdg_parent ) == PDG_Z ) { result = pdg_parent ; //std::cout << " ... Lepton is from W or Z , return code = " << result << std::endl; break; } - if ( abs( pdg_parent ) == PDGCode::PHOTON ) { + if ( abs( pdg_parent ) == PDG_PHOTON ) { result = pdg_parent ; //std::cout << " ... Lepton is from a virtual photon , return code = " << result << std::endl; break; } - if ( abs( pdg_parent ) == PDGCode::TAU ) { + if ( abs( pdg_parent ) == PDG_TAU ) { result = pdg_parent ; //std::cout << " ... Lepton is from a tau, return code = " << result << std::endl; break; } - if ( abs( pdg_parent ) == PDGCode::ELECTRON ) { // beam particle ? + if ( abs( pdg_parent ) == PDG_ELECTRON ) { // beam particle ? // beam particles should have generatorStatus = 4, // but that is not the case in files produced from Whizard + p6 if ( in.at(index).generatorStatus == 4 || ind.at ( in.at(index).parents_begin ) == 0 ) { @@ -735,7 +735,7 @@ int get_lepton_origin(const edm4hep::MCParticleData &p, } } - if ( pdg == PDGCode::ELECTRON && abs( pdg_parent ) == PDGCode::MUON ) { // mu -> e + if ( pdg == PDG_ELECTRON && abs( pdg_parent ) == PDG_MUON ) { // mu -> e result = pdg_parent; //std::cout << " ... Electron from a muon decay, return code = " << result << std::endl; break;