TopObjects/src/TtGenEvent.cc

0001 //
0002 //
0003
0004 #include "FWCore/Utilities/interface/EDMException.h"
0005 #include "CommonTools/CandUtils/interface/pdgIdUtils.h"
0006 #include "AnalysisDataFormats/TopObjects/interface/TtGenEvent.h"
0007 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0008
0009 /// default constructor from decaySubset and initSubset
0010 TtGenEvent::TtGenEvent(reco::GenParticleRefProd& decaySubset, reco::GenParticleRefProd& initSubset) {
0011   parts_ = decaySubset;
0012   initPartons_ = initSubset;
0013   if (top() && topBar())
0014     topPair_ = math::XYZTLorentzVector(top()->p4() + topBar()->p4());
0015 }
0016
0017 bool TtGenEvent::fromGluonFusion() const {
0018   const reco::GenParticleCollection& initPartsColl = *initPartons_;
0019   if (initPartsColl.size() == 2)
0020     if (initPartsColl[0].pdgId() == 21 && initPartsColl[1].pdgId() == 21)
0021       return true;
0022   return false;
0023 }
0024
0025 bool TtGenEvent::fromQuarkAnnihilation() const {
0026   const reco::GenParticleCollection& initPartsColl = *initPartons_;
0027   if (initPartsColl.size() == 2)
0028     if (std::abs(initPartsColl[0].pdgId()) < TopDecayID::tID && initPartsColl[0].pdgId() == -initPartsColl[1].pdgId())
0029       return true;
0030   return false;
0031 }
0032
0033 WDecay::LeptonType TtGenEvent::semiLeptonicChannel() const {
0034   WDecay::LeptonType type = WDecay::kNone;
0035   if (isSemiLeptonic() && singleLepton()) {
0036     if (std::abs(singleLepton()->pdgId()) == TopDecayID::elecID)
0037       type = WDecay::kElec;
0038     if (std::abs(singleLepton()->pdgId()) == TopDecayID::muonID)
0039       type = WDecay::kMuon;
0040     if (std::abs(singleLepton()->pdgId()) == TopDecayID::tauID)
0041       type = WDecay::kTau;
0042   }
0043   return type;
0044 }
0045
0046 std::pair<WDecay::LeptonType, WDecay::LeptonType> TtGenEvent::fullLeptonicChannel() const {
0047   WDecay::LeptonType typeA = WDecay::kNone, typeB = WDecay::kNone;
0048   if (isFullLeptonic()) {
0049     if (lepton()) {
0050       if (std::abs(lepton()->pdgId()) == TopDecayID::elecID)
0051         typeA = WDecay::kElec;
0052       if (std::abs(lepton()->pdgId()) == TopDecayID::muonID)
0053         typeA = WDecay::kMuon;
0054       if (std::abs(lepton()->pdgId()) == TopDecayID::tauID)
0055         typeA = WDecay::kTau;
0056     }
0057     if (leptonBar()) {
0058       if (std::abs(leptonBar()->pdgId()) == TopDecayID::elecID)
0059         typeB = WDecay::kElec;
0060       if (std::abs(leptonBar()->pdgId()) == TopDecayID::muonID)
0061         typeB = WDecay::kMuon;
0062       if (std::abs(leptonBar()->pdgId()) == TopDecayID::tauID)
0063         typeB = WDecay::kTau;
0064     }
0065   }
0066   return (std::pair<WDecay::LeptonType, WDecay::LeptonType>(typeA, typeB));
0067 }
0068
0069 const reco::GenParticle* TtGenEvent::lepton(bool excludeTauLeptons) const {
0070   const reco::GenParticle* cand = nullptr;
0071   const reco::GenParticleCollection& partsColl = *parts_;
0072   for (unsigned int i = 0; i < partsColl.size(); ++i) {
0073     if (reco::isLepton(partsColl[i]) && partsColl[i].mother() &&
0074         std::abs(partsColl[i].mother()->pdgId()) == TopDecayID::WID) {
0075       if (reco::flavour(partsColl[i]) > 0) {
0076         cand = &partsColl[i];
0077       }
0078     }
0079   }
0080   return cand;
0081 }
0082
0083 const reco::GenParticle* TtGenEvent::leptonBar(bool excludeTauLeptons) const {
0084   const reco::GenParticle* cand = nullptr;
0085   const reco::GenParticleCollection& partsColl = *parts_;
0086   for (unsigned int i = 0; i < partsColl.size(); ++i) {
0087     if (reco::isLepton(partsColl[i]) && partsColl[i].mother() &&
0088         std::abs(partsColl[i].mother()->pdgId()) == TopDecayID::WID) {
0089       if (reco::flavour(partsColl[i]) < 0) {
0090         cand = &partsColl[i];
0091       }
0092     }
0093   }
0094   return cand;
0095 }
0096
0097 const reco::GenParticle* TtGenEvent::singleLepton(bool excludeTauLeptons) const {
0098   const reco::GenParticle* cand = nullptr;
0099   if (isSemiLeptonic(excludeTauLeptons)) {
0100     const reco::GenParticleCollection& partsColl = *parts_;
0101     for (unsigned int i = 0; i < partsColl.size(); ++i) {
0102       if (reco::isLepton(partsColl[i]) && partsColl[i].mother() &&
0103           std::abs(partsColl[i].mother()->pdgId()) == TopDecayID::WID) {
0104         cand = &partsColl[i];
0105       }
0106     }
0107   }
0108   return cand;
0109 }
0110
0111 const reco::GenParticle* TtGenEvent::neutrino(bool excludeTauLeptons) const {
0112   const reco::GenParticle* cand = nullptr;
0113   const reco::GenParticleCollection& partsColl = *parts_;
0114   for (unsigned int i = 0; i < partsColl.size(); ++i) {
0115     if (reco::isNeutrino(partsColl[i]) && partsColl[i].mother() &&
0116         std::abs(partsColl[i].mother()->pdgId()) == TopDecayID::WID) {
0117       if (reco::flavour(partsColl[i]) > 0) {
0118         cand = &partsColl[i];
0119       }
0120     }
0121   }
0122   return cand;
0123 }
0124
0125 const reco::GenParticle* TtGenEvent::neutrinoBar(bool excludeTauLeptons) const {
0126   const reco::GenParticle* cand = nullptr;
0127   const reco::GenParticleCollection& partsColl = *parts_;
0128   for (unsigned int i = 0; i < partsColl.size(); ++i) {
0129     if (reco::isNeutrino(partsColl[i]) && partsColl[i].mother() &&
0130         std::abs(partsColl[i].mother()->pdgId()) == TopDecayID::WID) {
0131       if (reco::flavour(partsColl[i]) < 0) {
0132         cand = &partsColl[i];
0133       }
0134     }
0135   }
0136   return cand;
0137 }
0138
0139 const reco::GenParticle* TtGenEvent::singleNeutrino(bool excludeTauLeptons) const {
0140   const reco::GenParticle* cand = nullptr;
0141   if (isSemiLeptonic(excludeTauLeptons)) {
0142     const reco::GenParticleCollection& partsColl = *parts_;
0143     for (unsigned int i = 0; i < partsColl.size(); ++i) {
0144       if (reco::isNeutrino(partsColl[i]) && partsColl[i].mother() &&
0145           std::abs(partsColl[i].mother()->pdgId()) == TopDecayID::WID) {
0146         cand = &partsColl[i];
0147       }
0148     }
0149   }
0150   return cand;
0151 }
0152
0153 const reco::GenParticle* TtGenEvent::hadronicDecayQuark(bool invertFlavor) const {
0154   const reco::GenParticle* cand = nullptr;
0155   // catch W boson and check its daughters for a quark;
0156   // make sure the decay is semi-leptonic first; this
0157   // only makes sense if taus are not excluded from the
0158   // decision
0159   if (singleLepton(false)) {
0160     for (reco::GenParticleCollection::const_iterator w = parts_->begin(); w != parts_->end(); ++w) {
0161       if (std::abs(w->pdgId()) == TopDecayID::WID) {
0162         // make sure that the particle is a W daughter
0163         for (reco::GenParticle::const_iterator wd = w->begin(); wd != w->end(); ++wd) {
0164           // make sure that the parton is a quark
0165           if (std::abs(wd->pdgId()) < TopDecayID::tID) {
0166             if (invertFlavor ? reco::flavour(*wd) < 0 : reco::flavour(*wd) > 0) {
0167               cand = dynamic_cast<const reco::GenParticle*>(&(*wd));
0168               if (cand == nullptr) {
0169                 throw edm::Exception(edm::errors::InvalidReference, "Not a GenParticle");
0170               }
0171               break;
0172             }
0173           }
0174         }
0175       }
0176     }
0177   }
0178   return cand;
0179 }
0180
0181 const reco::GenParticle* TtGenEvent::hadronicDecayB(bool excludeTauLeptons) const {
0182   const reco::GenParticle* cand = nullptr;
0183   if (singleLepton(excludeTauLeptons)) {
0184     const reco::GenParticleCollection& partsColl = *parts_;
0185     const reco::GenParticle& singleLep = *singleLepton(excludeTauLeptons);
0186     for (unsigned int i = 0; i < partsColl.size(); ++i) {
0187       if (std::abs(partsColl[i].pdgId()) == TopDecayID::bID &&
0188           reco::flavour(singleLep) == reco::flavour(partsColl[i])) {
0189         cand = &partsColl[i];
0190       }
0191     }
0192   }
0193   return cand;
0194 }
0195
0196 const reco::GenParticle* TtGenEvent::hadronicDecayW(bool excludeTauLeptons) const {
0197   const reco::GenParticle* cand = nullptr;
0198   if (singleLepton(excludeTauLeptons)) {
0199     const reco::GenParticleCollection& partsColl = *parts_;
0200     const reco::GenParticle& singleLep = *singleLepton(excludeTauLeptons);
0201     for (unsigned int i = 0; i < partsColl.size(); ++i) {
0202       if (std::abs(partsColl[i].pdgId()) == TopDecayID::WID &&
0203           reco::flavour(singleLep) != -reco::flavour(partsColl[i])) {
0204         // PDG Id:13=mu- 24=W+ (+24)->(-13) (-24)->(+13) opposite sign
0205         cand = &partsColl[i];
0206       }
0207     }
0208   }
0209   return cand;
0210 }
0211
0212 const reco::GenParticle* TtGenEvent::hadronicDecayTop(bool excludeTauLeptons) const {
0213   const reco::GenParticle* cand = nullptr;
0214   if (singleLepton(excludeTauLeptons)) {
0215     const reco::GenParticleCollection& partsColl = *parts_;
0216     const reco::GenParticle& singleLep = *singleLepton(excludeTauLeptons);
0217     for (unsigned int i = 0; i < partsColl.size(); ++i) {
0218       if (std::abs(partsColl[i].pdgId()) == TopDecayID::tID &&
0219           reco::flavour(singleLep) == reco::flavour(partsColl[i])) {
0220         cand = &partsColl[i];
0221       }
0222     }
0223   }
0224   return cand;
0225 }
0226
0227 const reco::GenParticle* TtGenEvent::leptonicDecayB(bool excludeTauLeptons) const {
0228   const reco::GenParticle* cand = nullptr;
0229   if (singleLepton(excludeTauLeptons)) {
0230     const reco::GenParticleCollection& partsColl = *parts_;
0231     const reco::GenParticle& singleLep = *singleLepton(excludeTauLeptons);
0232     for (unsigned int i = 0; i < partsColl.size(); ++i) {
0233       if (std::abs(partsColl[i].pdgId()) == TopDecayID::bID &&
0234           reco::flavour(singleLep) != reco::flavour(partsColl[i])) {
0235         cand = &partsColl[i];
0236       }
0237     }
0238   }
0239   return cand;
0240 }
0241
0242 const reco::GenParticle* TtGenEvent::leptonicDecayW(bool excludeTauLeptons) const {
0243   const reco::GenParticle* cand = nullptr;
0244   if (singleLepton(excludeTauLeptons)) {
0245     const reco::GenParticleCollection& partsColl = *parts_;
0246     const reco::GenParticle& singleLep = *singleLepton(excludeTauLeptons);
0247     for (unsigned int i = 0; i < partsColl.size(); ++i) {
0248       if (std::abs(partsColl[i].pdgId()) == TopDecayID::WID &&
0249           reco::flavour(singleLep) == -reco::flavour(partsColl[i])) {
0250         // PDG Id:13=mu- 24=W+ (+24)->(-13) (-24)->(+13) opposite sign
0251         cand = &partsColl[i];
0252       }
0253     }
0254   }
0255   return cand;
0256 }
0257
0258 const reco::GenParticle* TtGenEvent::leptonicDecayTop(bool excludeTauLeptons) const {
0259   const reco::GenParticle* cand = nullptr;
0260   if (singleLepton(excludeTauLeptons)) {
0261     const reco::GenParticleCollection& partsColl = *parts_;
0262     const reco::GenParticle& singleLep = *singleLepton(excludeTauLeptons);
0263     for (unsigned int i = 0; i < partsColl.size(); ++i) {
0264       if (std::abs(partsColl[i].pdgId()) == TopDecayID::tID &&
0265           reco::flavour(singleLep) != reco::flavour(partsColl[i])) {
0266         cand = &partsColl[i];
0267       }
0268     }
0269   }
0270   return cand;
0271 }
0272
0273 std::vector<const reco::GenParticle*> TtGenEvent::leptonicDecayTopRadiation(bool excludeTauLeptons) const {
0274   if (leptonicDecayTop(excludeTauLeptons)) {
0275     return (leptonicDecayTop(excludeTauLeptons)->pdgId() > 0 ? radiatedGluons(TopDecayID::tID)
0276                                                              : radiatedGluons(-TopDecayID::tID));
0277   }
0278   std::vector<const reco::GenParticle*> rad;
0279   return (rad);
0280 }
0281
0282 std::vector<const reco::GenParticle*> TtGenEvent::hadronicDecayTopRadiation(bool excludeTauLeptons) const {
0283   if (hadronicDecayTop(excludeTauLeptons)) {
0284     return (hadronicDecayTop(excludeTauLeptons)->pdgId() > 0 ? radiatedGluons(TopDecayID::tID)
0285                                                              : radiatedGluons(-TopDecayID::tID));
0286   }
0287   std::vector<const reco::GenParticle*> rad;
0288   return (rad);
0289 }