CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
MCTruthHelper.h
Go to the documentation of this file.
1 #ifndef PhysicsTools_HepMCCandAlgos_GenParticlesHelper_h
2 #define PhysicsTools_HepMCCandAlgos_GenParticlesHelper_h
3 
5 #include "HepPDT/ParticleID.hh"
8 
9 #include <iostream>
10 #include <unordered_set>
11 
12 template <typename P>
14 public:
16  //these are robust, generator-independent functions for categorizing
17  //mainly final state particles, but also intermediate hadrons
18  //or radiating leptons
19 
20  //is particle prompt (not from hadron, muon, or tau decay)
21  bool isPrompt(const P &p) const;
22 
23  //is particle prompt and final state
24  bool isPromptFinalState(const P &p) const;
25 
26  //is particle a decayed hadron, muon, or tau (does not include resonance decays like W,Z,Higgs,top,etc)
27  //This flag is equivalent to status 2 in the current HepMC standard
28  //but older generators (pythia6, herwig6) predate this and use status 2 also for other intermediate
29  //particles/states
30  bool isDecayedLeptonHadron(const P &p) const;
31 
32  //is particle prompt and decayed
33  bool isPromptDecayed(const P &p) const;
34 
35  //this particle is a direct or indirect tau decay product
36  bool isTauDecayProduct(const P &p) const;
37 
38  //this particle is a direct or indirect decay product of a prompt tau
39  bool isPromptTauDecayProduct(const P &p) const;
40 
41  //this particle is a direct tau decay product
42  bool isDirectTauDecayProduct(const P &p) const;
43 
44  //this particle is a direct decay product from a prompt tau
45  bool isDirectPromptTauDecayProduct(const P &p) const;
46 
47  //this particle is a direct or indirect muon decay product
48  bool isMuonDecayProduct(const P &p) const;
49 
50  //this particle is a direct or indirect decay product of a prompt muon
51  bool isPromptMuonDecayProduct(const P &p) const;
52 
53  //this particle is a direct decay product from a hadron
54  bool isDirectHadronDecayProduct(const P &p) const;
55 
56  //is particle a hadron
57  bool isHadron(const P &p) const;
58 
60  //these are generator history-dependent functions for tagging particles
61  //associated with the hard process
62  //Currently implemented for Pythia 6 and Pythia 8 status codes and history
63 
64  //this particle is part of the hard process
65  bool isHardProcess(const P &p) const;
66 
67  //this particle is the direct descendant of a hard process particle of the same pdg id
68  bool fromHardProcess(const P &p) const;
69 
70  //this particle is the final state direct descendant of a hard process particle
71  bool fromHardProcessFinalState(const P &p) const;
72 
73  //this particle is the decayed direct descendant of a hard process particle
74  //such as a tau from the hard process
75  bool fromHardProcessDecayed(const P &p) const;
76 
77  //this particle is a direct or indirect decay product of a tau
78  //from the hard process
79  bool isHardProcessTauDecayProduct(const P &p) const;
80 
81  //this particle is a direct decay product of a tau
82  //from the hard process
83  bool isDirectHardProcessTauDecayProduct(const P &p) const;
84 
85  //this particle is the direct descendant of a hard process particle of the same pdg id
86  //For outgoing particles the kinematics are those before QCD or QED FSR
87  //This corresponds roughly to status code 3 in pythia 6
88  bool fromHardProcessBeforeFSR(const P &p) const;
89 
90  //this particle is the first copy of the particle in the chain with the same pdg id
91  bool isFirstCopy(const P &p) const;
92 
93  //this particle is the last copy of the particle in the chain with the same pdg id
94  //(and therefore is more likely, but not guaranteed, to carry the final physical momentum)
95  bool isLastCopy(const P &p) const;
96 
97  //this particle is the last copy of the particle in the chain with the same pdg id
98  //before QED or QCD FSR
99  //(and therefore is more likely, but not guaranteed, to carry the momentum after ISR)
100  //This flag only really makes sense for outgoing particles
101  bool isLastCopyBeforeFSR(const P &p) const;
102 
104  //These are utility functions used by the above
105 
106  //first mother in chain with a different pdg than the particle
107  const P *uniqueMother(const P &p) const;
108 
109  //return first copy of particle in chain (may be the particle itself)
110  const P *firstCopy(const P &p) const;
111 
112  //return last copy of particle in chain (may be the particle itself)
113  const P *lastCopy(const P &p) const;
114 
115  //return last copy of particle in chain before QED or QCD FSR (may be the particle itself)
116  const P *lastCopyBeforeFSR(const P &p) const;
117 
118  //return last copy of particle in chain before QED or QCD FSR, starting from itself (may be the particle itself)
119  const P *lastDaughterCopyBeforeFSR(const P &p) const;
120 
121  //return mother copy which is a hard process particle
122  const P *hardProcessMotherCopy(const P &p) const;
123 
124  //return previous copy of particle in chain (0 in case this is already the first copy)
125  const P *previousCopy(const P &p) const;
126 
127  //return next copy of particle in chain (0 in case this is already the last copy)
128  const P *nextCopy(const P &p) const;
129 
130  //return decayed mother (walk up the chain until found)
131  const P *findDecayedMother(const P &p) const;
132 
133  //return decayed mother matching requested abs(pdgid) (walk up the chain until found)
134  const P *findDecayedMother(const P &p, int abspdgid) const;
135 
137  //These are very basic utility functions to implement a common interface for reco::GenParticle
138  //and HepMC::GenParticle
139 
140  //pdgid
141  int pdgId(const reco::GenParticle &p) const;
142 
143  //pdgid
144  int pdgId(const HepMC::GenParticle &p) const;
145 
146  //abs(pdgid)
147  int absPdgId(const reco::GenParticle &p) const;
148 
149  //abs(pdgid)
150  int absPdgId(const HepMC::GenParticle &p) const;
151 
152  //number of mothers
153  unsigned int numberOfMothers(const reco::GenParticle &p) const;
154 
155  //number of mothers
156  unsigned int numberOfMothers(const HepMC::GenParticle &p) const;
157 
158  //mother
159  const reco::GenParticle *mother(const reco::GenParticle &p, unsigned int imoth = 0) const;
160 
161  //mother
162  const HepMC::GenParticle *mother(const HepMC::GenParticle &p, unsigned int imoth = 0) const;
163 
164  //number of daughters
165  unsigned int numberOfDaughters(const reco::GenParticle &p) const;
166 
167  //number of daughters
168  unsigned int numberOfDaughters(const HepMC::GenParticle &p) const;
169 
170  //ith daughter
171  const reco::GenParticle *daughter(const reco::GenParticle &p, unsigned int idau) const;
172 
173  //ith daughter
174  const HepMC::GenParticle *daughter(const HepMC::GenParticle &p, unsigned int idau) const;
175 
177  //Helper function to fill status flags
179 };
180 
182 //implementations
183 
185 template <typename P>
186 bool MCTruthHelper<P>::isPrompt(const P &p) const {
187  //particle from hadron/muon/tau decay -> not prompt
188  //checking all the way up the chain treats properly the radiated photon
189  //case as well
190  return findDecayedMother(p) == nullptr;
191 }
192 
194 template <typename P>
196  return p.status() == 1 && isPrompt(p);
197 }
198 
199 template <typename P>
201  return p.status() == 2 && (isHadron(p) || absPdgId(p) == 13 || absPdgId(p) == 15) && isLastCopy(p);
202 }
203 
204 template <typename P>
206  return isDecayedLeptonHadron(p) && isPrompt(p);
207 }
208 
210 template <typename P>
212  return findDecayedMother(p, 15) != nullptr;
213 }
214 
216 template <typename P>
218  const P *tau = findDecayedMother(p, 15);
219  return tau && isPrompt(*tau);
220 }
221 
223 template <typename P>
225  const P *tau = findDecayedMother(p, 15);
226  const P *dm = findDecayedMother(p);
227  return tau && tau == dm;
228 }
229 
231 template <typename P>
233  const P *tau = findDecayedMother(p, 15);
234  const P *dm = findDecayedMother(p);
235  return tau && tau == dm && isPrompt(*tau);
236 }
237 
239 template <typename P>
241  return findDecayedMother(p, 13) != 0;
242 }
243 
245 template <typename P>
247  const P *mu = findDecayedMother(p, 13);
248  return mu && isPrompt(*mu);
249 }
250 
252 template <typename P>
254  const P *um = uniqueMother(p);
255  return um && isHadron(*um) && isDecayedLeptonHadron(*um);
256 }
257 
259 template <typename P>
260 bool MCTruthHelper<P>::isHadron(const P &p) const {
261  HepPDT::ParticleID heppdtid(pdgId(p));
262  return heppdtid.isHadron();
263 }
264 
266 template <typename P>
267 bool MCTruthHelper<P>::isHardProcess(const P &p) const {
268  //status 3 in pythia6 means hard process;
269  if (p.status() == 3)
270  return true;
271 
272  //hard process codes for pythia8 are 21-29 inclusive (currently 21,22,23,24 are used)
273  if (p.status() > 20 && p.status() < 30)
274  return true;
275 
276  //if this is a final state or decayed particle,
277  //check if direct mother is a resonance decay in pythia8 but exclude FSR branchings
278  //(In pythia8 if a resonance decay product did not undergo any further branchings
279  //it will be directly stored as status 1 or 2 without any status 23 copy)
280  if (p.status() == 1 || p.status() == 2) {
281  const P *um = mother(p);
282  if (um) {
283  const P *firstcopy = firstCopy(*um);
284  bool fromResonance = firstcopy && firstcopy->status() == 22;
285 
286  const P *umNext = nextCopy(*um);
287  bool fsrBranching = umNext && umNext->status() > 50 && umNext->status() < 60;
288 
289  if (fromResonance && !fsrBranching)
290  return true;
291  }
292  }
293 
294  return false;
295 }
296 
298 template <typename P>
300  return hardProcessMotherCopy(p) != nullptr;
301 }
302 
304 template <typename P>
306  return p.status() == 1 && fromHardProcess(p);
307 }
308 
310 template <typename P>
312  return isDecayedLeptonHadron(p) && fromHardProcess(p);
313 }
314 
316 template <typename P>
318  const P *tau = findDecayedMother(p, 15);
319  return tau && fromHardProcessDecayed(*tau);
320 }
321 
323 template <typename P>
325  const P *tau = findDecayedMother(p, 15);
326  const P *dm = findDecayedMother(p);
327  return tau && tau == dm && fromHardProcess(*tau);
328 }
329 
330 template <typename P>
332  //pythia 6 documentation line roughly corresponds to this condition
333  if (p.status() == 3)
334  return true;
335 
336  //check hard process mother properties
337  const P *hpc = hardProcessMotherCopy(p);
338  if (!hpc)
339  return false;
340 
341  //for incoming partons in pythia8, more useful information is not
342  //easily available, so take only the incoming parton itself
343  if (hpc->status() == 21 && (&p) == hpc)
344  return true;
345 
346  //for intermediate particles in pythia 8, just take the last copy
347  if (hpc->status() == 22 && isLastCopy(p))
348  return true;
349 
350  //for outgoing particles in pythia 8, explicitly find the last copy
351  //before FSR starting from the hardProcess particle, and take only
352  //this one
353  if ((hpc->status() == 23 || hpc->status() == 1) && (&p) == lastDaughterCopyBeforeFSR(*hpc))
354  return true;
355 
356  //didn't satisfy any of the conditions
357  return false;
358 }
359 
361 template <typename P>
362 bool MCTruthHelper<P>::isFirstCopy(const P &p) const {
363  return &p == firstCopy(p);
364 }
365 
367 template <typename P>
368 bool MCTruthHelper<P>::isLastCopy(const P &p) const {
369  return &p == lastCopy(p);
370 }
371 
373 template <typename P>
375  return &p == lastCopyBeforeFSR(p);
376 }
377 
379 template <typename P>
380 const P *MCTruthHelper<P>::uniqueMother(const P &p) const {
381  const P *mo = &p;
382  std::unordered_set<const P *> dupCheck;
383  while (mo && pdgId(*mo) == pdgId(p)) {
384  dupCheck.insert(mo);
385  mo = mother(*mo);
386  if (dupCheck.count(mo))
387  return nullptr;
388  }
389  return mo;
390 }
391 
393 template <typename P>
394 const P *MCTruthHelper<P>::firstCopy(const P &p) const {
395  const P *pcopy = &p;
396  std::unordered_set<const P *> dupCheck;
397  while (previousCopy(*pcopy)) {
398  dupCheck.insert(pcopy);
399  pcopy = previousCopy(*pcopy);
400  if (dupCheck.count(pcopy))
401  return nullptr;
402  }
403  return pcopy;
404 }
405 
407 template <typename P>
408 const P *MCTruthHelper<P>::lastCopy(const P &p) const {
409  const P *pcopy = &p;
410  std::unordered_set<const P *> dupCheck;
411  while (nextCopy(*pcopy)) {
412  dupCheck.insert(pcopy);
413  pcopy = nextCopy(*pcopy);
414  if (dupCheck.count(pcopy))
415  return nullptr;
416  }
417  return pcopy;
418 }
419 
421 template <typename P>
422 const P *MCTruthHelper<P>::lastCopyBeforeFSR(const P &p) const {
423  //start with first copy and then walk down until there is FSR
424  const P *pcopy = firstCopy(p);
425  if (!pcopy)
426  return nullptr;
427  bool hasDaughterCopy = true;
428  std::unordered_set<const P *> dupCheck;
429  while (hasDaughterCopy) {
430  dupCheck.insert(pcopy);
431  hasDaughterCopy = false;
432  const unsigned int ndau = numberOfDaughters(*pcopy);
433  //look for FSR
434  for (unsigned int idau = 0; idau < ndau; ++idau) {
435  const P *dau = daughter(*pcopy, idau);
436  if (pdgId(*dau) == 21 || pdgId(*dau) == 22) {
437  //has fsr (or else decayed and is the last copy by construction)
438  return pcopy;
439  }
440  }
441  //look for daughter copy
442  for (unsigned int idau = 0; idau < ndau; ++idau) {
443  const P *dau = daughter(*pcopy, idau);
444  if (pdgId(*dau) == pdgId(p)) {
445  pcopy = dau;
446  hasDaughterCopy = true;
447  break;
448  }
449  }
450  if (dupCheck.count(pcopy))
451  return nullptr;
452  }
453  return pcopy;
454 }
455 
457 template <typename P>
459  //start with this particle and then walk down until there is FSR
460  const P *pcopy = &p;
461  bool hasDaughterCopy = true;
462  std::unordered_set<const P *> dupCheck;
463  while (hasDaughterCopy) {
464  dupCheck.insert(pcopy);
465  hasDaughterCopy = false;
466  const unsigned int ndau = numberOfDaughters(*pcopy);
467  //look for FSR
468  for (unsigned int idau = 0; idau < ndau; ++idau) {
469  const P *dau = daughter(*pcopy, idau);
470  if (pdgId(*dau) == 21 || pdgId(*dau) == 22) {
471  //has fsr (or else decayed and is the last copy by construction)
472  return pcopy;
473  }
474  }
475  //look for daughter copy
476  for (unsigned int idau = 0; idau < ndau; ++idau) {
477  const P *dau = daughter(*pcopy, idau);
478  if (pdgId(*dau) == pdgId(p)) {
479  pcopy = dau;
480  hasDaughterCopy = true;
481  break;
482  }
483  }
484  if (dupCheck.count(pcopy))
485  return nullptr;
486  }
487  return pcopy;
488 }
489 
491 template <typename P>
493  //is particle itself is hard process particle
494  if (isHardProcess(p))
495  return &p;
496 
497  //check if any other copies are hard process particles
498  const P *pcopy = &p;
499  std::unordered_set<const P *> dupCheck;
500  while (previousCopy(*pcopy)) {
501  dupCheck.insert(pcopy);
502  pcopy = previousCopy(*pcopy);
503  if (isHardProcess(*pcopy))
504  return pcopy;
505  if (dupCheck.count(pcopy))
506  break;
507  }
508  return nullptr;
509 }
510 
512 template <typename P>
513 const P *MCTruthHelper<P>::previousCopy(const P &p) const {
514  const unsigned int nmoth = numberOfMothers(p);
515  for (unsigned int imoth = 0; imoth < nmoth; ++imoth) {
516  const P *moth = mother(p, imoth);
517  if (pdgId(*moth) == pdgId(p)) {
518  return moth;
519  }
520  }
521 
522  return nullptr;
523 }
524 
526 template <typename P>
527 const P *MCTruthHelper<P>::nextCopy(const P &p) const {
528  const unsigned int ndau = numberOfDaughters(p);
529  for (unsigned int idau = 0; idau < ndau; ++idau) {
530  const P *dau = daughter(p, idau);
531  if (pdgId(*dau) == pdgId(p)) {
532  return dau;
533  }
534  }
535 
536  return nullptr;
537 }
538 
540 template <typename P>
541 const P *MCTruthHelper<P>::findDecayedMother(const P &p) const {
542  const P *mo = mother(p);
543  std::unordered_set<const P *> dupCheck;
544  while (mo && !isDecayedLeptonHadron(*mo)) {
545  dupCheck.insert(mo);
546  mo = mother(*mo);
547  if (dupCheck.count(mo))
548  return nullptr;
549  }
550  return mo;
551 }
552 
554 template <typename P>
555 const P *MCTruthHelper<P>::findDecayedMother(const P &p, int abspdgid) const {
556  const P *mo = mother(p);
557  std::unordered_set<const P *> dupCheck;
558  while (mo && (absPdgId(*mo) != abspdgid || !isDecayedLeptonHadron(*mo))) {
559  dupCheck.insert(mo);
560  mo = mother(*mo);
561  if (dupCheck.count(mo))
562  return nullptr;
563  }
564  return mo;
565 }
566 
568 template <typename P>
570  return p.pdgId();
571 }
572 
574 template <typename P>
576  return p.pdg_id();
577 }
578 
580 template <typename P>
582  return std::abs(p.pdgId());
583 }
584 
586 template <typename P>
588  return std::abs(p.pdg_id());
589 }
590 
592 template <typename P>
594  return p.numberOfMothers();
595 }
596 
598 template <typename P>
600  return p.production_vertex() ? p.production_vertex()->particles_in_size() : 0;
601 }
602 
604 template <typename P>
605 const reco::GenParticle *MCTruthHelper<P>::mother(const reco::GenParticle &p, unsigned int imoth) const {
606  return static_cast<const reco::GenParticle *>(p.mother(imoth));
607 }
608 
610 template <typename P>
611 const HepMC::GenParticle *MCTruthHelper<P>::mother(const HepMC::GenParticle &p, unsigned int imoth) const {
612  return p.production_vertex() && p.production_vertex()->particles_in_size()
613  ? *(p.production_vertex()->particles_in_const_begin() + imoth)
614  : nullptr;
615 }
616 
618 template <typename P>
620  return p.numberOfDaughters();
621 }
622 
624 template <typename P>
626  return p.end_vertex() ? p.end_vertex()->particles_out_size() : 0;
627 }
628 
630 template <typename P>
631 const reco::GenParticle *MCTruthHelper<P>::daughter(const reco::GenParticle &p, unsigned int idau) const {
632  return static_cast<const reco::GenParticle *>(p.daughter(idau));
633 }
634 
636 template <typename P>
637 const HepMC::GenParticle *MCTruthHelper<P>::daughter(const HepMC::GenParticle &p, unsigned int idau) const {
638  return *(p.end_vertex()->particles_out_const_begin() + idau);
639 }
640 
642 template <typename P>
644  statusFlags.setIsPrompt(isPrompt(p));
645  statusFlags.setIsDecayedLeptonHadron(isDecayedLeptonHadron(p));
646  statusFlags.setIsTauDecayProduct(isTauDecayProduct(p));
647  statusFlags.setIsPromptTauDecayProduct(isPromptTauDecayProduct(p));
648  statusFlags.setIsDirectTauDecayProduct(isDirectTauDecayProduct(p));
649  statusFlags.setIsDirectPromptTauDecayProduct(isDirectPromptTauDecayProduct(p));
650  statusFlags.setIsDirectHadronDecayProduct(isDirectHadronDecayProduct(p));
651  statusFlags.setIsHardProcess(isHardProcess(p));
652  statusFlags.setFromHardProcess(fromHardProcess(p));
653  statusFlags.setIsHardProcessTauDecayProduct(isHardProcessTauDecayProduct(p));
654  statusFlags.setIsDirectHardProcessTauDecayProduct(isDirectHardProcessTauDecayProduct(p));
655  statusFlags.setFromHardProcessBeforeFSR(fromHardProcessBeforeFSR(p));
656  statusFlags.setIsFirstCopy(isFirstCopy(p));
657  statusFlags.setIsLastCopy(isLastCopy(p));
658  statusFlags.setIsLastCopyBeforeFSR(isLastCopyBeforeFSR(p));
659 }
660 
661 #endif
const reco::GenParticle * mother(const reco::GenParticle &p, unsigned int imoth=0) const
const P * uniqueMother(const P &p) const
int absPdgId(const reco::GenParticle &p) const
bool isHadron(const P &p) const
const P * findDecayedMother(const P &p) const
bool isPromptMuonDecayProduct(const P &p) const
const reco::GenParticle * daughter(const reco::GenParticle &p, unsigned int idau) const
const P * lastCopy(const P &p) const
const P * lastCopyBeforeFSR(const P &p) const
bool isDirectHadronDecayProduct(const P &p) const
const P * hardProcessMotherCopy(const P &p) const
bool isLastCopy(const P &p) const
bool isPrompt(const P &p) const
bool isHardProcess(const P &p) const
bool isPromptFinalState(const P &p) const
bool isPromptDecayed(const P &p) const
const P * nextCopy(const P &p) const
bool isLastCopyBeforeFSR(const P &p) const
bool fromHardProcessDecayed(const P &p) const
bool isHardProcessTauDecayProduct(const P &p) const
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
const P * lastDaughterCopyBeforeFSR(const P &p) const
bool isFirstCopy(const P &p) const
int pdgId(const reco::GenParticle &p) const
bool isPromptTauDecayProduct(const P &p) const
bool isDecayedLeptonHadron(const P &p) const
std::pair< OmniClusterRef, TrackingParticleRef > P
const P * previousCopy(const P &p) const
bool fromHardProcessFinalState(const P &p) const
unsigned int numberOfMothers(const reco::GenParticle &p) const
bool isMuonDecayProduct(const P &p) const
bool fromHardProcessBeforeFSR(const P &p) const
void fillGenStatusFlags(const P &p, reco::GenStatusFlags &statusFlags) const
bool isDirectPromptTauDecayProduct(const P &p) const
bool isDirectTauDecayProduct(const P &p) const
unsigned int numberOfDaughters(const reco::GenParticle &p) const
bool isDirectHardProcessTauDecayProduct(const P &p) const
bool isTauDecayProduct(const P &p) const
bool fromHardProcess(const P &p) const
const P * firstCopy(const P &p) const