Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:23:34 

0001 // system include files
0002 #include <memory>
0003 #include <map>
0004 #include <algorithm>
0005 
0006 // user include files
0007 #include "FWCore/Framework/interface/Frameworkfwd.h"
0008 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0009 #include "FWCore/Utilities/interface/InputTag.h"
0010 
0011 #include "FWCore/Framework/interface/Event.h"
0012 #include "FWCore/Framework/interface/MakerMacros.h"
0013 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0014 
0015 #include "DataFormats/JetReco/interface/Jet.h"
0016 #include "DataFormats/JetReco/interface/GenJetCollection.h"
0017 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0018 
0019 #include "FWCore/ServiceRegistry/interface/Service.h"
0020 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0021 #include <TTree.h>
0022 #include <Math/VectorUtil.h>
0023 
0024 class matchGenHFHadrons : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0025 public:
0026   explicit matchGenHFHadrons(const edm::ParameterSet&);
0027 
0028 private:
0029   void analyze(const edm::Event& event, const edm::EventSetup& setup) override;
0030   void beginJob() override;
0031   void bookHadronBranches(TTree* tree, const int flavour);
0032   void clearEventVariables();
0033   void clearHadronVariables();
0034   void fillHadronTree(TTree* tree,
0035                       const int flavour,
0036                       const reco::GenJetCollection& genJets,
0037                       const std::vector<int>& genHadIndex,
0038                       const std::vector<int>& genHadJetIndex,
0039                       const std::vector<int>& genHadFlavour,
0040                       const std::vector<int>& genHadFromTopWeakDecay,
0041                       const std::vector<reco::GenParticle>& genHadPlusMothers,
0042                       const std::vector<std::vector<int> >& genHadPlusMothersIndices,
0043                       const std::vector<int>& genHadLeptonHadronIndex,
0044                       const std::vector<int>& genHadLeptonViaTau,
0045                       const std::vector<int>& genHadBHadronId = std::vector<int>(0));
0046   // A recursive function that scans through the particle chain to find a particle with specific pdgId or abs(pdgId)
0047   void findAncestorParticles(const int startId,
0048                              std::set<int>& resultIds,
0049                              const std::vector<reco::GenParticle>& genParticles,
0050                              const std::vector<std::vector<int> >& motherIndices,
0051                              const int endPdgId,
0052                              const bool endPdgIdIsAbs = false,
0053                              const int firstAnyLast = 0);
0054 
0055   // Jets configuration
0056   const double genJetPtMin_;
0057   const double genJetAbsEtaMax_;
0058 
0059   // Input tags
0060   const edm::EDGetTokenT<reco::GenJetCollection> genJetsToken_;
0061 
0062   const edm::EDGetTokenT<std::vector<int> > genBHadJetIndexToken_;
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   const edm::EDGetTokenT<std::vector<int> > genBHadLeptonHadronIndexToken_;
0069   const edm::EDGetTokenT<std::vector<int> > genBHadLeptonViaTauToken_;
0070 
0071   const edm::EDGetTokenT<std::vector<int> > genCHadJetIndexToken_;
0072   const edm::EDGetTokenT<std::vector<int> > genCHadFlavourToken_;
0073   const edm::EDGetTokenT<std::vector<int> > genCHadFromTopWeakDecayToken_;
0074   const edm::EDGetTokenT<std::vector<int> > genCHadBHadronIdToken_;
0075   const edm::EDGetTokenT<std::vector<reco::GenParticle> > genCHadPlusMothersToken_;
0076   const edm::EDGetTokenT<std::vector<std::vector<int> > > genCHadPlusMothersIndicesToken_;
0077   const edm::EDGetTokenT<std::vector<int> > genCHadIndexToken_;
0078   const edm::EDGetTokenT<std::vector<int> > genCHadLeptonHadronIndexToken_;
0079   const edm::EDGetTokenT<std::vector<int> > genCHadLeptonViaTauToken_;
0080 
0081   // Output to be written to the trees
0082   TTree* tree_events_;
0083   TTree* tree_bHadrons_;
0084   TTree* tree_cHadrons_;
0085 
0086   int nJets_;
0087 
0088   int nBJets_;
0089   int nBHadrons_;
0090   int nBHadronsTop_;
0091   int nBHadronsAdditional_;
0092   int nBHadronsPseudoadditional_;
0093 
0094   int nCJets_;
0095   int nCHadrons_;
0096   int nCHadronsAdditional_;
0097   int nCHadronsPseudoadditional_;
0098 
0099   int nBHadronsHiggs_;
0100   int additionalJetEventId_;
0101 
0102   // Additional performance information
0103   int hadron_flavour_;
0104   int hadron_nQuarksAll_;
0105   int hadron_nQuarksSameFlavour_;
0106   int hadron_fromTopWeakDecay_;
0107   int hadron_bHadronId_;
0108 
0109   double hadron_quark1_dR_;
0110   double hadron_quark2_dR_;
0111   double hadron_quark1_ptRatio_;
0112   double hadron_quark2_ptRatio_;
0113 
0114   int hadron_nLeptonsAll_;
0115   int hadron_nLeptonsViaTau_;
0116 
0117   double hadron_jetPtRatio_;
0118   double hadron_jetDR_;
0119 };
0120 
0121 bool sort_pairs_second(const std::pair<int, double>& a, const std::pair<int, double>& b) { return a.second < b.second; }
0122 
0123 matchGenHFHadrons::matchGenHFHadrons(const edm::ParameterSet& config)
0124     : genJetPtMin_(config.getParameter<double>("genJetPtMin")),
0125       genJetAbsEtaMax_(config.getParameter<double>("genJetAbsEtaMax")),
0126       genJetsToken_(consumes<reco::GenJetCollection>(config.getParameter<edm::InputTag>("genJets"))),
0127       genBHadJetIndexToken_(consumes<std::vector<int> >(config.getParameter<edm::InputTag>("genBHadJetIndex"))),
0128       genBHadFlavourToken_(consumes<std::vector<int> >(config.getParameter<edm::InputTag>("genBHadFlavour"))),
0129       genBHadFromTopWeakDecayToken_(
0130           consumes<std::vector<int> >(config.getParameter<edm::InputTag>("genBHadFromTopWeakDecay"))),
0131       genBHadPlusMothersToken_(
0132           consumes<std::vector<reco::GenParticle> >(config.getParameter<edm::InputTag>("genBHadPlusMothers"))),
0133       genBHadPlusMothersIndicesToken_(
0134           consumes<std::vector<std::vector<int> > >(config.getParameter<edm::InputTag>("genBHadPlusMothersIndices"))),
0135       genBHadIndexToken_(consumes<std::vector<int> >(config.getParameter<edm::InputTag>("genBHadIndex"))),
0136       genBHadLeptonHadronIndexToken_(
0137           consumes<std::vector<int> >(config.getParameter<edm::InputTag>("genBHadLeptonHadronIndex"))),
0138       genBHadLeptonViaTauToken_(consumes<std::vector<int> >(config.getParameter<edm::InputTag>("genBHadLeptonViaTau"))),
0139       genCHadJetIndexToken_(consumes<std::vector<int> >(config.getParameter<edm::InputTag>("genCHadJetIndex"))),
0140       genCHadFlavourToken_(consumes<std::vector<int> >(config.getParameter<edm::InputTag>("genCHadFlavour"))),
0141       genCHadFromTopWeakDecayToken_(
0142           consumes<std::vector<int> >(config.getParameter<edm::InputTag>("genCHadFromTopWeakDecay"))),
0143       genCHadBHadronIdToken_(consumes<std::vector<int> >(config.getParameter<edm::InputTag>("genCHadBHadronId"))),
0144       genCHadPlusMothersToken_(
0145           consumes<std::vector<reco::GenParticle> >(config.getParameter<edm::InputTag>("genCHadPlusMothers"))),
0146       genCHadPlusMothersIndicesToken_(
0147           consumes<std::vector<std::vector<int> > >(config.getParameter<edm::InputTag>("genCHadPlusMothersIndices"))),
0148       genCHadIndexToken_(consumes<std::vector<int> >(config.getParameter<edm::InputTag>("genCHadIndex"))),
0149       genCHadLeptonHadronIndexToken_(
0150           consumes<std::vector<int> >(config.getParameter<edm::InputTag>("genCHadLeptonHadronIndex"))),
0151       genCHadLeptonViaTauToken_(
0152           consumes<std::vector<int> >(config.getParameter<edm::InputTag>("genCHadLeptonViaTau"))) {
0153   usesResource(TFileService::kSharedResource);
0154 }
0155 
0156 void matchGenHFHadrons::analyze(const edm::Event& event, const edm::EventSetup& setup) {
0157   this->clearEventVariables();
0158 
0159   // Reading gen jets from the event
0160   edm::Handle<reco::GenJetCollection> genJets;
0161   event.getByToken(genJetsToken_, genJets);
0162 
0163   // Reading B hadrons related information
0164   edm::Handle<std::vector<int> > genBHadFlavour;
0165   event.getByToken(genBHadFlavourToken_, genBHadFlavour);
0166 
0167   edm::Handle<std::vector<int> > genBHadJetIndex;
0168   event.getByToken(genBHadJetIndexToken_, genBHadJetIndex);
0169 
0170   edm::Handle<std::vector<int> > genBHadFromTopWeakDecay;
0171   event.getByToken(genBHadFromTopWeakDecayToken_, genBHadFromTopWeakDecay);
0172 
0173   edm::Handle<std::vector<reco::GenParticle> > genBHadPlusMothers;
0174   event.getByToken(genBHadPlusMothersToken_, genBHadPlusMothers);
0175 
0176   edm::Handle<std::vector<std::vector<int> > > genBHadPlusMothersIndices;
0177   event.getByToken(genBHadPlusMothersIndicesToken_, genBHadPlusMothersIndices);
0178 
0179   edm::Handle<std::vector<int> > genBHadIndex;
0180   event.getByToken(genBHadIndexToken_, genBHadIndex);
0181 
0182   edm::Handle<std::vector<int> > genBHadLeptonHadronIndex;
0183   event.getByToken(genBHadLeptonHadronIndexToken_, genBHadLeptonHadronIndex);
0184 
0185   edm::Handle<std::vector<int> > genBHadLeptonViaTau;
0186   event.getByToken(genBHadLeptonViaTauToken_, genBHadLeptonViaTau);
0187 
0188   // Reading C hadrons related information
0189   edm::Handle<std::vector<int> > genCHadFlavour;
0190   event.getByToken(genCHadFlavourToken_, genCHadFlavour);
0191 
0192   edm::Handle<std::vector<int> > genCHadJetIndex;
0193   event.getByToken(genCHadJetIndexToken_, genCHadJetIndex);
0194 
0195   edm::Handle<std::vector<int> > genCHadFromTopWeakDecay;
0196   event.getByToken(genCHadFromTopWeakDecayToken_, genCHadFromTopWeakDecay);
0197 
0198   edm::Handle<std::vector<int> > genCHadBHadronId;
0199   event.getByToken(genCHadBHadronIdToken_, genCHadBHadronId);
0200 
0201   edm::Handle<std::vector<reco::GenParticle> > genCHadPlusMothers;
0202   event.getByToken(genCHadPlusMothersToken_, genCHadPlusMothers);
0203 
0204   edm::Handle<std::vector<std::vector<int> > > genCHadPlusMothersIndices;
0205   event.getByToken(genCHadPlusMothersIndicesToken_, genCHadPlusMothersIndices);
0206 
0207   edm::Handle<std::vector<int> > genCHadIndex;
0208   event.getByToken(genCHadIndexToken_, genCHadIndex);
0209 
0210   edm::Handle<std::vector<int> > genCHadLeptonHadronIndex;
0211   event.getByToken(genCHadLeptonHadronIndexToken_, genCHadLeptonHadronIndex);
0212 
0213   edm::Handle<std::vector<int> > genCHadLeptonViaTau;
0214   event.getByToken(genCHadLeptonViaTauToken_, genCHadLeptonViaTau);
0215 
0216   // Counting number of jets of different flavours in the event
0217   for (size_t jetId = 0; jetId < genJets->size(); ++jetId) {
0218     if (genJets->at(jetId).pt() < genJetPtMin_)
0219       continue;
0220     if (std::fabs(genJets->at(jetId).eta()) < genJetAbsEtaMax_)
0221       continue;
0222     nJets_++;
0223     if (std::find(genBHadJetIndex->begin(), genBHadJetIndex->end(), jetId) != genBHadJetIndex->end())
0224       nBJets_++;
0225     else if (std::find(genCHadJetIndex->begin(), genCHadJetIndex->end(), jetId) != genCHadJetIndex->end())
0226       nCJets_++;
0227   }
0228   // Counting number of different b hadrons
0229   for (size_t hadronId = 0; hadronId < genBHadFlavour->size(); ++hadronId) {
0230     const int flavour = genBHadFlavour->at(hadronId);
0231     const int flavourAbs = std::abs(flavour);
0232     // 0 - not from top decay; 1 - from top decay; 2 - not identified;
0233     const bool afterTopDecay = genBHadFromTopWeakDecay->at(hadronId) > 0 ? true : false;
0234     nBHadrons_++;
0235     if (flavourAbs == 25)
0236       nBHadronsHiggs_++;
0237     if (flavourAbs == 6)
0238       nBHadronsTop_++;
0239     if (flavourAbs != 6) {
0240       if (afterTopDecay)
0241         nBHadronsPseudoadditional_++;
0242       else
0243         nBHadronsAdditional_++;
0244     }
0245   }
0246   // Counting number of different c hadrons
0247   for (size_t hadronId = 0; hadronId < genCHadJetIndex->size(); ++hadronId) {
0248     // 0 - not from top decay; 1 - from top decay; 2 - not identified;
0249     const bool afterTopDecay = genCHadFromTopWeakDecay->at(hadronId) > 0 ? true : false;
0250     nCHadrons_++;
0251     // Skipping c hadrons that are coming from b hadrons
0252     if (genCHadBHadronId->at(hadronId) >= 0)
0253       continue;
0254 
0255     if (afterTopDecay)
0256       nCHadronsPseudoadditional_++;
0257     else
0258       nCHadronsAdditional_++;
0259   }
0260 
0261   // Map <jet index, number of specific hadrons in the jet>
0262   // B jets with b hadrons directly from top quark decay (including those not in acceptance)
0263   std::map<int, int> bJetFromTopIds_all;
0264   // B jets with b hadrons directly from top quark decay
0265   std::map<int, int> bJetFromTopIds;
0266   // B jets with b hadrons after top quark decay
0267   std::map<int, int> bJetAfterTopIds;
0268   // B jets with b hadrons before top quark decay chain
0269   std::map<int, int> bJetBeforeTopIds;
0270   // C jets with c hadrons before top quark decay chain
0271   std::map<int, int> cJetBeforeTopIds;
0272   // C jets with c hadrons after top quark decay
0273   std::map<int, int> cJetAfterTopIds;
0274 
0275   // Counting number of specific hadrons in each b jet
0276   for (size_t hadronId = 0; hadronId < genBHadIndex->size(); ++hadronId) {
0277     // Flavour of the hadron's origin
0278     const int flavour = genBHadFlavour->at(hadronId);
0279     // Skipping b hadrons from W-boson decays (in semileptonic or full-hadronic channels)
0280     if (std::abs(flavour) == 24)
0281       continue;
0282     // Whether hadron radiated before top quark decay
0283     const bool afterTopDecay = genBHadFromTopWeakDecay->at(hadronId) > 0 ? true : false;
0284     // Index of a jet associated to the hadron
0285     const int jetIndex = genBHadJetIndex->at(hadronId);
0286     // Skipping hadrons which have no associated jet
0287     if (jetIndex < 0)
0288       continue;
0289     // Jet from direct top quark decay [pdgId(top)=6]
0290     if (std::abs(flavour) == 6) {
0291       if (bJetFromTopIds_all.count(jetIndex) < 1)
0292         bJetFromTopIds_all[jetIndex] = 1;
0293       else
0294         bJetFromTopIds_all[jetIndex]++;
0295     }
0296     // Skipping if jet is not in acceptance
0297     if (genJets->at(jetIndex).pt() < genJetPtMin_)
0298       continue;
0299     if (std::fabs(genJets->at(jetIndex).eta()) > genJetAbsEtaMax_)
0300       continue;
0301     // Identifying jets with b hadrons not from top quark decay
0302     // Jet from direct top quark decay [pdgId(top)=6]
0303     if (std::abs(flavour) == 6) {
0304       if (bJetFromTopIds.count(jetIndex) < 1)
0305         bJetFromTopIds[jetIndex] = 1;
0306       else
0307         bJetFromTopIds[jetIndex]++;
0308     }
0309     // Skipping if jet is from top quark decay
0310     if (std::abs(flavour) == 6)
0311       continue;
0312     // Jet before top quark decay
0313     if (!afterTopDecay) {
0314       if (bJetBeforeTopIds.count(jetIndex) < 1)
0315         bJetBeforeTopIds[jetIndex] = 1;
0316       else
0317         bJetBeforeTopIds[jetIndex]++;
0318     }
0319     // Jet after top quark decay but not directly from top
0320     else if (afterTopDecay) {
0321       if (bJetAfterTopIds.count(jetIndex) < 1)
0322         bJetAfterTopIds[jetIndex] = 1;
0323       else
0324         bJetAfterTopIds[jetIndex]++;
0325     }
0326   }
0327 
0328   // Counting number of specific hadrons in each c jet
0329   for (size_t hadronId = 0; hadronId < genCHadJetIndex->size(); ++hadronId) {
0330     // Skipping c hadrons that are coming from b hadrons
0331     if (genCHadBHadronId->at(hadronId) >= 0)
0332       continue;
0333     // Flavour of the hadron's origin
0334     const int flavour = genCHadFlavour->at(hadronId);
0335     // Skipping c hadrons from W-boson decays (in semileptonic or full-hadronic channels)
0336     if (std::abs(flavour) == 24)
0337       continue;
0338     // Index of a jet associated to the hadron
0339     const int jetIndex = genCHadJetIndex->at(hadronId);
0340     // Whether hadron radiated after top quark decay
0341     const bool afterTopDecay = genCHadFromTopWeakDecay->at(hadronId) > 0 ? true : false;
0342     // Skipping hadrons which have no associated jet
0343     if (jetIndex < 0)
0344       continue;
0345     // Skipping if jet is not in acceptance
0346     if (genJets->at(jetIndex).pt() < genJetPtMin_)
0347       continue;
0348     if (std::fabs(genJets->at(jetIndex).eta()) > genJetAbsEtaMax_)
0349       continue;
0350     // Jet before top quark decay
0351     if (!afterTopDecay) {
0352       if (cJetBeforeTopIds.count(jetIndex) < 1)
0353         cJetBeforeTopIds[jetIndex] = 1;
0354       else
0355         cJetBeforeTopIds[jetIndex]++;
0356     }
0357     // Jet after top quark decay but not directly from top
0358     else if (afterTopDecay) {
0359       if (cJetAfterTopIds.count(jetIndex) < 1)
0360         cJetAfterTopIds[jetIndex] = 1;
0361       else
0362         cJetAfterTopIds[jetIndex]++;
0363     }
0364   }
0365 
0366   // Finding additional b jets (before top decay)
0367   std::vector<int> additionalBJetIds;
0368   for (std::map<int, int>::iterator it = bJetBeforeTopIds.begin(); it != bJetBeforeTopIds.end(); ++it) {
0369     const int jetId = it->first;
0370     // Skipping the jet if it contains a b hadron directly from top quark decay
0371     if (bJetFromTopIds.count(jetId) > 0)
0372       continue;
0373     additionalBJetIds.push_back(jetId);
0374   }
0375   // Finding pseudo-additional b jets (after top decay)
0376   std::vector<int> pseudoadditionalBJetIds;
0377   for (std::map<int, int>::iterator it = bJetAfterTopIds.begin(); it != bJetAfterTopIds.end(); ++it) {
0378     const int jetId = it->first;
0379     // Skipping the jet if it contains a b hadron directly from top quark decay
0380     if (bJetFromTopIds.count(jetId) > 0)
0381       continue;
0382     pseudoadditionalBJetIds.push_back(jetId);
0383   }
0384   // Finding additional c jets
0385   std::vector<int> additionalCJetIds;
0386   for (std::map<int, int>::iterator it = cJetBeforeTopIds.begin(); it != cJetBeforeTopIds.end(); ++it) {
0387     const int jetId = it->first;
0388     // Skipping the jet if it contains a b hadron directly from top quark decay
0389     if (bJetFromTopIds.count(jetId) > 0)
0390       continue;
0391     additionalCJetIds.push_back(jetId);
0392   }
0393   // Finding pseudo-additional c jets (after top decay)
0394   std::vector<int> pseudoadditionalCJetIds;
0395   for (std::map<int, int>::iterator it = cJetAfterTopIds.begin(); it != cJetAfterTopIds.end(); ++it) {
0396     const int jetId = it->first;
0397     // Skipping the jet if it contains a b hadron directly from top quark decay
0398     if (bJetFromTopIds.count(jetId) > 0)
0399       continue;
0400     pseudoadditionalCJetIds.push_back(jetId);
0401   }
0402 
0403   // Categorizing event based on number of additional b/c jets and number of corresponding hadrons in each of them
0404   // Exclusively for dileptonic final states of the ttbar system (may be unapplicable for other decay channels)
0405   additionalJetEventId_ = bJetFromTopIds.size() * 100;
0406   // tt + 1 additional b jet
0407   if (additionalBJetIds.size() == 1) {
0408     int nHadronsInJet = bJetBeforeTopIds[additionalBJetIds.at(0)];
0409     // tt + 1 additional b jet from 1 additional b hadron
0410     if (nHadronsInJet == 1)
0411       additionalJetEventId_ += 51;
0412     // tt + 1 additional b jet from >=2 additional b hadrons
0413     else
0414       additionalJetEventId_ += 52;
0415   }
0416   // tt + 2 additional b jets
0417   else if (additionalBJetIds.size() > 1) {
0418     int nHadronsInJet1 = bJetBeforeTopIds[additionalBJetIds.at(0)];
0419     int nHadronsInJet2 = bJetBeforeTopIds[additionalBJetIds.at(1)];
0420     // tt + 2 additional b jets each from 1 additional b hadron
0421     if (std::max(nHadronsInJet1, nHadronsInJet2) == 1)
0422       additionalJetEventId_ += 53;
0423     // tt + 2 additional b jets one of which from >=2 overlapping additional b hadrons
0424     else if (std::min(nHadronsInJet1, nHadronsInJet2) == 1 && std::max(nHadronsInJet1, nHadronsInJet2) > 1)
0425       additionalJetEventId_ += 54;
0426     // tt + 2 additional b jets each from >=2 additional b hadrons
0427     else if (std::min(nHadronsInJet1, nHadronsInJet2) > 1)
0428       additionalJetEventId_ += 55;
0429   }
0430   // tt + no additional b jets
0431   else if (additionalBJetIds.size() == 0) {
0432     // tt + >=1 pseudo-additional b jet with b hadrons after top quark decay
0433     if (pseudoadditionalBJetIds.size() > 0)
0434       additionalJetEventId_ += 56;
0435     // tt + 1 additional c jet
0436     else if (additionalCJetIds.size() == 1) {
0437       int nHadronsInJet = cJetBeforeTopIds[additionalCJetIds.at(0)];
0438       // tt + 1 additional c jet from 1 additional c hadron
0439       if (nHadronsInJet == 1)
0440         additionalJetEventId_ += 41;
0441       // tt + 1 additional c jet from >=2 overlapping additional c hadrons
0442       else
0443         additionalJetEventId_ += 42;
0444     }
0445     // tt + >=2 additional c jets
0446     else if (additionalCJetIds.size() > 1) {
0447       int nHadronsInJet1 = cJetBeforeTopIds[additionalCJetIds.at(0)];
0448       int nHadronsInJet2 = cJetBeforeTopIds[additionalCJetIds.at(1)];
0449       // tt + 2 additional c jets each from 1 additional c hadron
0450       if (std::max(nHadronsInJet1, nHadronsInJet2) == 1)
0451         additionalJetEventId_ += 43;
0452       // tt + 2 additional c jets one of which from >=2 overlapping additional c hadrons
0453       else if (std::min(nHadronsInJet1, nHadronsInJet2) == 1 && std::max(nHadronsInJet1, nHadronsInJet2) > 1)
0454         additionalJetEventId_ += 44;
0455       // tt + 2 additional c jets each from >=2 additional c hadrons
0456       else if (std::min(nHadronsInJet1, nHadronsInJet2) > 1)
0457         additionalJetEventId_ += 45;
0458     }
0459     // tt + no additional c jets
0460     else if (additionalCJetIds.size() == 0) {
0461       // tt + >=1 pseudo-additional c jet with c hadrons after top quark decay
0462       if (pseudoadditionalCJetIds.size() > 0)
0463         additionalJetEventId_ += 46;
0464       // tt + light jets
0465       else
0466         additionalJetEventId_ += 0;
0467     }
0468   }
0469 
0470   // Performance information filled in a tree with entry for each b hadron #######################################
0471   fillHadronTree(tree_bHadrons_,
0472                  5,
0473                  *genJets,
0474                  *genBHadIndex,
0475                  *genBHadJetIndex,
0476                  *genBHadFlavour,
0477                  *genBHadFromTopWeakDecay,
0478                  *genBHadPlusMothers,
0479                  *genBHadPlusMothersIndices,
0480                  *genBHadLeptonHadronIndex,
0481                  *genBHadLeptonViaTau);
0482 
0483   // Performance information filled in a tree with entry for each c hadron #######################################
0484   fillHadronTree(tree_cHadrons_,
0485                  4,
0486                  *genJets,
0487                  *genCHadIndex,
0488                  *genCHadJetIndex,
0489                  *genCHadFlavour,
0490                  *genCHadFromTopWeakDecay,
0491                  *genCHadPlusMothers,
0492                  *genCHadPlusMothersIndices,
0493                  *genCHadLeptonHadronIndex,
0494                  *genCHadLeptonViaTau,
0495                  *genCHadBHadronId);
0496 
0497   tree_events_->Fill();
0498 }
0499 
0500 void matchGenHFHadrons::findAncestorParticles(const int startId,
0501                                               std::set<int>& resultIds,
0502                                               const std::vector<reco::GenParticle>& genParticles,
0503                                               const std::vector<std::vector<int> >& motherIndices,
0504                                               const int endPdgId,
0505                                               const bool endPdgIdIsAbs,
0506                                               const int firstAnyLast) {
0507   if (startId < 0)
0508     return;
0509   int startPdgId = genParticles.at(startId).pdgId();
0510   if (endPdgIdIsAbs)
0511     startPdgId = std::abs(startPdgId);
0512 
0513   // No mothers stored for the requested particle
0514   if ((int)motherIndices.size() < startId + 1)
0515     return;
0516 
0517   // Getting ids of mothers of the starting particle
0518   const std::vector<int>& motherIds = motherIndices.at(startId);
0519 
0520   // In case the starting particle has pdgId as requested and it should not have same mothers
0521   if (startPdgId == endPdgId && firstAnyLast == 1) {
0522     int nSamePdgIdMothers = 0;
0523     // Counting how many mothers with the same pdgId the particle has
0524     for (size_t motherId = 0; motherId < motherIds.size(); ++motherId) {
0525       const int motherParticleId = motherIds.at(motherId);
0526       if (motherParticleId < 0)
0527         continue;
0528       int motherParticlePdgId = genParticles.at(motherParticleId).pdgId();
0529       if (endPdgIdIsAbs)
0530         motherParticlePdgId = std::abs(motherParticlePdgId);
0531       if (motherParticlePdgId == startPdgId)
0532         nSamePdgIdMothers++;
0533     }
0534     if (nSamePdgIdMothers < 1) {
0535       resultIds.insert(startId);
0536       return;
0537     }
0538   }
0539 
0540   // Looping over all mothers of the particle
0541   for (size_t motherId = 0; motherId < motherIds.size(); ++motherId) {
0542     const int motherParticleId = motherIds.at(motherId);
0543     if (motherParticleId < 0)
0544       continue;
0545     int motherParticlePdgId = genParticles.at(motherParticleId).pdgId();
0546     if (endPdgIdIsAbs)
0547       motherParticlePdgId = std::abs(motherParticlePdgId);
0548     // Checking whether pdgId of this particle matches the requested one
0549     bool isCandidate = false;
0550     bool inProperPlace = false;
0551     if (motherParticlePdgId == endPdgId)
0552       isCandidate = true;
0553     // Checking whether the candidate particle is first/last according to the request
0554     if (isCandidate) {
0555       if (firstAnyLast == -1 && startPdgId != motherParticlePdgId)
0556         inProperPlace = true;
0557       else if (firstAnyLast == 0)
0558         inProperPlace = true;
0559     }
0560     // Adding the mother if it has proper pdgId and place in chain
0561     if (isCandidate && inProperPlace) {
0562       resultIds.insert(motherParticleId);
0563       if (firstAnyLast < 0)
0564         continue;
0565     }
0566 
0567     // Checking mothers of this mother
0568     findAncestorParticles(
0569         motherParticleId, resultIds, genParticles, motherIndices, endPdgId, endPdgIdIsAbs, firstAnyLast);
0570   }
0571 }
0572 
0573 void matchGenHFHadrons::beginJob() {
0574   edm::Service<TFileService> fileService;
0575   if (!fileService)
0576     throw edm::Exception(edm::errors::Configuration, "TFileService is not registered in cfg file");
0577   tree_events_ = fileService->make<TTree>("tree_events", "tree_events");
0578   tree_bHadrons_ = fileService->make<TTree>("tree_bHadrons", "tree_bHadrons");
0579   tree_cHadrons_ = fileService->make<TTree>("tree_cHadrons", "tree_cHadrons");
0580 
0581   // Creating output branches
0582   tree_events_->Branch("nJets", &nJets_);
0583 
0584   tree_events_->Branch("nBJets", &nBJets_);
0585   tree_events_->Branch("nBHadrons", &nBHadrons_);
0586   tree_events_->Branch("nBHadronsTop", &nBHadronsTop_);
0587   tree_events_->Branch("nBHadronsAdditional", &nBHadronsAdditional_);
0588   tree_events_->Branch("nBHadronsPseudoadditional", &nBHadronsPseudoadditional_);
0589 
0590   tree_events_->Branch("nCJets", &nCJets_);
0591   tree_events_->Branch("nCHadrons", &nCHadrons_);
0592   tree_events_->Branch("nCHadronsAdditional", &nCHadronsAdditional_);
0593   tree_events_->Branch("nCHadronsPseudoadditional", &nCHadronsPseudoadditional_);
0594 
0595   tree_events_->Branch("nBHadronsHiggs", &nBHadronsHiggs_);
0596   tree_events_->Branch("additionalJetEventId", &additionalJetEventId_);
0597 
0598   bookHadronBranches(tree_bHadrons_, 5);
0599   bookHadronBranches(tree_cHadrons_, 4);
0600 }
0601 
0602 void matchGenHFHadrons::bookHadronBranches(TTree* tree, const int flavour) {
0603   tree->Branch("flavour", &hadron_flavour_);
0604   tree->Branch("nQuarksAll", &hadron_nQuarksAll_);
0605   tree->Branch("nQuarksSameFlavour", &hadron_nQuarksSameFlavour_);
0606   tree->Branch("fromTopWeakDecay", &hadron_fromTopWeakDecay_);
0607 
0608   tree->Branch("quark1_dR", &hadron_quark1_dR_);
0609   tree->Branch("quark2_dR", &hadron_quark2_dR_);
0610   tree->Branch("quark1_ptRatio", &hadron_quark1_ptRatio_);
0611   tree->Branch("quark2_ptRatio", &hadron_quark2_ptRatio_);
0612 
0613   tree->Branch("nLeptonsAll", &hadron_nLeptonsAll_);
0614   tree->Branch("nLeptonsViaTau", &hadron_nLeptonsViaTau_);
0615 
0616   tree->Branch("jetPtRatio", &hadron_jetPtRatio_);
0617   tree->Branch("jetDR", &hadron_jetDR_);
0618 
0619   if (flavour == 4)
0620     tree->Branch("bHadronId", &hadron_bHadronId_);
0621 }
0622 
0623 void matchGenHFHadrons::clearEventVariables() {
0624   nJets_ = 0;
0625 
0626   nBJets_ = 0;
0627   nBHadrons_ = 0;
0628   nBHadronsTop_ = 0;
0629   nBHadronsAdditional_ = 0;
0630   nBHadronsPseudoadditional_ = 0;
0631 
0632   nCJets_ = 0;
0633   nCHadrons_ = 0;
0634   nCHadronsAdditional_ = 0;
0635   nCHadronsPseudoadditional_ = 0;
0636 
0637   nBHadronsHiggs_ = 0;
0638   additionalJetEventId_ += -1;
0639 }
0640 
0641 void matchGenHFHadrons::clearHadronVariables() {
0642   hadron_flavour_ = 0;
0643   hadron_nQuarksAll_ = 0;
0644   hadron_nQuarksSameFlavour_ = 0;
0645   hadron_fromTopWeakDecay_ = -1;
0646   hadron_bHadronId_ = -1;
0647 
0648   hadron_quark1_dR_ = -0.1;
0649   hadron_quark2_dR_ = -0.1;
0650   hadron_quark1_ptRatio_ = -0.1;
0651   hadron_quark2_ptRatio_ = -0.1;
0652 
0653   hadron_nLeptonsAll_ = 0;
0654   hadron_nLeptonsViaTau_ = 0;
0655 
0656   hadron_jetPtRatio_ = -0.1;
0657   hadron_jetDR_ = -0.1;
0658 }
0659 
0660 void matchGenHFHadrons::fillHadronTree(TTree* tree,
0661                                        const int flavour,
0662                                        const reco::GenJetCollection& genJets,
0663                                        const std::vector<int>& genHadIndex,
0664                                        const std::vector<int>& genHadJetIndex,
0665                                        const std::vector<int>& genHadFlavour,
0666                                        const std::vector<int>& genHadFromTopWeakDecay,
0667                                        const std::vector<reco::GenParticle>& genHadPlusMothers,
0668                                        const std::vector<std::vector<int> >& genHadPlusMothersIndices,
0669                                        const std::vector<int>& genHadLeptonHadronIndex,
0670                                        const std::vector<int>& genHadLeptonViaTau,
0671                                        const std::vector<int>& genHadBHadronId) {
0672   for (size_t hadronId = 0; hadronId < genHadIndex.size(); ++hadronId) {
0673     clearHadronVariables();
0674 
0675     const int hadronParticleId = genHadIndex.at(hadronId);
0676     if (hadronParticleId < 0)
0677       continue;
0678 
0679     // Calculating hadron/jet pt ratio
0680     const int hadronJetId = genHadJetIndex.at(hadronId);
0681     if (hadronJetId >= 0) {
0682       hadron_jetPtRatio_ = genHadPlusMothers.at(genHadIndex.at(hadronId)).pt() / genJets.at(hadronJetId).pt();
0683       hadron_jetDR_ = ROOT::Math::VectorUtil::DeltaR(genHadPlusMothers.at(genHadIndex.at(hadronId)).polarP4(),
0684                                                      genJets.at(hadronJetId).polarP4());
0685     }
0686 
0687     if (genHadBHadronId.size() > 0)
0688       hadron_bHadronId_ = genHadBHadronId.at(hadronId);
0689     hadron_fromTopWeakDecay_ = genHadFromTopWeakDecay.at(hadronId);
0690 
0691     const int hadronParticlePdgId = genHadPlusMothers.at(hadronParticleId).pdgId();
0692     hadron_flavour_ = genHadFlavour.at(hadronId);
0693     // Finding all b quark candidates
0694     std::set<int> hadronLastQuarks;
0695     findAncestorParticles(
0696         hadronParticleId, hadronLastQuarks, genHadPlusMothers, genHadPlusMothersIndices, flavour, true, 1);
0697     hadron_nQuarksAll_ = hadronLastQuarks.size();
0698     // Identifying flavour of the proper candidate quark
0699     int quarkCandidateFlavourSign = hadronParticlePdgId / std::abs(hadronParticlePdgId);
0700     // Inverting flavour if this is for a b meson
0701     if (std::abs(hadronParticlePdgId) / 100 == 5 && std::abs(hadronParticlePdgId) % 1000 / 10 / 11 != 5)
0702       quarkCandidateFlavourSign *= -1;
0703     // Selecting candidates that have proper flavour sign
0704     std::vector<std::pair<int, double> > hadronLastQuarks_sameFlavour;
0705     for (std::set<int>::iterator it = hadronLastQuarks.begin(); it != hadronLastQuarks.end(); ++it) {
0706       const int quarkParticleId = *it;
0707       if (quarkParticleId < 0)
0708         continue;
0709       const int quarkParticlePdgId = genHadPlusMothers.at(quarkParticleId).pdgId();
0710       // Skipping quarks that have opposite flavour sign
0711       if (quarkParticlePdgId * quarkCandidateFlavourSign < 0)
0712         continue;
0713       double hadron_quark_dR = ROOT::Math::VectorUtil::DeltaR(genHadPlusMothers.at(hadronParticleId).polarP4(),
0714                                                               genHadPlusMothers.at(quarkParticleId).polarP4());
0715       hadronLastQuarks_sameFlavour.push_back(std::pair<int, double>(quarkParticleId, hadron_quark_dR));
0716     }
0717     hadron_nQuarksSameFlavour_ = hadronLastQuarks_sameFlavour.size();
0718 
0719     // Sorting quarks with proper flavour by their dR
0720     std::sort(hadronLastQuarks_sameFlavour.begin(), hadronLastQuarks_sameFlavour.end(), sort_pairs_second);
0721     const int hadronQuark1ParticleId =
0722         hadronLastQuarks_sameFlavour.size() > 0 ? hadronLastQuarks_sameFlavour.at(0).first : -1;
0723     const int hadronQuark2ParticleId =
0724         hadronLastQuarks_sameFlavour.size() > 1 ? hadronLastQuarks_sameFlavour.at(1).first : -1;
0725 
0726     if (hadronQuark1ParticleId >= 0) {
0727       hadron_quark1_dR_ = ROOT::Math::VectorUtil::DeltaR(genHadPlusMothers.at(hadronParticleId).polarP4(),
0728                                                          genHadPlusMothers.at(hadronQuark1ParticleId).polarP4());
0729       hadron_quark1_ptRatio_ =
0730           genHadPlusMothers.at(hadronParticleId).pt() / genHadPlusMothers.at(hadronQuark1ParticleId).pt();
0731     }
0732 
0733     if (hadronQuark2ParticleId >= 0) {
0734       hadron_quark2_dR_ = ROOT::Math::VectorUtil::DeltaR(genHadPlusMothers.at(hadronParticleId).polarP4(),
0735                                                          genHadPlusMothers.at(hadronQuark2ParticleId).polarP4());
0736       hadron_quark2_ptRatio_ =
0737           genHadPlusMothers.at(hadronParticleId).pt() / genHadPlusMothers.at(hadronQuark2ParticleId).pt();
0738     }
0739 
0740     // Counting number of leptons coming from the b hadron decay
0741     for (size_t leptonId = 0; leptonId < genHadLeptonHadronIndex.size(); ++leptonId) {
0742       const size_t leptonHadronId = genHadLeptonHadronIndex.at(leptonId);
0743       // Skipping the lepton id doesn't come from the current b hadron
0744       if (leptonHadronId != hadronId)
0745         continue;
0746       hadron_nLeptonsAll_++;
0747       if (genHadLeptonViaTau.at(leptonId) == 1)
0748         hadron_nLeptonsViaTau_++;
0749     }
0750 
0751     tree->Fill();
0752   }
0753 }
0754 
0755 DEFINE_FWK_MODULE(matchGenHFHadrons);
[ 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!