CMS 3D CMS Logo

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