TopSkimming/src/TtDecayChannelSelector.cc

0001 #include "FWCore/Utilities/interface/EDMException.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include "AnalysisDataFormats/TopObjects/interface/TopGenEvent.h"
0004 #include "TopQuarkAnalysis/TopSkimming/interface/TtDecayChannelSelector.h"
0005 #include "DataFormats/Candidate/interface/Candidate.h"
0006
0007 // static const string for status check in
0008 // TtDecayChannelSelector::search functions
0009 static const std::string kGenParticles = "genParticles";
0010
0011 // maximal number of possible leptonic decay
0012 // channels
0013 static const unsigned int kDecayChannels = 3;
0014
0015 TtDecayChannelSelector::TtDecayChannelSelector(const edm::ParameterSet& cfg)
0016     : invert_(cfg.getParameter<bool>("invert")),
0017       allowElectron_(false),
0018       allowMuon_(false),
0019       allow1Prong_(false),
0020       allow3Prong_(false) {
0021   // tau decays are not restricted if this PSet does not exist at all
0022   restrictTauDecays_ = cfg.existsAs<edm::ParameterSet>("restrictTauDecays");
0023   // determine allowed tau decays
0024   if (restrictTauDecays_) {
0025     edm::ParameterSet allowedTauDecays = cfg.getParameter<edm::ParameterSet>("restrictTauDecays");
0026     // tau decays are not restricted if none of the following parameters exists
0027     restrictTauDecays_ = (allowedTauDecays.existsAs<bool>("electron") || allowedTauDecays.existsAs<bool>("muon") ||
0028                           allowedTauDecays.existsAs<bool>("oneProng") || allowedTauDecays.existsAs<bool>("threeProng"));
0029     // specify the different possible restrictions of the tau decay channels
0030     allowElectron_ =
0031         (allowedTauDecays.existsAs<bool>("electron") ? allowedTauDecays.getParameter<bool>("electron") : false);
0032     allowMuon_ = (allowedTauDecays.existsAs<bool>("muon") ? allowedTauDecays.getParameter<bool>("muon") : false);
0033     allow1Prong_ =
0034         (allowedTauDecays.existsAs<bool>("oneProng") ? allowedTauDecays.getParameter<bool>("oneProng") : false);
0035     allow3Prong_ =
0036         (allowedTauDecays.existsAs<bool>("threeProng") ? allowedTauDecays.getParameter<bool>("threeProng") : false);
0037   }
0038   // allowed top decays PSet
0039   edm::ParameterSet allowedTopDecays = cfg.getParameter<edm::ParameterSet>("allowedTopDecays");
0040
0041   // fill decayBranchA_
0042   edm::ParameterSet decayBranchA = allowedTopDecays.getParameter<edm::ParameterSet>("decayBranchA");
0043   decayBranchA_.push_back(decayBranchA.getParameter<bool>("electron"));
0044   decayBranchA_.push_back(decayBranchA.getParameter<bool>("muon"));
0045   decayBranchA_.push_back(decayBranchA.getParameter<bool>("tau"));
0046
0047   // fill decay branchB_
0048   edm::ParameterSet decayBranchB = allowedTopDecays.getParameter<edm::ParameterSet>("decayBranchB");
0049   decayBranchB_.push_back(decayBranchB.getParameter<bool>("electron"));
0050   decayBranchB_.push_back(decayBranchB.getParameter<bool>("muon"));
0051   decayBranchB_.push_back(decayBranchB.getParameter<bool>("tau"));
0052
0053   // fill allowedDecays_
0054   for (unsigned int d = 0; d < kDecayChannels; ++d) {
0055     allowedDecays_.push_back(decayBranchA_[d] + decayBranchB_[d]);
0056   }
0057 }
0058
0059 bool TtDecayChannelSelector::operator()(const reco::GenParticleCollection& parts, std::string_view inputType) const {
0060   bool verbose =
0061       false;  // set this to true for debugging and add TtDecayChannelSelector category to the MessageLogger in your cfg file
0062   unsigned int iLep = 0;
0063   unsigned int iTop = 0, iBeauty = 0, iElec = 0, iMuon = 0, iTau = 0;
0064   for (reco::GenParticleCollection::const_iterator top = parts.begin(); top != parts.end(); ++top) {
0065     if (search(top, TopDecayID::tID, inputType)) {
0066       ++iTop;
0067       for (reco::GenParticle::const_iterator td = top->begin(); td != top->end(); ++td) {
0068         if (search(td, TopDecayID::bID, inputType)) {
0069           ++iBeauty;
0070         }
0071         if (search(td, TopDecayID::WID, inputType)) {
0072           for (reco::GenParticle::const_iterator wd = td->begin(); wd != td->end(); ++wd) {
0073             if (std::abs(wd->pdgId()) == TopDecayID::elecID) {
0074               ++iElec;
0075             }
0076             if (std::abs(wd->pdgId()) == TopDecayID::muonID) {
0077               ++iMuon;
0078             }
0079             if (std::abs(wd->pdgId()) == TopDecayID::tauID) {
0080               if (restrictTauDecays_) {
0081                 // count as iTau if it is leptonic, one-prong
0082                 // or three-prong and ignore increasing iLep
0083                 // though else
0084                 if (tauDecay(*wd)) {
0085                   ++iTau;
0086                 } else {
0087                   ++iLep;
0088                 }
0089               } else {
0090                 ++iTau;
0091               }
0092             }
0093           }
0094         }
0095       }
0096     }
0097   }
0098   if (verbose) {
0099     edm::LogVerbatim log("TtDecayChannelSelector");
0100     log << "----------------------"
0101         << "\n"
0102         << " iTop    : " << iTop << "\n"
0103         << " iBeauty : " << iBeauty << "\n"
0104         << " iElec   : " << iElec << "\n"
0105         << " iMuon   : " << iMuon << "\n"
0106         << " iTau    : " << iTau + iLep;
0107     if (restrictTauDecays_ && (iTau + iLep) > 0) {
0108       log << " (" << iTau << ")\n";
0109     } else {
0110       log << "\n";
0111     }
0112     log << "- - - - - - - - - - - "
0113         << "\n";
0114   }
0115   iLep += iElec + iMuon + iTau;
0116
0117   bool accept = false;
0118   unsigned int channel = decayChannel();
0119   if ((iTop == 2) && (iBeauty == 2)) {
0120     if (channel == iLep) {
0121       if (channel == 0) {
0122         // no lepton: accept without restriction we already
0123         // know that the number of leptons is correct
0124         accept = true;
0125       }
0126       if (channel == 1) {
0127         // one lepton: check that this one is allowed
0128         accept = (iElec && allowedDecays_[Elec]) || (iMuon && allowedDecays_[Muon]) || (iTau && allowedDecays_[Tau]);
0129       }
0130       if (channel == 2) {
0131         if (checkSum(allowedDecays_) == channel) {
0132           // no redundancy
0133           accept = (allowedDecays_[Elec] == (int)iElec) && (allowedDecays_[Muon] == (int)iMuon) &&
0134                    (allowedDecays_[Tau] == (int)iTau);
0135         } else {
0136           // reject events with wrong tau decays
0137           if (iElec + iMuon + iTau != channel) {
0138             accept = false;
0139           } else {
0140             if ((iElec == 2) || (iMuon == 2) || (iTau == 2)) {
0141               // same lepton twice: check that this is allowed.
0142               accept = (allowedDecays_[Elec] == (int)iElec) || (allowedDecays_[Muon] == (int)iMuon) ||
0143                        (allowedDecays_[Tau] == (int)iTau);
0144             } else {
0145               // two different leptons: look if there is a possible combination
0146               accept = (((iElec && decayBranchA_[Elec]) &&
0147                          ((iMuon && decayBranchB_[Muon]) || (iTau && decayBranchB_[Tau]))) ||
0148                         ((iMuon && decayBranchA_[Muon]) &&
0149                          ((iElec && decayBranchB_[Elec]) || (iTau && decayBranchB_[Tau]))) ||
0150                         ((iTau && decayBranchA_[Tau]) &&
0151                          ((iElec && decayBranchB_[Elec]) || (iMuon && decayBranchB_[Muon]))));
0152             }
0153           }
0154         }
0155       }
0156     }
0157     accept = ((!invert_ && accept) || (!(!invert_) && !accept));
0158   } else {
0159     edm::LogWarning("NoVtbDecay") << "Decay is not via Vtb";
0160   }
0161   if (verbose)
0162     edm::LogVerbatim("TtDecayChannelSelector") << " accept  : " << accept;
0163   return accept;
0164 }
0165
0166 bool TtDecayChannelSelector::search(reco::GenParticleCollection::const_iterator& part,
0167                                     int pdgId,
0168                                     std::string_view inputType) const {
0169   if (inputType == kGenParticles) {
0170     return (std::abs(part->pdgId()) == pdgId && part->status() == TopDecayID::unfrag) ? true : false;
0171   } else {
0172     return (std::abs(part->pdgId()) == pdgId) ? true : false;
0173   }
0174 }
0175
0176 bool TtDecayChannelSelector::search(reco::GenParticle::const_iterator& part,
0177                                     int pdgId,
0178                                     std::string_view inputType) const {
0179   if (inputType == kGenParticles) {
0180     return (std::abs(part->pdgId()) == pdgId && part->status() == TopDecayID::unfrag) ? true : false;
0181   } else {
0182     return (std::abs(part->pdgId()) == pdgId) ? true : false;
0183   }
0184 }
0185
0186 unsigned int TtDecayChannelSelector::countProngs(const reco::Candidate& part) const {
0187   // if stable, return 1 or 0
0188   if (part.status() == 1) {
0189     return (part.charge() != 0);
0190   }
0191   // if unstable, call recursively on daughters
0192   int prong = 0;
0193   for (reco::Candidate::const_iterator daughter = part.begin(); daughter != part.end(); ++daughter) {
0194     prong += countProngs(*daughter);
0195   }
0196   return prong;
0197 }
0198
0199 bool TtDecayChannelSelector::tauDecay(const reco::Candidate& tau) const {
0200   bool electronTau = false;
0201   bool muonTau = false;
0202   unsigned int nch = 0;
0203   // loop on tau decays, check for an elec
0204   // or muon and count charged particles
0205   for (reco::Candidate::const_iterator daughter = tau.begin(); daughter != tau.end(); ++daughter) {
0206     // if the tau daughter is again a tau, this means that the particle has
0207     // still to be propagated; in that case, return the result of the same
0208     // method applied on the daughter of the current particle
0209     if (daughter->pdgId() == tau.pdgId()) {
0210       return tauDecay(*daughter);
0211     }
0212     // check for electron from tau decay
0213     electronTau |= (std::abs(daughter->pdgId()) == TopDecayID::elecID);
0214     // check for muon from tau decay
0215     muonTau |= (std::abs(daughter->pdgId()) == TopDecayID::muonID);
0216     // count charged particles
0217     nch += countProngs(*daughter);
0218   }
0219   return ((allowElectron_ && electronTau) || (allowMuon_ && muonTau) ||
0220           (allow1Prong_ && !electronTau && !muonTau && nch == 1) ||
0221           (allow3Prong_ && !electronTau && !muonTau && nch == 3));
0222 }