diff --git a/PhysicsTools/JetMCAlgos/python/GenHFHadronMatcher_cfi.py b/PhysicsTools/JetMCAlgos/python/GenHFHadronMatcher_cfi.py index 89d8966270706..d4f5487e0db45 100644 --- a/PhysicsTools/JetMCAlgos/python/GenHFHadronMatcher_cfi.py +++ b/PhysicsTools/JetMCAlgos/python/GenHFHadronMatcher_cfi.py @@ -5,7 +5,7 @@ jetFlavourInfos = cms.InputTag("genJetFlavourInfos"), flavour = cms.int32(5), onlyJetClusteredHadrons = cms.bool(True), - noBBbarResonances = cms.bool(True), + noBBbarResonances = cms.bool(False), ) diff --git a/PhysicsTools/JetMCAlgos/test/BuildFile.xml b/PhysicsTools/JetMCAlgos/test/BuildFile.xml index 0aed79a8249d3..d5a58b05d2bc6 100644 --- a/PhysicsTools/JetMCAlgos/test/BuildFile.xml +++ b/PhysicsTools/JetMCAlgos/test/BuildFile.xml @@ -20,3 +20,11 @@ + + + + + + + + diff --git a/PhysicsTools/JetMCAlgos/test/genHFHadronMatcher.py b/PhysicsTools/JetMCAlgos/test/genHFHadronMatcher.py deleted file mode 100644 index ec143d9f94e3a..0000000000000 --- a/PhysicsTools/JetMCAlgos/test/genHFHadronMatcher.py +++ /dev/null @@ -1,115 +0,0 @@ -import FWCore.ParameterSet.Config as cms - -process = cms.Process("Analyzer") - - -## setting the format of the input files AOD/miniAOD -runOnAOD = True - -## enabling unscheduled mode for modules -process.options = cms.untracked.PSet( - wantSummary = cms.untracked.bool(True), - allowUnscheduled = cms.untracked.bool(True), -) - -## configure message logger -process.load("FWCore.MessageLogger.MessageLogger_cfi") -process.MessageLogger.cerr.threshold = 'INFO' -process.MessageLogger.cerr.FwkReport.reportEvery = 100 - -## define input -process.source = cms.Source("PoolSource", - fileNames = cms.untracked.vstring( - ## add your favourite file here (1000+ events in each file defined below) - ## Check runOnAOD parameter to match the format of input files - # AOD - '/store/relval/CMSSW_7_2_0_pre7/RelValTTbar_13/GEN-SIM-RECO/PU50ns_PRE_LS172_V12-v1/00000/1267B7ED-2F4E-E411-A0B9-0025905964A6.root', - - # miniAOD -# '/store/cmst3/user/gpetrucc/miniAOD/v1/TTbarH_HToBB_M-125_13TeV_pythia6_PU_S14_PAT.root', - - # other AOD samples -# '/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) - ) - -## define maximal number of events to loop over -process.maxEvents = cms.untracked.PSet( - input = cms.untracked.int32(1000) -) - -# Setting input particle collections to be used by the tools -genParticleCollection = '' -genJetInputParticleCollection = '' -genJetCollection = 'ak5GenJetsCustom' -if 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") - -# Producing own jets for testing purposes -from RecoJets.JetProducers.ak5GenJets_cfi import ak5GenJets -process.ak5GenJetsCustom = ak5GenJets.clone( - src = genJetInputParticleCollection, - rParam = cms.double(0.5), - 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.sequences.GenHFHadronMatching_cff import genJetFlavourPlusLeptonInfos -process.genJetFlavourPlusLeptonInfos = genJetFlavourPlusLeptonInfos.clone( - jets = genJetCollection, - rParam = cms.double(0.5), - jetAlgorithm = cms.string("AntiKt") -) - - -# Plugin for analysing B hadrons -# MUST use the same particle collection as in selectedHadronsAndPartons -from PhysicsTools.JetMCAlgos.sequences.GenHFHadronMatching_cff import matchGenBHadron -process.matchGenBHadron = matchGenBHadron.clone( - genParticles = genParticleCollection -) - -# Plugin for analysing C hadrons -# MUST use the same particle collection as in selectedHadronsAndPartons -from PhysicsTools.JetMCAlgos.sequences.GenHFHadronMatching_cff import matchGenCHadron -process.matchGenCHadron = matchGenCHadron.clone( - genParticles = genParticleCollection -) - - -## defining only the final modules to run: dependencies will be run automatically [allowUnscheduled = True] -process.p1 = cms.Path( - process.matchGenBHadron * - process.matchGenCHadron -) - -## module to store raw output from the processed modules into the ROOT file -process.out = cms.OutputModule("PoolOutputModule", - fileName = cms.untracked.string('genHFHadronMatcher_out.root'), - outputCommands = cms.untracked.vstring('drop *', 'keep *_matchGen*_*_*') - ) -process.outpath = cms.EndPath(process.out) diff --git a/PhysicsTools/JetMCAlgos/test/matchGenHFHadrons.cc b/PhysicsTools/JetMCAlgos/test/matchGenHFHadrons.cc new file mode 100644 index 0000000000000..2ce3f8210a3e1 --- /dev/null +++ b/PhysicsTools/JetMCAlgos/test/matchGenHFHadrons.cc @@ -0,0 +1,577 @@ +// system include files +#include +#include +#include + +// user include files +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/EDAnalyzer.h" +#include "FWCore/Utilities/interface/InputTag.h" + +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" + + +#include "DataFormats/JetReco/interface/Jet.h" +#include "DataFormats/JetReco/interface/GenJetCollection.h" +#include "DataFormats/HepMCCandidate/interface/GenParticle.h" + +#include "FWCore/ServiceRegistry/interface/Service.h" +#include "CommonTools/UtilAlgos/interface/TFileService.h" +#include +#include + + +class matchGenHFHadrons : public edm::EDAnalyzer { + public: + explicit matchGenHFHadrons(const edm::ParameterSet & ); + ~matchGenHFHadrons() {}; + + private: + virtual void analyze(const edm::Event& event, const edm::EventSetup& setup); + virtual void beginJob(); + void clearEventVariables(); + void clearBHadronVariables(); + // A recursive function that scans through the particle chain to find a particle with specific pdgId or abs(pdgId) + void findAncestorParticles(const int startId, std::set& resultIds, + const std::vector& genParticles, const std::vector >& motherIndices, + const int endPdgId, const bool endPdgIdIsAbs = false, const int firstAnyLast = 0); + + + // Jets 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_; + + // Output to be written to the trees + TTree* tree_events_; + TTree* tree_bHadrons_; + + int nJets_; + + int nBJets_; + int nBHadrons_; + int nBHadronsTop_; + int nBHadronsAdditional_; + int nBHadronsPseudoadditional_; + + int nCJets_; + int nCHadrons_; + int nCHadronsAdditional_; + int nCHadronsPseudoadditional_; + + int nBHadronsHiggs_; + int additionalJetEventId_; + + // Additional performance information + int bHadron_flavour_; + int bHadron_nQuarksAll_; + int bHadron_nQuarksSameFlavour_; + + double bHadron_quark1_dR_; + double bHadron_quark2_dR_; + double bHadron_quark1_ptRatio_; + double bHadron_quark2_ptRatio_; + + int bHadron_nLeptonsAll_; + int bHadron_nLeptonsViaTau_; + + double bHadron_jetPtRatio_; +}; + +bool sort_pairs_second(const std::pair& a, const std::pair& b) +{ + return a.second < b.second; +} + +matchGenHFHadrons::matchGenHFHadrons ( const edm::ParameterSet& config): + genJetPtMin_(config.getParameter("genJetPtMin")), + genJetAbsEtaMax_(config.getParameter("genJetAbsEtaMax")), + genJetsToken_(consumes(config.getParameter("genJets"))), + genBHadJetIndexToken_(consumes >(config.getParameter("genBHadJetIndex"))), + genBHadFlavourToken_(consumes >(config.getParameter("genBHadFlavour"))), + genBHadFromTopWeakDecayToken_(consumes >(config.getParameter("genBHadFromTopWeakDecay"))), + genBHadPlusMothersToken_(consumes >(config.getParameter("genBHadPlusMothers"))), + genBHadPlusMothersIndicesToken_(consumes > >(config.getParameter("genBHadPlusMothersIndices"))), + genBHadIndexToken_(consumes >(config.getParameter("genBHadIndex"))), + genBHadLeptonHadronIndexToken_(consumes >(config.getParameter("genBHadLeptonHadronIndex"))), + genBHadLeptonViaTauToken_(consumes >(config.getParameter("genBHadLeptonViaTau"))), + genCHadJetIndexToken_(consumes >(config.getParameter("genCHadJetIndex"))), + genCHadFlavourToken_(consumes >(config.getParameter("genCHadFlavour"))), + genCHadFromTopWeakDecayToken_(consumes >(config.getParameter("genCHadFromTopWeakDecay"))), + genCHadBHadronIdToken_(consumes >(config.getParameter("genCHadBHadronId"))) +{ + +} + + +void matchGenHFHadrons::analyze(const edm::Event& event, const edm::EventSetup& setup) +{ + this->clearEventVariables(); + + // Reading gen jets from the event + edm::Handle genJets; + event.getByToken(genJetsToken_, genJets); + + // Reading B hadrons related information + edm::Handle > genBHadFlavour; + event.getByToken(genBHadFlavourToken_, genBHadFlavour); + + edm::Handle > genBHadJetIndex; + event.getByToken(genBHadJetIndexToken_, genBHadJetIndex); + + edm::Handle > genBHadFromTopWeakDecay; + event.getByToken(genBHadFromTopWeakDecayToken_, genBHadFromTopWeakDecay); + + edm::Handle > genBHadPlusMothers; + event.getByToken(genBHadPlusMothersToken_, genBHadPlusMothers); + + edm::Handle > > genBHadPlusMothersIndices; + event.getByToken(genBHadPlusMothersIndicesToken_, genBHadPlusMothersIndices); + + edm::Handle > genBHadIndex; + event.getByToken(genBHadIndexToken_, genBHadIndex); + + edm::Handle > genBHadLeptonHadronIndex; + event.getByToken(genBHadLeptonHadronIndexToken_, genBHadLeptonHadronIndex); + + edm::Handle > genBHadLeptonViaTau; + event.getByToken(genBHadLeptonViaTauToken_, genBHadLeptonViaTau); + + // Reading C hadrons related information + edm::Handle > genCHadFlavour; + event.getByToken(genCHadFlavourToken_, genCHadFlavour); + + edm::Handle > genCHadJetIndex; + event.getByToken(genCHadJetIndexToken_, genCHadJetIndex); + + edm::Handle > genCHadFromTopWeakDecay; + event.getByToken(genCHadFromTopWeakDecayToken_, genCHadFromTopWeakDecay); + + edm::Handle > genCHadBHadronId; + event.getByToken(genCHadBHadronIdToken_, genCHadBHadronId); + + + // Counting number of jets of different flavours in the event + for(size_t jetId = 0; jetId < genJets->size(); ++jetId) { + if(genJets->at(jetId).pt() < genJetPtMin_) continue; + if(std::fabs(genJets->at(jetId).eta()) < genJetAbsEtaMax_) continue; + nJets_++; + if(std::find(genBHadJetIndex->begin(), genBHadJetIndex->end(), jetId) != genBHadJetIndex->end()) nBJets_++; + else if(std::find(genCHadJetIndex->begin(), genCHadJetIndex->end(), jetId) != genCHadJetIndex->end()) nCJets_++; + } + // Counting number of different b hadrons + for(size_t hadronId = 0; hadronId < genBHadFlavour->size(); ++hadronId) { + const int flavour = genBHadFlavour->at(hadronId); + const int flavourAbs = std::abs(flavour); + const bool afterTop = genBHadFromTopWeakDecay->at(hadronId) == 1 ? true : false; + nBHadrons_++; + if(flavourAbs==25) nBHadronsHiggs_++; + if(flavourAbs==6) nBHadronsTop_++; + if(flavourAbs!=6) { + if(afterTop) nBHadronsPseudoadditional_++; + else nBHadronsAdditional_++; + } + } + // Counting number of different c hadrons + for(size_t hadronId = 0; hadronId < genCHadJetIndex->size(); ++hadronId) { + const bool afterTop = genCHadFromTopWeakDecay->at(hadronId) == 1 ? true : false; + nCHadrons_++; + // Skipping c hadrons that are coming from b hadrons + if(genCHadBHadronId->at(hadronId) >= 0) continue; + + if(afterTop) nCHadronsPseudoadditional_++; + else nCHadronsAdditional_++; + } + + // Map + // B jets with b hadrons directly from top quark decay + std::map bJetFromTopIds_all; + // B jets with b hadrons directly from top quark decay + std::map bJetFromTopIds; + // B jets with b hadrons after top quark decay + std::map bJetAfterTopIds; + // B jets with b hadrons before top quark decay chain + std::map bJetBeforeTopIds; + // C jets with c hadrons before top quark decay chain + std::map cJetBeforeTopIds; + // C jets with c hadrons after top quark decay + std::map cJetAfterTopIds; + + // Counting number of specific hadrons in each b jet + for(size_t hadronId = 0; hadronId < genBHadIndex->size(); ++hadronId) { + // Flavour of the hadron's origin + const int flavour = genBHadFlavour->at(hadronId); + // Whether hadron radiated before top quark decay + const bool fromTopDecay = genBHadFromTopWeakDecay->at(hadronId); + // Index of a jet associated to the hadron + const int jetIndex = genBHadJetIndex->at(hadronId); + // Skipping hadrons which have no associated jet + if(jetIndex < 0) continue; + // Jet from direct top quark decay [pdgId(top)=6] + if(std::abs(flavour) == 6) { + if(bJetFromTopIds_all.count(jetIndex) < 1) bJetFromTopIds_all[jetIndex] = 1; + else bJetFromTopIds_all[jetIndex]++; + } + // Skipping if jet is not in acceptance + if(genJets->at(jetIndex).pt() < genJetPtMin_) continue; + if(std::fabs(genJets->at(jetIndex).eta()) > genJetAbsEtaMax_) continue; + // Identifying jets with b hadrons not from top quark decay + // Jet from direct top quark decay [pdgId(top)=6] + if(std::abs(flavour) == 6) { + if(bJetFromTopIds.count(jetIndex) < 1) bJetFromTopIds[jetIndex] = 1; + else bJetFromTopIds[jetIndex]++; + } + // Skipping if jet is from top quark decay + if(std::abs(flavour) == 6) continue; + // Jet before top quark decay + if(!fromTopDecay) { + if(bJetBeforeTopIds.count(jetIndex) < 1) bJetBeforeTopIds[jetIndex] = 1; + else bJetBeforeTopIds[jetIndex]++; + } + // Jet after top quark decay but not directly from top + else if(fromTopDecay) { + if(bJetAfterTopIds.count(jetIndex) < 1) bJetAfterTopIds[jetIndex] = 1; + else bJetAfterTopIds[jetIndex]++; + } + } + + // Counting number of specific hadrons in each c jet + for(size_t hadronId = 0; hadronId < genCHadJetIndex->size(); ++hadronId) { + // Skipping c hadrons that are coming from b hadrons + if(genCHadBHadronId->at(hadronId) >= 0) continue; + // Index of a jet associated to the hadron + const int jetIndex = genCHadJetIndex->at(hadronId); + // Whether hadron radiated before top quark decay + const bool fromTopDecay = genCHadFromTopWeakDecay->at(hadronId); + // Skipping hadrons which have no associated jet + if(jetIndex < 0) continue; + // Skipping if jet is not in acceptance + if(genJets->at(jetIndex).pt() < genJetPtMin_) continue; + if(std::fabs(genJets->at(jetIndex).eta()) > genJetAbsEtaMax_) continue; + // Jet before top quark decay + if(!fromTopDecay) { + if(cJetBeforeTopIds.count(jetIndex) < 1) cJetBeforeTopIds[jetIndex] = 1; + else cJetBeforeTopIds[jetIndex]++; + } + // Jet after top quark decay but not directly from top + else if(fromTopDecay) { + if(cJetAfterTopIds.count(jetIndex) < 1) cJetAfterTopIds[jetIndex] = 1; + else cJetAfterTopIds[jetIndex]++; + } + } + + // Finding additional b jets (before top decay) + std::vector additionalBJetIds; + for(std::map::iterator it = bJetBeforeTopIds.begin(); it != bJetBeforeTopIds.end(); ++it) { + const int jetId = it->first; + // Skipping the jet if it contains a b hadron directly from top quark decay + if(bJetFromTopIds.count(jetId) > 0) continue; + additionalBJetIds.push_back(jetId); + } + // Finding pseudo-additional b jets (after top decay) + std::vector pseudoadditionalBJetIds; + for(std::map::iterator it = bJetAfterTopIds.begin(); it != bJetAfterTopIds.end(); ++it) { + const int jetId = it->first; + // Skipping the jet if it contains a b hadron directly from top quark decay + if(bJetAfterTopIds.count(jetId) > 0) continue; + pseudoadditionalBJetIds.push_back(jetId); + } + // Finding additional c jets + std::vector additionalCJetIds; + for(std::map::iterator it = cJetBeforeTopIds.begin(); it != cJetBeforeTopIds.end(); ++it) { + const int jetId = it->first; + additionalCJetIds.push_back(jetId); + } + // Finding pseudo-additional c jets (after top decay) + std::vector pseudoadditionalCJetIds; + for(std::map::iterator it = cJetAfterTopIds.begin(); it != cJetAfterTopIds.end(); ++it) { + const int jetId = it->first; + // Skipping the jet if it contains a b hadron directly from top quark decay + if(cJetAfterTopIds.count(jetId) > 0) continue; + pseudoadditionalCJetIds.push_back(jetId); + } + + // Categorizing event based on number of additional b/c jets + // and number of corresponding hadrons in each of them + additionalJetEventId_ = bJetFromTopIds.size()*100; + // tt + 1 additional b jet + if (additionalBJetIds.size() == 1) { + int nHadronsInJet = bJetBeforeTopIds[additionalBJetIds.at(0)]; + // 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 (additionalBJetIds.size() > 1) { + int nHadronsInJet1 = bJetBeforeTopIds[additionalBJetIds.at(0)]; + int nHadronsInJet2 = bJetBeforeTopIds[additionalBJetIds.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 overlapping 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::max(nHadronsInJet1, nHadronsInJet2) == 1) additionalJetEventId_ = 55; + } + // tt + no additional b jets + else if(additionalBJetIds.size() == 0) { + // tt + >=1 pseudo-additional b jet with b hadrons after top quark decay + if(pseudoadditionalBJetIds.size() > 0) additionalJetEventId_ = 56; + // tt + 1 additional c jet + else if(additionalCJetIds.size() == 1) { + int nHadronsInJet = cJetBeforeTopIds[additionalCJetIds.at(0)]; + // tt + 1 additional c jet from 1 additional c hadron + if(nHadronsInJet == 1) additionalJetEventId_ = 41; + // tt + 1 additional c jet from >=2 overlapping additional c hadrons + else additionalJetEventId_ = 42; + } + // tt + >=2 additional c jets + else if(additionalCJetIds.size() > 1) { + int nHadronsInJet1 = cJetBeforeTopIds[additionalCJetIds.at(0)]; + int nHadronsInJet2 = cJetBeforeTopIds[additionalCJetIds.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 overlapping 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::max(nHadronsInJet1, nHadronsInJet2) == 1) additionalJetEventId_ = 45; + } + // tt + no additional c jets + else if(additionalCJetIds.size() == 0) { + // tt + >=1 pseudo-additional c jet with c hadrons after top quark decay + if(pseudoadditionalCJetIds.size() > 0) additionalJetEventId_ = 46; + // tt + light jets + else additionalJetEventId_ = 0; + } + } + + + // Performance information filled in a tree with entry for each b hadron ####################################### + // Looping over all b hadrons + for(size_t hadronId = 0; hadronId < genBHadIndex->size(); ++hadronId) { + clearBHadronVariables(); + + const int hadronParticleId = genBHadIndex->at(hadronId); + if(hadronParticleId < 0) continue; + + // Calculating hadron/jet pt ratio + const int hadronJetId = genBHadJetIndex->at(hadronId); + if(hadronJetId >= 0) { + bHadron_jetPtRatio_ = genBHadPlusMothers->at(genBHadIndex->at(hadronId)).pt() / genJets->at(hadronJetId).pt(); + } + + const int hadronParticlePdgId = genBHadPlusMothers->at(hadronParticleId).pdgId(); + bHadron_flavour_ = genBHadFlavour->at(hadronId); + // Finding all b quark candidates + std::set hadronLastQuarks; + findAncestorParticles(hadronParticleId, hadronLastQuarks, *genBHadPlusMothers, *genBHadPlusMothersIndices, 5, true, 1); + bHadron_nQuarksAll_ = hadronLastQuarks.size(); + // Identifying flavour of the proper candidate quark + int quarkCandidateFlavourSign = hadronParticlePdgId > 0 ? 1 : -1; + // Inverting flavour if this is a b meson + if(hadronParticlePdgId / 1000 == 0) quarkCandidateFlavourSign *= -1; + // Selecting candidates that have proper flavour sign + std::vector > hadronLastQuarks_sameFlavour; + for(std::set::iterator it=hadronLastQuarks.begin(); it!=hadronLastQuarks.end(); ++it) { + const int quarkParticleId = *it; + if(quarkParticleId < 0) continue; + const int quarkParticlePdgId = genBHadPlusMothers->at(quarkParticleId).pdgId(); + // Skipping quarks that have opposite flavour sign + if(quarkParticlePdgId * quarkCandidateFlavourSign < 0) continue; + double hadron_quark_dR = ROOT::Math::VectorUtil::DeltaR(genBHadPlusMothers->at(hadronParticleId).polarP4(), genBHadPlusMothers->at(quarkParticleId).polarP4()); + hadronLastQuarks_sameFlavour.push_back( std::pair(quarkParticleId, hadron_quark_dR) ); + } + bHadron_nQuarksSameFlavour_ = hadronLastQuarks_sameFlavour.size(); + + // Sorting quarks with proper flavour by their dR + std::sort(hadronLastQuarks_sameFlavour.begin(), hadronLastQuarks_sameFlavour.end(), sort_pairs_second); + const int hadronQuark1ParticleId = hadronLastQuarks_sameFlavour.size() > 0 ? hadronLastQuarks_sameFlavour.at(0).first : -1; + const int hadronQuark2ParticleId = hadronLastQuarks_sameFlavour.size() > 1 ? hadronLastQuarks_sameFlavour.at(1).first : -1; + + if(hadronQuark1ParticleId >= 0) { + bHadron_quark1_dR_ = ROOT::Math::VectorUtil::DeltaR(genBHadPlusMothers->at(hadronParticleId).polarP4(), genBHadPlusMothers->at(hadronQuark1ParticleId).polarP4()); + bHadron_quark1_ptRatio_ = genBHadPlusMothers->at(hadronParticleId).pt() / genBHadPlusMothers->at(hadronQuark1ParticleId).pt(); + } + + if(hadronQuark2ParticleId >= 0) { + bHadron_quark2_dR_ = ROOT::Math::VectorUtil::DeltaR(genBHadPlusMothers->at(hadronParticleId).polarP4(), genBHadPlusMothers->at(hadronQuark2ParticleId).polarP4()); + bHadron_quark2_ptRatio_ = genBHadPlusMothers->at(hadronParticleId).pt() / genBHadPlusMothers->at(hadronQuark2ParticleId).pt(); + } + + // Counting number of leptons coming from the b hadron decay + for(size_t leptonId = 0; leptonId < genBHadLeptonHadronIndex->size(); ++leptonId) { + const size_t leptonHadronId = genBHadLeptonHadronIndex->at(leptonId); + // Skipping the lepton id doesn't come from the current b hadron + if(leptonHadronId != hadronId) continue; + bHadron_nLeptonsAll_++; + if(genBHadLeptonViaTau->at(leptonId) == 1) bHadron_nLeptonsViaTau_++; + } + + + tree_bHadrons_->Fill(); + } + + tree_events_->Fill(); +} + + +void matchGenHFHadrons::findAncestorParticles(const int startId, std::set& resultIds, + const std::vector& genParticles, const std::vector >& motherIndices, + const int endPdgId, const bool endPdgIdIsAbs, const int firstAnyLast) +{ + if(startId < 0) return; + int startPdgId = genParticles.at(startId).pdgId(); + if(endPdgIdIsAbs) startPdgId = std::abs(startPdgId); + + // No mothers stored for the requested particle + if((int)motherIndices.size() < startId+1) return; + + // Getting ids of mothers of the starting particle + const std::vector& motherIds = motherIndices.at(startId); + + // In case the starting particle has pdgId as requested and it should not have same mothers + if(startPdgId == endPdgId && firstAnyLast == 1) { + int nSamePdgIdMothers = 0; + // Counting how many mothers with the same pdgId the particle has + for(size_t motherId = 0; motherId < motherIds.size(); ++motherId) { + const int motherParticleId = motherIds.at(motherId); + if(motherParticleId < 0) continue; + int motherParticlePdgId = genParticles.at(motherParticleId).pdgId(); + if(endPdgIdIsAbs) motherParticlePdgId = std::abs(motherParticlePdgId); + if(motherParticlePdgId == startPdgId) nSamePdgIdMothers++; + } + if(nSamePdgIdMothers < 1) { + resultIds.insert(startId); + return; + } + } + + // Looping over all mothers of the particle + for(size_t motherId = 0; motherId < motherIds.size(); ++motherId) { + const int motherParticleId = motherIds.at(motherId); + if(motherParticleId < 0) continue; + int motherParticlePdgId = genParticles.at(motherParticleId).pdgId(); + if(endPdgIdIsAbs) motherParticlePdgId = std::abs(motherParticlePdgId); + // Checking whether pdgId of this particle matches the requested one + bool isCandidate = false; + bool inProperPlace = false; + if(motherParticlePdgId == endPdgId) isCandidate = true; + // Checking whether the candidate particle is first/last according to the request + if(isCandidate) { + if(firstAnyLast == -1 && startPdgId != motherParticlePdgId) inProperPlace = true; + else if(firstAnyLast == 0) inProperPlace = true; + } + // Adding the mother if it has proper pdgId and place in chain + if(isCandidate && inProperPlace) { + resultIds.insert(motherParticleId); + if(firstAnyLast < 0) continue; + } + + // Checking mothers of this mother + findAncestorParticles(motherParticleId, resultIds, genParticles, motherIndices, endPdgId, endPdgIdIsAbs, firstAnyLast); + } + + +} + + +void matchGenHFHadrons::beginJob() +{ + edm::Service fileService; + if(!fileService) throw edm::Exception(edm::errors::Configuration, "TFileService is not registered in cfg file"); + tree_events_ = fileService->make("tree_events", "tree_events"); + tree_bHadrons_ = fileService->make("tree_bHadrons", "tree_bHadrons"); + + // Creating output branches + tree_events_->Branch("nJets", &nJets_); + + tree_events_->Branch("nBJets", &nBJets_); + tree_events_->Branch("nBHadrons", &nBHadrons_); + tree_events_->Branch("nBHadronsTop", &nBHadronsTop_); + tree_events_->Branch("nBHadronsAdditional", &nBHadronsAdditional_); + tree_events_->Branch("nBHadronsPseudoadditional", &nBHadronsPseudoadditional_); + + tree_events_->Branch("nCJets", &nCJets_); + tree_events_->Branch("nCHadrons", &nCHadrons_); + tree_events_->Branch("nCHadronsAdditional", &nCHadronsAdditional_); + tree_events_->Branch("nCHadronsPseudoadditional", &nCHadronsPseudoadditional_); + + tree_events_->Branch("nBHadronsHiggs", &nBHadronsHiggs_); + tree_events_->Branch("additionalJetEventId", &additionalJetEventId_); + + tree_bHadrons_->Branch("bHadron_flavour", &bHadron_flavour_); + tree_bHadrons_->Branch("bHadron_nQuarksAll", &bHadron_nQuarksAll_); + tree_bHadrons_->Branch("bHadron_nQuarksSameFlavour", &bHadron_nQuarksSameFlavour_); + + tree_bHadrons_->Branch("bHadron_quark1_dR", &bHadron_quark1_dR_); + tree_bHadrons_->Branch("bHadron_quark2_dR", &bHadron_quark2_dR_); + tree_bHadrons_->Branch("bHadron_quark1_ptRatio", &bHadron_quark1_ptRatio_); + tree_bHadrons_->Branch("bHadron_quark2_ptRatio", &bHadron_quark2_ptRatio_); + + tree_bHadrons_->Branch("bHadron_nLeptonsAll", &bHadron_nLeptonsAll_); + tree_bHadrons_->Branch("bHadron_nLeptonsViaTau", &bHadron_nLeptonsViaTau_); + + tree_bHadrons_->Branch("bHadron_jetPtRatio", &bHadron_jetPtRatio_); +} + + +void matchGenHFHadrons::clearEventVariables() +{ + nJets_ = 0; + + nBJets_ = 0; + nBHadrons_ = 0; + nBHadronsTop_ = 0; + nBHadronsAdditional_ = 0; + nBHadronsPseudoadditional_ = 0; + + nCJets_ = 0; + nCHadrons_ = 0; + nCHadronsAdditional_ = 0; + nCHadronsPseudoadditional_ = 0; + + nBHadronsHiggs_ = 0; + additionalJetEventId_ = -1; +} + + +void matchGenHFHadrons::clearBHadronVariables() +{ + bHadron_flavour_ = 0; + bHadron_nQuarksAll_ = 0; + bHadron_nQuarksSameFlavour_ = 0; + + bHadron_quark1_dR_ = -0.1; + bHadron_quark2_dR_ = -0.1; + bHadron_quark1_ptRatio_ = -0.1; + bHadron_quark2_ptRatio_ = -0.1; + + + bHadron_nLeptonsAll_ = 0; + bHadron_nLeptonsViaTau_ = 0; + + bHadron_jetPtRatio_ = -0.1; +} + + + + +DEFINE_FWK_MODULE( matchGenHFHadrons ); diff --git a/PhysicsTools/JetMCAlgos/test/matchGenHFHadrons.py b/PhysicsTools/JetMCAlgos/test/matchGenHFHadrons.py new file mode 100644 index 0000000000000..7598b208af78b --- /dev/null +++ b/PhysicsTools/JetMCAlgos/test/matchGenHFHadrons.py @@ -0,0 +1,156 @@ +import FWCore.ParameterSet.Config as cms +import FWCore.ParameterSet.VarParsing as VarParsing +import sys + +process = cms.Process("Analyzer") + +## defining command line options +options = VarParsing.VarParsing ('standard') +options.register('runOnAOD', True, VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.bool, "run on AOD") + +## processing provided 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]) + +## enabling unscheduled mode for modules +process.options = cms.untracked.PSet( + wantSummary = cms.untracked.bool(True), + allowUnscheduled = cms.untracked.bool(True), +) + +## 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 (1000+ events in each file defined below) + '/store/relval/CMSSW_7_2_0_pre7/RelValTTbar_13/GEN-SIM-RECO/PU50ns_PRE_LS172_V12-v1/00000/1267B7ED-2F4E-E411-A0B9-0025905964A6.root', + ## other AOD samples +# '/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 (1000+ events in each file defined below) + '/store/cmst3/user/gpetrucc/miniAOD/v1/TTbarH_HToBB_M-125_13TeV_pythia6_PU_S14_PAT.root', + ), + skipEvents = cms.untracked.uint32(0) + ) + +## define maximal number of events to loop over +process.maxEvents = cms.untracked.PSet( + input = cms.untracked.int32(1000) +) + +# Setting 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") + +# Producing own jets for testing purposes +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.sequences.GenHFHadronMatching_cff import genJetFlavourPlusLeptonInfos +process.genJetFlavourPlusLeptonInfos = genJetFlavourPlusLeptonInfos.clone( + jets = genJetCollection, + rParam = cms.double(0.4), + jetAlgorithm = cms.string("AntiKt") +) + + +# Plugin for analysing B hadrons +# MUST use the same particle collection as in selectedHadronsAndPartons +from PhysicsTools.JetMCAlgos.sequences.GenHFHadronMatching_cff import matchGenBHadron +process.matchGenBHadron = matchGenBHadron.clone( + genParticles = genParticleCollection +) + +# Plugin for analysing C hadrons +# MUST use the same particle collection as in selectedHadronsAndPartons +from PhysicsTools.JetMCAlgos.sequences.GenHFHadronMatching_cff import matchGenCHadron +process.matchGenCHadron = matchGenCHadron.clone( + genParticles = genParticleCollection +) + + +## configuring the testing analyzer that produces output tree +process.matchGenHFHadrons = cms.EDAnalyzer("matchGenHFHadrons", + # phase space of jets to be stored + genJetPtMin = cms.double(20), + genJetAbsEtaMax = cms.double(2.4), + # input tags holding information about matching + genJets = cms.InputTag(genJetCollection), + 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"), +) + +## setting up output root file +process.TFileService = cms.Service("TFileService", + fileName = cms.string("matchGenHFHadrons_trees.root") +) + + + +## defining only the final modules to run: dependencies will be run automatically [allowUnscheduled = True] +process.p1 = cms.Path( + process.matchGenHFHadrons +) + +## module to store raw output from the processed modules into the ROOT file +process.out = cms.OutputModule("PoolOutputModule", + fileName = cms.untracked.string('matchGenHFHadrons_out.root'), + outputCommands = cms.untracked.vstring('drop *', 'keep *_matchGen*_*_*') + ) +process.outpath = cms.EndPath(process.out)