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   ~ttHFGenFilter() override;
0048
0049   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0050
0051 private:
0052   void beginStream(edm::StreamID) override;
0053   bool filter(edm::Event&, const edm::EventSetup&) override;
0054   void endStream() override;
0055
0056   virtual bool HasAdditionalBHadron(const std::vector<int>&,
0057                                     const std::vector<int>&,
0058                                     const std::vector<reco::GenParticle>&,
0059                                     std::vector<const reco::Candidate*>&);
0060   virtual bool analyzeMothersRecursive(const reco::Candidate*, std::vector<const reco::Candidate*>& AllTopMothers);
0061   virtual std::vector<const reco::Candidate*> GetTops(const std::vector<reco::GenParticle>&,
0062                                                       std::vector<const reco::Candidate*>& AllTopMothers);
0063   virtual void FindAllTopMothers(const reco::Candidate* particle, std::vector<const reco::Candidate*>&);
0064
0065   const edm::EDGetTokenT<reco::GenParticleCollection> genParticlesToken_;
0066   const edm::EDGetTokenT<std::vector<int> > genBHadFlavourToken_;
0067   const edm::EDGetTokenT<std::vector<int> > genBHadFromTopWeakDecayToken_;
0068   const edm::EDGetTokenT<std::vector<reco::GenParticle> > genBHadPlusMothersToken_;
0069   const edm::EDGetTokenT<std::vector<std::vector<int> > > genBHadPlusMothersIndicesToken_;
0070   const edm::EDGetTokenT<std::vector<int> > genBHadIndexToken_;
0071   bool OnlyHardProcessBHadrons_;
0072   bool taggingMode_;
0073 };
0074
0075 //
0076 // constants, enums and typedefs
0077 //
0078
0079 //
0080 // static data member definitions
0081 //
0082
0083 //
0084 // constructors and destructor
0085 //
0086 ttHFGenFilter::ttHFGenFilter(const edm::ParameterSet& iConfig)
0087     : genParticlesToken_(consumes<reco::GenParticleCollection>(iConfig.getParameter<edm::InputTag>("genParticles"))),
0088       genBHadFlavourToken_(consumes<std::vector<int> >(iConfig.getParameter<edm::InputTag>("genBHadFlavour"))),
0089       genBHadFromTopWeakDecayToken_(
0090           consumes<std::vector<int> >(iConfig.getParameter<edm::InputTag>("genBHadFromTopWeakDecay"))),
0091       genBHadPlusMothersToken_(
0092           consumes<std::vector<reco::GenParticle> >(iConfig.getParameter<edm::InputTag>("genBHadPlusMothers"))),
0093       genBHadPlusMothersIndicesToken_(
0094           consumes<std::vector<std::vector<int> > >(iConfig.getParameter<edm::InputTag>("genBHadPlusMothersIndices"))),
0095       genBHadIndexToken_(consumes<std::vector<int> >(iConfig.getParameter<edm::InputTag>("genBHadIndex"))) {
0096   //now do what ever initialization is needed
0097   OnlyHardProcessBHadrons_ = iConfig.getParameter<bool>("OnlyHardProcessBHadrons");
0098   taggingMode_ = iConfig.getParameter<bool>("taggingMode");
0099
0100   produces<bool>();
0101 }
0102
0103 ttHFGenFilter::~ttHFGenFilter() {
0104   // do anything here that needs to be done at destruction time
0105   // (e.g. close files, deallocate resources etc.)
0106 }
0107
0108 //
0109 // member functions
0110 //
0111
0112 // ------------ method called on each new Event  ------------
0113 bool ttHFGenFilter::filter(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0114   using namespace std;
0115
0116   //get GenParticleCollection
0117   edm::Handle<reco::GenParticleCollection> genParticles;
0118   iEvent.getByToken(genParticlesToken_, genParticles);
0119
0120   //get Information on B-Hadrons
0121   // the information is calculated in GenHFHadronMatcher
0122
0123   edm::Handle<std::vector<int> > genBHadFlavour;
0124   iEvent.getByToken(genBHadFlavourToken_, genBHadFlavour);
0125
0126   edm::Handle<std::vector<int> > genBHadFromTopWeakDecay;
0127   iEvent.getByToken(genBHadFromTopWeakDecayToken_, genBHadFromTopWeakDecay);
0128
0129   edm::Handle<std::vector<reco::GenParticle> > genBHadPlusMothers;
0130   iEvent.getByToken(genBHadPlusMothersToken_, genBHadPlusMothers);
0131
0132   edm::Handle<std::vector<std::vector<int> > > genBHadPlusMothersIndices;
0133   iEvent.getByToken(genBHadPlusMothersIndicesToken_, genBHadPlusMothersIndices);
0134
0135   edm::Handle<std::vector<int> > genBHadIndex;
0136   iEvent.getByToken(genBHadIndexToken_, genBHadIndex);
0137
0138   //check whether the event has additional b-hadrons not coming from top/topbar decay
0139
0140   std::vector<const reco::Candidate*> AllTopMothers;
0141   std::vector<const reco::Candidate*> Tops = GetTops(*genParticles, AllTopMothers);
0142
0143   bool pass = HasAdditionalBHadron(*genBHadIndex, *genBHadFlavour, *genBHadPlusMothers, AllTopMothers);
0144
0145   iEvent.put(std::make_unique<bool>(pass));
0146
0147   return taggingMode_ || pass;
0148 }
0149
0150 bool ttHFGenFilter::HasAdditionalBHadron(const std::vector<int>& genBHadIndex,
0151                                          const std::vector<int>& genBHadFlavour,
0152                                          const std::vector<reco::GenParticle>& genBHadPlusMothers,
0153                                          std::vector<const reco::Candidate*>& AllTopMothers) {
0154   for (unsigned int i = 0; i < genBHadIndex.size(); i++) {
0155     const reco::GenParticle* bhadron = genBHadIndex[i] >= 0 && genBHadIndex[i] < int(genBHadPlusMothers.size())
0156                                            ? &(genBHadPlusMothers[genBHadIndex[i]])
0157                                            : nullptr;
0158     int motherflav = genBHadFlavour[i];
0159     bool from_tth = (abs(motherflav) == 6 || abs(motherflav) == 25);  //b-hadron comes from top or higgs decay
0160     bool fromhp = false;
0161
0162     if (!OnlyHardProcessBHadrons_) {
0163       if (bhadron != nullptr && !from_tth) {
0164         return true;
0165       }
0166     }
0167
0168     if (OnlyHardProcessBHadrons_) {
0169       if (bhadron != nullptr && !from_tth) {
0170         //std::cout << "PT: " << bhadron->pt() << " , eta: " << bhadron->eta() << std::endl;
0171
0172         fromhp = analyzeMothersRecursive(bhadron, AllTopMothers);
0173         if (fromhp) {
0174           return true;
0175         }
0176       }
0177       if (i == genBHadIndex.size() - 1) {
0178         return false;
0179       }
0180     }
0181   }
0182
0183   return false;
0184 }
0185
0186 //recursive function, that loops over all chain particles  until it finds the protons
0187 bool ttHFGenFilter::analyzeMothersRecursive(const reco::Candidate* particle,
0188                                             std::vector<const reco::Candidate*>& AllTopMothers) {
0189   if (particle->status() > 20 && particle->status() < 30) {  //particle comes from hardest process in event
0190     return true;
0191   }
0192   for (unsigned int k = 0; k < AllTopMothers.size(); k++) {
0193     if (particle == AllTopMothers[k]) {  //particle comes from  ISR
0194       return true;
0195     }
0196   }
0197   bool isFromHardProcess = false;
0198   for (unsigned int i = 0; i < particle->numberOfMothers(); i++) {
0199     const reco::Candidate* mother = particle->mother(i);
0200     isFromHardProcess = analyzeMothersRecursive(mother, AllTopMothers);
0201     if (isFromHardProcess) {
0202       return true;
0203     }
0204   }
0205   return isFromHardProcess;
0206 }
0207
0208 std::vector<const reco::Candidate*> ttHFGenFilter::GetTops(const std::vector<reco::GenParticle>& genParticles,
0209                                                            std::vector<const reco::Candidate*>& AllTopMothers) {
0210   //  std::vector<const reco::Candidate*> tops;
0211   const reco::GenParticle* firstTop = nullptr;
0212   const reco::GenParticle* firstTopBar = nullptr;
0213   bool foundTop = false;
0214   bool foundTopBar = false;
0215   std::vector<const reco::GenParticle*> Tops;
0216
0217   //loop over all genParticles and find  a Top and AntiTop quark
0218   //then find all mothers of the Tops
0219   //this is then used to  if the b-hadron is an inital state radiation particle
0220
0221   for (reco::GenParticleCollection::const_iterator i_particle = genParticles.begin(); i_particle != genParticles.end();
0222        ++i_particle) {
0223     const reco::GenParticle* thisParticle = &*i_particle;
0224     if (!foundTop && thisParticle->pdgId() == 6) {
0225       firstTop = thisParticle;
0226       for (unsigned int i = 0; i < firstTop->numberOfMothers(); i++) {
0227         FindAllTopMothers(thisParticle->mother(i), AllTopMothers);
0228       }
0229       foundTop = true;
0230     } else if (!foundTopBar && thisParticle->pdgId() == -6) {
0231       firstTopBar = thisParticle;
0232       for (unsigned int i = 0; i < firstTopBar->numberOfMothers(); i++) {
0233         FindAllTopMothers(firstTopBar->mother(i), AllTopMothers);
0234       }
0235       foundTopBar = true;
0236     }
0237     if (foundTopBar && foundTop) {
0238       //tops.push_back(firstTop);
0239       //tops.push_back(firstTopBar);
0240       //return tops;
0241       break;
0242     }
0243   }
0244   return AllTopMothers;
0245 }
0246
0247 // finds all mothers of "particle", that are not tops or protons
0248
0249 void ttHFGenFilter::FindAllTopMothers(const reco::Candidate* particle,
0250                                       std::vector<const reco::Candidate*>& AllTopMothers) {
0251   // std::cout << "particle mother: " << particle->mother(0) << std::endl;
0252   for (unsigned int i = 0; i < particle->numberOfMothers(); i++) {
0253     if (abs(particle->mother(i)->pdgId()) != 6 && particle->mother(i)->pdgId() != 2212) {
0254       AllTopMothers.push_back(particle->mother(i));
0255       if (particle->mother(i)->pdgId() != 2212 || particle->mother(i)->numberOfMothers() > 1) {
0256         //std::cout << "Size of vector in loop = " << AllTopMothers.size() << std::endl;
0257         FindAllTopMothers(particle->mother(i), AllTopMothers);
0258       }
0259     }
0260   }
0261 }
0262
0263 // ------------ method called once each stream before processing any runs, lumis or events  ------------
0264 void ttHFGenFilter::beginStream(edm::StreamID) {}
0265
0266 // ------------ method called once each stream after processing all runs, lumis and events  ------------
0267 void ttHFGenFilter::endStream() {}
0268
0269 // ------------ method called when starting to processes a run  ------------
0270 /*
0271 void
0272 ttHFGenFilter::beginRun(edm::Run const&, edm::EventSetup const&)
0273 {
0274 }
0275 */
0276
0277 // ------------ method called when ending the processing of a run  ------------
0278 /*
0279 void
0280 ttHFGenFilter::endRun(edm::Run const&, edm::EventSetup const&)
0281 {
0282 }
0283 */
0284
0285 // ------------ method called when starting to processes a luminosity block  ------------
0286 /*
0287 void
0288 ttHFGenFilter::beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
0289 {
0290 }
0291 */
0292
0293 // ------------ method called when ending the processing of a luminosity block  ------------
0294 /*
0295 void
0296 ttHFGenFilter::endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
0297 {
0298 }
0299 */
0300
0301 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0302 void ttHFGenFilter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0303   //The following says we do not know what parameters are allowed so do no validation
0304   // Please change this to state exactly what you do use, even if it is no parameters
0305   edm::ParameterSetDescription desc;
0306   desc.setUnknown();
0307   descriptions.addDefault(desc);
0308 }
0309 //define this as a plug-in
0310 DEFINE_FWK_MODULE(ttHFGenFilter);