Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:31:25 

0001 // -*- C++ -*-
0002 //
0003 // Package:    TopQuarkAnalysis/TopTools
0004 // Class:      GenTtbarCategorizer
0005 //
0006 /**\class GenTtbarCategorizer GenTtbarCategorizer.cc TopQuarkAnalysis/TopTools/plugins/GenTtbarCategorizer.cc
0007 
0008  Description: Categorization of different tt+xx processes, returning unique ID for each process as e.g. tt+bb, tt+b, tt+2b, tt+cc, ...
0009 
0010  Implementation:
0011      
0012      The classification scheme returns an ID per event, and works as follows:
0013      
0014      All jets in the following need to be in the acceptance as given by the config parameters |eta|, pt.
0015      A c jet must contain at least one c hadron and should contain no b hadrons
0016      
0017      First, b jets from top are identified, i.e. jets containing a b hadron from t->b decay
0018      They are encoded in the ID as numberOfBjetsFromTop*100, i.e.
0019      0xx: no b jets from top in acceptance
0020      1xx: 1 b jet from top in acceptance
0021      2xx: both b jets from top in acceptance
0022      
0023      Then, b jets from W are identified, i.e. jets containing a b hadron from W->b decay
0024      They are encoded in the ID as numberOfBjetsFromW*1000, i.e.
0025      0xxx: no b jets from W in acceptance
0026      1xxx: 1 b jet from W in acceptance
0027      2xxx: 2 b jets from W in acceptance
0028      
0029      Then, c jets from W are identified, i.e. jets containing a c hadron from W->c decay, but no b hadrons
0030      They are encoded in the ID as numberOfCjetsFromW*10000, i.e.
0031      0xxxx: no c jets from W in acceptance
0032      1xxxx: 1 c jet from W in acceptance
0033      2xxxx: 2 c jets from W in acceptance
0034      
0035      From the remaining jets, the ID is formed based on the additional b jets (IDs x5x) and c jets (IDs x4x) in the following order:
0036      x55: at least 2 additional b jets with at least two of them having >= 2 b hadrons in each
0037      x54: at least 2 additional b jets with one of them having >= 2 b hadrons, the others having =1 b hadron
0038      x53: at least 2 additional b jets with all having =1 b hadron
0039      x52: exactly 1 additional b jet having >=2 b hadrons
0040      x51: exactly 1 additional b jet having =1 b hadron
0041      x45: at least 2 additional c jets with at least two of them having >= 2 c hadrons in each
0042      x44: at least 2 additional c jets with one of them having >= 2 c hadrons, the others having =1 c hadron
0043      x43: at least 2 additional c jets with all having =1 c hadron
0044      x42: exactly 1 additional c jet having >=2 c hadrons
0045      x41: exactly 1 additional c jet having =1 c hadron
0046      x00: No additional b or c jet, i.e. only light flavour jets or no additional jets
0047 */
0048 //
0049 // Original Author:  Johannes Hauk, Nazar Bartosik
0050 //         Created:  Sun, 14 Jun 2015 19:42:58 GMT
0051 //
0052 //
0053 
0054 // system include files
0055 #include <memory>
0056 #include <algorithm>
0057 #include <functional>
0058 
0059 // user include files
0060 #include "FWCore/Framework/interface/Frameworkfwd.h"
0061 #include "FWCore/Framework/interface/global/EDProducer.h"
0062 #include "FWCore/Framework/interface/Event.h"
0063 #include "FWCore/Framework/interface/MakerMacros.h"
0064 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0065 #include "FWCore/Utilities/interface/InputTag.h"
0066 #include "DataFormats/JetReco/interface/GenJetCollection.h"
0067 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0068 
0069 //
0070 // class declaration
0071 //
0072 
0073 class GenTtbarCategorizer : public edm::global::EDProducer<> {
0074 public:
0075   explicit GenTtbarCategorizer(const edm::ParameterSet&);
0076 
0077   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0078 
0079 private:
0080   void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0081 
0082   std::vector<int> nHadronsOrderedJetIndices(const std::map<int, int>& m_jetIndex) const;
0083 
0084   // ----------member data ---------------------------
0085 
0086   // Jet configuration
0087   const double genJetPtMin_;
0088   const double genJetAbsEtaMax_;
0089 
0090   // Input tags
0091   const edm::EDGetTokenT<reco::GenJetCollection> genJetsToken_;
0092 
0093   const edm::EDGetTokenT<std::vector<int> > genBHadJetIndexToken_;
0094   const edm::EDGetTokenT<std::vector<int> > genBHadFlavourToken_;
0095   const edm::EDGetTokenT<std::vector<int> > genBHadFromTopWeakDecayToken_;
0096   const edm::EDGetTokenT<std::vector<reco::GenParticle> > genBHadPlusMothersToken_;
0097   const edm::EDGetTokenT<std::vector<std::vector<int> > > genBHadPlusMothersIndicesToken_;
0098   const edm::EDGetTokenT<std::vector<int> > genBHadIndexToken_;
0099   const edm::EDGetTokenT<std::vector<int> > genBHadLeptonHadronIndexToken_;
0100   const edm::EDGetTokenT<std::vector<int> > genBHadLeptonViaTauToken_;
0101 
0102   const edm::EDGetTokenT<std::vector<int> > genCHadJetIndexToken_;
0103   const edm::EDGetTokenT<std::vector<int> > genCHadFlavourToken_;
0104   const edm::EDGetTokenT<std::vector<int> > genCHadFromTopWeakDecayToken_;
0105   const edm::EDGetTokenT<std::vector<int> > genCHadBHadronIdToken_;
0106 };
0107 
0108 //
0109 // constants, enums and typedefs
0110 //
0111 
0112 //
0113 // static data member definitions
0114 //
0115 
0116 //
0117 // constructors and destructor
0118 //
0119 GenTtbarCategorizer::GenTtbarCategorizer(const edm::ParameterSet& iConfig)
0120     : genJetPtMin_(iConfig.getParameter<double>("genJetPtMin")),
0121       genJetAbsEtaMax_(iConfig.getParameter<double>("genJetAbsEtaMax")),
0122       genJetsToken_(consumes<reco::GenJetCollection>(iConfig.getParameter<edm::InputTag>("genJets"))),
0123       genBHadJetIndexToken_(consumes<std::vector<int> >(iConfig.getParameter<edm::InputTag>("genBHadJetIndex"))),
0124       genBHadFlavourToken_(consumes<std::vector<int> >(iConfig.getParameter<edm::InputTag>("genBHadFlavour"))),
0125       genBHadFromTopWeakDecayToken_(
0126           consumes<std::vector<int> >(iConfig.getParameter<edm::InputTag>("genBHadFromTopWeakDecay"))),
0127       genBHadPlusMothersToken_(
0128           consumes<std::vector<reco::GenParticle> >(iConfig.getParameter<edm::InputTag>("genBHadPlusMothers"))),
0129       genBHadPlusMothersIndicesToken_(
0130           consumes<std::vector<std::vector<int> > >(iConfig.getParameter<edm::InputTag>("genBHadPlusMothersIndices"))),
0131       genBHadIndexToken_(consumes<std::vector<int> >(iConfig.getParameter<edm::InputTag>("genBHadIndex"))),
0132       genBHadLeptonHadronIndexToken_(
0133           consumes<std::vector<int> >(iConfig.getParameter<edm::InputTag>("genBHadLeptonHadronIndex"))),
0134       genBHadLeptonViaTauToken_(
0135           consumes<std::vector<int> >(iConfig.getParameter<edm::InputTag>("genBHadLeptonViaTau"))),
0136       genCHadJetIndexToken_(consumes<std::vector<int> >(iConfig.getParameter<edm::InputTag>("genCHadJetIndex"))),
0137       genCHadFlavourToken_(consumes<std::vector<int> >(iConfig.getParameter<edm::InputTag>("genCHadFlavour"))),
0138       genCHadFromTopWeakDecayToken_(
0139           consumes<std::vector<int> >(iConfig.getParameter<edm::InputTag>("genCHadFromTopWeakDecay"))),
0140       genCHadBHadronIdToken_(consumes<std::vector<int> >(iConfig.getParameter<edm::InputTag>("genCHadBHadronId"))) {
0141   produces<int>("genTtbarId");
0142 }
0143 
0144 //
0145 // member functions
0146 //
0147 
0148 // ------------ method called to produce the data  ------------
0149 void GenTtbarCategorizer::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0150   // Access gen jets
0151   edm::Handle<reco::GenJetCollection> genJets;
0152   iEvent.getByToken(genJetsToken_, genJets);
0153 
0154   // Access B hadrons information
0155   edm::Handle<std::vector<int> > genBHadFlavour;
0156   iEvent.getByToken(genBHadFlavourToken_, genBHadFlavour);
0157 
0158   edm::Handle<std::vector<int> > genBHadJetIndex;
0159   iEvent.getByToken(genBHadJetIndexToken_, genBHadJetIndex);
0160 
0161   edm::Handle<std::vector<int> > genBHadFromTopWeakDecay;
0162   iEvent.getByToken(genBHadFromTopWeakDecayToken_, genBHadFromTopWeakDecay);
0163 
0164   edm::Handle<std::vector<reco::GenParticle> > genBHadPlusMothers;
0165   iEvent.getByToken(genBHadPlusMothersToken_, genBHadPlusMothers);
0166 
0167   edm::Handle<std::vector<std::vector<int> > > genBHadPlusMothersIndices;
0168   iEvent.getByToken(genBHadPlusMothersIndicesToken_, genBHadPlusMothersIndices);
0169 
0170   edm::Handle<std::vector<int> > genBHadIndex;
0171   iEvent.getByToken(genBHadIndexToken_, genBHadIndex);
0172 
0173   edm::Handle<std::vector<int> > genBHadLeptonHadronIndex;
0174   iEvent.getByToken(genBHadLeptonHadronIndexToken_, genBHadLeptonHadronIndex);
0175 
0176   edm::Handle<std::vector<int> > genBHadLeptonViaTau;
0177   iEvent.getByToken(genBHadLeptonViaTauToken_, genBHadLeptonViaTau);
0178 
0179   // Access C hadrons information
0180   edm::Handle<std::vector<int> > genCHadFlavour;
0181   iEvent.getByToken(genCHadFlavourToken_, genCHadFlavour);
0182 
0183   edm::Handle<std::vector<int> > genCHadJetIndex;
0184   iEvent.getByToken(genCHadJetIndexToken_, genCHadJetIndex);
0185 
0186   edm::Handle<std::vector<int> > genCHadFromTopWeakDecay;
0187   iEvent.getByToken(genCHadFromTopWeakDecayToken_, genCHadFromTopWeakDecay);
0188 
0189   edm::Handle<std::vector<int> > genCHadBHadronId;
0190   iEvent.getByToken(genCHadBHadronIdToken_, genCHadBHadronId);
0191 
0192   // Map <jet index, number of specific hadrons in jet>
0193   // B jets with b hadrons directly from t->b decay
0194   std::map<int, int> bJetFromTopIds;
0195   // B jets with b hadrons from W->b decay
0196   std::map<int, int> bJetFromWIds;
0197   // C jets with c hadrons from W->c decay
0198   std::map<int, int> cJetFromWIds;
0199   // B jets with b hadrons before top quark decay chain
0200   std::map<int, int> bJetAdditionalIds;
0201   // C jets with c hadrons before top quark decay chain
0202   std::map<int, int> cJetAdditionalIds;
0203 
0204   // Count number of specific b hadrons in each jet
0205   for (size_t hadronId = 0; hadronId < genBHadIndex->size(); ++hadronId) {
0206     // Index of jet associated to the hadron
0207     const int jetIndex = genBHadJetIndex->at(hadronId);
0208     // Skip hadrons which have no associated jet
0209     if (jetIndex < 0)
0210       continue;
0211     // Skip if jet is not in acceptance
0212     if (genJets->at(jetIndex).pt() < genJetPtMin_)
0213       continue;
0214     if (std::fabs(genJets->at(jetIndex).eta()) > genJetAbsEtaMax_)
0215       continue;
0216     // Flavour of the hadron's origin
0217     const int flavour = genBHadFlavour->at(hadronId);
0218     // Jet from t->b decay [pdgId(top)=6]
0219     if (std::abs(flavour) == 6) {
0220       if (bJetFromTopIds.count(jetIndex) < 1)
0221         bJetFromTopIds[jetIndex] = 1;
0222       else
0223         bJetFromTopIds[jetIndex]++;
0224       continue;
0225     }
0226     // Jet from W->b decay [pdgId(W)=24]
0227     if (std::abs(flavour) == 24) {
0228       if (bJetFromWIds.count(jetIndex) < 1)
0229         bJetFromWIds[jetIndex] = 1;
0230       else
0231         bJetFromWIds[jetIndex]++;
0232       continue;
0233     }
0234     // Identify jets with b hadrons not from top-quark or W-boson decay
0235     if (bJetAdditionalIds.count(jetIndex) < 1)
0236       bJetAdditionalIds[jetIndex] = 1;
0237     else
0238       bJetAdditionalIds[jetIndex]++;
0239   }
0240 
0241   // Cleaning up b jets from W->b decays
0242   for (std::map<int, int>::iterator it = bJetFromWIds.begin(); it != bJetFromWIds.end();) {
0243     // Cannot be a b jet from t->b decay
0244     if (bJetFromTopIds.count(it->first) > 0)
0245       bJetFromWIds.erase(it++);
0246     else
0247       ++it;
0248   }
0249 
0250   // Cleaning up additional b jets
0251   for (std::map<int, int>::iterator it = bJetAdditionalIds.begin(); it != bJetAdditionalIds.end();) {
0252     // Cannot be a b jet from t->b decay
0253     if (bJetFromTopIds.count(it->first) > 0)
0254       bJetAdditionalIds.erase(it++);
0255     // Cannot be a b jet from W->b decay
0256     else if (bJetFromWIds.count(it->first) > 0)
0257       bJetAdditionalIds.erase(it++);
0258     else
0259       ++it;
0260   }
0261 
0262   // Count number of specific c hadrons in each c jet
0263   for (size_t hadronId = 0; hadronId < genCHadJetIndex->size(); ++hadronId) {
0264     // Index of jet associated to the hadron
0265     const int jetIndex = genCHadJetIndex->at(hadronId);
0266     // Skip hadrons which have no associated jet
0267     if (jetIndex < 0)
0268       continue;
0269     // Skip c hadrons that are coming from b hadrons
0270     if (genCHadBHadronId->at(hadronId) >= 0)
0271       continue;
0272     // Skip if jet is not in acceptance
0273     if (genJets->at(jetIndex).pt() < genJetPtMin_)
0274       continue;
0275     if (std::fabs(genJets->at(jetIndex).eta()) > genJetAbsEtaMax_)
0276       continue;
0277     // Skip if jet contains a b hadron
0278     if (bJetFromTopIds.count(jetIndex) > 0)
0279       continue;
0280     if (bJetFromWIds.count(jetIndex) > 0)
0281       continue;
0282     if (bJetAdditionalIds.count(jetIndex) > 0)
0283       continue;
0284     // Flavour of the hadron's origin
0285     const int flavour = genCHadFlavour->at(hadronId);
0286     // Jet from W->c decay [pdgId(W)=24]
0287     if (std::abs(flavour) == 24) {
0288       if (cJetFromWIds.count(jetIndex) < 1)
0289         cJetFromWIds[jetIndex] = 1;
0290       else
0291         cJetFromWIds[jetIndex]++;
0292       continue;
0293     }
0294     // Identify jets with c hadrons not from W-boson decay
0295     if (cJetAdditionalIds.count(jetIndex) < 1)
0296       cJetAdditionalIds[jetIndex] = 1;
0297     else
0298       cJetAdditionalIds[jetIndex]++;
0299   }
0300 
0301   // Cleaning up additional c jets
0302   for (std::map<int, int>::iterator it = cJetAdditionalIds.begin(); it != cJetAdditionalIds.end();) {
0303     // Cannot be a c jet from W->c decay
0304     if (cJetFromWIds.count(it->first) > 0)
0305       cJetAdditionalIds.erase(it++);
0306     else
0307       ++it;
0308   }
0309 
0310   // Categorize event based on number of additional b/c jets
0311   // and number of corresponding hadrons in each of them
0312   int additionalJetEventId = bJetFromTopIds.size() * 100 + bJetFromWIds.size() * 1000 + cJetFromWIds.size() * 10000;
0313   // tt + 1 additional b jet
0314   if (bJetAdditionalIds.size() == 1) {
0315     const int nHadronsInJet = bJetAdditionalIds.begin()->second;
0316     // tt + 1 additional b jet from 1 additional b hadron
0317     if (nHadronsInJet == 1)
0318       additionalJetEventId += 51;
0319     // tt + 1 additional b jet from >=2 additional b hadrons
0320     else
0321       additionalJetEventId += 52;
0322   }
0323   // tt + >=2 additional b jets
0324   else if (bJetAdditionalIds.size() > 1) {
0325     // Check first two additional b jets (rare cases could have more)
0326     const std::vector<int> orderedJetIndices = nHadronsOrderedJetIndices(bJetAdditionalIds);
0327     int nHadronsInJet1 = bJetAdditionalIds.at(orderedJetIndices.at(0));
0328     int nHadronsInJet2 = bJetAdditionalIds.at(orderedJetIndices.at(1));
0329     // tt + >=2 additional b jets each from 1 additional b hadron
0330     if (std::max(nHadronsInJet1, nHadronsInJet2) == 1)
0331       additionalJetEventId += 53;
0332     // tt + >=2 additional b jets one of which from >=2 additional b hadrons
0333     else if (std::min(nHadronsInJet1, nHadronsInJet2) == 1 && std::max(nHadronsInJet1, nHadronsInJet2) > 1)
0334       additionalJetEventId += 54;
0335     // tt + >=2 additional b jets each from >=2 additional b hadrons
0336     else if (std::min(nHadronsInJet1, nHadronsInJet2) > 1)
0337       additionalJetEventId += 55;
0338   }
0339   // tt + no additional b jets
0340   else {
0341     // tt + 1 additional c jet
0342     if (cJetAdditionalIds.size() == 1) {
0343       const int nHadronsInJet = cJetAdditionalIds.begin()->second;
0344       // tt + 1 additional c jet from 1 additional c hadron
0345       if (nHadronsInJet == 1)
0346         additionalJetEventId += 41;
0347       // tt + 1 additional c jet from >=2 additional c hadrons
0348       else
0349         additionalJetEventId += 42;
0350     }
0351     // tt + >=2 additional c jets
0352     else if (cJetAdditionalIds.size() > 1) {
0353       // Check first two additional c jets (rare cases could have more)
0354       const std::vector<int> orderedJetIndices = nHadronsOrderedJetIndices(cJetAdditionalIds);
0355       int nHadronsInJet1 = cJetAdditionalIds.at(orderedJetIndices.at(0));
0356       int nHadronsInJet2 = cJetAdditionalIds.at(orderedJetIndices.at(1));
0357       // tt + >=2 additional c jets each from 1 additional c hadron
0358       if (std::max(nHadronsInJet1, nHadronsInJet2) == 1)
0359         additionalJetEventId += 43;
0360       // tt + >=2 additional c jets one of which from >=2 additional c hadrons
0361       else if (std::min(nHadronsInJet1, nHadronsInJet2) == 1 && std::max(nHadronsInJet1, nHadronsInJet2) > 1)
0362         additionalJetEventId += 44;
0363       // tt + >=2 additional c jets each from >=2 additional c hadrons
0364       else if (std::min(nHadronsInJet1, nHadronsInJet2) > 1)
0365         additionalJetEventId += 45;
0366     }
0367     // tt + no additional c jets
0368     else {
0369       // tt + light jets
0370       additionalJetEventId += 0;
0371     }
0372   }
0373 
0374   std::unique_ptr<int> ttbarId(new int);
0375   *ttbarId = additionalJetEventId;
0376   iEvent.put(std::move(ttbarId), "genTtbarId");
0377 }
0378 
0379 // ------------ method returns a vector of jet indices from the given map, sorted by N hadrons in descending order  ------------
0380 std::vector<int> GenTtbarCategorizer::nHadronsOrderedJetIndices(const std::map<int, int>& m_jetIndex) const {
0381   const int nElements = m_jetIndex.size();
0382   std::vector<std::pair<int, int> > v_jetNhadIndexPair;
0383   v_jetNhadIndexPair.reserve(nElements);
0384   for (std::map<int, int>::const_iterator it = m_jetIndex.begin(); it != m_jetIndex.end(); ++it) {
0385     const int jetIndex = it->first;
0386     const int nHadrons = it->second;
0387     v_jetNhadIndexPair.push_back(std::pair<int, int>(nHadrons, jetIndex));
0388   }
0389   // Sorting the vector of pairs by their key value
0390   std::sort(v_jetNhadIndexPair.begin(), v_jetNhadIndexPair.end(), std::greater<std::pair<int, int> >());
0391   // Building the vector of indices in the proper order
0392   std::vector<int> v_orderedJetIndices;
0393   v_orderedJetIndices.reserve(nElements);
0394   for (std::vector<std::pair<int, int> >::const_iterator it = v_jetNhadIndexPair.begin();
0395        it != v_jetNhadIndexPair.end();
0396        ++it) {
0397     v_orderedJetIndices.push_back(it->second);
0398   }
0399 
0400   return v_orderedJetIndices;
0401 }
0402 
0403 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0404 void GenTtbarCategorizer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0405   //The following says we do not know what parameters are allowed so do no validation
0406   // Please change this to state exactly what you do use, even if it is no parameters
0407   edm::ParameterSetDescription desc;
0408 
0409   desc.add<double>("genJetPtMin", 20.);
0410   desc.add<double>("genJetAbsEtaMax", 2.4);
0411   desc.add<edm::InputTag>("genJets", edm::InputTag("ak4GenJets"));
0412   desc.add<edm::InputTag>("genBHadJetIndex", edm::InputTag("matchGenBHadron", "genBHadJetIndex"));
0413   desc.add<edm::InputTag>("genBHadFlavour", edm::InputTag("matchGenBHadron", "genBHadFlavour"));
0414   desc.add<edm::InputTag>("genBHadFromTopWeakDecay", edm::InputTag("matchGenBHadron", "genBHadFromTopWeakDecay"));
0415   desc.add<edm::InputTag>("genBHadPlusMothers", edm::InputTag("matchGenBHadron", "genBHadPlusMothers"));
0416   desc.add<edm::InputTag>("genBHadPlusMothersIndices", edm::InputTag("matchGenBHadron", "genBHadPlusMothersIndices"));
0417   desc.add<edm::InputTag>("genBHadIndex", edm::InputTag("matchGenBHadron", "genBHadIndex"));
0418   desc.add<edm::InputTag>("genBHadLeptonHadronIndex", edm::InputTag("matchGenBHadron", "genBHadLeptonHadronIndex"));
0419   desc.add<edm::InputTag>("genBHadLeptonViaTau", edm::InputTag("matchGenBHadron", "genBHadLeptonViaTau"));
0420 
0421   desc.add<edm::InputTag>("genCHadJetIndex", edm::InputTag("matchGenCHadron", "genCHadJetIndex"));
0422   desc.add<edm::InputTag>("genCHadFlavour", edm::InputTag("matchGenCHadron", "genCHadFlavour"));
0423   desc.add<edm::InputTag>("genCHadFromTopWeakDecay", edm::InputTag("matchGenCHadron", "genCHadFromTopWeakDecay"));
0424   desc.add<edm::InputTag>("genCHadBHadronId", edm::InputTag("matchGenCHadron", "genCHadBHadronId"));
0425 
0426   descriptions.add("categorizeGenTtbar", desc);
0427 }
0428 
0429 //define this as a plug-in
0430 DEFINE_FWK_MODULE(GenTtbarCategorizer);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!