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/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

@@ -5,6 +5,11 @@
 #include "nano/external/interface/TopTriggerSF.h"
 //#include "nano/external/interface/TTbarModeDefs.h"
 
+#include "nano/external/interface/JetCorrectionUncertainty.h"
+#include "nano/external/interface/JetCorrectorParameters.h"
+#include "nano/external/interface/JetResolution.h"
+#include <TRandom3.h>
+
 class topObjectSelection : public nanoBase
 {
 protected:
@@ -14,6 +19,14 @@ class topObjectSelection : public nanoBase
 
   Float_t b_maxBDiscr_nonb;
 
+  UInt_t m_unFlag;
+
+  JetCorrectionUncertainty *jecUnc;
+
+  JMENano::JetResolution jetResObj;
+  JMENano::JetResolutionScaleFactor jetResSFObj;
+  TRandom3 *rndEngine;
+
 public:
   // YOU MUST SET UP ALL IN THE BELOW!!!
   // (SetCutValues() will force you to do it)
@@ -58,6 +71,9 @@ class topObjectSelection : public nanoBase
   Float_t  cut_BJetConeSizeOverlap;
   Float_t *cut_BJetTypeBTag; // For example, set it as cut_BJetTypeBTag = Jet_btagCSVV2;
   Float_t  cut_BJetBTagCut;
+
+  Float_t m_fDRcone_JER;
+  Float_t m_fResFactorMathcer;
 
 public: 
   // Tip: If you want to use your own additional cut with the existing cut, 
@@ -73,8 +89,22 @@ 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;};
+
+  // 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);
+
+  Int_t GetMatchGenJet(UInt_t nIdxJet, Float_t fResolution);
 
 public:
   std::vector<TParticle> muonSelection();
@@ -87,15 +117,33 @@ 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() {}
+  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() {
+    if ( jecUnc != NULL ) delete jecUnc;
+    if ( rndEngine != NULL ) delete rndEngine;
+  }
 
   // In this function you need to set all the cut conditions in the above
   // If you do not set this function up (so that you didn't set the cuts), the compiler will deny your code, 
   // so you can be noticed that you forgot the setting up.
   // And you don't need to run this function indivisually; it will be run in the creator of this class.
   virtual int SetCutValues() = 0;
+
+  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/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,36 @@
 
 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), jecUnc(NULL), rndEngine(NULL)
+{
+  //if ( ( m_unFlag & ( OptFlag_JES_Up | OptFlag_JES_Dn | OptFlag_JER_Up | OptFlag_JER_Dn ) ) != 0 )
+  if ( m_isMC ) {
+    std::string env = getenv("CMSSW_BASE");
+
+    std::string strPathJetResSFObj = env + "/src/nano/analysis/data/jer/"
+      "Summer16_25nsV1_MC_SF_AK4PFchs.txt";
+    jetResSFObj = JMENano::JetResolutionScaleFactor(strPathJetResSFObj.c_str());
+
+    std::string strPathJetResObj = env + "/src/nano/analysis/data/jer/"
+      "Summer16_25nsV1_MC_PtResolution_AK4PFchs.txt";
+    jetResObj = JMENano::JetResolution(strPathJetResObj.c_str());
+
+    rndEngine = new TRandom3(12345);
+
+    if ( ( m_unFlag & ( OptFlag_JES_Up | OptFlag_JES_Dn ) ) != 0 ) {
+      std::string strPathJecUnc = env + "/src/nano/analysis/data/jer/"
+        //"Summer16_23Sep2016V4_MC_Uncertainty_AK4PFchs.txt";
+        "Summer16_23Sep2016V4_MC_UncertaintySources_AK4PFchs.txt";
+
+      JetCorrectorParameters JetCorPar(strPathJecUnc, "Total");
+      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
+}
 
 vector<TParticle> topObjectSelection::elecSelection() {
   vector<TParticle> elecs; 
@@ -126,17 +153,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);
+
+    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 +217,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);
+
+    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 +240,95 @@ vector<TParticle> topObjectSelection::bjetSelection() {
   }
   return bjets;
 }
+
+
+// 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.
+void topObjectSelection::GetJetMassPt(UInt_t nIdx, 
+  Float_t &fJetMass, Float_t &fJetPt, Float_t &fJetEta, Float_t &fJetPhi) 
+{
+  Float_t fCorrFactor = 1.0;
+
+  fJetMass = Jet_mass[ nIdx ];
+  fJetPt   = Jet_pt[ nIdx ];
+  fJetEta  = Jet_eta[ nIdx ];
+  fJetPhi  = Jet_phi[ nIdx ];
+
+  //if ( m_isMC && ( m_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 = jetResObj.getResolution(jetPars); // Note: this is relative resolution.
+    const double cJER = jetResSFObj.getScaleFactor(jetPars, 
+      ( ( m_unFlag & ( OptFlag_JER_Up | OptFlag_JER_Dn ) ) == 0 ? Variation::NOMINAL : 
+        ( ( m_unFlag & OptFlag_JER_Up ) != 0 ? Variation::UP : Variation::DOWN ) ));
+
+    // We need corresponding genJet
+    //Int_t nIdxGen = Jet_genJetIdx[ nIdx ];
+    Int_t nIdxGen = GetMatchGenJet(nIdx, fJetPt * jetRes);
+    const double genJetPt = GenJet_pt[ nIdxGen ];
+
+    // JER (nominal and up and down) - apply scaling method if matched genJet is found,
+    //       apply gaussian smearing method if unmatched
+    /*if ( nIdxGen >= 0 && //deltaR(genJet->p4(), jet.p4()) < 0.2 && // From CATTool
+         std::abs(genJetPt - fJetPt) < jetRes * 3 * fJetPt ) 
+      fCorrFactor = std::max(0., (genJetPt + ( fJetPt - genJetPt ) * cJER) / fJetPt);*/
+    if ( nIdxGen >= 0 ) {
+      fCorrFactor = ( genJetPt + ( fJetPt - genJetPt ) * cJER ) / fJetPt;
+    } else {
+      const double smear = rndEngine->Gaus(0, 1);
+      fCorrFactor = ( cJER <= 1 ? 1 : 1 + smear * jetRes * sqrt(cJER * cJER - 1) );
+    }
+
+    if ( ( m_unFlag & ( OptFlag_JES_Up | OptFlag_JES_Dn ) ) != 0 ) { // JES
+      // The evaluator needs pT and eta of the current jet
+      jecUnc->setJetPt(fCorrFactor * fJetPt);
+      jecUnc->setJetEta(fJetEta);
+
+      if ( ( m_unFlag & OptFlag_JES_Up ) != 0 ) {
+        fCorrFactor *= 1 + jecUnc->getUncertainty(true);
+      } else {
+        fCorrFactor *= 1 - jecUnc->getUncertainty(true);
+      }
+    }
+  }
+
+  fJetMass *= fCorrFactor;
+  fJetPt   *= fCorrFactor;
+}
+
+
+Int_t topObjectSelection::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 );
+}
+
+