diff --git a/TreeFillers/BuildFile.xml b/TreeFillers/BuildFile.xml index 40fda67..16b67fe 100644 --- a/TreeFillers/BuildFile.xml +++ b/TreeFillers/BuildFile.xml @@ -8,6 +8,8 @@ + + diff --git a/TreeFillers/interface/FatJetFiller.h b/TreeFillers/interface/FatJetFiller.h index f41247d..72e6347 100644 --- a/TreeFillers/interface/FatJetFiller.h +++ b/TreeFillers/interface/FatJetFiller.h @@ -23,10 +23,12 @@ class FatJetFiller : public BaseFiller { size i_eta = 0; size i_phi = 0; size i_mass = 0; + size i_toRawFact = 0; size i_csv = 0; size i_id = 0; size i_hadronFlavor = 0; size i_partonFlavor = 0; + size i_JECUnc = 0; size i_genIDX = 0; size i_gen_pt = 0; size i_gen_eta = 0; @@ -45,6 +47,7 @@ class FatJetFiller : public BaseFiller { size i_sj1_raw_pt = 0; size i_sj1_raw_mass = 0; size i_sj1_csv = 0; + size i_sj1_JECUnc = 0; size i_sj1_hadronFlavor = 0; size i_sj1_partonFlavor = 0; @@ -54,12 +57,16 @@ class FatJetFiller : public BaseFiller { size i_sj2_mass = 0; size i_sj2_raw_pt = 0; size i_sj2_raw_mass = 0; + size i_sj2_JECUnc = 0; size i_sj2_csv = 0; size i_sj2_hadronFlavor = 0; size i_sj2_partonFlavor = 0; bool isRealData = false; bool fillGenJets = false; + std::string jetType = ""; + std::string subjetType = ""; + std::string jetDef ; std::string subjetDef; edm::EDGetTokenT > token_jets; @@ -68,6 +75,12 @@ class FatJetFiller : public BaseFiller { edm::Handle > han_jets; edm::Handle han_genJets; + + edm::ESHandle jetCorParameters; + std::unique_ptr jetCorUnc; + + edm::ESHandle subjetCorParameters; + std::unique_ptr subjetCorUnc; }; } diff --git a/TreeFillers/interface/JetFiller.h b/TreeFillers/interface/JetFiller.h index 01a043f..26a7f3d 100644 --- a/TreeFillers/interface/JetFiller.h +++ b/TreeFillers/interface/JetFiller.h @@ -5,6 +5,9 @@ #include "DataFormats/PatCandidates/interface/Jet.h" #include "DataFormats/JetReco/interface/GenJet.h" #include "DataFormats/JetReco/interface/GenJetCollection.h" +#include "CondFormats/JetMETObjects/interface/JetCorrectionUncertainty.h" +#include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h" +#include "FWCore/Framework/interface/ESHandle.h" namespace AnaTM{ class JetFiller : public BaseFiller { @@ -25,10 +28,14 @@ class JetFiller : public BaseFiller { size i_eta = 0; size i_phi = 0; size i_mass = 0; + size i_toRawFact = 0; + size i_metUnc_rawPx = 0; + size i_metUnc_rawPy = 0; size i_csv = 0; size i_id = 0; size i_hadronFlavor = 0; size i_partonFlavor = 0; + size i_JECUnc = 0; size i_genIDX = 0; size i_gen_pt = 0; size i_gen_eta = 0; @@ -36,6 +43,7 @@ class JetFiller : public BaseFiller { size i_gen_mass = 0; bool isRealData = false; + std::string jetType = ""; bool fillGenJets = false; edm::EDGetTokenT > token_jets; edm::EDGetTokenT token_genJets; @@ -44,6 +52,9 @@ class JetFiller : public BaseFiller { edm::Handle > han_jets; edm::Handle han_genJets; + edm::ESHandle jetCorParameters; + std::unique_ptr jetCorUnc; + }; } diff --git a/TreeFillers/python/fatJetFiller_cff.py b/TreeFillers/python/fatJetFiller_cff.py index 1718124..a9d1046 100644 --- a/TreeFillers/python/fatJetFiller_cff.py +++ b/TreeFillers/python/fatJetFiller_cff.py @@ -9,6 +9,8 @@ genjets = cms.InputTag('ak8GenJetsNoNu'), jetDef = cms.string('AK8PuppiNoLep'), subjetDef = cms.string('SoftDrop'), + jetType = cms.string('AK8PFPuppi'), + subjetType = cms.string('AK4PFPuppi'), ) ak8PuppiFatJetFiller = cms.PSet( ignore = cms.bool(False), @@ -19,4 +21,6 @@ genjets = cms.InputTag('ak8GenJetsNoNu'), jetDef = cms.string('AK8Puppi'), subjetDef = cms.string('SoftDrop'), + jetType = cms.string('AK8PFPuppi'), + subjetType = cms.string('AK4PFPuppi'), ) \ No newline at end of file diff --git a/TreeFillers/python/jetFiller_cff.py b/TreeFillers/python/jetFiller_cff.py index c4f9a62..0eb3b7a 100644 --- a/TreeFillers/python/jetFiller_cff.py +++ b/TreeFillers/python/jetFiller_cff.py @@ -7,14 +7,17 @@ minJetPT = cms.double(10), jets = cms.InputTag('selectedUpdatedPatJetsAK4PFCHS'), genjets = cms.InputTag('slimmedGenJets'), + jetType = cms.string('AK4PFchs'), ) ak4PuppiJetFiller = ak4JetFiller.clone( branchName = cms.string("ak4PuppiJet"), fillGenJets = cms.bool(False), jets = cms.InputTag('selectedUpdatedPatJetsAK4PFPuppi'), + jetType = cms.string('AK4PFPuppi'), ) ak4PuppiNoLepJetFiller = ak4JetFiller.clone( branchName = cms.string("ak4PuppiNoLepJet"), fillGenJets = cms.bool(False), jets = cms.InputTag('selectedPatJetsAK4PFPuppiNoLep'), + jetType = cms.string('AK4PFPuppi'), ) \ No newline at end of file diff --git a/TreeFillers/src/FatJetFiller.cc b/TreeFillers/src/FatJetFiller.cc index 1553823..d8900af 100644 --- a/TreeFillers/src/FatJetFiller.cc +++ b/TreeFillers/src/FatJetFiller.cc @@ -2,6 +2,7 @@ #include "AnalysisTreeMaker/TreeFillers/interface/JetFiller.h" #include "AnalysisTreeMaker/TreeFillers/interface/FatJetFiller.h" #include "AnalysisTreeMaker/TreeFillers/interface/FillerConstants.h" +#include "JetMETCorrections/Objects/interface/JetCorrectionsRecord.h" using ASTypes::size8; using ASTypes::int8; @@ -13,6 +14,9 @@ FatJetFiller::FatJetFiller(const edm::ParameterSet& fullParamSet, const std::str isRealData(isRealData) { if(ignore()) return; + jetType =cfg.getParameter("jetType"); + subjetType =cfg.getParameter("subjetType"); + fillGenJets =cfg.getParameter("fillGenJets"); token_jets =cc.consumes >(cfg.getParameter("jets")); minJetPT = cfg.getParameter("minJetPT"); @@ -26,12 +30,14 @@ FatJetFiller::FatJetFiller(const edm::ParameterSet& fullParamSet, const std::str i_eta = data.addMulti(branchName,"eta" , 0); i_phi = data.addMulti(branchName,"phi" , 0); i_mass = data.addMulti(branchName,"mass" , 0); + i_toRawFact = data.addMulti(branchName,"toRawFact" , 0); i_csv = data.addMulti(branchName,"csv" , 0); i_id = data.addMulti(branchName,"id" , 0); if(!isRealData){ i_hadronFlavor = data.addMulti(branchName,"hadronFlavor" , 0); i_partonFlavor = data.addMulti(branchName,"partonFlavor" , 0); + i_JECUnc = data.addMulti(branchName,"JECUnc" , 0); if(fillGenJets){ i_genIDX = data.addMulti(branchName,"genIDX" , 0); i_gen_pt = data.addMulti(branchName,"gen_pt" , 0); @@ -53,6 +59,7 @@ FatJetFiller::FatJetFiller(const edm::ParameterSet& fullParamSet, const std::str i_sj1_raw_mass = data.addMulti(branchName,"sj1_raw_mass" , 0); i_sj1_csv = data.addMulti(branchName,"sj1_csv" , 0); if(!isRealData){ + i_sj1_JECUnc = data.addMulti(branchName,"sj1_JECUnc" , 0); i_sj1_hadronFlavor = data.addMulti(branchName,"sj1_hadronFlavor" , 0); i_sj1_partonFlavor = data.addMulti(branchName,"sj1_partonFlavor" , 0); } @@ -64,6 +71,7 @@ FatJetFiller::FatJetFiller(const edm::ParameterSet& fullParamSet, const std::str i_sj2_raw_mass = data.addMulti(branchName,"sj2_raw_mass" , 0); i_sj2_csv = data.addMulti(branchName,"sj2_csv" , 0); if(!isRealData){ + i_sj2_JECUnc = data.addMulti(branchName,"sj2_JECUnc" , 0); i_sj2_hadronFlavor = data.addMulti(branchName,"sj2_hadronFlavor" , 0); i_sj2_partonFlavor = data.addMulti(branchName,"sj2_partonFlavor" , 0); } @@ -77,6 +85,17 @@ void FatJetFiller::load(const edm::Event& iEvent, const edm::EventSetup& iSetup) if(!isRealData && fillGenJets) iEvent.getByToken(token_genJets ,han_genJets ); + + if(!isRealData){ + iSetup.get().get(jetType.c_str(),jetCorParameters); + JetCorrectorParameters const & JetCorPar = (*jetCorParameters)["Uncertainty"]; + jetCorUnc.reset(new JetCorrectionUncertainty(JetCorPar)); + + iSetup.get().get(subjetType.c_str(),subjetCorParameters); + JetCorrectorParameters const & subjetCorPar = (*subjetCorParameters)["Uncertainty"]; + subjetCorUnc.reset(new JetCorrectionUncertainty(subjetCorPar)); + } + loadedStatus = true; }; std::vector FatJetFiller::processGenJets(){ @@ -110,6 +129,8 @@ void FatJetFiller::fill(){ data.fillMulti(i_eta ,float(jet.eta() )); data.fillMulti(i_phi ,float(jet.phi() )); data.fillMulti(i_mass ,float(jet.mass())); + data.fillMulti(i_toRawFact, float(jet.jecFactor("Uncorrected"))); + data.fillMulti(i_csv , float(jet.bDiscriminator("pfCombinedInclusiveSecondaryVertexV2BJetTags"))); @@ -119,11 +140,17 @@ void FatJetFiller::fill(){ data.fillMulti(i_id ,idStat); - if(!isRealData && fillGenJets){ + if(!isRealData){ data.fillMulti(i_hadronFlavor ,ASTypes::convertTo(jet.hadronFlavour(),"JetFiller::hadronFlavor") ); data.fillMulti(i_partonFlavor ,ASTypes::convertTo(jet.partonFlavour(),"JetFiller::partonFlavor") ); - auto genRef = jet.genJetFwdRef().backRef(); - data.fillMulti(i_genIDX , genRef.isNull() ? size8(255) : genIndicies[genRef.key()] ); + jetCorUnc->setJetEta(jet.eta()); + jetCorUnc->setJetPt(jet.pt()); // here you must use the CORRECTED jet pt + data.fillMulti(i_JECUnc ,float(jetCorUnc->getUncertainty(true))); + + if(fillGenJets){ + auto genRef = jet.genJetFwdRef().backRef(); + data.fillMulti(i_genIDX , genRef.isNull() ? size8(255) : genIndicies[genRef.key()] ); + } } data.fillMulti(i_bbt , @@ -144,6 +171,9 @@ void FatJetFiller::fill(){ if(!isRealData){ data.fillMulti(i_sj1_hadronFlavor,ASTypes::convertTo(subjets[0]->hadronFlavour(),"FatJetFiller::hadronFlavor")); data.fillMulti(i_sj1_partonFlavor,ASTypes::convertTo(subjets[0]->partonFlavour(),"FatJetFiller::partonFlavor")); + subjetCorUnc->setJetEta(subjets[0]->eta()); + subjetCorUnc->setJetPt(subjets[0]->pt()); // here you must use the CORRECTED jet pt + data.fillMulti(i_sj1_JECUnc ,float(subjetCorUnc->getUncertainty(true))); } } else { data.fillMulti(i_sj1_pt ,float(0)); @@ -156,6 +186,8 @@ void FatJetFiller::fill(){ if(!isRealData){ data.fillMulti(i_sj1_hadronFlavor,size8(0)); data.fillMulti(i_sj1_partonFlavor,size8(0)); + data.fillMulti(i_sj1_JECUnc ,float(0)); + } } if(subjets.size() >1){ @@ -169,6 +201,9 @@ void FatJetFiller::fill(){ if(!isRealData){ data.fillMulti(i_sj2_hadronFlavor,ASTypes::convertTo(subjets[1]->hadronFlavour(),"FatJetFiller::hadronFlavor")); data.fillMulti(i_sj2_partonFlavor,ASTypes::convertTo(subjets[1]->partonFlavour(),"FatJetFiller::partonFlavor")); + subjetCorUnc->setJetEta(subjets[1]->eta()); + subjetCorUnc->setJetPt(subjets[1]->pt()); // here you must use the CORRECTED jet pt + data.fillMulti(i_sj2_JECUnc ,float(subjetCorUnc->getUncertainty(true))); } } else { data.fillMulti(i_sj2_pt ,float(0)); @@ -181,6 +216,7 @@ void FatJetFiller::fill(){ if(!isRealData){ data.fillMulti(i_sj2_hadronFlavor,size8(0)); data.fillMulti(i_sj2_partonFlavor,size8(0)); + data.fillMulti(i_sj1_JECUnc ,float(0)); } } diff --git a/TreeFillers/src/JetFiller.cc b/TreeFillers/src/JetFiller.cc index 70a0318..33004d8 100644 --- a/TreeFillers/src/JetFiller.cc +++ b/TreeFillers/src/JetFiller.cc @@ -1,6 +1,7 @@ #include "AnalysisTreeMaker/TreeFillers/interface/JetFiller.h" #include "AnalysisTreeMaker/TreeFillers/interface/FillerConstants.h" +#include "JetMETCorrections/Objects/interface/JetCorrectionsRecord.h" using ASTypes::size8; using ASTypes::int8; namespace AnaTM{ @@ -10,6 +11,7 @@ JetFiller::JetFiller(const edm::ParameterSet& fullParamSet, const std::string& p ,isRealData(isRealData) { if(ignore()) return; + jetType =cfg.getParameter("jetType"); fillGenJets =cfg.getParameter("fillGenJets"); token_jets =cc.consumes >(cfg.getParameter("jets")); minJetPT = cfg.getParameter("minJetPT"); @@ -20,12 +22,16 @@ JetFiller::JetFiller(const edm::ParameterSet& fullParamSet, const std::string& p i_eta = data.addMulti(branchName,"eta" , 0); i_phi = data.addMulti(branchName,"phi" , 0); i_mass = data.addMulti(branchName,"mass" , 0); + i_toRawFact = data.addMulti(branchName,"toRawFact" , 0); + i_metUnc_rawPx = data.addMulti(branchName,"metUnc_rawPx" , 0); + i_metUnc_rawPy = data.addMulti(branchName,"metUnc_rawPy" , 0); i_csv = data.addMulti(branchName,"csv" , 0); i_id = data.addMulti(branchName,"id" , 0); if(!isRealData){ i_hadronFlavor = data.addMulti(branchName,"hadronFlavor" , 0); i_partonFlavor = data.addMulti(branchName,"partonFlavor" , 0); + i_JECUnc = data.addMulti(branchName,"JECUnc" , 0); if(fillGenJets){ i_genIDX = data.addMulti(branchName,"genIDX" , 0); i_gen_pt = data.addMulti(branchName,"gen_pt" , 0); @@ -43,6 +49,12 @@ void JetFiller::load(const edm::Event& iEvent, const edm::EventSetup& iSetup) { if(!isRealData && fillGenJets) iEvent.getByToken(token_genJets ,han_genJets ); + if(!isRealData){ + iSetup.get().get(jetType.c_str(),jetCorParameters); + JetCorrectorParameters const & JetCorPar = (*jetCorParameters)["Uncertainty"]; + jetCorUnc.reset(new JetCorrectionUncertainty(JetCorPar)); + } + loadedStatus = true; }; @@ -138,6 +150,28 @@ void JetFiller::fill(){ data.fillMulti(i_eta ,float(jet.eta() )); data.fillMulti(i_phi ,float(jet.phi() )); data.fillMulti(i_mass ,float(jet.mass())); + const float rawFactor = jet.jecFactor("Uncorrected"); + data.fillMulti(i_toRawFact ,rawFactor); + + + reco::LeafCandidate::LorentzVector raw_p4 = jet.p4() * rawFactor; + const float emf = ( jet.neutralEmEnergy() + jet.chargedEmEnergy() )/raw_p4.E(); + if(emf >0.9){ + data.fillMulti(i_metUnc_rawPx ,0.0); + data.fillMulti(i_metUnc_rawPy ,0.0); + } else { + for ( unsigned int idau = 0; idau < jet.numberOfDaughters(); ++idau) { + const auto * pfcand = jet.daughter(idau); + if(pfcand->isGlobalMuon() || pfcand->isStandAloneMuon()){ + raw_p4 -= pfcand->p4(); + } + } + data.fillMulti(i_metUnc_rawPx ,float(raw_p4.px())); + data.fillMulti(i_metUnc_rawPy ,float(raw_p4.py())); + } + + + data.fillMulti(i_csv , float(jet.bDiscriminator("pfCombinedInclusiveSecondaryVertexV2BJetTags"))); @@ -152,6 +186,11 @@ void JetFiller::fill(){ if(!isRealData){ data.fillMulti(i_hadronFlavor ,ASTypes::convertTo(jet.hadronFlavour(),"JetFiller::hadronFlavor") ); data.fillMulti(i_partonFlavor ,ASTypes::convertTo(jet.partonFlavour(),"JetFiller::partonFlavor") ); + + jetCorUnc->setJetEta(jet.eta()); + jetCorUnc->setJetPt(jet.pt()); // here you must use the CORRECTED jet pt + data.fillMulti(i_JECUnc ,float(jetCorUnc->getUncertainty(true))); + if(fillGenJets){ auto genRef = jet.genJetFwdRef().backRef(); data.fillMulti(i_genIDX , genRef.isNull() ? size8(255) : genIndicies[genRef.key()] );