Merging jet uncertainty evaluator

analysis/BuildFile.xml

@@ -1,4 +1,6 @@
 <use name="nano/external" />
+<use name="JetMETCorrections/Objects"/>
+<use name="JetMETCorrections/Modules"/>
 <use name="root"/>
 <use name="rootxml"/>
   <lib name="Eve"/>

analysis/interface/nanoBase.h

@@ -16,6 +16,11 @@
 #include "nano/external/interface/MuonScaleFactorEvaluator.h"
 #include "nano/external/interface/ElecScaleFactorEvaluator.h"
 
+#include "nano/external/interface/JetCorrectionUncertainty.h"
+#include "nano/external/interface/JetCorrectorParameters.h"
+#include "nano/external/interface/JetResolution.h"
+#include <TRandom3.h>
+
 #include "nano/external/interface/BTagCalibrationStandalone.h"
 //#include "nano/external/interface/TopTriggerSF.h"
 //#include "nano/external/interface/TTbarModeDefs.h"
@@ -32,6 +37,17 @@ class nanoBase : public Events
   nanoBase(TTree *tree, Bool_t isMC=false) : nanoBase(tree,0,0,isMC) {}
   virtual ~nanoBase();
   virtual void Loop() = 0;
+
+  // In uncertainty study we need to switch the kinematic variables of jets
+  // The following variables are for this switch
+  // In the topObjectSelection.cc these variables are used instead of Jet_pt, Jet_mass, and so on.
+  // In default, these are same as the original ones, but when a user wants to study systematic uncertainty 
+  // so that he/she needs to switch them to the evaluated ones, 
+  // just touching them in anlalyser class will be okay, and this is for it.
+  virtual void GetJetMassPt(UInt_t nIdx, 
+    Float_t &fJetMass, Float_t &fJetPt, Float_t &fJetEta, Float_t &fJetPhi, UInt_t unFlag);
+
+  Int_t GetMatchGenJet(UInt_t nIdxJet, Float_t fResolution);
 
   // Output Variables
   TFile* m_output;
@@ -44,8 +60,29 @@ class nanoBase : public Events
   MuonScaleFactorEvaluator m_muonSF;
   ElecScaleFactorEvaluator m_elecSF;
   BTagCalibrationReader m_btagSF, m_btagSF_up, m_btagSF_dn;
+  JetCorrectionUncertainty *m_jecUnc;
+  JMENano::JetResolution m_jetResObj;
+  JMENano::JetResolutionScaleFactor m_jetResSFObj;
+  TRandom3 *m_rndEngine;
+
+  Float_t m_fDRcone_JER;
+  Float_t m_fResFactorMathcer;
 
   Bool_t m_isMC;
+
+  enum {
+    OptBit_JER_Up = 0, 
+    OptBit_JER_Dn, 
+    OptBit_JES_Up, 
+    OptBit_JES_Dn
+  };
+
+  enum {
+    OptFlag_JER_Up = ( 1 << OptBit_JER_Up ), 
+    OptFlag_JER_Dn = ( 1 << OptBit_JER_Dn ), 
+    OptFlag_JES_Up = ( 1 << OptBit_JES_Up ), 
+    OptFlag_JES_Dn = ( 1 << OptBit_JES_Dn )
+  };
 };
 
 #endif

analysis/interface/topEventSelectionDL.h

@@ -52,8 +52,8 @@ class topEventSelectionDL : public topObjectSelection
 
   virtual int EventSelection();
 
-  topEventSelectionDL(TTree *tree=0, TTree *had=0, TTree *hadTruth=0, Bool_t isMC = false, Bool_t dl = true, Bool_t sle = false, Bool_t slm = false);
-  topEventSelectionDL(TTree *tree=0, Bool_t isMC=false, Bool_t dl=true, Bool_t sle=false, Bool_t slm=false) : topEventSelectionDL(tree, 0, 0, isMC, dl, sle, slm) {}
+  topEventSelectionDL(TTree *tree=0, TTree *had=0, TTree *hadTruth=0, Bool_t isMC = false, Bool_t dl = true, Bool_t sle = false, Bool_t slm = false, UInt_t unFlag = 0);
+  topEventSelectionDL(TTree *tree=0, Bool_t isMC=false, Bool_t dl=true, Bool_t sle=false, Bool_t slm=false, UInt_t unFlag = 0) : topEventSelectionDL(tree, 0, 0, isMC, dl, sle, slm, unFlag) {}
   ~topEventSelectionDL();  
   virtual void Loop() = 0;
 

analysis/interface/topEventSelectionSL.h

@@ -51,8 +51,8 @@ class topEventSelectionSL : public topObjectSelection
   virtual Bool_t TrigForEl() {return HLT_Ele27_WPTight_Gsf;}
   virtual Bool_t TrigForMu() {return HLT_IsoTkMu24 || HLT_IsoMu24;}
 
-  topEventSelectionSL(TTree *tree=0, TTree *had=0, TTree *hadTruth=0, Bool_t isMC = false, Bool_t sle = false, Bool_t slm = false);
-  topEventSelectionSL(TTree *tree=0, Bool_t isMC=false, Bool_t sle=false, Bool_t slm=false) : topEventSelectionSL(tree, 0, 0, isMC, sle, slm) {}
+  topEventSelectionSL(TTree *tree=0, TTree *had=0, TTree *hadTruth=0, Bool_t isMC = false, Bool_t sle = false, Bool_t slm = false, UInt_t unFlag = 0);
+  topEventSelectionSL(TTree *tree=0, Bool_t isMC=false, Bool_t sle=false, Bool_t slm=false, UInt_t unFlag = 0) : topEventSelectionSL(tree, 0, 0, isMC, sle, slm, unFlag) {}
   ~topEventSelectionSL();  
   virtual void Loop() = 0;
 };

analysis/interface/topObjectSelection.h

@@ -14,6 +14,8 @@ class topObjectSelection : public nanoBase
 
   Float_t b_maxBDiscr_nonb;
 
+  UInt_t m_unFlag;
+
 public:
   // YOU MUST SET UP ALL IN THE BELOW!!!
   // (SetCutValues() will force you to do it)
@@ -73,8 +75,11 @@ class topObjectSelection : public nanoBase
     return Muon_isPFcand[ nIdx ] && ( Muon_globalMu[ nIdx ] || Muon_trackerMu[ nIdx ] );
   };
   virtual bool additionalConditionForGenJet(UInt_t nIdx) {return true;};
-  virtual bool additionalConditionForJet(UInt_t nIdx) {return true;};
-  virtual bool additionalConditionForBJet(UInt_t nIdx) {return true;};
+
+  virtual bool additionalConditionForJet(UInt_t nIdx, Float_t &fJetPt, Float_t &fJetEta, Float_t &fJetPhi, Float_t &fJetMass) 
+    {return true;};
+  virtual bool additionalConditionForBJet(UInt_t nIdx, Float_t &fJetPt, Float_t &fJetEta, Float_t &fJetPhi, Float_t &fJetMass) 
+    {return true;};
 
 public:
   std::vector<TParticle> muonSelection();
@@ -87,8 +92,9 @@ class topObjectSelection : public nanoBase
 
   std::vector<TParticle> genJetSelection();
 
-  topObjectSelection(TTree *tree=0, TTree *had=0, TTree *hadTruth=0, Bool_t isMC = false);
-  topObjectSelection(TTree *tree=0, Bool_t isMC=false) : topObjectSelection(tree, 0, 0, isMC) {}
+  topObjectSelection(TTree *tree=0, TTree *had=0, TTree *hadTruth=0, Bool_t isMC = false, UInt_t unFlag = 0);
+  topObjectSelection(TTree *tree=0, Bool_t isMC=false, UInt_t unFlag = 0) : 
+    topObjectSelection(tree, 0, 0, isMC, unFlag) {}
   ~topObjectSelection() {}
 
   // In this function you need to set all the cut conditions in the above

analysis/src/nanoBase.cc

@@ -33,6 +33,39 @@ nanoBase::nanoBase(TTree *tree, TTree *had, TTree *hadTruth, Bool_t isMC) :
   m_btagSF.load(calib, BTagEntry::FLAV_B, "comb");
   m_btagSF_up.load(calib, BTagEntry::FLAV_B, "comb");
   m_btagSF_dn.load(calib, BTagEntry::FLAV_B, "comb");
+
+  m_jecUnc = NULL;
+  m_rndEngine = NULL;
+
+  if (isMC) {
+    std::string env = getenv("CMSSW_BASE");
+
+    std::string strPathJetResSFObj = env+"/src/nano/analysis/data/jer/"
+      "Summer16_25nsV1_MC_SF_AK4PFchs.txt";
+    std::string strPathJetResObj = env+"/src/nano/analysis/data/jer/"
+      "Summer16_25nsV1_MC_PtResolution_AK4PFchs.txt";
+    std::string strPathJecUnc = env + "/src/nano/analysis/data/jer/"
+      "Summer16_23Sep2016V4_MC_UncertaintySources_AK4PFchs.txt";
+
+    if (!exists_test(strPathJetResSFObj) || !exists_test(strPathJetResObj) || 
+        !exists_test(strPathJecUnc))
+    {
+      std::cout << "Missing data file, run getFiles and try again" << std::endl;
+      exit(50);
+    }
+
+    m_jetResSFObj = JMENano::JetResolutionScaleFactor(strPathJetResSFObj.c_str());
+    m_jetResObj = JMENano::JetResolution(strPathJetResObj.c_str());
+
+    m_rndEngine = new TRandom3(12345);
+
+
+    JetCorrectorParameters JetCorPar(strPathJecUnc, "Total");
+    m_jecUnc = new JetCorrectionUncertainty(JetCorPar);
+  }
+
+  m_fDRcone_JER = 0.4; // For AK4 jets
+  m_fResFactorMathcer = 3; // According to https://twiki.cern.ch/twiki/bin/viewauth/CMS/JetResolution
 }
 
 nanoBase::~nanoBase()
@@ -41,6 +74,90 @@ nanoBase::~nanoBase()
     m_output->Write();
     m_output->Close();
   }
+
+  if (m_jecUnc != NULL) delete m_jecUnc;
+  if (m_rndEngine != NULL) delete m_rndEngine;
+}
+
+// In uncertainty study we need to switch the kinematic variables of jets
+// In default, these are same as the original ones, but when a user wants to study systematic uncertainty 
+// so that he/she needs to switch them to the evaluated ones, 
+// just touching them in anlalyser class will be okay, and this is for it.
+void nanoBase::GetJetMassPt(UInt_t nIdx, 
+  Float_t &fJetMass, Float_t &fJetPt, Float_t &fJetEta, Float_t &fJetPhi, UInt_t unFlag) 
+{
+  Float_t fCorrFactor = 1.0;
+
+  fJetMass = Jet_mass[nIdx];
+  fJetPt   = Jet_pt[nIdx];
+  fJetEta  = Jet_eta[nIdx];
+  fJetPhi  = Jet_phi[nIdx];
+
+  //if ( m_isMC && ( unFlag & ( OptFlag_JER_Up | OptFlag_JER_Dn | OptFlag_JES_Up | OptFlag_JES_Dn ) ) != 0 )
+  if (m_isMC) {
+    // Evaluating the central part cJER of the factor
+    JMENano::JetParameters jetPars = {{JMENano::Binning::JetPt, fJetPt},
+                                  {JMENano::Binning::JetEta, fJetEta},
+                                  {JMENano::Binning::Rho, fixedGridRhoFastjetAll}};
+
+    const double jetRes = m_jetResObj.getResolution(jetPars); // Note: this is relative resolution.
+    const double cJER = m_jetResSFObj.getScaleFactor(jetPars, 
+      ((unFlag & (OptFlag_JER_Up | OptFlag_JER_Dn)) == 0 ? Variation::NOMINAL : 
+        ((unFlag & OptFlag_JER_Up) != 0 ? Variation::UP : Variation::DOWN)));
+
+    // We need corresponding genJet
+    Int_t nIdxGen = GetMatchGenJet(nIdx, fJetPt * jetRes);
+    const double genJetPt = GenJet_pt[nIdxGen];
+
+    if ( nIdxGen >= 0 ) {
+      fCorrFactor = (genJetPt+(fJetPt-genJetPt)*cJER)/fJetPt;
+    } else {
+      const double smear = m_rndEngine->Gaus(0, 1);
+      fCorrFactor = (cJER <= 1 ? 1 : 1+smear*jetRes*sqrt(cJER*cJER-1));
+    }
+
+    if ((unFlag & (OptFlag_JES_Up | OptFlag_JES_Dn)) != 0) { // JES
+      // The evaluator needs pT and eta of the current jet
+      m_jecUnc->setJetPt(fCorrFactor*fJetPt);
+      m_jecUnc->setJetEta(fJetEta);
+
+      if ((unFlag & OptFlag_JES_Up) != 0) {
+        fCorrFactor *= 1+m_jecUnc->getUncertainty(true);
+      } else {
+        fCorrFactor *= 1-m_jecUnc->getUncertainty(true);
+      }
+    }
+  }
+
+  fJetMass *= fCorrFactor;
+  fJetPt   *= fCorrFactor;
+}
+
+
+Int_t nanoBase::GetMatchGenJet(UInt_t nIdxJet, Float_t fResolution) {
+  UInt_t i;
+
+  double dEta, dPhi, dR;
+  double dRFound = m_fDRcone_JER;
+  UInt_t nIdxFound = 999;
+
+  for (i = 0; i < nGenJet; i++) {
+    dEta = Jet_eta[nIdxJet]-GenJet_eta[i];
+    dPhi = std::abs(Jet_phi[nIdxJet]-GenJet_phi[i]);
+    if (dPhi > (double)M_PI) dPhi -= (double)(2*M_PI);
+
+    dR = std::sqrt(dEta*dEta+dPhi*dPhi);
+
+    if (dR >= (double)m_fDRcone_JER*0.5) continue;
+    if (dRFound > dR) {
+      if (std::abs(Jet_pt[nIdxJet] - GenJet_pt[i]) >= m_fResFactorMathcer*fResolution) continue;
+
+      dRFound = dR;
+      nIdxFound = i;
+    }
+  }
+
+  return (dRFound < 0.75*(double)m_fDRcone_JER ? (Int_t)nIdxFound : -1);
 }
 
 void nanoBase::Loop()

analysis/src/topEventSelectionDL.cc

@@ -2,8 +2,8 @@
 
 using std::vector;
 
-topEventSelectionDL::topEventSelectionDL(TTree *tree, TTree *had, TTree *hadTruth, Bool_t isMC, Bool_t dl, Bool_t sle, Bool_t slm) :
-  topObjectSelection(tree, had, hadTruth, isMC),
+topEventSelectionDL::topEventSelectionDL(TTree *tree, TTree *had, TTree *hadTruth, Bool_t isMC, Bool_t dl, Bool_t sle, Bool_t slm, UInt_t unFlag) :
+  topObjectSelection(tree, had, hadTruth, isMC, unFlag),
   m_isDL(dl), m_isSL_e(sle), m_isSL_m(slm) 
 {
   SetCutValues();

analysis/src/topEventSelectionSL.cc

@@ -2,8 +2,8 @@
 
 using std::vector;
 
-topEventSelectionSL::topEventSelectionSL(TTree *tree, TTree *had, TTree *hadTruth, Bool_t isMC, Bool_t sle, Bool_t slm) :
-  topObjectSelection(tree, had, hadTruth, isMC),
+topEventSelectionSL::topEventSelectionSL(TTree *tree, TTree *had, TTree *hadTruth, Bool_t isMC, Bool_t sle, Bool_t slm, UInt_t unFlag) :
+  topObjectSelection(tree, had, hadTruth, isMC, unFlag),
   h_nevents(0),
   h_genweights(0),
   h_cutFlow(0),

analysis/src/topObjectSelection.cc

@@ -2,9 +2,10 @@
 
 using std::vector;
 
-topObjectSelection::topObjectSelection(TTree *tree, TTree *had, TTree *hadTruth, Bool_t isMC) :
-  nanoBase(tree, had, hadTruth, isMC)
-{}
+topObjectSelection::topObjectSelection(TTree *tree, TTree *had, TTree *hadTruth, Bool_t isMC, UInt_t unFlag) : 
+  nanoBase(tree, had, hadTruth, isMC), m_unFlag(unFlag)
+{
+}
 
 vector<TParticle> topObjectSelection::elecSelection() {
   vector<TParticle> elecs; 
@@ -126,17 +127,20 @@ vector<TParticle> topObjectSelection::jetSelection() {
   double csvWgtHF_dn = 1.0, csvWgtLF_dn = 1.0, csvWgtC_dn = 1.0, csvWgtTotal_dn = 1.0;
   vector<TParticle> jets;
   for (UInt_t i = 0; i < nJet; ++i){
-    if (Jet_pt[i] < cut_JetPt) continue;
-    if (std::abs(Jet_eta[i]) > cut_JetEta) continue;
+    Float_t fJetMass, fJetPt, fJetEta, fJetPhi;
+    GetJetMassPt(i, fJetMass, fJetPt, fJetEta, fJetPhi, m_unFlag);
+
+    if (fJetPt < cut_JetPt) continue;
+    if (std::abs(fJetEta) > cut_JetEta) continue;
     if (Jet_jetId[i] < cut_JetID) continue;
     TLorentzVector mom;
-    mom.SetPtEtaPhiM(Jet_pt[i], Jet_eta[i], Jet_phi[i], Jet_mass[i]);
+    mom.SetPtEtaPhiM(fJetPt, fJetEta, fJetPhi, fJetMass);
     bool hasOverLap = false;
     for (auto lep : recoleps){
         if (mom.TLorentzVector::DeltaR(lep) < cut_JetConeSizeOverlap) hasOverLap = true;
     }
     if (hasOverLap) continue;
-    if ( !additionalConditionForJet(i) ) continue;
+    if ( !additionalConditionForJet(i, fJetPt, fJetEta, fJetPhi, fJetMass) ) continue;
 
     auto jet = TParticle();
     jet.SetMomentum(mom);
@@ -187,18 +191,21 @@ vector<TParticle> topObjectSelection::jetSelection() {
 vector<TParticle> topObjectSelection::bjetSelection() {
   vector<TParticle> bjets;
   for (UInt_t i = 0; i < nJet; ++i ) {
-    if (Jet_pt[i] < cut_BJetPt) continue;
-    if (std::abs(Jet_eta[i]) > cut_BJetEta) continue;
+    Float_t fJetMass, fJetPt, fJetEta, fJetPhi;
+    GetJetMassPt(i, fJetMass, fJetPt, fJetEta, fJetPhi, m_unFlag);
+
+    if (fJetPt < cut_BJetPt) continue;
+    if (std::abs(fJetEta) > cut_BJetEta) continue;
     if (Jet_jetId[i] < cut_BJetID) continue;
     if (cut_BJetTypeBTag[i] < cut_BJetBTagCut) continue;
     TLorentzVector mom;
-    mom.SetPtEtaPhiM(Jet_pt[i], Jet_eta[i], Jet_phi[i], Jet_mass[i]);
+    mom.SetPtEtaPhiM(fJetPt, fJetEta, fJetPhi, fJetMass);
     bool hasOverLap = false;
     for (auto lep : recoleps) {
       if (mom.TLorentzVector::DeltaR(lep) < cut_BJetConeSizeOverlap) hasOverLap = true;
     }
     if (hasOverLap) continue;
-    if ( !additionalConditionForBJet(i) ) continue;
+    if ( !additionalConditionForBJet(i, fJetPt, fJetEta, fJetPhi, fJetMass) ) continue;
 
     auto bjet = TParticle();
     bjet.SetMomentum(mom);
@@ -207,3 +214,5 @@ vector<TParticle> topObjectSelection::bjetSelection() {
   }
   return bjets;
 }
+
+