JetMCAlgos/plugins/ttHFGenFilter.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    ttHFGenFilter/ttHFGenFilter
0004 // Class:      ttHFGenFilter
0005 //
0006 /**\class ttHFGenFilter ttHFGenFilter.cc ttHFGenFilter/ttHFGenFilter/plugins/ttHFGenFilter.cc
0007
0008  Description: Filters ttbar events, where additional b-hadrons come from ISR or hard process
0009
0010  Implementation:
0011      The ttHFGenFilter is a edm::Filter, that returns true, if a tt+jets event contains a b-hadron, that does not come from top decay or when the b-hadron comes from ISR.
0012     The filter uses the GenHFHadronMatcher to classify, whether the b-hadron is an "additional b-hadron". This is done in the ttHFGenFilter::HasAdditionalBHadron() function.
0013     To classify, whether the b-hadron comes from ISR:
0014         first all mothers of the top-topbar pair are found
0015         then the mother-chain of all b-hadrons is checked if it contains at least one mother, that is also in the mother-chain of the top-topbar pair
0016
0017 */
0018 //
0019 // Original Author:  Andrej Saibel
0020 //         Created:  Tue, 05 Jul 2016 09:36:09 GMT
0021 //
0022 // VERY IMPORTANT: when running this code, you should make sure, that the GenHFHadronMatcher runs in onlyJetClusteredHadrons = cms.bool(False) mode
0023 //                 If you want to filter all events with additional b-hadrons use: OnlyHardProcessBHadrons = cms.bool(False)
0024
0025 // system include files
0026 #include <memory>
0027
0028 // user include files
0029 #include "FWCore/Framework/interface/Frameworkfwd.h"
0030 #include "FWCore/Framework/interface/stream/EDFilter.h"
0031
0032 #include "FWCore/Framework/interface/Event.h"
0033 #include "FWCore/Framework/interface/MakerMacros.h"
0034 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0035
0036 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0037 #include "FWCore/Utilities/interface/StreamID.h"
0038 #include "FWCore/Utilities/interface/InputTag.h"
0039
0040 //
0041 // class declaration
0042 //
0043
0044 class ttHFGenFilter : public edm::stream::EDFilter<> {
0045 public:
0046   explicit ttHFGenFilter(const edm::ParameterSet&);
0047
0048   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0049
0050 private:
0051   bool filter(edm::Event&, const edm::EventSetup&) override;
0052
0053   virtual bool HasAdditionalBHadron(const std::vector<int>&,
0054                                     const std::vector<int>&,
0055                                     const std::vector<reco::GenParticle>&,
0056                                     std::vector<const reco::Candidate*>&);
0057   virtual bool analyzeMothersRecursive(const reco::Candidate*, std::vector<const reco::Candidate*>& AllTopMothers);
0058   virtual std::vector<const reco::Candidate*> GetTops(const std::vector<reco::GenParticle>&,
0059                                                       std::vector<const reco::Candidate*>& AllTopMothers);
0060   virtual void FindAllTopMothers(const reco::Candidate* particle, std::vector<const reco::Candidate*>&);
0061
0062   const edm::EDGetTokenT<reco::GenParticleCollection> genParticlesToken_;
0063   const edm::EDGetTokenT<std::vector<int> > genBHadFlavourToken_;
0064   const edm::EDGetTokenT<std::vector<int> > genBHadFromTopWeakDecayToken_;
0065   const edm::EDGetTokenT<std::vector<reco::GenParticle> > genBHadPlusMothersToken_;
0066   const edm::EDGetTokenT<std::vector<std::vector<int> > > genBHadPlusMothersIndicesToken_;
0067   const edm::EDGetTokenT<std::vector<int> > genBHadIndexToken_;
0068   bool OnlyHardProcessBHadrons_;
0069   bool taggingMode_;
0070 };
0071
0072 //
0073 // constants, enums and typedefs
0074 //
0075
0076 //
0077 // static data member definitions
0078 //
0079
0080 //
0081 // constructors and destructor
0082 //
0083 ttHFGenFilter::ttHFGenFilter(const edm::ParameterSet& iConfig)
0084     : genParticlesToken_(consumes<reco::GenParticleCollection>(iConfig.getParameter<edm::InputTag>("genParticles"))),
0085       genBHadFlavourToken_(consumes<std::vector<int> >(iConfig.getParameter<edm::InputTag>("genBHadFlavour"))),
0086       genBHadFromTopWeakDecayToken_(
0087           consumes<std::vector<int> >(iConfig.getParameter<edm::InputTag>("genBHadFromTopWeakDecay"))),
0088       genBHadPlusMothersToken_(
0089           consumes<std::vector<reco::GenParticle> >(iConfig.getParameter<edm::InputTag>("genBHadPlusMothers"))),
0090       genBHadPlusMothersIndicesToken_(
0091           consumes<std::vector<std::vector<int> > >(iConfig.getParameter<edm::InputTag>("genBHadPlusMothersIndices"))),
0092       genBHadIndexToken_(consumes<std::vector<int> >(iConfig.getParameter<edm::InputTag>("genBHadIndex"))) {
0093   //now do what ever initialization is needed
0094   OnlyHardProcessBHadrons_ = iConfig.getParameter<bool>("OnlyHardProcessBHadrons");
0095   taggingMode_ = iConfig.getParameter<bool>("taggingMode");
0096
0097   produces<bool>();
0098 }
0099
0100 //
0101 // member functions
0102 //
0103
0104 // ------------ method called on each new Event  ------------
0105 bool ttHFGenFilter::filter(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0106   using namespace std;
0107
0108   //get GenParticleCollection
0109   edm::Handle<reco::GenParticleCollection> genParticles;
0110   iEvent.getByToken(genParticlesToken_, genParticles);
0111
0112   //get Information on B-Hadrons
0113   // the information is calculated in GenHFHadronMatcher
0114
0115   edm::Handle<std::vector<int> > genBHadFlavour;
0116   iEvent.getByToken(genBHadFlavourToken_, genBHadFlavour);
0117
0118   edm::Handle<std::vector<int> > genBHadFromTopWeakDecay;
0119   iEvent.getByToken(genBHadFromTopWeakDecayToken_, genBHadFromTopWeakDecay);
0120
0121   edm::Handle<std::vector<reco::GenParticle> > genBHadPlusMothers;
0122   iEvent.getByToken(genBHadPlusMothersToken_, genBHadPlusMothers);
0123
0124   edm::Handle<std::vector<std::vector<int> > > genBHadPlusMothersIndices;
0125   iEvent.getByToken(genBHadPlusMothersIndicesToken_, genBHadPlusMothersIndices);
0126
0127   edm::Handle<std::vector<int> > genBHadIndex;
0128   iEvent.getByToken(genBHadIndexToken_, genBHadIndex);
0129
0130   //check whether the event has additional b-hadrons not coming from top/topbar decay
0131
0132   std::vector<const reco::Candidate*> AllTopMothers;
0133   std::vector<const reco::Candidate*> Tops = GetTops(*genParticles, AllTopMothers);
0134
0135   bool pass = HasAdditionalBHadron(*genBHadIndex, *genBHadFlavour, *genBHadPlusMothers, AllTopMothers);
0136
0137   iEvent.put(std::make_unique<bool>(pass));
0138
0139   return taggingMode_ || pass;
0140 }
0141
0142 bool ttHFGenFilter::HasAdditionalBHadron(const std::vector<int>& genBHadIndex,
0143                                          const std::vector<int>& genBHadFlavour,
0144                                          const std::vector<reco::GenParticle>& genBHadPlusMothers,
0145                                          std::vector<const reco::Candidate*>& AllTopMothers) {
0146   for (unsigned int i = 0; i < genBHadIndex.size(); i++) {
0147     const reco::GenParticle* bhadron = genBHadIndex[i] >= 0 && genBHadIndex[i] < int(genBHadPlusMothers.size())
0148                                            ? &(genBHadPlusMothers[genBHadIndex[i]])
0149                                            : nullptr;
0150     int motherflav = genBHadFlavour[i];
0151     bool from_tth = (abs(motherflav) == 6 || abs(motherflav) == 25);  //b-hadron comes from top or higgs decay
0152     bool fromhp = false;
0153
0154     if (!OnlyHardProcessBHadrons_) {
0155       if (bhadron != nullptr && !from_tth) {
0156         return true;
0157       }
0158     }
0159
0160     if (OnlyHardProcessBHadrons_) {
0161       if (bhadron != nullptr && !from_tth) {
0162         //std::cout << "PT: " << bhadron->pt() << " , eta: " << bhadron->eta() << std::endl;
0163
0164         fromhp = analyzeMothersRecursive(bhadron, AllTopMothers);
0165         if (fromhp) {
0166           return true;
0167         }
0168       }
0169       if (i == genBHadIndex.size() - 1) {
0170         return false;
0171       }
0172     }
0173   }
0174
0175   return false;
0176 }
0177
0178 //recursive function, that loops over all chain particles  until it finds the protons
0179 bool ttHFGenFilter::analyzeMothersRecursive(const reco::Candidate* particle,
0180                                             std::vector<const reco::Candidate*>& AllTopMothers) {
0181   if (particle->status() > 20 && particle->status() < 30) {  //particle comes from hardest process in event
0182     return true;
0183   }
0184   for (unsigned int k = 0; k < AllTopMothers.size(); k++) {
0185     if (particle == AllTopMothers[k]) {  //particle comes from  ISR
0186       return true;
0187     }
0188   }
0189   bool isFromHardProcess = false;
0190   for (unsigned int i = 0; i < particle->numberOfMothers(); i++) {
0191     const reco::Candidate* mother = particle->mother(i);
0192     isFromHardProcess = analyzeMothersRecursive(mother, AllTopMothers);
0193     if (isFromHardProcess) {
0194       return true;
0195     }
0196   }
0197   return isFromHardProcess;
0198 }
0199
0200 std::vector<const reco::Candidate*> ttHFGenFilter::GetTops(const std::vector<reco::GenParticle>& genParticles,
0201                                                            std::vector<const reco::Candidate*>& AllTopMothers) {
0202   //  std::vector<const reco::Candidate*> tops;
0203   const reco::GenParticle* firstTop = nullptr;
0204   const reco::GenParticle* firstTopBar = nullptr;
0205   bool foundTop = false;
0206   bool foundTopBar = false;
0207   std::vector<const reco::GenParticle*> Tops;
0208
0209   //loop over all genParticles and find  a Top and AntiTop quark
0210   //then find all mothers of the Tops
0211   //this is then used to  if the b-hadron is an inital state radiation particle
0212
0213   for (reco::GenParticleCollection::const_iterator i_particle = genParticles.begin(); i_particle != genParticles.end();
0214        ++i_particle) {
0215     const reco::GenParticle* thisParticle = &*i_particle;
0216     if (!foundTop && thisParticle->pdgId() == 6) {
0217       firstTop = thisParticle;
0218       for (unsigned int i = 0; i < firstTop->numberOfMothers(); i++) {
0219         FindAllTopMothers(thisParticle->mother(i), AllTopMothers);
0220       }
0221       foundTop = true;
0222     } else if (!foundTopBar && thisParticle->pdgId() == -6) {
0223       firstTopBar = thisParticle;
0224       for (unsigned int i = 0; i < firstTopBar->numberOfMothers(); i++) {
0225         FindAllTopMothers(firstTopBar->mother(i), AllTopMothers);
0226       }
0227       foundTopBar = true;
0228     }
0229     if (foundTopBar && foundTop) {
0230       //tops.push_back(firstTop);
0231       //tops.push_back(firstTopBar);
0232       //return tops;
0233       break;
0234     }
0235   }
0236   return AllTopMothers;
0237 }
0238
0239 // finds all mothers of "particle", that are not tops or protons
0240
0241 void ttHFGenFilter::FindAllTopMothers(const reco::Candidate* particle,
0242                                       std::vector<const reco::Candidate*>& AllTopMothers) {
0243   // std::cout << "particle mother: " << particle->mother(0) << std::endl;
0244   for (unsigned int i = 0; i < particle->numberOfMothers(); i++) {
0245     if (abs(particle->mother(i)->pdgId()) != 6 && particle->mother(i)->pdgId() != 2212) {
0246       AllTopMothers.push_back(particle->mother(i));
0247       if (particle->mother(i)->pdgId() != 2212 || particle->mother(i)->numberOfMothers() > 1) {
0248         //std::cout << "Size of vector in loop = " << AllTopMothers.size() << std::endl;
0249         FindAllTopMothers(particle->mother(i), AllTopMothers);
0250       }
0251     }
0252   }
0253 }
0254
0255 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0256 void ttHFGenFilter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0257   //The following says we do not know what parameters are allowed so do no validation
0258   // Please change this to state exactly what you do use, even if it is no parameters
0259   edm::ParameterSetDescription desc;
0260   desc.setUnknown();
0261   descriptions.addDefault(desc);
0262 }
0263 //define this as a plug-in
0264 DEFINE_FWK_MODULE(ttHFGenFilter);