diff --git a/analyzers/dataframe/FCCAnalyses/MCParticleCodes.h b/analyzers/dataframe/FCCAnalyses/MCParticleCodes.h new file mode 100644 index 00000000000..6b0a7b4eadb --- /dev/null +++ b/analyzers/dataframe/FCCAnalyses/MCParticleCodes.h @@ -0,0 +1,42 @@ + +#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 PDGCode : int { + + PDG_UNKNOWN = 0, + + // Leptons + PDG_ELECTRON = 11, + PDG_MUON = 13, + PDG_TAU = 15, + PDG_E_NEUTRINO = 12, + PDG_MU_NEUTRINO = 14, + PDG_TAU_NEUTRINO = 16, + + // Partons + 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 + PDG_PHOTON = 22, + PDG_Z = 23, + PDG_W = 24, + PDG_HIGGS = 25, + + // ... + }; + +} // namespace FCCAnalyses + + +#endif diff --git a/analyzers/dataframe/src/JetTaggingUtils.cc b/analyzers/dataframe/src/JetTaggingUtils.cc index 7ed529c18a3..fc560342e47 100644 --- a/analyzers/dataframe/src/JetTaggingUtils.cc +++ b/analyzers/dataframe/src/JetTaggingUtils.cc @@ -1,4 +1,5 @@ #include "FCCAnalyses/JetTaggingUtils.h" +#include "FCCAnalyses/MCParticleCodes.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) > PDG_QUARK_BOTTOM && parton.PDG != PDG_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] == 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 != 21) { + + } else if (parton.PDG != PDG_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) == PDG_QUARK_BOTTOM && gRandom->Uniform() <= efficiency) result[j] = 1; - if (in.at(j) == 4 && gRandom->Uniform() <= mistag_c) + + if (in.at(j) == PDG_QUARK_CHARM && gRandom->Uniform() <= mistag_c) result[j] = 1; - if (in.at(j) < 4 && gRandom->Uniform() <= mistag_l) + + if (in.at(j) < PDG_QUARK_CHARM && gRandom->Uniform() <= mistag_l) result[j] = 1; - if (in.at(j) == 21 && gRandom->Uniform() <= mistag_g) + + if (in.at(j) == PDG_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) == PDG_QUARK_CHARM && gRandom->Uniform() <= efficiency) result[j] = 1; - if (in.at(j) == 5 && gRandom->Uniform() <= mistag_b) + + if (in.at(j) == PDG_QUARK_BOTTOM && gRandom->Uniform() <= mistag_b) result[j] = 1; - if (in.at(j) < 4 && gRandom->Uniform() <= mistag_l) + + if (in.at(j) < PDG_QUARK_CHARM && gRandom->Uniform() <= mistag_l) result[j] = 1; - if (in.at(j) == 21 && gRandom->Uniform() <= mistag_g) + + if (in.at(j) == PDG_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) < PDG_QUARK_CHARM && gRandom->Uniform() <= efficiency) result[j] = 1; - if (in.at(j) == 5 && gRandom->Uniform() <= mistag_b) + + if (in.at(j) == PDG_QUARK_BOTTOM && gRandom->Uniform() <= mistag_b) result[j] = 1; - if (in.at(j) == 4 && gRandom->Uniform() <= mistag_c) + + if (in.at(j) == PDG_QUARK_CHARM && gRandom->Uniform() <= mistag_c) result[j] = 1; - if (in.at(j) == 21 && gRandom->Uniform() <= mistag_g) + + if (in.at(j) == PDG_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) == PDG_GLUON && gRandom->Uniform() <= efficiency) result[j] = 1; - if (in.at(j) == 5 && gRandom->Uniform() <= mistag_b) + + if (in.at(j) == PDG_QUARK_BOTTOM && gRandom->Uniform() <= mistag_b) result[j] = 1; - if (in.at(j) == 4 && gRandom->Uniform() <= mistag_c) + + if (in.at(j) == PDG_QUARK_CHARM && gRandom->Uniform() <= mistag_c) result[j] = 1; - if (in.at(j) < 4 && 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 009d0d4c255..b8f59977608 100644 --- a/analyzers/dataframe/src/MCParticle.cc +++ b/analyzers/dataframe/src/MCParticle.cc @@ -1,4 +1,5 @@ #include "FCCAnalyses/MCParticle.h" +#include "FCCAnalyses/MCParticleCodes.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 != PDG_ELECTRON && pdg != PDG_MUON && pdg != PDG_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 ) == 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 ) == 22 ) { + 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 ) == 15 ) { + 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 ) == 11 ) { // 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 ) { @@ -734,7 +735,7 @@ int get_lepton_origin(const edm4hep::MCParticleData &p, } } - if ( pdg == 11 && abs( pdg_parent ) == 13 ) { // 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;