Support for JECs

brentstone · Mar 27, 2018 · 470ff00 · 470ff00
1 parent b2633fa
commit 470ff00
Show file tree

Hide file tree

Showing 7 changed files with 111 additions and 3 deletions.
diff --git a/TreeFillers/BuildFile.xml b/TreeFillers/BuildFile.xml
@@ -8,6 +8,8 @@
 <use   name="DataFormats/PatCandidates"/>
 <use   name="SimDataFormats/PileupSummaryInfo"/>
 <use   name="SimDataFormats/GeneratorProducts"/>
+<use   name="CondFormats/JetMETObjects"/>
+<use   name="JetMETCorrections/Objects"/>
 <export>
   <lib   name="1"/>
 </export>

diff --git 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<std::vector<pat::Jet> >          token_jets;
@@ -68,6 +75,12 @@ class FatJetFiller : public BaseFiller {
 
     edm::Handle<std::vector<pat::Jet> >          han_jets;
     edm::Handle<reco::GenJetCollection>          han_genJets;
+
+    edm::ESHandle<JetCorrectorParametersCollection> jetCorParameters;
+    std::unique_ptr<JetCorrectionUncertainty> jetCorUnc;
+
+    edm::ESHandle<JetCorrectorParametersCollection> subjetCorParameters;
+    std::unique_ptr<JetCorrectionUncertainty> subjetCorUnc;
 };
 }
 

diff --git 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,17 +28,22 @@ 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;
 	size i_gen_phi            = 0;
 	size i_gen_mass           = 0;
 
     bool isRealData           = false;
+    std::string jetType           = "";
     bool fillGenJets          = false;
     edm::EDGetTokenT<std::vector<pat::Jet> >          token_jets;
     edm::EDGetTokenT<reco::GenJetCollection>          token_genJets;
@@ -44,6 +52,9 @@ class JetFiller : public BaseFiller {
     edm::Handle<std::vector<pat::Jet> >          han_jets;
     edm::Handle<reco::GenJetCollection>          han_genJets;
 
+    edm::ESHandle<JetCorrectorParametersCollection> jetCorParameters;
+    std::unique_ptr<JetCorrectionUncertainty> jetCorUnc;
+
 };
 }
 

diff --git 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'),
 )
diff --git 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'),
 )
diff --git 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<std::string>("jetType");
+    subjetType   =cfg.getParameter<std::string>("subjetType");
+
     fillGenJets  =cfg.getParameter<bool>("fillGenJets");
     token_jets   =cc.consumes<std::vector<pat::Jet> >(cfg.getParameter<edm::InputTag>("jets"));
     minJetPT     = cfg.getParameter<double>("minJetPT");
@@ -26,12 +30,14 @@ FatJetFiller::FatJetFiller(const edm::ParameterSet& fullParamSet, const std::str
     i_eta            = data.addMulti<float>(branchName,"eta"                   , 0);
     i_phi            = data.addMulti<float>(branchName,"phi"                   , 0);
     i_mass           = data.addMulti<float>(branchName,"mass"                  , 0);
+    i_toRawFact      = data.addMulti<float>(branchName,"toRawFact"             , 0);
     i_csv            = data.addMulti<float>(branchName,"csv"                   , 0);
     i_id             = data.addMulti<size8>(branchName,"id"                    , 0);
 
     if(!isRealData){
         i_hadronFlavor   = data.addMulti<int8>(branchName,"hadronFlavor"          , 0);
         i_partonFlavor   = data.addMulti<int8>(branchName,"partonFlavor"          , 0);
+        i_JECUnc         = data.addMulti<float>(branchName,"JECUnc"               , 0);
         if(fillGenJets){
             i_genIDX         = data.addMulti<size8>(branchName,"genIDX"                , 0);
             i_gen_pt         = data.addMulti<float>(branchName,"gen_pt"                , 0);
@@ -53,6 +59,7 @@ FatJetFiller::FatJetFiller(const edm::ParameterSet& fullParamSet, const std::str
     i_sj1_raw_mass       = data.addMulti<float>(branchName,"sj1_raw_mass"                        , 0);
     i_sj1_csv            = data.addMulti<float>(branchName,"sj1_csv"                             , 0);
     if(!isRealData){
+        i_sj1_JECUnc         = data.addMulti<float>(branchName,"sj1_JECUnc"               , 0);
         i_sj1_hadronFlavor   = data.addMulti<int8>(branchName,"sj1_hadronFlavor"                    , 0);
         i_sj1_partonFlavor   = data.addMulti<int8>(branchName,"sj1_partonFlavor"                    , 0);
     }
@@ -64,6 +71,7 @@ FatJetFiller::FatJetFiller(const edm::ParameterSet& fullParamSet, const std::str
     i_sj2_raw_mass       = data.addMulti<float>(branchName,"sj2_raw_mass"                        , 0);
     i_sj2_csv            = data.addMulti<float>(branchName,"sj2_csv"                             , 0);
     if(!isRealData){
+        i_sj2_JECUnc         = data.addMulti<float>(branchName,"sj2_JECUnc"               , 0);
         i_sj2_hadronFlavor   = data.addMulti<int8>(branchName,"sj2_hadronFlavor"                    , 0);
         i_sj2_partonFlavor   = data.addMulti<int8>(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<JetCorrectionsRecord>().get(jetType.c_str(),jetCorParameters);
+        JetCorrectorParameters const & JetCorPar = (*jetCorParameters)["Uncertainty"];
+        jetCorUnc.reset(new JetCorrectionUncertainty(JetCorPar));
+
+        iSetup.get<JetCorrectionsRecord>().get(subjetType.c_str(),subjetCorParameters);
+        JetCorrectorParameters const & subjetCorPar = (*subjetCorParameters)["Uncertainty"];
+        subjetCorUnc.reset(new JetCorrectionUncertainty(subjetCorPar));
+    }
+
     loadedStatus = true;
 };
 std::vector<size8> 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<int8>(jet.hadronFlavour(),"JetFiller::hadronFlavor") );
             data.fillMulti(i_partonFlavor ,ASTypes::convertTo<int8>(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<int8>(subjets[0]->hadronFlavour(),"FatJetFiller::hadronFlavor"));
                 data.fillMulti(i_sj1_partonFlavor,ASTypes::convertTo<int8>(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<int8>(subjets[1]->hadronFlavour(),"FatJetFiller::hadronFlavor"));
                 data.fillMulti(i_sj2_partonFlavor,ASTypes::convertTo<int8>(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
@@ -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<std::string>("jetType");
     fillGenJets  =cfg.getParameter<bool>("fillGenJets");
     token_jets   =cc.consumes<std::vector<pat::Jet> >(cfg.getParameter<edm::InputTag>("jets"));
     minJetPT     = cfg.getParameter<double>("minJetPT");
@@ -20,12 +22,16 @@ JetFiller::JetFiller(const edm::ParameterSet& fullParamSet, const std::string& p
     i_eta            = data.addMulti<float>(branchName,"eta"                   , 0);
     i_phi            = data.addMulti<float>(branchName,"phi"                   , 0);
     i_mass           = data.addMulti<float>(branchName,"mass"                  , 0);
+    i_toRawFact      = data.addMulti<float>(branchName,"toRawFact"             , 0);
+    i_metUnc_rawPx   = data.addMulti<float>(branchName,"metUnc_rawPx"          , 0);
+    i_metUnc_rawPy   = data.addMulti<float>(branchName,"metUnc_rawPy"          , 0);
     i_csv            = data.addMulti<float>(branchName,"csv"                   , 0);
     i_id             = data.addMulti<size8>(branchName,"id"                    , 0);
 
     if(!isRealData){
         i_hadronFlavor   = data.addMulti<int8>(branchName,"hadronFlavor"          , 0);
         i_partonFlavor   = data.addMulti<int8>(branchName,"partonFlavor"          , 0);
+        i_JECUnc         = data.addMulti<float>(branchName,"JECUnc"                , 0);
         if(fillGenJets){
             i_genIDX         = data.addMulti<size8>(branchName,"genIDX"                , 0);
             i_gen_pt         = data.addMulti<float>(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<JetCorrectionsRecord>().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<int8>(jet.hadronFlavour(),"JetFiller::hadronFlavor") );
             data.fillMulti(i_partonFlavor ,ASTypes::convertTo<int8>(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()] );