From 38d9cde14c2f4d13018e338dd89c704b07d5ad8c Mon Sep 17 00:00:00 2001 From: Nazar Bartosik Date: Tue, 4 Aug 2015 15:36:04 +0200 Subject: [PATCH 1/3] Added plugin for generated tt+X event categorisation --- .../TopTools/plugins/GenTtbarCategorizer.cc | 434 ++++++++++++++++++ .../python/GenTtbarCategorizer_cfi.py | 24 + TopQuarkAnalysis/TopTools/test/BuildFile.xml | 8 + .../TopTools/test/TestGenTtbarCategories.cc | 185 ++++++++ .../TopTools/test/testGenTtbarCategories.py | 142 ++++++ 5 files changed, 793 insertions(+) create mode 100644 TopQuarkAnalysis/TopTools/plugins/GenTtbarCategorizer.cc create mode 100644 TopQuarkAnalysis/TopTools/python/GenTtbarCategorizer_cfi.py create mode 100644 TopQuarkAnalysis/TopTools/test/TestGenTtbarCategories.cc create mode 100644 TopQuarkAnalysis/TopTools/test/testGenTtbarCategories.py diff --git a/TopQuarkAnalysis/TopTools/plugins/GenTtbarCategorizer.cc b/TopQuarkAnalysis/TopTools/plugins/GenTtbarCategorizer.cc new file mode 100644 index 0000000000000..92631909541ab --- /dev/null +++ b/TopQuarkAnalysis/TopTools/plugins/GenTtbarCategorizer.cc @@ -0,0 +1,434 @@ +// -*- C++ -*- +// +// Package: TopQuarkAnalysis/TopTools +// Class: GenTtbarCategorizer +// +/**\class GenTtbarCategorizer GenTtbarCategorizer.cc TopQuarkAnalysis/TopTools/plugins/GenTtbarCategorizer.cc + + Description: Categorization of different tt+xx processes, returning unique ID for each process as e.g. tt+bb, tt+b, tt+2b, tt+cc, ... + + Implementation: + + The classification scheme returns an ID per event, and works as follows: + + All jets in the following need to be in the acceptance as given by the config parameters |eta|, pt. + A c jet must contain at least one c hadron and should contain no b hadrons + + First, b jets from top are identified, i.e. jets containing a b hadron from t->b decay + They are encoded in the ID as numberOfBjetsFromTop*100, i.e. + 0xx: no b jets from top in acceptance + 1xx: 1 b jet from top in acceptance + 2xx: both b jets from top in acceptance + + Then, b jets from W are identified, i.e. jets containing a b hadron from W->b decay + They are encoded in the ID as numberOfBjetsFromW*1000, i.e. + 0xxx: no b jets from W in acceptance + 1xxx: 1 b jet from W in acceptance + 2xxx: 2 b jets from W in acceptance + + Then, c jets from W are identified, i.e. jets containing a c hadron from W->c decay, but no b hadrons + They are encoded in the ID as numberOfCjetsFromW*10000, i.e. + 0xxxx: no c jets from W in acceptance + 1xxxx: 1 c jet from W in acceptance + 2xxxx: 2 c jets from W in acceptance + + 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: + x55: at least 2 additional b jets with at least two of them having >= 2 b hadrons in each + x54: at least 2 additional b jets with one of them having >= 2 b hadrons, the others having =1 b hadron + x53: at least 2 additional b jets with all having =1 b hadron + x52: exactly 1 additional b jet having >=2 b hadrons + x51: exactly 1 additional b jet having =1 b hadron + x45: at least 2 additional c jets with at least two of them having >= 2 c hadrons in each + x44: at least 2 additional c jets with one of them having >= 2 c hadrons, the others having =1 c hadron + x43: at least 2 additional c jets with all having =1 c hadron + x42: exactly 1 additional c jet having >=2 c hadrons + x41: exactly 1 additional c jet having =1 c hadron + x00: No additional b or c jet, i.e. only light flavour jets or no additional jets +*/ +// +// Original Author: Johannes Hauk, Nazar Bartosik +// Created: Sun, 14 Jun 2015 19:42:58 GMT +// +// + + +// system include files +#include +#include +#include + +// user include files +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/EDProducer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Utilities/interface/InputTag.h" +#include "DataFormats/JetReco/interface/GenJetCollection.h" +#include "DataFormats/HepMCCandidate/interface/GenParticle.h" + + +// +// class declaration +// + +class GenTtbarCategorizer : public edm::EDProducer { + public: + explicit GenTtbarCategorizer(const edm::ParameterSet&); + ~GenTtbarCategorizer(); + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + + private: + virtual void beginJob() override; + virtual void produce(edm::Event&, const edm::EventSetup&) override; + virtual void endJob() override; + + //virtual void beginRun(edm::Run const&, edm::EventSetup const&) override; + //virtual void endRun(edm::Run const&, edm::EventSetup const&) override; + //virtual void beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) override; + //virtual void endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) override; + + std::vector nHadronsOrderedJetIndices(const std::map& m_jetIndex); + + // ----------member data --------------------------- + + // Jet configuration + const double genJetPtMin_; + const double genJetAbsEtaMax_; + + // Input tags + const edm::EDGetTokenT genJetsToken_; + + const edm::EDGetTokenT > genBHadJetIndexToken_; + const edm::EDGetTokenT > genBHadFlavourToken_; + const edm::EDGetTokenT > genBHadFromTopWeakDecayToken_; + const edm::EDGetTokenT > genBHadPlusMothersToken_; + const edm::EDGetTokenT > > genBHadPlusMothersIndicesToken_; + const edm::EDGetTokenT > genBHadIndexToken_; + const edm::EDGetTokenT > genBHadLeptonHadronIndexToken_; + const edm::EDGetTokenT > genBHadLeptonViaTauToken_; + + const edm::EDGetTokenT > genCHadJetIndexToken_; + const edm::EDGetTokenT > genCHadFlavourToken_; + const edm::EDGetTokenT > genCHadFromTopWeakDecayToken_; + const edm::EDGetTokenT > genCHadBHadronIdToken_; + + +}; + +// +// constants, enums and typedefs +// + + +// +// static data member definitions +// + +// +// constructors and destructor +// +GenTtbarCategorizer::GenTtbarCategorizer(const edm::ParameterSet& iConfig): +genJetPtMin_(iConfig.getParameter("genJetPtMin")), +genJetAbsEtaMax_(iConfig.getParameter("genJetAbsEtaMax")), +genJetsToken_(consumes(iConfig.getParameter("genJets"))), +genBHadJetIndexToken_(consumes >(iConfig.getParameter("genBHadJetIndex"))), +genBHadFlavourToken_(consumes >(iConfig.getParameter("genBHadFlavour"))), +genBHadFromTopWeakDecayToken_(consumes >(iConfig.getParameter("genBHadFromTopWeakDecay"))), +genBHadPlusMothersToken_(consumes >(iConfig.getParameter("genBHadPlusMothers"))), +genBHadPlusMothersIndicesToken_(consumes > >(iConfig.getParameter("genBHadPlusMothersIndices"))), +genBHadIndexToken_(consumes >(iConfig.getParameter("genBHadIndex"))), +genBHadLeptonHadronIndexToken_(consumes >(iConfig.getParameter("genBHadLeptonHadronIndex"))), +genBHadLeptonViaTauToken_(consumes >(iConfig.getParameter("genBHadLeptonViaTau"))), +genCHadJetIndexToken_(consumes >(iConfig.getParameter("genCHadJetIndex"))), +genCHadFlavourToken_(consumes >(iConfig.getParameter("genCHadFlavour"))), +genCHadFromTopWeakDecayToken_(consumes >(iConfig.getParameter("genCHadFromTopWeakDecay"))), +genCHadBHadronIdToken_(consumes >(iConfig.getParameter("genCHadBHadronId"))) +{ + produces("genTtbarId"); +} + + +GenTtbarCategorizer::~GenTtbarCategorizer() +{} + + +// +// member functions +// + +// ------------ method called to produce the data ------------ +void +GenTtbarCategorizer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) +{ + // Access gen jets + edm::Handle genJets; + iEvent.getByToken(genJetsToken_, genJets); + + + // Access B hadrons information + edm::Handle > genBHadFlavour; + iEvent.getByToken(genBHadFlavourToken_, genBHadFlavour); + + edm::Handle > genBHadJetIndex; + iEvent.getByToken(genBHadJetIndexToken_, genBHadJetIndex); + + edm::Handle > genBHadFromTopWeakDecay; + iEvent.getByToken(genBHadFromTopWeakDecayToken_, genBHadFromTopWeakDecay); + + edm::Handle > genBHadPlusMothers; + iEvent.getByToken(genBHadPlusMothersToken_, genBHadPlusMothers); + + edm::Handle > > genBHadPlusMothersIndices; + iEvent.getByToken(genBHadPlusMothersIndicesToken_, genBHadPlusMothersIndices); + + edm::Handle > genBHadIndex; + iEvent.getByToken(genBHadIndexToken_, genBHadIndex); + + edm::Handle > genBHadLeptonHadronIndex; + iEvent.getByToken(genBHadLeptonHadronIndexToken_, genBHadLeptonHadronIndex); + + edm::Handle > genBHadLeptonViaTau; + iEvent.getByToken(genBHadLeptonViaTauToken_, genBHadLeptonViaTau); + + + // Access C hadrons information + edm::Handle > genCHadFlavour; + iEvent.getByToken(genCHadFlavourToken_, genCHadFlavour); + + edm::Handle > genCHadJetIndex; + iEvent.getByToken(genCHadJetIndexToken_, genCHadJetIndex); + + edm::Handle > genCHadFromTopWeakDecay; + iEvent.getByToken(genCHadFromTopWeakDecayToken_, genCHadFromTopWeakDecay); + + edm::Handle > genCHadBHadronId; + iEvent.getByToken(genCHadBHadronIdToken_, genCHadBHadronId); + + + // Map + // B jets with b hadrons directly from t->b decay + std::map bJetFromTopIds; + // B jets with b hadrons from W->b decay + std::map bJetFromWIds; + // C jets with c hadrons from W->c decay + std::map cJetFromWIds; + // B jets with b hadrons before top quark decay chain + std::map bJetAdditionalIds; + // C jets with c hadrons before top quark decay chain + std::map cJetAdditionalIds; + + + // Count number of specific b hadrons in each jet + for(size_t hadronId = 0; hadronId < genBHadIndex->size(); ++hadronId) { + // Index of jet associated to the hadron + const int jetIndex = genBHadJetIndex->at(hadronId); + // Skip hadrons which have no associated jet + if(jetIndex < 0) continue; + // Skip if jet is not in acceptance + if(genJets->at(jetIndex).pt() < genJetPtMin_) continue; + if(std::fabs(genJets->at(jetIndex).eta()) > genJetAbsEtaMax_) continue; + // Flavour of the hadron's origin + const int flavour = genBHadFlavour->at(hadronId); + // Jet from t->b decay [pdgId(top)=6] + if(std::abs(flavour) == 6) { + if(bJetFromTopIds.count(jetIndex) < 1) bJetFromTopIds[jetIndex] = 1; + else bJetFromTopIds[jetIndex]++; + continue; + } + // Jet from W->b decay [pdgId(W)=24] + if(std::abs(flavour) == 24) { + if(bJetFromWIds.count(jetIndex) < 1) bJetFromWIds[jetIndex] = 1; + else bJetFromWIds[jetIndex]++; + continue; + } + // Identify jets with b hadrons not from top-quark or W-boson decay + if(bJetAdditionalIds.count(jetIndex) < 1) bJetAdditionalIds[jetIndex] = 1; + else bJetAdditionalIds[jetIndex]++; + } + + // Cleaning up b jets from W->b decays + for(std::map::iterator it = bJetFromWIds.begin(); it != bJetFromWIds.end(); ) { + // Cannot be a b jet from t->b decay + if(bJetFromTopIds.count(it->first) > 0) bJetFromWIds.erase(it++); + else ++it; + } + + // Cleaning up additional b jets + for(std::map::iterator it = bJetAdditionalIds.begin(); it != bJetAdditionalIds.end(); ) { + // Cannot be a b jet from t->b decay + if(bJetFromTopIds.count(it->first) > 0) bJetAdditionalIds.erase(it++); + // Cannot be a b jet from W->b decay + else if(bJetFromWIds.count(it->first) > 0) bJetAdditionalIds.erase(it++); + else ++it; + } + + // Count number of specific c hadrons in each c jet + for(size_t hadronId = 0; hadronId < genCHadJetIndex->size(); ++hadronId) { + // Index of jet associated to the hadron + const int jetIndex = genCHadJetIndex->at(hadronId); + // Skip hadrons which have no associated jet + if(jetIndex < 0) continue; + // Skip c hadrons that are coming from b hadrons + if(genCHadBHadronId->at(hadronId) >= 0) continue; + // Skip if jet is not in acceptance + if(genJets->at(jetIndex).pt() < genJetPtMin_) continue; + if(std::fabs(genJets->at(jetIndex).eta()) > genJetAbsEtaMax_) continue; + // Skip if jet contains a b hadron + if(bJetFromTopIds.count(jetIndex) > 0) continue; + if(bJetFromWIds.count(jetIndex) > 0) continue; + if(bJetAdditionalIds.count(jetIndex) > 0) continue; + // Flavour of the hadron's origin + const int flavour = genCHadFlavour->at(hadronId); + // Jet from W->c decay [pdgId(W)=24] + if(std::abs(flavour) == 24) { + if(cJetFromWIds.count(jetIndex) < 1) cJetFromWIds[jetIndex] = 1; + else cJetFromWIds[jetIndex]++; + continue; + } + // Identify jets with c hadrons not from W-boson decay + if(cJetAdditionalIds.count(jetIndex) < 1) cJetAdditionalIds[jetIndex] = 1; + else cJetAdditionalIds[jetIndex]++; + } + + // Cleaning up additional c jets + for(std::map::iterator it = cJetAdditionalIds.begin(); it != cJetAdditionalIds.end(); ) { + // Cannot be a c jet from W->c decay + if(cJetFromWIds.count(it->first) > 0) cJetAdditionalIds.erase(it++); + else ++it; + } + + // Categorize event based on number of additional b/c jets + // and number of corresponding hadrons in each of them + int additionalJetEventId = bJetFromTopIds.size()*100 + bJetFromWIds.size()*1000 + cJetFromWIds.size()*10000; + // tt + 1 additional b jet + if(bJetAdditionalIds.size() == 1){ + const int nHadronsInJet = bJetAdditionalIds.begin()->second; + // tt + 1 additional b jet from 1 additional b hadron + if(nHadronsInJet == 1) additionalJetEventId += 51; + // tt + 1 additional b jet from >=2 additional b hadrons + else additionalJetEventId += 52; + } + // tt + >=2 additional b jets + else if(bJetAdditionalIds.size() > 1){ + // Check first two additional b jets (rare cases could have more) + const std::vector orderedJetIndices = nHadronsOrderedJetIndices(bJetAdditionalIds); + int nHadronsInJet1 = bJetAdditionalIds.at(orderedJetIndices.at(0)); + int nHadronsInJet2 = bJetAdditionalIds.at(orderedJetIndices.at(1)); + // tt + >=2 additional b jets each from 1 additional b hadron + if(std::max(nHadronsInJet1, nHadronsInJet2) == 1) additionalJetEventId += 53; + // tt + >=2 additional b jets one of which from >=2 additional b hadrons + else if(std::min(nHadronsInJet1, nHadronsInJet2) == 1 && std::max(nHadronsInJet1, nHadronsInJet2) > 1) additionalJetEventId += 54; + // tt + >=2 additional b jets each from >=2 additional b hadrons + else if(std::min(nHadronsInJet1, nHadronsInJet2) > 1) additionalJetEventId += 55; + } + // tt + no additional b jets + else{ + // tt + 1 additional c jet + if(cJetAdditionalIds.size() == 1){ + const int nHadronsInJet = cJetAdditionalIds.begin()->second; + // tt + 1 additional c jet from 1 additional c hadron + if(nHadronsInJet == 1) additionalJetEventId += 41; + // tt + 1 additional c jet from >=2 additional c hadrons + else additionalJetEventId += 42; + } + // tt + >=2 additional c jets + else if(cJetAdditionalIds.size() > 1){ + // Check first two additional c jets (rare cases could have more) + const std::vector orderedJetIndices = nHadronsOrderedJetIndices(cJetAdditionalIds); + int nHadronsInJet1 = cJetAdditionalIds.at(orderedJetIndices.at(0)); + int nHadronsInJet2 = cJetAdditionalIds.at(orderedJetIndices.at(1)); + // tt + >=2 additional c jets each from 1 additional c hadron + if(std::max(nHadronsInJet1, nHadronsInJet2) == 1) additionalJetEventId += 43; + // tt + >=2 additional c jets one of which from >=2 additional c hadrons + else if(std::min(nHadronsInJet1, nHadronsInJet2) == 1 && std::max(nHadronsInJet1, nHadronsInJet2) > 1) additionalJetEventId += 44; + // tt + >=2 additional c jets each from >=2 additional c hadrons + else if(std::min(nHadronsInJet1, nHadronsInJet2) > 1) additionalJetEventId += 45; + } + // tt + no additional c jets + else{ + // tt + light jets + additionalJetEventId += 0; + } + } + + std::auto_ptr ttbarId(new int); + *ttbarId = additionalJetEventId; + iEvent.put(ttbarId, "genTtbarId"); +} + +// ------------ method called once each job just before starting event loop ------------ +void +GenTtbarCategorizer::beginJob() +{} + +// ------------ method called once each job just after ending the event loop ------------ +void +GenTtbarCategorizer::endJob() +{} + +// ------------ method called when starting to processes a run ------------ +/* +void +GenTtbarCategorizer::beginRun(edm::Run const&, edm::EventSetup const&) +{ +} +*/ + +// ------------ method called when ending the processing of a run ------------ +/* +void +GenTtbarCategorizer::endRun(edm::Run const&, edm::EventSetup const&) +{ +} +*/ + +// ------------ method called when starting to processes a luminosity block ------------ +/* +void +GenTtbarCategorizer::beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) +{ +} +*/ + +// ------------ method called when ending the processing of a luminosity block ------------ +/* +void +GenTtbarCategorizer::endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) +{ +} +*/ + +// ------------ method returns a vector of jet indices from the given map, sorted by N hadrons in descending order ------------ +std::vector GenTtbarCategorizer::nHadronsOrderedJetIndices(const std::map& m_jetIndex) +{ + std::vector > v_jetNhadIndexPair; + for(std::map::const_iterator it = m_jetIndex.begin(); it != m_jetIndex.end(); ++it) { + const int jetIndex = it->first; + const int nHadrons = it->second; + v_jetNhadIndexPair.push_back( std::pair(nHadrons, jetIndex) ); + } + // Sorting the vector of pairs by their key value + std::sort(v_jetNhadIndexPair.begin(), v_jetNhadIndexPair.end(), std::greater >()); + // Building the vector of indices in the proper order + std::vector orderedJetIndices; + for(std::vector >::const_iterator it = v_jetNhadIndexPair.begin(); it != v_jetNhadIndexPair.end(); ++it) { + orderedJetIndices.push_back(it->second); + } + + return orderedJetIndices; +} + +// ------------ method fills 'descriptions' with the allowed parameters for the module ------------ +void +GenTtbarCategorizer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + //The following says we do not know what parameters are allowed so do no validation + // Please change this to state exactly what you do use, even if it is no parameters + edm::ParameterSetDescription desc; + desc.setUnknown(); + descriptions.addDefault(desc); +} + +//define this as a plug-in +DEFINE_FWK_MODULE(GenTtbarCategorizer); diff --git a/TopQuarkAnalysis/TopTools/python/GenTtbarCategorizer_cfi.py b/TopQuarkAnalysis/TopTools/python/GenTtbarCategorizer_cfi.py new file mode 100644 index 0000000000000..cd761472a86c3 --- /dev/null +++ b/TopQuarkAnalysis/TopTools/python/GenTtbarCategorizer_cfi.py @@ -0,0 +1,24 @@ +import FWCore.ParameterSet.Config as cms + +categorizeGenTtbar = cms.EDProducer("GenTtbarCategorizer", + # Phase space of additional jets + genJetPtMin = cms.double(20.), + genJetAbsEtaMax = cms.double(2.4), + # Input tags holding information about b/c hadron matching + genJets = cms.InputTag("ak4GenJets"), + genBHadJetIndex = cms.InputTag("matchGenBHadron", "genBHadJetIndex"), + genBHadFlavour = cms.InputTag("matchGenBHadron", "genBHadFlavour"), + genBHadFromTopWeakDecay = cms.InputTag("matchGenBHadron", "genBHadFromTopWeakDecay"), + genBHadPlusMothers = cms.InputTag("matchGenBHadron", "genBHadPlusMothers"), + genBHadPlusMothersIndices = cms.InputTag("matchGenBHadron", "genBHadPlusMothersIndices"), + genBHadIndex = cms.InputTag("matchGenBHadron", "genBHadIndex"), + genBHadLeptonHadronIndex = cms.InputTag("matchGenBHadron", "genBHadLeptonHadronIndex"), + genBHadLeptonViaTau = cms.InputTag("matchGenBHadron", "genBHadLeptonViaTau"), + genCHadJetIndex = cms.InputTag("matchGenCHadron", "genCHadJetIndex"), + genCHadFlavour = cms.InputTag("matchGenCHadron", "genCHadFlavour"), + genCHadFromTopWeakDecay = cms.InputTag("matchGenCHadron", "genCHadFromTopWeakDecay"), + genCHadBHadronId = cms.InputTag("matchGenCHadron", "genCHadBHadronId"), +) + + + diff --git a/TopQuarkAnalysis/TopTools/test/BuildFile.xml b/TopQuarkAnalysis/TopTools/test/BuildFile.xml index ad0843ec15962..8ae1fb544eab9 100644 --- a/TopQuarkAnalysis/TopTools/test/BuildFile.xml +++ b/TopQuarkAnalysis/TopTools/test/BuildFile.xml @@ -4,3 +4,11 @@ + + + + + + + + diff --git a/TopQuarkAnalysis/TopTools/test/TestGenTtbarCategories.cc b/TopQuarkAnalysis/TopTools/test/TestGenTtbarCategories.cc new file mode 100644 index 0000000000000..28a88e611380e --- /dev/null +++ b/TopQuarkAnalysis/TopTools/test/TestGenTtbarCategories.cc @@ -0,0 +1,185 @@ +// -*- C++ -*- +// +// Package: TopQuarkAnalysis/TopTools +// Class: TestGenTtbarCategories +// +/**\class TestGenTtbarCategories TestGenTtbarCategories.cc PhysicsTools/JetMCAlgos/test/TestGenTtbarCategories.cc + + Description: Analyzer for testing corresponding producer GenTtbarCategorizer + + Implementation: + [Notes on implementation] +*/ +// +// Original Author: Johannes Hauk, Nazar Bartosik +// Created: Sun, 14 Jun 2015 21:00:23 GMT +// +// + + +// system include files +#include +#include + +// user include files +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/EDAnalyzer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Utilities/interface/InputTag.h" +#include "FWCore/ServiceRegistry/interface/Service.h" +#include "CommonTools/UtilAlgos/interface/TFileService.h" + +#include + +// +// class declaration +// + +class TestGenTtbarCategories : public edm::EDAnalyzer { + public: + explicit TestGenTtbarCategories(const edm::ParameterSet&); + ~TestGenTtbarCategories(); + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + + + private: + virtual void beginJob() override; + virtual void analyze(const edm::Event&, const edm::EventSetup&) override; + virtual void endJob() override; + + //virtual void beginRun(edm::Run const&, edm::EventSetup const&) override; + //virtual void endRun(edm::Run const&, edm::EventSetup const&) override; + //virtual void beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) override; + //virtual void endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) override; + + // ----------member data --------------------------- + + // Input tags + const edm::EDGetTokenT genTtbarIdToken_; + + // Variables to fill + int ttbarId_; + int ttbarAdditionalJetId_; + int nBjetsFromTop_; + int nBjetsFromW_; + int nCjetsFromW_; + + // Tree to be filled + TTree* tree_; + +}; + +// +// constants, enums and typedefs +// + +// +// static data member definitions +// + +// +// constructors and destructor +// +TestGenTtbarCategories::TestGenTtbarCategories(const edm::ParameterSet& iConfig): +genTtbarIdToken_(consumes(iConfig.getParameter("genTtbarId"))) +{} + + +TestGenTtbarCategories::~TestGenTtbarCategories() +{} + + +// +// member functions +// + +// ------------ method called for each event ------------ +void +TestGenTtbarCategories::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) +{ + edm::Handle genTtbarId; + iEvent.getByToken(genTtbarIdToken_, genTtbarId); + + // ID including information about b/c jets in acceptance from t->b/W->b/W->c decays as well as additional ones + ttbarId_ = *genTtbarId; + + // ID based only on additional b/c jets + ttbarAdditionalJetId_ = ttbarId_%100; + + // Number of b/c jets from t->b or W->b/c decays + nBjetsFromTop_ = ttbarId_%1000/100; + nBjetsFromW_ = ttbarId_%10000/1000; + nCjetsFromW_ = ttbarId_%100000/10000; + + // Filling the tree + tree_->Fill(); +} + + +// ------------ method called once each job just before starting event loop ------------ +void +TestGenTtbarCategories::beginJob() +{ + edm::Service fileService; + if(!fileService) throw edm::Exception(edm::errors::Configuration, "TFileService is not registered in cfg file"); + + tree_ = fileService->make("tree", "tree"); + tree_->Branch("ttbarId", &ttbarId_); + tree_->Branch("ttbarAdditionalJetId", &ttbarAdditionalJetId_); + tree_->Branch("nBjetsFromTop", &nBjetsFromTop_); + tree_->Branch("nBjetsFromW", &nBjetsFromW_); + tree_->Branch("nCjetsFromW", &nCjetsFromW_); +} + +// ------------ method called once each job just after ending the event loop ------------ +void +TestGenTtbarCategories::endJob() +{} + +// ------------ method called when starting to processes a run ------------ +/* +void +TestGenTtbarCategories::beginRun(edm::Run const&, edm::EventSetup const&) +{ +} +*/ + +// ------------ method called when ending the processing of a run ------------ +/* +void +TestGenTtbarCategories::endRun(edm::Run const&, edm::EventSetup const&) +{ +} +*/ + +// ------------ method called when starting to processes a luminosity block ------------ +/* +void +TestGenTtbarCategories::beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) +{ +} +*/ + +// ------------ method called when ending the processing of a luminosity block ------------ +/* +void +TestGenTtbarCategories::endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) +{ +} +*/ + +// ------------ method fills 'descriptions' with the allowed parameters for the module ------------ +void +TestGenTtbarCategories::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + //The following says we do not know what parameters are allowed so do no validation + // Please change this to state exactly what you do use, even if it is no parameters + edm::ParameterSetDescription desc; + desc.setUnknown(); + descriptions.addDefault(desc); +} + +//define this as a plug-in +DEFINE_FWK_MODULE(TestGenTtbarCategories); diff --git a/TopQuarkAnalysis/TopTools/test/testGenTtbarCategories.py b/TopQuarkAnalysis/TopTools/test/testGenTtbarCategories.py new file mode 100644 index 0000000000000..c66af0b2c7f70 --- /dev/null +++ b/TopQuarkAnalysis/TopTools/test/testGenTtbarCategories.py @@ -0,0 +1,142 @@ +import FWCore.ParameterSet.Config as cms +import FWCore.ParameterSet.VarParsing as VarParsing +import sys + +## Define the process +process = cms.Process("Analyzer") +process.options = cms.untracked.PSet( + wantSummary = cms.untracked.bool(True), + allowUnscheduled = cms.untracked.bool(True), +) + +## Set up command line options +options = VarParsing.VarParsing ('standard') +options.register('runOnAOD', True, VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.bool, "run on AOD") + +## Get and parse command line options +if( hasattr(sys, "argv") ): + for args in sys.argv : + arg = args.split(',') + for val in arg: + val = val.split('=') + if(len(val)==2): + setattr(options,val[0], val[1]) + +## Configure message logger +process.load("FWCore.MessageLogger.MessageLogger_cfi") +process.MessageLogger.cerr.threshold = 'INFO' +process.MessageLogger.cerr.FwkReport.reportEvery = 100 + +## Define input +if options.runOnAOD: + process.source = cms.Source("PoolSource", + fileNames = cms.untracked.vstring( + # add your favourite AOD files here +# '/store/relval/CMSSW_7_2_0_pre7/RelValTTbar_13/GEN-SIM-RECO/PU50ns_PRE_LS172_V12-v1/00000/1267B7ED-2F4E-E411-A0B9-0025905964A6.root', + '/store/mc/Spring14dr/TTJets_MSDecaysCKM_central_Tune4C_13TeV-madgraph-tauola/AODSIM/PU_S14_POSTLS170_V6-v1/00000/00120F7A-84F5-E311-9FBE-002618943910.root', +# '/store/mc/Spring14dr/TT_Tune4C_13TeV-pythia8-tauola/AODSIM/Flat20to50_POSTLS170_V5-v1/00000/023E1847-ADDC-E311-91A2-003048FFD754.root', +# '/store/mc/Spring14dr/TTbarH_HToBB_M-125_13TeV_pythia6/AODSIM/PU20bx25_POSTLS170_V5-v1/00000/1CAB7E58-0BD0-E311-B688-00266CFFBC3C.root', +# '/store/mc/Spring14dr/TTbarH_M-125_13TeV_amcatnlo-pythia8-tauola/AODSIM/PU20bx25_POSTLS170_V5-v1/00000/0E3D08A9-C610-E411-A862-0025B3E0657E.root', + ), + skipEvents = cms.untracked.uint32(0) + ) +else: + process.source = cms.Source("PoolSource", + fileNames = cms.untracked.vstring( + # add your favourite miniAOD files here +# '/store/mc/Phys14DR/TTJets_MSDecaysCKM_central_Tune4C_13TeV-madgraph-tauola/MINIAODSIM/PU20bx25_PHYS14_25_V1-v1/00000/FE26BEB8-D575-E411-A13E-00266CF2AE10.root', +# '/store/mc/RunIISpring15DR74/TTJets_TuneCUETP8M1_13TeV-amcatnloFXFX-pythia8/MINIAODSIM/Asympt25ns_MCRUN2_74_V9-v1/00000/FE3EF0E9-9F02-E511-BA6F-549F35AE4FA2.root', + '/store/mc/RunIISpring15DR74/TTJets_TuneCUETP8M1_13TeV-amcatnloFXFX-scaledown-pythia8/MINIAODSIM/Asympt50ns_MCRUN2_74_V9A-v2/50000/FEA9ADAC-FB0B-E511-9B8C-00266CF9AB9C.root', +# '/store/mc/RunIISpring15DR74/TT_TuneZ2star_13TeV-powheg-pythia6-tauola/MINIAODSIM/Asympt25ns_MCRUN2_74_V9-v1/50000/C6F89D57-640A-E511-8D41-549F35AD8BC9.root', +# '/store/mc/RunIISpring15DR74/TT_TuneZ2star_13TeV-powheg-pythia6-tauola/MINIAODSIM/Asympt50ns_MCRUN2_74_V9A-v3/00000/D4BE2E58-CE08-E511-9247-0025907DC9CC.root', + ), + skipEvents = cms.untracked.uint32(0) + ) + +## Define maximum number of events to loop over +process.maxEvents = cms.untracked.PSet( + input = cms.untracked.int32(1000) +) + +## Set input particle collections to be used by the tools +genParticleCollection = '' +genJetInputParticleCollection = '' +genJetCollection = 'ak4GenJetsCustom' +if options.runOnAOD: + genParticleCollection = 'genParticles' + genJetInputParticleCollection = 'genParticlesForJets' + ## producing a subset of genParticles to be used for jet clustering in AOD + from RecoJets.Configuration.GenJetParticles_cff import genParticlesForJets + process.genParticlesForJets = genParticlesForJets.clone() +else: + genParticleCollection = 'prunedGenParticles' + genJetInputParticleCollection = 'packedGenParticles' + +## Supplies PDG ID to real name resolution of MC particles +process.load("SimGeneral.HepPDTESSource.pythiapdt_cfi") + +## Produce own jets (re-clustering in miniAOD needed at present to avoid crash) +from RecoJets.JetProducers.ak4GenJets_cfi import ak4GenJets +process.ak4GenJetsCustom = ak4GenJets.clone( + src = genJetInputParticleCollection, + rParam = cms.double(0.4), + jetAlgorithm = cms.string("AntiKt") +) + +## Ghost particle collection used for Hadron-Jet association +# MUST use proper input particle collection +from PhysicsTools.JetMCAlgos.HadronAndPartonSelector_cfi import selectedHadronsAndPartons +process.selectedHadronsAndPartons = selectedHadronsAndPartons.clone( + particles = genParticleCollection +) + +## Input particle collection for matching to gen jets (partons + leptons) +# MUST use use proper input jet collection: the jets to which hadrons should be associated +# rParam and jetAlgorithm MUST match those used for jets to be associated with hadrons +# More details on the tool: https://twiki.cern.ch/twiki/bin/view/CMSPublic/SWGuideBTagMCTools#New_jet_flavour_definition +from PhysicsTools.JetMCAlgos.AK4PFJetsMCFlavourInfos_cfi import ak4JetFlavourInfos +process.genJetFlavourInfos = ak4JetFlavourInfos.clone( + jets = genJetCollection, +) + + +## Plugin for analysing B hadrons +# MUST use the same particle collection as in selectedHadronsAndPartons +from PhysicsTools.JetMCAlgos.GenHFHadronMatcher_cff import matchGenBHadron +process.matchGenBHadron = matchGenBHadron.clone( + genParticles = genParticleCollection, + jetFlavourInfos = "genJetFlavourInfos" +) + +## Plugin for analysing C hadrons +# MUST use the same particle collection as in selectedHadronsAndPartons +from PhysicsTools.JetMCAlgos.GenHFHadronMatcher_cff import matchGenCHadron +process.matchGenCHadron = matchGenCHadron.clone( + genParticles = genParticleCollection, + jetFlavourInfos = "genJetFlavourInfos" +) + +## Producer for ttbar categorisation ID +# MUST use same genJetCollection as used for tools above +from TopQuarkAnalysis.TopTools.GenTtbarCategorizer_cfi import categorizeGenTtbar +process.categorizeGenTtbar = categorizeGenTtbar.clone( + genJetPtMin = 20., + genJetAbsEtaMax = 2.4, + genJets = genJetCollection, +) + +## Configure test analyzer +process.testGenTtbarCategories = cms.EDAnalyzer("TestGenTtbarCategories", + genTtbarId = cms.InputTag("categorizeGenTtbar", "genTtbarId") +) + +## Output root file +process.TFileService = cms.Service("TFileService", + fileName = cms.string("genTtbarId.root") +) + +## Path +process.p1 = cms.Path( + process.testGenTtbarCategories +) + From f97cf64232c336e9f50bb8d935f8d26e90922bf1 Mon Sep 17 00:00:00 2001 From: Nazar Bartosik Date: Tue, 18 Aug 2015 13:18:43 +0200 Subject: [PATCH 2/3] Optimised vector filling performance with reserved size if known --- TopQuarkAnalysis/TopTools/plugins/GenTtbarCategorizer.cc | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/TopQuarkAnalysis/TopTools/plugins/GenTtbarCategorizer.cc b/TopQuarkAnalysis/TopTools/plugins/GenTtbarCategorizer.cc index 92631909541ab..ef7a56f83fafc 100644 --- a/TopQuarkAnalysis/TopTools/plugins/GenTtbarCategorizer.cc +++ b/TopQuarkAnalysis/TopTools/plugins/GenTtbarCategorizer.cc @@ -403,7 +403,9 @@ GenTtbarCategorizer::endLuminosityBlock(edm::LuminosityBlock const&, edm::EventS // ------------ method returns a vector of jet indices from the given map, sorted by N hadrons in descending order ------------ std::vector GenTtbarCategorizer::nHadronsOrderedJetIndices(const std::map& m_jetIndex) { + const int nElements = m_jetIndex.size(); std::vector > v_jetNhadIndexPair; + v_jetNhadIndexPair.reserve(nElements); for(std::map::const_iterator it = m_jetIndex.begin(); it != m_jetIndex.end(); ++it) { const int jetIndex = it->first; const int nHadrons = it->second; @@ -412,12 +414,13 @@ std::vector GenTtbarCategorizer::nHadronsOrderedJetIndices(const std::map >()); // Building the vector of indices in the proper order - std::vector orderedJetIndices; + std::vector v_orderedJetIndices; + v_orderedJetIndices.reserve(nElements); for(std::vector >::const_iterator it = v_jetNhadIndexPair.begin(); it != v_jetNhadIndexPair.end(); ++it) { - orderedJetIndices.push_back(it->second); + v_orderedJetIndices.push_back(it->second); } - return orderedJetIndices; + return v_orderedJetIndices; } // ------------ method fills 'descriptions' with the allowed parameters for the module ------------ From 9b75e9071408c13038e2c45a477e51fe291ea46d Mon Sep 17 00:00:00 2001 From: Nazar Bartosik Date: Wed, 2 Sep 2015 13:30:35 +0200 Subject: [PATCH 3/3] No jet reclustering in miniAOD for genTtbar event categorisation. Exclusded leptons from jets in AOD. --- .../TopTools/test/testGenTtbarCategories.py | 25 +++++++++---------- 1 file changed, 12 insertions(+), 13 deletions(-) diff --git a/TopQuarkAnalysis/TopTools/test/testGenTtbarCategories.py b/TopQuarkAnalysis/TopTools/test/testGenTtbarCategories.py index c66af0b2c7f70..d51bce417f3be 100644 --- a/TopQuarkAnalysis/TopTools/test/testGenTtbarCategories.py +++ b/TopQuarkAnalysis/TopTools/test/testGenTtbarCategories.py @@ -60,28 +60,27 @@ ## Set input particle collections to be used by the tools genParticleCollection = '' -genJetInputParticleCollection = '' -genJetCollection = 'ak4GenJetsCustom' +genJetCollection = '' if options.runOnAOD: genParticleCollection = 'genParticles' - genJetInputParticleCollection = 'genParticlesForJets' + genJetCollection = 'ak4GenJetsCustom' ## producing a subset of genParticles to be used for jet clustering in AOD - from RecoJets.Configuration.GenJetParticles_cff import genParticlesForJets - process.genParticlesForJets = genParticlesForJets.clone() + from RecoJets.Configuration.GenJetParticles_cff import genParticlesForJetsNoNu + process.genParticlesForJetsCustom = genParticlesForJetsNoNu.clone() + ## Produce own jets (re-clustering in miniAOD needed at present to avoid crash) + from RecoJets.JetProducers.ak4GenJets_cfi import ak4GenJets + process.ak4GenJetsCustom = ak4GenJets.clone( + src = 'genParticlesForJetsCustom', + rParam = cms.double(0.4), + jetAlgorithm = cms.string("AntiKt") + ) else: genParticleCollection = 'prunedGenParticles' - genJetInputParticleCollection = 'packedGenParticles' + genJetCollection = 'slimmedGenJets' ## Supplies PDG ID to real name resolution of MC particles process.load("SimGeneral.HepPDTESSource.pythiapdt_cfi") -## Produce own jets (re-clustering in miniAOD needed at present to avoid crash) -from RecoJets.JetProducers.ak4GenJets_cfi import ak4GenJets -process.ak4GenJetsCustom = ak4GenJets.clone( - src = genJetInputParticleCollection, - rParam = cms.double(0.4), - jetAlgorithm = cms.string("AntiKt") -) ## Ghost particle collection used for Hadron-Jet association # MUST use proper input particle collection