HepMCCandAlgos/interface/MCTruthHelper.h

0001 #ifndef PhysicsTools_HepMCCandAlgos_GenParticlesHelper_h
0002 #define PhysicsTools_HepMCCandAlgos_GenParticlesHelper_h
0003
0004 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0005 #include "HepPDT/ParticleID.hh"
0006 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0007 #include "DataFormats/HepMCCandidate/interface/GenStatusFlags.h"
0008
0009 #include "HepMC3/GenEvent.h"
0010 #include "HepMC3/GenParticle.h"
0011 #include "HepMC3/GenVertex.h"
0012
0013 #include <iostream>
0014 #include <unordered_set>
0015
0016 template <typename P>
0017 class MCTruthHelper {
0018 public:
0019   /////////////////////////////////////////////////////////////////////////////
0020   //these are robust, generator-independent functions for categorizing
0021   //mainly final state particles, but also intermediate hadrons
0022   //or radiating leptons
0023
0024   //is particle prompt (not from hadron, muon, or tau decay)
0025   bool isPrompt(const P &p) const;
0026
0027   //is particle prompt and final state
0028   bool isPromptFinalState(const P &p) const;
0029
0030   //is particle a decayed hadron, muon, or tau (does not include resonance decays like W,Z,Higgs,top,etc)
0031   //This flag is equivalent to status 2 in the current HepMC standard
0032   //but older generators (pythia6, herwig6) predate this and use status 2 also for other intermediate
0033   //particles/states
0034   bool isDecayedLeptonHadron(const P &p) const;
0035
0036   //is particle prompt and decayed
0037   bool isPromptDecayed(const P &p) const;
0038
0039   //this particle is a direct or indirect tau decay product
0040   bool isTauDecayProduct(const P &p) const;
0041
0042   //this particle is a direct or indirect decay product of a prompt tau
0043   bool isPromptTauDecayProduct(const P &p) const;
0044
0045   //this particle is a direct tau decay product
0046   bool isDirectTauDecayProduct(const P &p) const;
0047
0048   //this particle is a direct decay product from a prompt tau
0049   bool isDirectPromptTauDecayProduct(const P &p) const;
0050
0051   //this particle is a direct or indirect muon decay product
0052   bool isMuonDecayProduct(const P &p) const;
0053
0054   //this particle is a direct or indirect decay product of a prompt muon
0055   bool isPromptMuonDecayProduct(const P &p) const;
0056
0057   //this particle is a direct decay product from a hadron
0058   bool isDirectHadronDecayProduct(const P &p) const;
0059
0060   //is particle a hadron
0061   bool isHadron(const P &p) const;
0062
0063   /////////////////////////////////////////////////////////////////////////////
0064   //these are generator history-dependent functions for tagging particles
0065   //associated with the hard process
0066   //Currently implemented for Pythia 6 and Pythia 8 status codes and history
0067
0068   //this particle is part of the hard process
0069   bool isHardProcess(const P &p) const;
0070
0071   //this particle is the direct descendant of a hard process particle of the same pdg id
0072   bool fromHardProcess(const P &p) const;
0073
0074   //this particle is the final state direct descendant of a hard process particle
0075   bool fromHardProcessFinalState(const P &p) const;
0076
0077   //this particle is the decayed direct descendant of a hard process particle
0078   //such as a tau from the hard process
0079   bool fromHardProcessDecayed(const P &p) const;
0080
0081   //this particle is a direct or indirect decay product of a tau
0082   //from the hard process
0083   bool isHardProcessTauDecayProduct(const P &p) const;
0084
0085   //this particle is a direct decay product of a tau
0086   //from the hard process
0087   bool isDirectHardProcessTauDecayProduct(const P &p) const;
0088
0089   //this particle is the direct descendant of a hard process particle of the same pdg id
0090   //For outgoing particles the kinematics are those before QCD or QED FSR
0091   //This corresponds roughly to status code 3 in pythia 6
0092   bool fromHardProcessBeforeFSR(const P &p) const;
0093
0094   //this particle is the first copy of the particle in the chain with the same pdg id
0095   bool isFirstCopy(const P &p) const;
0096
0097   //this particle is the last copy of the particle in the chain with the same pdg id
0098   //(and therefore is more likely, but not guaranteed, to carry the final physical momentum)
0099   bool isLastCopy(const P &p) const;
0100
0101   //this particle is the last copy of the particle in the chain with the same pdg id
0102   //before QED or QCD FSR
0103   //(and therefore is more likely, but not guaranteed, to carry the momentum after ISR)
0104   //This flag only really makes sense for outgoing particles
0105   bool isLastCopyBeforeFSR(const P &p) const;
0106
0107   /////////////////////////////////////////////////////////////////////////////
0108   //These are utility functions used by the above
0109
0110   //first mother in chain with a different pdg than the particle
0111   const P *uniqueMother(const P &p) const;
0112
0113   //return first copy of particle in chain (may be the particle itself)
0114   const P *firstCopy(const P &p) const;
0115
0116   //return last copy of particle in chain (may be the particle itself)
0117   const P *lastCopy(const P &p) const;
0118
0119   //return last copy of particle in chain before QED or QCD FSR (may be the particle itself)
0120   const P *lastCopyBeforeFSR(const P &p) const;
0121
0122   //return last copy of particle in chain before QED or QCD FSR, starting from itself (may be the particle itself)
0123   const P *lastDaughterCopyBeforeFSR(const P &p) const;
0124
0125   //return mother copy which is a hard process particle
0126   const P *hardProcessMotherCopy(const P &p) const;
0127
0128   //return previous copy of particle in chain (0 in case this is already the first copy)
0129   const P *previousCopy(const P &p) const;
0130
0131   //return next copy of particle in chain (0 in case this is already the last copy)
0132   const P *nextCopy(const P &p) const;
0133
0134   //return decayed mother (walk up the chain until found)
0135   const P *findDecayedMother(const P &p) const;
0136
0137   //return decayed mother matching requested abs(pdgid) (walk up the chain until found)
0138   const P *findDecayedMother(const P &p, int abspdgid) const;
0139
0140   /////////////////////////////////////////////////////////////////////////////
0141   //These are very basic utility functions to implement a common interface for reco::GenParticle
0142   //and HepMC::GenParticle
0143
0144   //pdgid
0145   int pdgId(const reco::GenParticle &p) const;
0146
0147   //pdgid
0148   int pdgId(const HepMC::GenParticle &p) const;
0149
0150   //pdgid
0151   int pdgId(const HepMC3::GenParticle &p) const;
0152
0153   //abs(pdgid)
0154   int absPdgId(const reco::GenParticle &p) const;
0155
0156   //abs(pdgid)
0157   int absPdgId(const HepMC::GenParticle &p) const;
0158
0159   //abs(pdgid)
0160   int absPdgId(const HepMC3::GenParticle &p) const;
0161
0162   //number of mothers
0163   unsigned int numberOfMothers(const reco::GenParticle &p) const;
0164
0165   //number of mothers
0166   unsigned int numberOfMothers(const HepMC::GenParticle &p) const;
0167
0168   //number of mothers
0169   unsigned int numberOfMothers(const HepMC3::GenParticle &p) const;
0170
0171   //mother
0172   const reco::GenParticle *mother(const reco::GenParticle &p, unsigned int imoth = 0) const;
0173
0174   //mother
0175   const HepMC::GenParticle *mother(const HepMC::GenParticle &p, unsigned int imoth = 0) const;
0176
0177   //mother
0178   const HepMC3::GenParticle *mother(const HepMC3::GenParticle &p, unsigned int imoth = 0) const;
0179
0180   //number of daughters
0181   unsigned int numberOfDaughters(const reco::GenParticle &p) const;
0182
0183   //number of daughters
0184   unsigned int numberOfDaughters(const HepMC::GenParticle &p) const;
0185
0186   //number of daughters
0187   unsigned int numberOfDaughters(const HepMC3::GenParticle &p) const;
0188
0189   //ith daughter
0190   const reco::GenParticle *daughter(const reco::GenParticle &p, unsigned int idau) const;
0191
0192   //ith daughter
0193   const HepMC::GenParticle *daughter(const HepMC::GenParticle &p, unsigned int idau) const;
0194
0195   //ith daughter
0196   const HepMC3::GenParticle *daughter(const HepMC3::GenParticle &p, unsigned int idau) const;
0197
0198   /////////////////////////////////////////////////////////////////////////////
0199   //Helper function to fill status flags
0200   void fillGenStatusFlags(const P &p, reco::GenStatusFlags &statusFlags) const;
0201 };
0202
0203 /////////////////////////////////////////////////////////////////////////////
0204 //implementations
0205
0206 /////////////////////////////////////////////////////////////////////////////
0207 template <typename P>
0208 bool MCTruthHelper<P>::isPrompt(const P &p) const {
0209   //particle from hadron/muon/tau decay -> not prompt
0210   //checking all the way up the chain treats properly the radiated photon
0211   //case as well
0212   return findDecayedMother(p) == nullptr;
0213 }
0214
0215 /////////////////////////////////////////////////////////////////////////////
0216 template <typename P>
0217 bool MCTruthHelper<P>::isPromptFinalState(const P &p) const {
0218   return p.status() == 1 && isPrompt(p);
0219 }
0220
0221 template <typename P>
0222 bool MCTruthHelper<P>::isDecayedLeptonHadron(const P &p) const {
0223   return p.status() == 2 && (isHadron(p) || absPdgId(p) == 13 || absPdgId(p) == 15) && isLastCopy(p);
0224 }
0225
0226 template <typename P>
0227 bool MCTruthHelper<P>::isPromptDecayed(const P &p) const {
0228   return isDecayedLeptonHadron(p) && isPrompt(p);
0229 }
0230
0231 /////////////////////////////////////////////////////////////////////////////
0232 template <typename P>
0233 bool MCTruthHelper<P>::isTauDecayProduct(const P &p) const {
0234   return findDecayedMother(p, 15) != nullptr;
0235 }
0236
0237 /////////////////////////////////////////////////////////////////////////////
0238 template <typename P>
0239 bool MCTruthHelper<P>::isPromptTauDecayProduct(const P &p) const {
0240   const P *tau = findDecayedMother(p, 15);
0241   return tau && isPrompt(*tau);
0242 }
0243
0244 /////////////////////////////////////////////////////////////////////////////
0245 template <typename P>
0246 bool MCTruthHelper<P>::isDirectTauDecayProduct(const P &p) const {
0247   const P *tau = findDecayedMother(p, 15);
0248   const P *dm = findDecayedMother(p);
0249   return tau && tau == dm;
0250 }
0251
0252 /////////////////////////////////////////////////////////////////////////////
0253 template <typename P>
0254 bool MCTruthHelper<P>::isDirectPromptTauDecayProduct(const P &p) const {
0255   const P *tau = findDecayedMother(p, 15);
0256   const P *dm = findDecayedMother(p);
0257   return tau && tau == dm && isPrompt(*tau);
0258 }
0259
0260 /////////////////////////////////////////////////////////////////////////////
0261 template <typename P>
0262 bool MCTruthHelper<P>::isMuonDecayProduct(const P &p) const {
0263   return findDecayedMother(p, 13) != 0;
0264 }
0265
0266 /////////////////////////////////////////////////////////////////////////////
0267 template <typename P>
0268 bool MCTruthHelper<P>::isPromptMuonDecayProduct(const P &p) const {
0269   const P *mu = findDecayedMother(p, 13);
0270   return mu && isPrompt(*mu);
0271 }
0272
0273 /////////////////////////////////////////////////////////////////////////////
0274 template <typename P>
0275 bool MCTruthHelper<P>::isDirectHadronDecayProduct(const P &p) const {
0276   const P *um = uniqueMother(p);
0277   return um && isHadron(*um) && isDecayedLeptonHadron(*um);
0278 }
0279
0280 /////////////////////////////////////////////////////////////////////////////
0281 template <typename P>
0282 bool MCTruthHelper<P>::isHadron(const P &p) const {
0283   HepPDT::ParticleID heppdtid(pdgId(p));
0284   return heppdtid.isHadron();
0285 }
0286
0287 /////////////////////////////////////////////////////////////////////////////
0288 template <typename P>
0289 bool MCTruthHelper<P>::isHardProcess(const P &p) const {
0290   //status 3 in pythia6 means hard process;
0291   if (p.status() == 3)
0292     return true;
0293
0294   //hard process codes for pythia8 are 21-29 inclusive (currently 21,22,23,24 are used)
0295   if (p.status() > 20 && p.status() < 30)
0296     return true;
0297
0298   //if this is a final state or decayed particle,
0299   //check if direct mother is a resonance decay in pythia8 but exclude FSR branchings
0300   //(In pythia8 if a resonance decay product did not undergo any further branchings
0301   //it will be directly stored as status 1 or 2 without any status 23 copy)
0302   if (p.status() == 1 || p.status() == 2) {
0303     const P *um = mother(p);
0304     if (um) {
0305       const P *firstcopy = firstCopy(*um);
0306       bool fromResonance = firstcopy && firstcopy->status() == 22;
0307
0308       const P *umNext = nextCopy(*um);
0309       bool fsrBranching = umNext && umNext->status() > 50 && umNext->status() < 60;
0310
0311       if (fromResonance && !fsrBranching)
0312         return true;
0313     }
0314   }
0315
0316   return false;
0317 }
0318
0319 /////////////////////////////////////////////////////////////////////////////
0320 template <typename P>
0321 bool MCTruthHelper<P>::fromHardProcess(const P &p) const {
0322   return hardProcessMotherCopy(p) != nullptr;
0323 }
0324
0325 /////////////////////////////////////////////////////////////////////////////
0326 template <typename P>
0327 bool MCTruthHelper<P>::fromHardProcessFinalState(const P &p) const {
0328   return p.status() == 1 && fromHardProcess(p);
0329 }
0330
0331 /////////////////////////////////////////////////////////////////////////////
0332 template <typename P>
0333 bool MCTruthHelper<P>::fromHardProcessDecayed(const P &p) const {
0334   return isDecayedLeptonHadron(p) && fromHardProcess(p);
0335 }
0336
0337 /////////////////////////////////////////////////////////////////////////////
0338 template <typename P>
0339 bool MCTruthHelper<P>::isHardProcessTauDecayProduct(const P &p) const {
0340   const P *tau = findDecayedMother(p, 15);
0341   return tau && fromHardProcessDecayed(*tau);
0342 }
0343
0344 /////////////////////////////////////////////////////////////////////////////
0345 template <typename P>
0346 bool MCTruthHelper<P>::isDirectHardProcessTauDecayProduct(const P &p) const {
0347   const P *tau = findDecayedMother(p, 15);
0348   const P *dm = findDecayedMother(p);
0349   return tau && tau == dm && fromHardProcess(*tau);
0350 }
0351
0352 template <typename P>
0353 bool MCTruthHelper<P>::fromHardProcessBeforeFSR(const P &p) const {
0354   //pythia 6 documentation line roughly corresponds to this condition
0355   if (p.status() == 3)
0356     return true;
0357
0358   //check hard process mother properties
0359   const P *hpc = hardProcessMotherCopy(p);
0360   if (!hpc)
0361     return false;
0362
0363   //for incoming partons in pythia8, more useful information is not
0364   //easily available, so take only the incoming parton itself
0365   if (hpc->status() == 21 && (&p) == hpc)
0366     return true;
0367
0368   //for intermediate particles in pythia 8, just take the last copy
0369   if (hpc->status() == 22 && isLastCopy(p))
0370     return true;
0371
0372   //for outgoing particles in pythia 8, explicitly find the last copy
0373   //before FSR starting from the hardProcess particle, and take only
0374   //this one
0375   if ((hpc->status() == 23 || hpc->status() == 1) && (&p) == lastDaughterCopyBeforeFSR(*hpc))
0376     return true;
0377
0378   //didn't satisfy any of the conditions
0379   return false;
0380 }
0381
0382 /////////////////////////////////////////////////////////////////////////////
0383 template <typename P>
0384 bool MCTruthHelper<P>::isFirstCopy(const P &p) const {
0385   return &p == firstCopy(p);
0386 }
0387
0388 /////////////////////////////////////////////////////////////////////////////
0389 template <typename P>
0390 bool MCTruthHelper<P>::isLastCopy(const P &p) const {
0391   return &p == lastCopy(p);
0392 }
0393
0394 /////////////////////////////////////////////////////////////////////////////
0395 template <typename P>
0396 bool MCTruthHelper<P>::isLastCopyBeforeFSR(const P &p) const {
0397   return &p == lastCopyBeforeFSR(p);
0398 }
0399
0400 /////////////////////////////////////////////////////////////////////////////
0401 template <typename P>
0402 const P *MCTruthHelper<P>::uniqueMother(const P &p) const {
0403   const P *mo = &p;
0404   std::unordered_set<const P *> dupCheck;
0405   while (mo && pdgId(*mo) == pdgId(p)) {
0406     dupCheck.insert(mo);
0407     mo = mother(*mo);
0408     if (dupCheck.count(mo))
0409       return nullptr;
0410   }
0411   return mo;
0412 }
0413
0414 /////////////////////////////////////////////////////////////////////////////
0415 template <typename P>
0416 const P *MCTruthHelper<P>::firstCopy(const P &p) const {
0417   const P *pcopy = &p;
0418   std::unordered_set<const P *> dupCheck;
0419   while (previousCopy(*pcopy)) {
0420     dupCheck.insert(pcopy);
0421     pcopy = previousCopy(*pcopy);
0422     if (dupCheck.count(pcopy))
0423       return nullptr;
0424   }
0425   return pcopy;
0426 }
0427
0428 /////////////////////////////////////////////////////////////////////////////
0429 template <typename P>
0430 const P *MCTruthHelper<P>::lastCopy(const P &p) const {
0431   const P *pcopy = &p;
0432   std::unordered_set<const P *> dupCheck;
0433   while (nextCopy(*pcopy)) {
0434     dupCheck.insert(pcopy);
0435     pcopy = nextCopy(*pcopy);
0436     if (dupCheck.count(pcopy))
0437       return nullptr;
0438   }
0439   return pcopy;
0440 }
0441
0442 /////////////////////////////////////////////////////////////////////////////
0443 template <typename P>
0444 const P *MCTruthHelper<P>::lastCopyBeforeFSR(const P &p) const {
0445   //start with first copy and then walk down until there is FSR
0446   const P *pcopy = firstCopy(p);
0447   if (!pcopy)
0448     return nullptr;
0449   bool hasDaughterCopy = true;
0450   std::unordered_set<const P *> dupCheck;
0451   while (hasDaughterCopy) {
0452     dupCheck.insert(pcopy);
0453     hasDaughterCopy = false;
0454     const unsigned int ndau = numberOfDaughters(*pcopy);
0455     //look for FSR
0456     for (unsigned int idau = 0; idau < ndau; ++idau) {
0457       const P *dau = daughter(*pcopy, idau);
0458       if (pdgId(*dau) == 21 || pdgId(*dau) == 22) {
0459         //has fsr (or else decayed and is the last copy by construction)
0460         return pcopy;
0461       }
0462     }
0463     //look for daughter copy
0464     for (unsigned int idau = 0; idau < ndau; ++idau) {
0465       const P *dau = daughter(*pcopy, idau);
0466       if (pdgId(*dau) == pdgId(p)) {
0467         pcopy = dau;
0468         hasDaughterCopy = true;
0469         break;
0470       }
0471     }
0472     if (dupCheck.count(pcopy))
0473       return nullptr;
0474   }
0475   return pcopy;
0476 }
0477
0478 /////////////////////////////////////////////////////////////////////////////
0479 template <typename P>
0480 const P *MCTruthHelper<P>::lastDaughterCopyBeforeFSR(const P &p) const {
0481   //start with this particle and then walk down until there is FSR
0482   const P *pcopy = &p;
0483   bool hasDaughterCopy = true;
0484   std::unordered_set<const P *> dupCheck;
0485   while (hasDaughterCopy) {
0486     dupCheck.insert(pcopy);
0487     hasDaughterCopy = false;
0488     const unsigned int ndau = numberOfDaughters(*pcopy);
0489     //look for FSR
0490     for (unsigned int idau = 0; idau < ndau; ++idau) {
0491       const P *dau = daughter(*pcopy, idau);
0492       if (pdgId(*dau) == 21 || pdgId(*dau) == 22) {
0493         //has fsr (or else decayed and is the last copy by construction)
0494         return pcopy;
0495       }
0496     }
0497     //look for daughter copy
0498     for (unsigned int idau = 0; idau < ndau; ++idau) {
0499       const P *dau = daughter(*pcopy, idau);
0500       if (pdgId(*dau) == pdgId(p)) {
0501         pcopy = dau;
0502         hasDaughterCopy = true;
0503         break;
0504       }
0505     }
0506     if (dupCheck.count(pcopy))
0507       return nullptr;
0508   }
0509   return pcopy;
0510 }
0511
0512 /////////////////////////////////////////////////////////////////////////////
0513 template <typename P>
0514 const P *MCTruthHelper<P>::hardProcessMotherCopy(const P &p) const {
0515   //is particle itself is hard process particle
0516   if (isHardProcess(p))
0517     return &p;
0518
0519   //check if any other copies are hard process particles
0520   const P *pcopy = &p;
0521   std::unordered_set<const P *> dupCheck;
0522   while (previousCopy(*pcopy)) {
0523     dupCheck.insert(pcopy);
0524     pcopy = previousCopy(*pcopy);
0525     if (isHardProcess(*pcopy))
0526       return pcopy;
0527     if (dupCheck.count(pcopy))
0528       break;
0529   }
0530   return nullptr;
0531 }
0532
0533 /////////////////////////////////////////////////////////////////////////////
0534 template <typename P>
0535 const P *MCTruthHelper<P>::previousCopy(const P &p) const {
0536   const unsigned int nmoth = numberOfMothers(p);
0537   for (unsigned int imoth = 0; imoth < nmoth; ++imoth) {
0538     const P *moth = mother(p, imoth);
0539     if (pdgId(*moth) == pdgId(p)) {
0540       return moth;
0541     }
0542   }
0543
0544   return nullptr;
0545 }
0546
0547 /////////////////////////////////////////////////////////////////////////////
0548 template <typename P>
0549 const P *MCTruthHelper<P>::nextCopy(const P &p) const {
0550   const unsigned int ndau = numberOfDaughters(p);
0551   for (unsigned int idau = 0; idau < ndau; ++idau) {
0552     const P *dau = daughter(p, idau);
0553     if (pdgId(*dau) == pdgId(p)) {
0554       return dau;
0555     }
0556   }
0557
0558   return nullptr;
0559 }
0560
0561 /////////////////////////////////////////////////////////////////////////////
0562 template <typename P>
0563 const P *MCTruthHelper<P>::findDecayedMother(const P &p) const {
0564   const P *mo = mother(p);
0565   std::unordered_set<const P *> dupCheck;
0566   while (mo && !isDecayedLeptonHadron(*mo)) {
0567     dupCheck.insert(mo);
0568     mo = mother(*mo);
0569     if (dupCheck.count(mo))
0570       return nullptr;
0571   }
0572   return mo;
0573 }
0574
0575 /////////////////////////////////////////////////////////////////////////////
0576 template <typename P>
0577 const P *MCTruthHelper<P>::findDecayedMother(const P &p, int abspdgid) const {
0578   const P *mo = mother(p);
0579   std::unordered_set<const P *> dupCheck;
0580   while (mo && (absPdgId(*mo) != abspdgid || !isDecayedLeptonHadron(*mo))) {
0581     dupCheck.insert(mo);
0582     mo = mother(*mo);
0583     if (dupCheck.count(mo))
0584       return nullptr;
0585   }
0586   return mo;
0587 }
0588
0589 //////////////////////////////////////////////////////////////
0590 template <typename P>
0591 int MCTruthHelper<P>::pdgId(const reco::GenParticle &p) const {
0592   return p.pdgId();
0593 }
0594
0595 //////////////////////////////////////////////////////////////
0596 template <typename P>
0597 int MCTruthHelper<P>::pdgId(const HepMC::GenParticle &p) const {
0598   return p.pdg_id();
0599 }
0600
0601 //////////////////////////////////////////////////////////////
0602 template <typename P>
0603 int MCTruthHelper<P>::pdgId(const HepMC3::GenParticle &p) const {
0604   return p.pid();
0605 }
0606
0607 //////////////////////////////////////////////////////////////
0608 template <typename P>
0609 int MCTruthHelper<P>::absPdgId(const reco::GenParticle &p) const {
0610   return std::abs(p.pdgId());
0611 }
0612
0613 //////////////////////////////////////////////////////////////
0614 template <typename P>
0615 int MCTruthHelper<P>::absPdgId(const HepMC::GenParticle &p) const {
0616   return std::abs(p.pdg_id());
0617 }
0618
0619 //////////////////////////////////////////////////////////////
0620 template <typename P>
0621 int MCTruthHelper<P>::absPdgId(const HepMC3::GenParticle &p) const {
0622   return std::abs(p.pid());
0623 }
0624
0625 /////////////////////////////////////////////////////////////////////////////
0626 template <typename P>
0627 unsigned int MCTruthHelper<P>::numberOfMothers(const reco::GenParticle &p) const {
0628   return p.numberOfMothers();
0629 }
0630
0631 /////////////////////////////////////////////////////////////////////////////
0632 template <typename P>
0633 unsigned int MCTruthHelper<P>::numberOfMothers(const HepMC::GenParticle &p) const {
0634   return p.production_vertex() ? p.production_vertex()->particles_in_size() : 0;
0635 }
0636
0637 /////////////////////////////////////////////////////////////////////////////
0638 template <typename P>
0639 unsigned int MCTruthHelper<P>::numberOfMothers(const HepMC3::GenParticle &p) const {
0640   return p.production_vertex() ? p.production_vertex()->particles_in_size() : 0;
0641 }
0642
0643 /////////////////////////////////////////////////////////////////////////////
0644 template <typename P>
0645 const reco::GenParticle *MCTruthHelper<P>::mother(const reco::GenParticle &p, unsigned int imoth) const {
0646   return static_cast<const reco::GenParticle *>(p.mother(imoth));
0647 }
0648
0649 /////////////////////////////////////////////////////////////////////////////
0650 template <typename P>
0651 const HepMC::GenParticle *MCTruthHelper<P>::mother(const HepMC::GenParticle &p, unsigned int imoth) const {
0652   return p.production_vertex() && p.production_vertex()->particles_in_size()
0653              ? *(p.production_vertex()->particles_in_const_begin() + imoth)
0654              : nullptr;
0655 }
0656
0657 /////////////////////////////////////////////////////////////////////////////
0658 template <typename P>
0659 const HepMC3::GenParticle *MCTruthHelper<P>::mother(const HepMC3::GenParticle &p, unsigned int imoth) const {
0660   if (numberOfMothers(p) > 0) {
0661     unsigned int imothi = 0;
0662     for (const HepMC3::ConstGenParticlePtr &mother : (p.production_vertex())->particles_in()) {
0663       if (imothi == imoth)
0664         return mother.get();
0665       imothi++;
0666     }
0667   }
0668   return nullptr;
0669 }
0670
0671 /////////////////////////////////////////////////////////////////////////////
0672 template <typename P>
0673 unsigned int MCTruthHelper<P>::numberOfDaughters(const reco::GenParticle &p) const {
0674   return p.numberOfDaughters();
0675 }
0676
0677 /////////////////////////////////////////////////////////////////////////////
0678 template <typename P>
0679 unsigned int MCTruthHelper<P>::numberOfDaughters(const HepMC::GenParticle &p) const {
0680   return p.end_vertex() ? p.end_vertex()->particles_out_size() : 0;
0681 }
0682
0683 /////////////////////////////////////////////////////////////////////////////
0684 template <typename P>
0685 unsigned int MCTruthHelper<P>::numberOfDaughters(const HepMC3::GenParticle &p) const {
0686   return p.end_vertex() ? p.end_vertex()->particles_out_size() : 0;
0687 }
0688
0689 /////////////////////////////////////////////////////////////////////////////
0690 template <typename P>
0691 const reco::GenParticle *MCTruthHelper<P>::daughter(const reco::GenParticle &p, unsigned int idau) const {
0692   return static_cast<const reco::GenParticle *>(p.daughter(idau));
0693 }
0694
0695 /////////////////////////////////////////////////////////////////////////////
0696 template <typename P>
0697 const HepMC::GenParticle *MCTruthHelper<P>::daughter(const HepMC::GenParticle &p, unsigned int idau) const {
0698   return *(p.end_vertex()->particles_out_const_begin() + idau);
0699 }
0700
0701 /////////////////////////////////////////////////////////////////////////////
0702 template <typename P>
0703 const HepMC3::GenParticle *MCTruthHelper<P>::daughter(const HepMC3::GenParticle &p, unsigned int idau) const {
0704   if (numberOfDaughters(p) > 0) {
0705     unsigned int idaui = 0;
0706     for (const HepMC3::ConstGenParticlePtr &daughter : (p.end_vertex())->particles_out()) {
0707       if (idaui == idau)
0708         return daughter.get();
0709       idaui++;
0710     }
0711   }
0712   return nullptr;
0713 }
0714
0715 /////////////////////////////////////////////////////////////////////////////
0716 template <typename P>
0717 void MCTruthHelper<P>::fillGenStatusFlags(const P &p, reco::GenStatusFlags &statusFlags) const {
0718   statusFlags.setIsPrompt(isPrompt(p));
0719   statusFlags.setIsDecayedLeptonHadron(isDecayedLeptonHadron(p));
0720   statusFlags.setIsTauDecayProduct(isTauDecayProduct(p));
0721   statusFlags.setIsPromptTauDecayProduct(isPromptTauDecayProduct(p));
0722   statusFlags.setIsDirectTauDecayProduct(isDirectTauDecayProduct(p));
0723   statusFlags.setIsDirectPromptTauDecayProduct(isDirectPromptTauDecayProduct(p));
0724   statusFlags.setIsDirectHadronDecayProduct(isDirectHadronDecayProduct(p));
0725   statusFlags.setIsHardProcess(isHardProcess(p));
0726   statusFlags.setFromHardProcess(fromHardProcess(p));
0727   statusFlags.setIsHardProcessTauDecayProduct(isHardProcessTauDecayProduct(p));
0728   statusFlags.setIsDirectHardProcessTauDecayProduct(isDirectHardProcessTauDecayProduct(p));
0729   statusFlags.setFromHardProcessBeforeFSR(fromHardProcessBeforeFSR(p));
0730   statusFlags.setIsFirstCopy(isFirstCopy(p));
0731   statusFlags.setIsLastCopy(isLastCopy(p));
0732   statusFlags.setIsLastCopyBeforeFSR(isLastCopyBeforeFSR(p));
0733 }
0734
0735 #endif