TopObjects/src/TopGenEvent.cc

0001 #include "FWCore/Utilities/interface/EDMException.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003
0004 #include "CommonTools/CandUtils/interface/pdgIdUtils.h"
0005 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0006 #include "AnalysisDataFormats/TopObjects/interface/TopGenEvent.h"
0007
0008 /// default contructor
0009 TopGenEvent::TopGenEvent(reco::GenParticleRefProd& decaySubset, reco::GenParticleRefProd& initSubset) {
0010   parts_ = decaySubset;
0011   initPartons_ = initSubset;
0012 }
0013
0014 const reco::GenParticle* TopGenEvent::candidate(int id, unsigned int parentId) const {
0015   const reco::GenParticle* cand = nullptr;
0016   const reco::GenParticleCollection& partsColl = *parts_;
0017   for (unsigned int i = 0; i < partsColl.size(); ++i) {
0018     if (partsColl[i].pdgId() == id) {
0019       if (parentId == 0 ? true : partsColl[i].mother() && std::abs(partsColl[i].mother()->pdgId()) == (int)parentId) {
0020         cand = &partsColl[i];
0021       }
0022     }
0023   }
0024   return cand;
0025 }
0026
0027 void TopGenEvent::print() const {
0028   edm::LogVerbatim log("TopGenEvent");
0029   log << "\n"
0030       << "--------------------------------------\n"
0031       << "- Dump TopGenEvent Content           -\n"
0032       << "--------------------------------------\n";
0033   for (reco::GenParticleCollection::const_iterator part = parts_->begin(); part < parts_->end(); ++part) {
0034     log << "pdgId:" << std::setw(5) << part->pdgId() << ", "
0035         << "mass:" << std::setw(11) << part->p4().mass() << ", "
0036         << "energy:" << std::setw(11) << part->energy() << ", "
0037         << "pt:" << std::setw(11) << part->pt() << "\n";
0038   }
0039 }
0040
0041 int TopGenEvent::numberOfLeptons(bool fromWBoson) const {
0042   int lep = 0;
0043   const reco::GenParticleCollection& partsColl = *parts_;
0044   for (unsigned int i = 0; i < partsColl.size(); ++i) {
0045     if (reco::isLepton(partsColl[i])) {
0046       if (fromWBoson) {
0047         if (partsColl[i].mother() && std::abs(partsColl[i].mother()->pdgId()) == TopDecayID::WID) {
0048           ++lep;
0049         }
0050       } else {
0051         ++lep;
0052       }
0053     }
0054   }
0055   return lep;
0056 }
0057
0058 int TopGenEvent::numberOfLeptons(WDecay::LeptonType typeRestriction, bool fromWBoson) const {
0059   int leptonType = -1;
0060   switch (typeRestriction) {
0061       // resolve whether or not there is
0062       // any restriction in lepton types
0063     case WDecay::kElec:
0064       leptonType = TopDecayID::elecID;
0065       break;
0066     case WDecay::kMuon:
0067       leptonType = TopDecayID::muonID;
0068       break;
0069     case WDecay::kTau:
0070       leptonType = TopDecayID::tauID;
0071       break;
0072     case WDecay::kNone:
0073       break;
0074   }
0075   int lep = 0;
0076   const reco::GenParticleCollection& partsColl = *parts_;
0077   for (unsigned int i = 0; i < partsColl.size(); ++i) {
0078     if (fromWBoson) {
0079       // restrict to particles originating from the W boson
0080       if (!(partsColl[i].mother() && std::abs(partsColl[i].mother()->pdgId()) == TopDecayID::WID)) {
0081         continue;
0082       }
0083     }
0084     if (leptonType > 0) {
0085       // in case of lepton type restriction
0086       if (std::abs(partsColl[i].pdgId()) == leptonType) {
0087         ++lep;
0088       }
0089     } else {
0090       // take any lepton type into account else
0091       if (reco::isLepton(partsColl[i])) {
0092         ++lep;
0093       }
0094     }
0095   }
0096   return lep;
0097 }
0098
0099 int TopGenEvent::numberOfBQuarks(bool fromTopQuark) const {
0100   int bq = 0;
0101   const reco::GenParticleCollection& partsColl = *parts_;
0102   for (unsigned int i = 0; i < partsColl.size(); ++i) {
0103     //depend if radiation qqbar are included or not
0104     if (std::abs(partsColl[i].pdgId()) == TopDecayID::bID) {
0105       if (fromTopQuark) {
0106         if (partsColl[i].mother() && std::abs(partsColl[i].mother()->pdgId()) == TopDecayID::tID) {
0107           ++bq;
0108         }
0109       } else {
0110         ++bq;
0111       }
0112     }
0113   }
0114   return bq;
0115 }
0116
0117 std::vector<const reco::GenParticle*> TopGenEvent::topSisters() const {
0118   std::vector<const reco::GenParticle*> sisters;
0119   for (reco::GenParticleCollection::const_iterator part = parts_->begin(); part < parts_->end(); ++part) {
0120     if (part->numberOfMothers() == 0 && std::abs(part->pdgId()) != TopDecayID::tID) {
0121       // choose top sister which do not have a
0122       // mother and are whether top nor anti-top
0123       if (dynamic_cast<const reco::GenParticle*>(&(*part)) == nullptr) {
0124         throw edm::Exception(edm::errors::InvalidReference, "Not a GenParticle");
0125       }
0126       sisters.push_back(part->clone());
0127     }
0128   }
0129   return sisters;
0130 }
0131
0132 const reco::GenParticle* TopGenEvent::daughterQuarkOfTop(bool invertCharge) const {
0133   const reco::GenParticle* cand = nullptr;
0134   for (reco::GenParticleCollection::const_iterator top = parts_->begin(); top < parts_->end(); ++top) {
0135     if (top->pdgId() == (invertCharge ? -TopDecayID::tID : TopDecayID::tID)) {
0136       for (reco::GenParticle::const_iterator quark = top->begin(); quark < top->end(); ++quark) {
0137         if (std::abs(quark->pdgId()) <= TopDecayID::bID) {
0138           cand = dynamic_cast<const reco::GenParticle*>(&(*quark));
0139           if (cand == nullptr) {
0140             throw edm::Exception(edm::errors::InvalidReference, "Not a GenParticle");
0141           }
0142           break;
0143         }
0144       }
0145     }
0146   }
0147   return cand;
0148 }
0149
0150 const reco::GenParticle* TopGenEvent::daughterQuarkOfWPlus(bool invertQuarkCharge, bool invertBosonCharge) const {
0151   const reco::GenParticle* cand = nullptr;
0152   const reco::GenParticleCollection& partsColl = *parts_;
0153   for (unsigned int i = 0; i < partsColl.size(); ++i) {
0154     if (partsColl[i].mother() &&
0155         partsColl[i].mother()->pdgId() == (invertBosonCharge ? -TopDecayID::WID : TopDecayID::WID) &&
0156         std::abs(partsColl[i].pdgId()) <= TopDecayID::bID &&
0157         (invertQuarkCharge ? reco::flavour(partsColl[i]) < 0 : reco::flavour(partsColl[i]) > 0)) {
0158       cand = &partsColl[i];
0159     }
0160   }
0161   return cand;
0162 }
0163
0164 std::vector<const reco::GenParticle*> TopGenEvent::lightQuarks(bool includingBQuarks) const {
0165   std::vector<const reco::GenParticle*> lightQuarks;
0166   for (reco::GenParticleCollection::const_iterator part = parts_->begin(); part < parts_->end(); ++part) {
0167     if ((includingBQuarks && std::abs(part->pdgId()) == TopDecayID::bID) || std::abs(part->pdgId()) < TopDecayID::bID) {
0168       if (dynamic_cast<const reco::GenParticle*>(&(*part)) == nullptr) {
0169         throw edm::Exception(edm::errors::InvalidReference, "Not a GenParticle");
0170       }
0171       lightQuarks.push_back(part->clone());
0172     }
0173   }
0174   return lightQuarks;
0175 }
0176
0177 std::vector<const reco::GenParticle*> TopGenEvent::radiatedGluons(int pdgId) const {
0178   std::vector<const reco::GenParticle*> rads;
0179   for (reco::GenParticleCollection::const_iterator part = parts_->begin(); part < parts_->end(); ++part) {
0180     if (part->mother() && part->mother()->pdgId() == pdgId) {
0181       if (part->pdgId() == TopDecayID::glueID) {
0182         if (dynamic_cast<const reco::GenParticle*>(&(*part)) == nullptr) {
0183           throw edm::Exception(edm::errors::InvalidReference, "Not a GenParticle");
0184         }
0185       }
0186       rads.push_back(part->clone());
0187     }
0188   }
0189   return rads;
0190 }