diff --git a/DeepNtuplizer/interface/helpers.h b/DeepNtuplizer/interface/helpers.h index 41f9860e16858..1965b41cbd00e 100644 --- a/DeepNtuplizer/interface/helpers.h +++ b/DeepNtuplizer/interface/helpers.h @@ -7,8 +7,8 @@ #include "DataFormats/PatCandidates/interface/Jet.h" namespace deep_ntuples { - enum JetFlavor {UNDEFINED, G, L, C, B, BB}; - JetFlavor jet_flavour(const pat::Jet& jet); + enum JetFlavor {UNDEFINED, G, L, C, B, BB, LeptonicB, LeptonicB_C}; + JetFlavor jet_flavour(const pat::Jet& jet, std::vector neutrinosLepB, std::vector neutrinosLepB_C); } #endif //DEEPNTUPLES_DEEPNTUPLIZER_INTERFACE_HELPERS_H_ diff --git a/DeepNtuplizer/interface/ntuple_JetInfo.h b/DeepNtuplizer/interface/ntuple_JetInfo.h index 664de3d50f077..9693be984257d 100644 --- a/DeepNtuplizer/interface/ntuple_JetInfo.h +++ b/DeepNtuplizer/interface/ntuple_JetInfo.h @@ -56,7 +56,10 @@ class ntuple_JetInfo: public ntuple_content{ genJetMatchWithNuToken_ = genJetMatchWithNuToken; } - + void setGenParticlesToken( + edm::EDGetTokenT genParticlesToken) { + genParticlesToken_ = genParticlesToken; + } //private: @@ -80,12 +83,19 @@ class ntuple_JetInfo: public ntuple_content{ edm::EDGetTokenT > genJetMatchReclusterToken_; edm::EDGetTokenT > genJetMatchWithNuToken_; + edm::EDGetTokenT genParticlesToken_; + edm::Handle > genJetMatchRecluster; edm::Handle > genJetMatchWithNu; + edm::Handle genParticlesHandle; + TRandom3 TRandom_; float gluonReduction_; + std::vector < reco::GenParticle> neutrinosLepB; + std::vector < reco::GenParticle> neutrinosLepB_C; + // labels (MC truth) // regressions pt, Deta, Dphi float gen_pt_; @@ -95,6 +105,8 @@ class ntuple_JetInfo: public ntuple_content{ int isC_; int isUDS_; int isG_; + int isLeptonicB_; + int isLeptonicB_C_; // global variables float npv_; diff --git a/DeepNtuplizer/plugins/DeepNtuplizer.cc b/DeepNtuplizer/plugins/DeepNtuplizer.cc index badc1fc340bc0..88636fd2325b3 100644 --- a/DeepNtuplizer/plugins/DeepNtuplizer.cc +++ b/DeepNtuplizer/plugins/DeepNtuplizer.cc @@ -120,8 +120,12 @@ DeepNtuplizer::DeepNtuplizer(const edm::ParameterSet& iConfig): jetinfo->setGenJetMatchWithNuToken( consumes >( iConfig.getParameter( "genJetMatchWithNu" ))); - addModule(jetinfo); + jetinfo->setGenParticlesToken( + consumes( + iConfig.getParameter("pruned"))); + + addModule(jetinfo); ntuple_pfCands * pfcands = new ntuple_pfCands(); addModule(pfcands); diff --git a/DeepNtuplizer/python/DeepNtuplizer_cfi.py b/DeepNtuplizer/python/DeepNtuplizer_cfi.py index 125d30e544f9d..30ea51aada3ca 100644 --- a/DeepNtuplizer/python/DeepNtuplizer_cfi.py +++ b/DeepNtuplizer/python/DeepNtuplizer_cfi.py @@ -7,6 +7,7 @@ SVs = cms.InputTag("slimmedSecondaryVertices"), genJetMatchWithNu = cms.InputTag("patGenJetMatchWithNu"), genJetMatchRecluster = cms.InputTag("patGenJetMatchRecluster"), + pruned = cms.InputTag("prunedGenParticles"), jetPtMin = cms.double(20.0), jetPtMax = cms.double(100), jetAbsEtaMin = cms.double(0.0), diff --git a/DeepNtuplizer/src/helpers.cc b/DeepNtuplizer/src/helpers.cc index 53f3b313647ce..9e7bd2c60b9d5 100644 --- a/DeepNtuplizer/src/helpers.cc +++ b/DeepNtuplizer/src/helpers.cc @@ -1,17 +1,29 @@ #include "../interface/helpers.h" namespace deep_ntuples { - JetFlavor jet_flavour(const pat::Jet& jet) { + JetFlavor jet_flavour(const pat::Jet& jet, std::vector neutrinosLepB, std::vector neutrinosLepB_C) { int hflav = abs(jet.hadronFlavour()); int pflav = abs(jet.partonFlavour()); size_t nbs = jet.jetFlavourInfo().getbHadrons().size(); size_t ncs = jet.jetFlavourInfo().getcHadrons().size(); if(hflav == 5) { //B jet - if(nbs > 1) return JetFlavor::BB; - else if(nbs == 1) return JetFlavor::B; - else return JetFlavor::UNDEFINED; - } + if(nbs > 1) return JetFlavor::BB; + else if(nbs == 1) { + for (std::vector::iterator it = neutrinosLepB.begin(); it != neutrinosLepB.end(); ++it){ + if(reco::deltaR(it->eta(),it->phi(),jet.eta(),jet.phi()) < 0.4) { + return JetFlavor::LeptonicB; + } + } + for (std::vector::iterator it = neutrinosLepB_C.begin(); it != neutrinosLepB_C.end(); ++it){ + if(reco::deltaR(it->eta(),it->phi(),jet.eta(),jet.phi()) < 0.4) { + return JetFlavor::LeptonicB_C; + } + } + return JetFlavor::B; + } + else return JetFlavor::UNDEFINED; + } else if(hflav == 4) { //C jet return JetFlavor::C; } diff --git a/DeepNtuplizer/src/ntuple_JetInfo.cc b/DeepNtuplizer/src/ntuple_JetInfo.cc index 06185f6a892ea..59cc5d35f36e1 100644 --- a/DeepNtuplizer/src/ntuple_JetInfo.cc +++ b/DeepNtuplizer/src/ntuple_JetInfo.cc @@ -16,152 +16,179 @@ using namespace std; void ntuple_JetInfo::getInput(const edm::ParameterSet& iConfig){ - gluonReduction_=(iConfig.getParameter("gluonReduction")); - jetPtMin_=(iConfig.getParameter("jetPtMin")); - jetPtMax_=(iConfig.getParameter("jetPtMax")); - jetAbsEtaMin_=(iConfig.getParameter("jetAbsEtaMin")); - jetAbsEtaMax_=(iConfig.getParameter("jetAbsEtaMax")); - - vector disc_names = iConfig.getParameter >("bDiscriminators"); - for(auto& name : disc_names) { - discriminators_[name] = 0.; - } + gluonReduction_=(iConfig.getParameter("gluonReduction")); + jetPtMin_=(iConfig.getParameter("jetPtMin")); + jetPtMax_=(iConfig.getParameter("jetPtMax")); + jetAbsEtaMin_=(iConfig.getParameter("jetAbsEtaMin")); + jetAbsEtaMax_=(iConfig.getParameter("jetAbsEtaMax")); + + vector disc_names = iConfig.getParameter >("bDiscriminators"); + for(auto& name : disc_names) { + discriminators_[name] = 0.; + } } void ntuple_JetInfo::initBranches(TTree* tree){ - //more general event info, here applied per jet - addBranch(tree,"npv" ,&npv_ ,"npv/f" ); - addBranch(tree,"event_no" ,&event_no_ ,"npv/i" ); - addBranch(tree,"jet_no" ,&jet_no_ ,"npv/i" ); - - - // truthe labels - addBranch(tree,"gen_pt" ,&gen_pt_ ,"gen_pt_/f" ); - addBranch(tree,"Delta_gen_pt" ,&Delta_gen_pt_,"Delta_gen_pt_/f" ); - addBranch(tree,"isB",&isB_, "isB_/i"); - addBranch(tree,"isC",&isC_, "isC_/i"); - addBranch(tree,"isUDS",&isUDS_, "isUDS_/i"); - addBranch(tree,"isG",&isG_, "isG_/i"); - - - // jet variables - //b=tree->Branch("jet_pt", &jet_pt_); - addBranch(tree,"jet_pt", &jet_pt_); - - addBranch(tree,"jet_corr_pt", &jet_corr_pt_); - addBranch(tree,"jet_eta", &jet_eta_); - - // quark gluon - addBranch(tree,"jet_qgl", &jet_qgl_); // qg tagger from jmar - addBranch(tree,"QG_ptD", &QG_ptD_); // momentum fraction per jet constituent - addBranch(tree,"QG_axis2", &QG_axis2_); // jet shape i.e. gluon are wider than quarks - addBranch(tree,"QG_mult", &QG_mult_); // multiplicity i.e. total num of PFcands reconstructed - // in the jet - - - addBranch(tree,"gen_pt_Recluster" ,&gen_pt_Recluster_ ,"gen_pt_Recluster_/f" ); - addBranch(tree,"gen_pt_WithNu" ,&gen_pt_WithNu_ ,"gen_pt_WithNu_/f" ); - addBranch(tree,"Delta_gen_pt_Recluster" ,&Delta_gen_pt_Recluster_ ,"Delta_gen_pt_Recluster_/f" ); - addBranch(tree,"Delta_gen_pt_WithNu" ,&Delta_gen_pt_WithNu_ ,"Delta_gen_pt_WithNu_/f" ); - - if(1) // discriminators might need to be filled differently. FIXME - for(auto& entry : discriminators_) { - string better_name(entry.first); - std::replace(better_name.begin(), better_name.end(), ':', '_'); - addBranch(tree,better_name.c_str(), &entry.second, (better_name+"/F").c_str()); - } + //more general event info, here applied per jet + addBranch(tree,"npv" ,&npv_ ,"npv/f" ); + addBranch(tree,"event_no" ,&event_no_ ,"event_no/i" ); + addBranch(tree,"jet_no" ,&jet_no_ ,"jet_no/i" ); + + + // truthe labels + addBranch(tree,"gen_pt" ,&gen_pt_ ,"gen_pt_/f" ); + addBranch(tree,"Delta_gen_pt" ,&Delta_gen_pt_,"Delta_gen_pt_/f" ); + addBranch(tree,"isB",&isB_, "isB_/i"); + addBranch(tree,"isC",&isC_, "isC_/i"); + addBranch(tree,"isUDS",&isUDS_, "isUDS_/i"); + addBranch(tree,"isG",&isG_, "isG_/i"); + addBranch(tree,"isLeptonicB",&isLeptonicB_, "isLeptonicB_/i"); + addBranch(tree,"isLeptonicB_C",&isLeptonicB_C_, "isLeptonicB_C_/i"); + + // jet variables + //b=tree->Branch("jet_pt", &jet_pt_); + addBranch(tree,"jet_pt", &jet_pt_); + + addBranch(tree,"jet_corr_pt", &jet_corr_pt_); + addBranch(tree,"jet_eta", &jet_eta_); + + // quark gluon + addBranch(tree,"jet_qgl", &jet_qgl_); // qg tagger from jmar + addBranch(tree,"QG_ptD", &QG_ptD_); // momentum fraction per jet constituent + addBranch(tree,"QG_axis2", &QG_axis2_); // jet shape i.e. gluon are wider than quarks + addBranch(tree,"QG_mult", &QG_mult_); // multiplicity i.e. total num of PFcands reconstructed + // in the jet + + + addBranch(tree,"gen_pt_Recluster" ,&gen_pt_Recluster_ ,"gen_pt_Recluster_/f" ); + addBranch(tree,"gen_pt_WithNu" ,&gen_pt_WithNu_ ,"gen_pt_WithNu_/f" ); + addBranch(tree,"Delta_gen_pt_Recluster" ,&Delta_gen_pt_Recluster_ ,"Delta_gen_pt_Recluster_/f" ); + addBranch(tree,"Delta_gen_pt_WithNu" ,&Delta_gen_pt_WithNu_ ,"Delta_gen_pt_WithNu_/f" ); + + if(1) // discriminators might need to be filled differently. FIXME + for(auto& entry : discriminators_) { + string better_name(entry.first); + std::replace(better_name.begin(), better_name.end(), ':', '_'); + addBranch(tree,better_name.c_str(), &entry.second, (better_name+"/F").c_str()); + } } void ntuple_JetInfo::readEvent(const edm::Event& iEvent){ - iEvent.getByToken(qglToken_, qglHandle); - iEvent.getByToken(ptDToken_, ptDHandle); - iEvent.getByToken(axis2Token_, axis2Handle); - iEvent.getByToken(multToken_, multHandle); - - - iEvent.getByToken(genJetMatchReclusterToken_, genJetMatchRecluster); - iEvent.getByToken(genJetMatchWithNuToken_, genJetMatchWithNu); - - //technically a branch fill but per event, therefore here - event_no_=iEvent.id().event(); + iEvent.getByToken(qglToken_, qglHandle); + iEvent.getByToken(ptDToken_, ptDHandle); + iEvent.getByToken(axis2Token_, axis2Handle); + iEvent.getByToken(multToken_, multHandle); + + + iEvent.getByToken(genJetMatchReclusterToken_, genJetMatchRecluster); + iEvent.getByToken(genJetMatchWithNuToken_, genJetMatchWithNu); + + iEvent.getByToken(genParticlesToken_, genParticlesHandle); + + neutrinosLepB.clear(); + neutrinosLepB_C.clear(); + + for (const reco::Candidate &genC : *genParticlesHandle) { + const reco::GenParticle &gen = static_cast< const reco::GenParticle &>(genC); + if(abs(gen.pdgId())==12||abs(gen.pdgId())==14||abs(gen.pdgId())==16) { + const reco::GenParticle* mother = static_cast< const reco::GenParticle*> (gen.mother()); + if(mother!=NULL) { + if((abs(mother->pdgId())>500&&abs(mother->pdgId())<600)||(abs(mother->pdgId())>5000&&abs(mother->pdgId())<6000)) { + neutrinosLepB.emplace_back(gen); + } + if((abs(mother->pdgId())>400&&abs(mother->pdgId())<500)||(abs(mother->pdgId())>4000&&abs(mother->pdgId())<5000)) { + neutrinosLepB_C.emplace_back(gen); + } + } + else { + std::cout << "No mother" << std::endl; + } + } + } + //technically a branch fill but per event, therefore here + event_no_=iEvent.id().event(); } //use either of these functions bool ntuple_JetInfo::fillBranches(const pat::Jet & jet, const size_t& jetidx, const edm::View * coll){ - if(!coll) - throw std::runtime_error("ntuple_JetInfo::fillBranches: no jet collection"); + if(!coll) + throw std::runtime_error("ntuple_JetInfo::fillBranches: no jet collection"); + + /// cuts /// - /// cuts /// + // some cuts to contrin training region + if ( jet.pt() < jetPtMin_ || jet.pt() > jetPtMax_ ) return false; // apply jet pT cut + if ( jet.eta() < fabs(jetAbsEtaMin_) ||jet.eta() > fabs(jetAbsEtaMax_) ) return false; // apply jet eta cut - // some cuts to contrin training region - if ( jet.pt() < jetPtMin_ || jet.pt() > jetPtMax_ ) return false; // apply jet pT cut - if ( jet.eta() < fabs(jetAbsEtaMin_) ||jet.eta() > fabs(jetAbsEtaMax_) ) return false; // apply jet eta cut + // often we have way to many gluons that we do not need. This randomply reduces the gluons + if (gluonReduction_>0 && jet.partonFlavour()==21) + if(TRandom_.Uniform()>gluonReduction_) return false; - // often we have way to many gluons that we do not need. This randomply reduces the gluons - if (gluonReduction_>0 && jet.partonFlavour()==21) - if(TRandom_.Uniform()>gluonReduction_) return false; + if(jet.genJet()==NULL)return false; - if(jet.genJet()==NULL)return false; //branch fills - for(auto& entry : discriminators_) { - entry.second = catchInfs(jet.bDiscriminator(entry.first),-0.1); - } + for(auto& entry : discriminators_) { + entry.second = catchInfs(jet.bDiscriminator(entry.first),-0.1); + } + - npv_ = vertices()->size(); - jet_no_=jetidx; + npv_ = vertices()->size(); + jet_no_=jetidx; - const auto jetRef = reco::CandidatePtr(coll->ptrs().at( jetidx)); - jet_qgl_ = (*qglHandle)[jetRef]; - QG_ptD_ = (*ptDHandle)[jetRef]; - QG_axis2_ = (*axis2Handle)[jetRef]; - QG_mult_ = (*multHandle)[jetRef]; + const auto jetRef = reco::CandidatePtr(coll->ptrs().at( jetidx)); + jet_qgl_ = (*qglHandle)[jetRef]; + QG_ptD_ = (*ptDHandle)[jetRef]; + QG_axis2_ = (*axis2Handle)[jetRef]; + QG_mult_ = (*multHandle)[jetRef]; - //std::vector > p= coll->ptrs(); + //std::vector > p= coll->ptrs(); - isB_=0; isC_=0; isUDS_=0; isG_=0; - switch(deep_ntuples::jet_flavour(jet)) { - case deep_ntuples::JetFlavor::L: isUDS_=1; break; - case deep_ntuples::JetFlavor::BB: - case deep_ntuples::JetFlavor::B: isB_=1; break; - case deep_ntuples::JetFlavor::C: isC_=1; break; - case deep_ntuples::JetFlavor::G: isG_=1; break; - default : break; - } + isB_=0; isC_=0; isUDS_=0; isG_=0, isLeptonicB_=0, isLeptonicB_C_=0; + switch(deep_ntuples::jet_flavour(jet, neutrinosLepB, neutrinosLepB_C)) { + case deep_ntuples::JetFlavor::L: isUDS_=1; break; + case deep_ntuples::JetFlavor::BB: + case deep_ntuples::JetFlavor::B: isB_=1; break; + case deep_ntuples::JetFlavor::C: isC_=1; break; + case deep_ntuples::JetFlavor::G: isG_=1; break; + case deep_ntuples::JetFlavor::LeptonicB: isLeptonicB_=1; break; + case deep_ntuples::JetFlavor::LeptonicB_C: isLeptonicB_C_=1; break; + default : break; + } - if(!isB_ && !isC_ && !isUDS_ && !isG_) return false; + if(!isB_ && !isC_ && !isUDS_ && !isG_ && !isLeptonicB_ && !isLeptonicB_C_) return false; - pat::JetCollection h; + pat::JetCollection h; - jet_pt_ = jet.correctedJet("Uncorrected").pt(); - jet_eta_ = jet.eta(); - jet_corr_pt_ = jet.pt(); + jet_pt_ = jet.correctedJet("Uncorrected").pt(); + jet_eta_ = jet.eta(); + jet_corr_pt_ = jet.pt(); - gen_pt_ = jet.genJet()->pt(); - Delta_gen_pt_ = jet.genJet()->pt()- jet_pt_; + gen_pt_ = jet.genJet()->pt(); + Delta_gen_pt_ = jet.genJet()->pt()- jet_pt_; - const edm::RefToBase patJetRef = coll->refAt(jetidx); - reco::GenJetRef genjetRecluster = (*genJetMatchRecluster)[patJetRef]; + const edm::RefToBase patJetRef = coll->refAt(jetidx); + reco::GenJetRef genjetRecluster = (*genJetMatchRecluster)[patJetRef]; - gen_pt_Recluster_ = 0.; - if (genjetRecluster.isNonnull() && genjetRecluster.isAvailable()) { - gen_pt_Recluster_ = genjetRecluster->pt(); - } - reco::GenJetRef genjetWithNu = (*genJetMatchWithNu)[patJetRef]; + gen_pt_Recluster_ = 0.; + if (genjetRecluster.isNonnull() && genjetRecluster.isAvailable()) { + gen_pt_Recluster_ = genjetRecluster->pt(); + } + reco::GenJetRef genjetWithNu = (*genJetMatchWithNu)[patJetRef]; - gen_pt_WithNu_ = 0.; - if (genjetWithNu.isNonnull() && genjetWithNu.isAvailable()) { - gen_pt_WithNu_ = genjetWithNu->pt(); - } + gen_pt_WithNu_ = 0.; + if (genjetWithNu.isNonnull() && genjetWithNu.isAvailable()) { + gen_pt_WithNu_ = genjetWithNu->pt(); + } - Delta_gen_pt_Recluster_=gen_pt_Recluster_-jet.pt(); - Delta_gen_pt_WithNu_=gen_pt_WithNu_-jet.pt(); + Delta_gen_pt_Recluster_=gen_pt_Recluster_-jet.pt(); + Delta_gen_pt_WithNu_=gen_pt_WithNu_-jet.pt(); - return true; + return true; } diff --git a/DeepNtuplizer/test/run_DeepNtuplizer.py b/DeepNtuplizer/test/run_DeepNtuplizer.py index 065a38dea06a3..be01264a4a16f 100644 --- a/DeepNtuplizer/test/run_DeepNtuplizer.py +++ b/DeepNtuplizer/test/run_DeepNtuplizer.py @@ -59,7 +59,7 @@ if options.nJobs > 1: print ("running over these files:") print (process.source.fileNames) -process.source.fileNames = ['file:/uscms/home/verzetti/nobackup/CMSSW_8_0_25/src/DeepNTuples/copy_numEvent100.root'] +#process.source.fileNames = ['file:/uscms/home/verzetti/nobackup/CMSSW_8_0_25/src/DeepNTuples/copy_numEvent100.root'] process.source.skipEvents = cms.untracked.uint32(options.skipEvents) process.maxEvents = cms.untracked.PSet( @@ -73,13 +73,18 @@ 'deepNNTagInfos', ] bTagDiscriminators = [ - 'pfCombinedInclusiveSecondaryVertexV2BJetTags', + 'softPFMuonBJetTags', + 'softPFElectronBJetTags', + 'pfJetBProbabilityBJetTags', + 'pfJetProbabilityBJetTags', + 'pfCombinedInclusiveSecondaryVertexV2BJetTags', 'deepFlavourJetTags:probudsg', #to be fixed with new names 'deepFlavourJetTags:probb', 'deepFlavourJetTags:probc', 'deepFlavourJetTags:probbb', 'deepFlavourJetTags:probcc', ] + jetCorrectionsAK4 = ('AK4PFchs', ['L1FastJet', 'L2Relative', 'L3Absolute'], 'None') from PhysicsTools.PatAlgos.tools.jetTools import updateJetCollection