File indexing completed on 2024-04-06 12:23:24
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 <iostream>
0010 #include <unordered_set>
0011
0012 template <typename P>
0013 class MCTruthHelper {
0014 public:
0015
0016
0017
0018
0019
0020
0021 bool isPrompt(const P &p) const;
0022
0023
0024 bool isPromptFinalState(const P &p) const;
0025
0026
0027
0028
0029
0030 bool isDecayedLeptonHadron(const P &p) const;
0031
0032
0033 bool isPromptDecayed(const P &p) const;
0034
0035
0036 bool isTauDecayProduct(const P &p) const;
0037
0038
0039 bool isPromptTauDecayProduct(const P &p) const;
0040
0041
0042 bool isDirectTauDecayProduct(const P &p) const;
0043
0044
0045 bool isDirectPromptTauDecayProduct(const P &p) const;
0046
0047
0048 bool isMuonDecayProduct(const P &p) const;
0049
0050
0051 bool isPromptMuonDecayProduct(const P &p) const;
0052
0053
0054 bool isDirectHadronDecayProduct(const P &p) const;
0055
0056
0057 bool isHadron(const P &p) const;
0058
0059
0060
0061
0062
0063
0064
0065 bool isHardProcess(const P &p) const;
0066
0067
0068 bool fromHardProcess(const P &p) const;
0069
0070
0071 bool fromHardProcessFinalState(const P &p) const;
0072
0073
0074
0075 bool fromHardProcessDecayed(const P &p) const;
0076
0077
0078
0079 bool isHardProcessTauDecayProduct(const P &p) const;
0080
0081
0082
0083 bool isDirectHardProcessTauDecayProduct(const P &p) const;
0084
0085
0086
0087
0088 bool fromHardProcessBeforeFSR(const P &p) const;
0089
0090
0091 bool isFirstCopy(const P &p) const;
0092
0093
0094
0095 bool isLastCopy(const P &p) const;
0096
0097
0098
0099
0100
0101 bool isLastCopyBeforeFSR(const P &p) const;
0102
0103
0104
0105
0106
0107 const P *uniqueMother(const P &p) const;
0108
0109
0110 const P *firstCopy(const P &p) const;
0111
0112
0113 const P *lastCopy(const P &p) const;
0114
0115
0116 const P *lastCopyBeforeFSR(const P &p) const;
0117
0118
0119 const P *lastDaughterCopyBeforeFSR(const P &p) const;
0120
0121
0122 const P *hardProcessMotherCopy(const P &p) const;
0123
0124
0125 const P *previousCopy(const P &p) const;
0126
0127
0128 const P *nextCopy(const P &p) const;
0129
0130
0131 const P *findDecayedMother(const P &p) const;
0132
0133
0134 const P *findDecayedMother(const P &p, int abspdgid) const;
0135
0136
0137
0138
0139
0140
0141 int pdgId(const reco::GenParticle &p) const;
0142
0143
0144 int pdgId(const HepMC::GenParticle &p) const;
0145
0146
0147 int absPdgId(const reco::GenParticle &p) const;
0148
0149
0150 int absPdgId(const HepMC::GenParticle &p) const;
0151
0152
0153 unsigned int numberOfMothers(const reco::GenParticle &p) const;
0154
0155
0156 unsigned int numberOfMothers(const HepMC::GenParticle &p) const;
0157
0158
0159 const reco::GenParticle *mother(const reco::GenParticle &p, unsigned int imoth = 0) const;
0160
0161
0162 const HepMC::GenParticle *mother(const HepMC::GenParticle &p, unsigned int imoth = 0) const;
0163
0164
0165 unsigned int numberOfDaughters(const reco::GenParticle &p) const;
0166
0167
0168 unsigned int numberOfDaughters(const HepMC::GenParticle &p) const;
0169
0170
0171 const reco::GenParticle *daughter(const reco::GenParticle &p, unsigned int idau) const;
0172
0173
0174 const HepMC::GenParticle *daughter(const HepMC::GenParticle &p, unsigned int idau) const;
0175
0176
0177
0178 void fillGenStatusFlags(const P &p, reco::GenStatusFlags &statusFlags) const;
0179 };
0180
0181
0182
0183
0184
0185 template <typename P>
0186 bool MCTruthHelper<P>::isPrompt(const P &p) const {
0187
0188
0189
0190 return findDecayedMother(p) == nullptr;
0191 }
0192
0193
0194 template <typename P>
0195 bool MCTruthHelper<P>::isPromptFinalState(const P &p) const {
0196 return p.status() == 1 && isPrompt(p);
0197 }
0198
0199 template <typename P>
0200 bool MCTruthHelper<P>::isDecayedLeptonHadron(const P &p) const {
0201 return p.status() == 2 && (isHadron(p) || absPdgId(p) == 13 || absPdgId(p) == 15) && isLastCopy(p);
0202 }
0203
0204 template <typename P>
0205 bool MCTruthHelper<P>::isPromptDecayed(const P &p) const {
0206 return isDecayedLeptonHadron(p) && isPrompt(p);
0207 }
0208
0209
0210 template <typename P>
0211 bool MCTruthHelper<P>::isTauDecayProduct(const P &p) const {
0212 return findDecayedMother(p, 15) != nullptr;
0213 }
0214
0215
0216 template <typename P>
0217 bool MCTruthHelper<P>::isPromptTauDecayProduct(const P &p) const {
0218 const P *tau = findDecayedMother(p, 15);
0219 return tau && isPrompt(*tau);
0220 }
0221
0222
0223 template <typename P>
0224 bool MCTruthHelper<P>::isDirectTauDecayProduct(const P &p) const {
0225 const P *tau = findDecayedMother(p, 15);
0226 const P *dm = findDecayedMother(p);
0227 return tau && tau == dm;
0228 }
0229
0230
0231 template <typename P>
0232 bool MCTruthHelper<P>::isDirectPromptTauDecayProduct(const P &p) const {
0233 const P *tau = findDecayedMother(p, 15);
0234 const P *dm = findDecayedMother(p);
0235 return tau && tau == dm && isPrompt(*tau);
0236 }
0237
0238
0239 template <typename P>
0240 bool MCTruthHelper<P>::isMuonDecayProduct(const P &p) const {
0241 return findDecayedMother(p, 13) != 0;
0242 }
0243
0244
0245 template <typename P>
0246 bool MCTruthHelper<P>::isPromptMuonDecayProduct(const P &p) const {
0247 const P *mu = findDecayedMother(p, 13);
0248 return mu && isPrompt(*mu);
0249 }
0250
0251
0252 template <typename P>
0253 bool MCTruthHelper<P>::isDirectHadronDecayProduct(const P &p) const {
0254 const P *um = uniqueMother(p);
0255 return um && isHadron(*um) && isDecayedLeptonHadron(*um);
0256 }
0257
0258
0259 template <typename P>
0260 bool MCTruthHelper<P>::isHadron(const P &p) const {
0261 HepPDT::ParticleID heppdtid(pdgId(p));
0262 return heppdtid.isHadron();
0263 }
0264
0265
0266 template <typename P>
0267 bool MCTruthHelper<P>::isHardProcess(const P &p) const {
0268
0269 if (p.status() == 3)
0270 return true;
0271
0272
0273 if (p.status() > 20 && p.status() < 30)
0274 return true;
0275
0276
0277
0278
0279
0280 if (p.status() == 1 || p.status() == 2) {
0281 const P *um = mother(p);
0282 if (um) {
0283 const P *firstcopy = firstCopy(*um);
0284 bool fromResonance = firstcopy && firstcopy->status() == 22;
0285
0286 const P *umNext = nextCopy(*um);
0287 bool fsrBranching = umNext && umNext->status() > 50 && umNext->status() < 60;
0288
0289 if (fromResonance && !fsrBranching)
0290 return true;
0291 }
0292 }
0293
0294 return false;
0295 }
0296
0297
0298 template <typename P>
0299 bool MCTruthHelper<P>::fromHardProcess(const P &p) const {
0300 return hardProcessMotherCopy(p) != nullptr;
0301 }
0302
0303
0304 template <typename P>
0305 bool MCTruthHelper<P>::fromHardProcessFinalState(const P &p) const {
0306 return p.status() == 1 && fromHardProcess(p);
0307 }
0308
0309
0310 template <typename P>
0311 bool MCTruthHelper<P>::fromHardProcessDecayed(const P &p) const {
0312 return isDecayedLeptonHadron(p) && fromHardProcess(p);
0313 }
0314
0315
0316 template <typename P>
0317 bool MCTruthHelper<P>::isHardProcessTauDecayProduct(const P &p) const {
0318 const P *tau = findDecayedMother(p, 15);
0319 return tau && fromHardProcessDecayed(*tau);
0320 }
0321
0322
0323 template <typename P>
0324 bool MCTruthHelper<P>::isDirectHardProcessTauDecayProduct(const P &p) const {
0325 const P *tau = findDecayedMother(p, 15);
0326 const P *dm = findDecayedMother(p);
0327 return tau && tau == dm && fromHardProcess(*tau);
0328 }
0329
0330 template <typename P>
0331 bool MCTruthHelper<P>::fromHardProcessBeforeFSR(const P &p) const {
0332
0333 if (p.status() == 3)
0334 return true;
0335
0336
0337 const P *hpc = hardProcessMotherCopy(p);
0338 if (!hpc)
0339 return false;
0340
0341
0342
0343 if (hpc->status() == 21 && (&p) == hpc)
0344 return true;
0345
0346
0347 if (hpc->status() == 22 && isLastCopy(p))
0348 return true;
0349
0350
0351
0352
0353 if ((hpc->status() == 23 || hpc->status() == 1) && (&p) == lastDaughterCopyBeforeFSR(*hpc))
0354 return true;
0355
0356
0357 return false;
0358 }
0359
0360
0361 template <typename P>
0362 bool MCTruthHelper<P>::isFirstCopy(const P &p) const {
0363 return &p == firstCopy(p);
0364 }
0365
0366
0367 template <typename P>
0368 bool MCTruthHelper<P>::isLastCopy(const P &p) const {
0369 return &p == lastCopy(p);
0370 }
0371
0372
0373 template <typename P>
0374 bool MCTruthHelper<P>::isLastCopyBeforeFSR(const P &p) const {
0375 return &p == lastCopyBeforeFSR(p);
0376 }
0377
0378
0379 template <typename P>
0380 const P *MCTruthHelper<P>::uniqueMother(const P &p) const {
0381 const P *mo = &p;
0382 std::unordered_set<const P *> dupCheck;
0383 while (mo && pdgId(*mo) == pdgId(p)) {
0384 dupCheck.insert(mo);
0385 mo = mother(*mo);
0386 if (dupCheck.count(mo))
0387 return nullptr;
0388 }
0389 return mo;
0390 }
0391
0392
0393 template <typename P>
0394 const P *MCTruthHelper<P>::firstCopy(const P &p) const {
0395 const P *pcopy = &p;
0396 std::unordered_set<const P *> dupCheck;
0397 while (previousCopy(*pcopy)) {
0398 dupCheck.insert(pcopy);
0399 pcopy = previousCopy(*pcopy);
0400 if (dupCheck.count(pcopy))
0401 return nullptr;
0402 }
0403 return pcopy;
0404 }
0405
0406
0407 template <typename P>
0408 const P *MCTruthHelper<P>::lastCopy(const P &p) const {
0409 const P *pcopy = &p;
0410 std::unordered_set<const P *> dupCheck;
0411 while (nextCopy(*pcopy)) {
0412 dupCheck.insert(pcopy);
0413 pcopy = nextCopy(*pcopy);
0414 if (dupCheck.count(pcopy))
0415 return nullptr;
0416 }
0417 return pcopy;
0418 }
0419
0420
0421 template <typename P>
0422 const P *MCTruthHelper<P>::lastCopyBeforeFSR(const P &p) const {
0423
0424 const P *pcopy = firstCopy(p);
0425 if (!pcopy)
0426 return nullptr;
0427 bool hasDaughterCopy = true;
0428 std::unordered_set<const P *> dupCheck;
0429 while (hasDaughterCopy) {
0430 dupCheck.insert(pcopy);
0431 hasDaughterCopy = false;
0432 const unsigned int ndau = numberOfDaughters(*pcopy);
0433
0434 for (unsigned int idau = 0; idau < ndau; ++idau) {
0435 const P *dau = daughter(*pcopy, idau);
0436 if (pdgId(*dau) == 21 || pdgId(*dau) == 22) {
0437
0438 return pcopy;
0439 }
0440 }
0441
0442 for (unsigned int idau = 0; idau < ndau; ++idau) {
0443 const P *dau = daughter(*pcopy, idau);
0444 if (pdgId(*dau) == pdgId(p)) {
0445 pcopy = dau;
0446 hasDaughterCopy = true;
0447 break;
0448 }
0449 }
0450 if (dupCheck.count(pcopy))
0451 return nullptr;
0452 }
0453 return pcopy;
0454 }
0455
0456
0457 template <typename P>
0458 const P *MCTruthHelper<P>::lastDaughterCopyBeforeFSR(const P &p) const {
0459
0460 const P *pcopy = &p;
0461 bool hasDaughterCopy = true;
0462 std::unordered_set<const P *> dupCheck;
0463 while (hasDaughterCopy) {
0464 dupCheck.insert(pcopy);
0465 hasDaughterCopy = false;
0466 const unsigned int ndau = numberOfDaughters(*pcopy);
0467
0468 for (unsigned int idau = 0; idau < ndau; ++idau) {
0469 const P *dau = daughter(*pcopy, idau);
0470 if (pdgId(*dau) == 21 || pdgId(*dau) == 22) {
0471
0472 return pcopy;
0473 }
0474 }
0475
0476 for (unsigned int idau = 0; idau < ndau; ++idau) {
0477 const P *dau = daughter(*pcopy, idau);
0478 if (pdgId(*dau) == pdgId(p)) {
0479 pcopy = dau;
0480 hasDaughterCopy = true;
0481 break;
0482 }
0483 }
0484 if (dupCheck.count(pcopy))
0485 return nullptr;
0486 }
0487 return pcopy;
0488 }
0489
0490
0491 template <typename P>
0492 const P *MCTruthHelper<P>::hardProcessMotherCopy(const P &p) const {
0493
0494 if (isHardProcess(p))
0495 return &p;
0496
0497
0498 const P *pcopy = &p;
0499 std::unordered_set<const P *> dupCheck;
0500 while (previousCopy(*pcopy)) {
0501 dupCheck.insert(pcopy);
0502 pcopy = previousCopy(*pcopy);
0503 if (isHardProcess(*pcopy))
0504 return pcopy;
0505 if (dupCheck.count(pcopy))
0506 break;
0507 }
0508 return nullptr;
0509 }
0510
0511
0512 template <typename P>
0513 const P *MCTruthHelper<P>::previousCopy(const P &p) const {
0514 const unsigned int nmoth = numberOfMothers(p);
0515 for (unsigned int imoth = 0; imoth < nmoth; ++imoth) {
0516 const P *moth = mother(p, imoth);
0517 if (pdgId(*moth) == pdgId(p)) {
0518 return moth;
0519 }
0520 }
0521
0522 return nullptr;
0523 }
0524
0525
0526 template <typename P>
0527 const P *MCTruthHelper<P>::nextCopy(const P &p) const {
0528 const unsigned int ndau = numberOfDaughters(p);
0529 for (unsigned int idau = 0; idau < ndau; ++idau) {
0530 const P *dau = daughter(p, idau);
0531 if (pdgId(*dau) == pdgId(p)) {
0532 return dau;
0533 }
0534 }
0535
0536 return nullptr;
0537 }
0538
0539
0540 template <typename P>
0541 const P *MCTruthHelper<P>::findDecayedMother(const P &p) const {
0542 const P *mo = mother(p);
0543 std::unordered_set<const P *> dupCheck;
0544 while (mo && !isDecayedLeptonHadron(*mo)) {
0545 dupCheck.insert(mo);
0546 mo = mother(*mo);
0547 if (dupCheck.count(mo))
0548 return nullptr;
0549 }
0550 return mo;
0551 }
0552
0553
0554 template <typename P>
0555 const P *MCTruthHelper<P>::findDecayedMother(const P &p, int abspdgid) const {
0556 const P *mo = mother(p);
0557 std::unordered_set<const P *> dupCheck;
0558 while (mo && (absPdgId(*mo) != abspdgid || !isDecayedLeptonHadron(*mo))) {
0559 dupCheck.insert(mo);
0560 mo = mother(*mo);
0561 if (dupCheck.count(mo))
0562 return nullptr;
0563 }
0564 return mo;
0565 }
0566
0567
0568 template <typename P>
0569 int MCTruthHelper<P>::pdgId(const reco::GenParticle &p) const {
0570 return p.pdgId();
0571 }
0572
0573
0574 template <typename P>
0575 int MCTruthHelper<P>::pdgId(const HepMC::GenParticle &p) const {
0576 return p.pdg_id();
0577 }
0578
0579
0580 template <typename P>
0581 int MCTruthHelper<P>::absPdgId(const reco::GenParticle &p) const {
0582 return std::abs(p.pdgId());
0583 }
0584
0585
0586 template <typename P>
0587 int MCTruthHelper<P>::absPdgId(const HepMC::GenParticle &p) const {
0588 return std::abs(p.pdg_id());
0589 }
0590
0591
0592 template <typename P>
0593 unsigned int MCTruthHelper<P>::numberOfMothers(const reco::GenParticle &p) const {
0594 return p.numberOfMothers();
0595 }
0596
0597
0598 template <typename P>
0599 unsigned int MCTruthHelper<P>::numberOfMothers(const HepMC::GenParticle &p) const {
0600 return p.production_vertex() ? p.production_vertex()->particles_in_size() : 0;
0601 }
0602
0603
0604 template <typename P>
0605 const reco::GenParticle *MCTruthHelper<P>::mother(const reco::GenParticle &p, unsigned int imoth) const {
0606 return static_cast<const reco::GenParticle *>(p.mother(imoth));
0607 }
0608
0609
0610 template <typename P>
0611 const HepMC::GenParticle *MCTruthHelper<P>::mother(const HepMC::GenParticle &p, unsigned int imoth) const {
0612 return p.production_vertex() && p.production_vertex()->particles_in_size()
0613 ? *(p.production_vertex()->particles_in_const_begin() + imoth)
0614 : nullptr;
0615 }
0616
0617
0618 template <typename P>
0619 unsigned int MCTruthHelper<P>::numberOfDaughters(const reco::GenParticle &p) const {
0620 return p.numberOfDaughters();
0621 }
0622
0623
0624 template <typename P>
0625 unsigned int MCTruthHelper<P>::numberOfDaughters(const HepMC::GenParticle &p) const {
0626 return p.end_vertex() ? p.end_vertex()->particles_out_size() : 0;
0627 }
0628
0629
0630 template <typename P>
0631 const reco::GenParticle *MCTruthHelper<P>::daughter(const reco::GenParticle &p, unsigned int idau) const {
0632 return static_cast<const reco::GenParticle *>(p.daughter(idau));
0633 }
0634
0635
0636 template <typename P>
0637 const HepMC::GenParticle *MCTruthHelper<P>::daughter(const HepMC::GenParticle &p, unsigned int idau) const {
0638 return *(p.end_vertex()->particles_out_const_begin() + idau);
0639 }
0640
0641
0642 template <typename P>
0643 void MCTruthHelper<P>::fillGenStatusFlags(const P &p, reco::GenStatusFlags &statusFlags) const {
0644 statusFlags.setIsPrompt(isPrompt(p));
0645 statusFlags.setIsDecayedLeptonHadron(isDecayedLeptonHadron(p));
0646 statusFlags.setIsTauDecayProduct(isTauDecayProduct(p));
0647 statusFlags.setIsPromptTauDecayProduct(isPromptTauDecayProduct(p));
0648 statusFlags.setIsDirectTauDecayProduct(isDirectTauDecayProduct(p));
0649 statusFlags.setIsDirectPromptTauDecayProduct(isDirectPromptTauDecayProduct(p));
0650 statusFlags.setIsDirectHadronDecayProduct(isDirectHadronDecayProduct(p));
0651 statusFlags.setIsHardProcess(isHardProcess(p));
0652 statusFlags.setFromHardProcess(fromHardProcess(p));
0653 statusFlags.setIsHardProcessTauDecayProduct(isHardProcessTauDecayProduct(p));
0654 statusFlags.setIsDirectHardProcessTauDecayProduct(isDirectHardProcessTauDecayProduct(p));
0655 statusFlags.setFromHardProcessBeforeFSR(fromHardProcessBeforeFSR(p));
0656 statusFlags.setIsFirstCopy(isFirstCopy(p));
0657 statusFlags.setIsLastCopy(isLastCopy(p));
0658 statusFlags.setIsLastCopyBeforeFSR(isLastCopyBeforeFSR(p));
0659 }
0660
0661 #endif