'Soft' codes for lepton trigger & ID SF

analysis/interface/nanoBase.h

@@ -9,6 +9,8 @@
 
 #include <TString.h>
 
+#include <string>
+
 #include "nano/external/interface/pileUpTool.h"
 #include "nano/external/interface/RoccoR.h"
 #include "nano/external/interface/lumiTool.h"
@@ -19,6 +21,7 @@
 #include "nano/external/interface/BTagCalibrationStandalone.h"
 //#include "nano/external/interface/TopTriggerSF.h"
 //#include "nano/external/interface/TTbarModeDefs.h"
+#include "nano/external/interface/computePtEtaTable.h"
 
 #include "TMVA/Tools.h"
 #include "TMVA/Reader.h"
@@ -45,6 +48,12 @@ class nanoBase : public Events
   ElecScaleFactorEvaluator m_elecSF;
   BTagCalibrationReader m_btagSF, m_btagSF_up, m_btagSF_dn;
 
+  std::string m_strTrigSFEl, m_strTrigSFMu;
+  std::string m_strLeptonSFEl, m_strLeptonSFMu;
+
+  computePtEtaTable m_tableTrigSFEl, m_tableTrigSFMu;
+  computePtEtaTable m_tableLeptonSFEl, m_tableLeptonSFMu;
+
   Bool_t m_isMC;
 };
 

analysis/src/nanoBase.cc

@@ -33,6 +33,25 @@ 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");
+
+  string strBaseData = env+"/src/nano/analysis/data/scaleFactor/2016/";
+  m_strTrigSFEl = strBaseData+"HLT_Ele32_eta2p1_WPTight_Gsf_FullRunRange.root";
+  m_strTrigSFMu = strBaseData+"EfficienciesAndSF_Run2016.root";
+  m_strLeptonSFEl = strBaseData+"2016LegacyReReco_ElectronTight.root";
+  m_strLeptonSFMu = strBaseData+"Run2016_SF_ID.root";
+
+  if (!exists_test(m_strTrigSFEl)   || !exists_test(m_strTrigSFMu) || 
+      !exists_test(m_strLeptonSFEl) || !exists_test(m_strLeptonSFMu))
+  {
+    std::cout << "Missing data file, run getFiles and try again" << std::endl;
+    exit(50);
+  }
+
+  m_tableTrigSFEl.LoadData(m_strTrigSFEl, "SF");
+  m_tableTrigSFMu.LoadData(m_strTrigSFMu, "IsoMu24_OR_IsoTkMu24_PtEtaBins/abseta_pt_ratio");
+
+  m_tableLeptonSFEl.LoadData(m_strLeptonSFEl, "EGamma_SF2D");
+  m_tableLeptonSFMu.LoadData(m_strLeptonSFMu, "NUM_TightID_DEN_genTracks_eta_pt");
 }
 
 nanoBase::~nanoBase()

analysis/src/topEventSelectionSL.cc

@@ -145,26 +145,27 @@ int topEventSelectionSL::EventSelection()
     recolep = muons[0];
     b_channel = CH_MU;
 
-    b_mueffweight    = muonSF_.getScaleFactor(recolep, 13, 0);
-    b_mueffweight_up = muonSF_.getScaleFactor(recolep, 13, 1);
-    b_mueffweight_dn = muonSF_.getScaleFactor(recolep, 13, -1);
+    b_mueffweight    = m_tableLeptonSFMu.getFactor(recolep.Pt(), recolep.Eta());
+    b_mueffweight_up = m_tableLeptonSFMu.getFactor(recolep.Pt(), recolep.Eta(),  1);
+    b_mueffweight_dn = m_tableLeptonSFMu.getFactor(recolep.Pt(), recolep.Eta(), -1);
   } else if (elecs.size() == 1) {
     recolep = elecs[0];
     b_channel = CH_EL;
 
-    b_eleffweight    = elecSF_.getScaleFactor(recolep, 11, 0);
-    b_eleffweight_up = elecSF_.getScaleFactor(recolep, 11, 1);
-    b_eleffweight_dn = elecSF_.getScaleFactor(recolep, 11, -1);
+    b_eleffweight    = m_tableLeptonSFEl.getFactor(recolep.Pt(), recolep.Eta());
+    b_eleffweight_up = m_tableLeptonSFEl.getFactor(recolep.Pt(), recolep.Eta(),  1);
+    b_eleffweight_dn = m_tableLeptonSFEl.getFactor(recolep.Pt(), recolep.Eta(), -1);
   }
 
   recolep.Momentum(b_lep);
   b_lep_pid = recolep.GetPdgCode();
   recoleps.push_back(b_lep);
 
   b_tri = b_tri_up = b_tri_dn = 0;
-  b_tri = 1; //computeTrigSF(recolep1, recolep2);
-  b_tri_up = 1; //computeTrigSF(recolep1, recolep2, 1);
-  b_tri_dn = 1; //computeTrigSF(recolep1, recolep2, -1);
+  computePtEtaTable &tableTrigSF = (std::abs(b_lep_pid) == 11 ? m_tableTrigSFEl : m_tableTrigSFMu);
+  b_tri    = tableTrigSF.getFactor(recolep.Pt(), recolep.Eta());
+  b_tri_up = tableTrigSF.getFactor(recolep.Pt(), recolep.Eta(),  1);
+  b_tri_dn = tableTrigSF.getFactor(recolep.Pt(), recolep.Eta(), -1);
 
   // Veto Leptons
 

external/interface/computePtEtaTable.h

@@ -0,0 +1,34 @@
+#ifndef computePtEtaTable_H
+#define computePtEtaTable_H
+
+
+#include <TFile.h>
+#include <TH2F.h>
+
+
+class computePtEtaTable {
+public:
+  bool m_bValid;
+
+  std::vector<std::vector<Float_t>> m_listVal;
+  std::vector<std::vector<Float_t>> m_listErr;
+
+  std::vector<Float_t> m_listPtBin;
+  std::vector<Float_t> m_listEtaBin;
+
+  bool m_bUseAbsEta;
+
+public:
+  computePtEtaTable() {};
+  ~computePtEtaTable() {};
+
+  bool isValid() {return m_bValid;};
+
+  int LoadData(std::string strPath, std::string strHistName);
+  double getFactor(Float_t fPt, Float_t fEta, int direction=0);
+};
+
+
+#endif // computePtEtaTable_H
+
+

external/scripts/mergeSFroots.py

@@ -0,0 +1,133 @@
+#!/usr/bin/python3
+
+
+import ROOT
+import os, sys, ctypes, array
+
+
+strOut = sys.argv[ 1 ] # The output file name
+strNameHist = sys.argv[ 2 ] # The name of histogram
+nOffsetArg = 3
+
+listSrc = []
+listLumi = []
+
+for i in range(( len(sys.argv) - nOffsetArg) // 2): 
+  listSrc.append(sys.argv[ 2 * i + nOffsetArg ])
+  listLumi.append(float(sys.argv[ 2 * i + nOffsetArg + 1 ]))
+
+fInvSumLumi = 1.0 / sum(listLumi)
+for i in range(len(listLumi)): listLumi[ i ] *= fInvSumLumi # What we actually want is relative lumis
+
+# Extracting binning info
+listBinX = []
+listBinY = []
+
+fFirst = ROOT.TFile.Open(listSrc[ 0 ])
+hFirst = fFirst.Get(strNameHist)
+strTitle = hFirst.GetTitle()
+
+hBinX = hFirst.ProjectionX()
+hBinY = hFirst.ProjectionY()
+
+nNumBinX = hBinX.GetNbinsX()
+nNumBinY = hBinY.GetNbinsX()
+
+for i in range(nNumBinX + 1): listBinX.append(hBinX.GetBinLowEdge(i + 1))
+for i in range(nNumBinY + 1): listBinY.append(hBinY.GetBinLowEdge(i + 1))
+
+fFirst.Close()
+
+# Preparing the output file and directory
+fNew = ROOT.TFile.Open(strOut, "RECREATE")
+
+# Some of source contains histogram into a directory. It will be also duplicated
+strDir = "/".join(strNameHist.split("/")[ : -1 ])
+if strDir != "": 
+  fNew.mkdir(strDir)
+  fNew.cd(strDir)
+
+# There will be data
+hNew = ROOT.TH2F(strNameHist.split("/")[ -1 ], strTitle, 
+  nNumBinX, array.array("d", listBinX), nNumBinY, array.array("d", listBinY))
+
+# The main part: Reading all source and combining the info from them
+# I used the following fomula: 
+#   SF  = sum_i ( SF_i * (rel. lumi)_i )
+#   err = sqrt( sum_i ( err_i * (rel. lumi)_i )^2 )
+# First, SF and err^2 will be calculated in stored in hNew. After that, sqrt(err^2) will be stored.
+fSrc = None
+strPathCurr = ""
+class NonConsistenceException(Exception): pass
+
+try: 
+  for nIdxSrc in range(len(listSrc)): 
+    strPathCurr = listSrc[ nIdxSrc ]
+    fSrc = ROOT.TFile.Open(strPathCurr)
+    hSrc = fSrc.Get(strNameHist)
+
+    if hSrc == None: raise NonConsistenceException
+
+    # Btw, for preventing a mistake, I added a check whether histograms are consistence or not
+    # The rule is easy: Does these histograms have same binning?
+    if nIdxSrc != 0: # The reference is the first histo, where we already got the bin infos of this
+      hBinXCurr = hSrc.ProjectionX()
+      hBinYCurr = hSrc.ProjectionY()
+
+      if nNumBinX != hBinXCurr.GetNbinsX() or nNumBinY != hBinYCurr.GetNbinsX(): 
+        raise NonConsistenceException
+
+      for i in range(nNumBinX + 1): 
+        if listBinX[ i ] != hBinXCurr.GetBinLowEdge(i + 1): 
+          raise NonConsistenceException
+
+      for i in range(nNumBinY + 1): 
+        if listBinY[ i ] != hBinYCurr.GetBinLowEdge(i + 1): 
+          raise NonConsistenceException
+
+    # Now, no problem. Time to calculate!
+    for i in range(nNumBinX): 
+      for j in range(nNumBinY): 
+        fCX = 0.5 * ( listBinX[ i ] + listBinX[ i + 1 ] )
+        fCY = 0.5 * ( listBinY[ j ] + listBinY[ j + 1 ] )
+
+        nIdxX = ctypes.c_int()
+        nIdxY = ctypes.c_int()
+        nIdxZ = ctypes.c_int()
+
+        # There is a strange issue that the bin indices don't correspond to the displayed bins, 
+        # which means, e.g., if I call hSrc.GetBinContent(1, 3) 
+        # (see HLT_Ele32_eta2p1_WPTight_Gsf_FullRunRange.root)
+        # then this method returns the content in (1, 1)-bin.
+        # Even this phenomonon does not occur always; only some of histograms has this issue.
+        # This complicated way to get the bin indices is for avoiding this problem.
+        hSrc.GetBinXYZ(hSrc.FindBin(fCX, fCY), nIdxX, nIdxY, nIdxZ)
+
+        fSF = hNew.GetBinContent(i + 1, j + 1)
+        fErr = hNew.GetBinError(i + 1, j + 1)
+
+        fSF  += hSrc.GetBinContent(nIdxX.value, nIdxY.value) * listLumi[ nIdxSrc ]
+        fErr += ( hSrc.GetBinError(nIdxX.value, nIdxY.value) * listLumi[ nIdxSrc ] ) ** 2
+
+        hNew.SetBinContent(i + 1, j + 1, fSF)
+        hNew.SetBinError(i + 1, j + 1, fErr)
+
+    fSrc.Close()
+
+except NonConsistenceException:
+  sys.stderr.write("FATAL ERROR: Non-consistence of histograms; %s\n"%strPathCurr)
+
+  fSrc.Close() # In here, one of source file had been opened but not closed
+  fNew.Close()
+
+  sys.exit(1)
+
+# As mentioned, storing sqrt(err^2)
+for nIdxX in range(nNumBinX): 
+  for nIdxY in range(nNumBinY): 
+    hNew.SetBinError(nIdxX + 1, nIdxY + 1, ( hNew.GetBinError(nIdxX + 1, nIdxY + 1) ) ** 0.5)
+
+fNew.Write()
+fNew.Close()
+
+

external/src/computePtEtaTable.cc

@@ -0,0 +1,102 @@
+#include "nano/external/interface/computePtEtaTable.h"
+
+#include <iostream>
+
+
+int computePtEtaTable::LoadData(std::string strPath, std::string strHistName) {
+  Int_t i, j;
+
+  TFile *fData;
+  TH2 *hData;
+
+  bool bFlipped;
+
+  TH1D *hEta, *hPt, *hTmp;
+  Int_t nNEta, nNPt, nTmp;
+
+  Float_t fBinCX, fBinCY;
+  Int_t nIdxX, nIdxY, nIdxZ;
+
+  m_bValid = false;
+
+  // Opening!
+  fData = TFile::Open(strPath.c_str());
+  if ( fData == NULL ) return -1;
+
+  hData = (TH2 *)fData->Get(strHistName.c_str());
+  m_bValid = true;
+
+  // To extract the binning info
+  hEta = hData->ProjectionX();
+  hPt  = hData->ProjectionY();
+
+  nNEta = hEta->GetNbinsX();
+  nNPt  = hPt->GetNbinsX();
+
+  bFlipped = false;
+
+  // Maybe... some of data is flipped; eta along y and pt along x
+  if (std::abs(hEta->GetBinLowEdge(nNEta+1)) > std::abs(hPt->GetBinLowEdge(nNPt+1))) {
+    hTmp = hEta; hEta = hPt; hPt = hTmp;
+    nTmp = nNEta; nNEta = nNPt; nNPt = nTmp;
+    bFlipped = true;
+  }
+
+  // Some of them has eta region [-2.4, 2.4], while others have [0, 2.4]
+  m_bUseAbsEta = hEta->GetBinLowEdge(1) > -0.1;
+
+  // Keeping the binning info
+  for (i = 1; i <= nNEta+1; i++) m_listEtaBin.push_back(hEta->GetBinLowEdge(i));
+  for (i = 1; i <= nNPt +1; i++) m_listPtBin.push_back(hPt->GetBinLowEdge(i));
+
+  // Making the store of weight and its uncertainty
+  for (i = 1; i <= nNEta; i++) {
+    m_listVal.emplace_back();
+    m_listErr.emplace_back();
+
+    for (j = 1; j <= nNPt; j++) {
+      if (!bFlipped) {
+        fBinCX = 0.5*(m_listEtaBin[i-1]+m_listEtaBin[i]);
+        fBinCY = 0.5*(m_listPtBin[j-1]+m_listPtBin[j]);
+      } else {
+        fBinCX = 0.5*(m_listPtBin[j-1]+m_listPtBin[j]);
+        fBinCY = 0.5*(m_listEtaBin[i-1]+m_listEtaBin[i]);
+      }
+
+      // There is a strange issue that the bin indices don't correspond to the displayed bins, 
+      // which means, e.g., if I call hSrc.GetBinContent(1, 3) 
+      // (see HLT_Ele32_eta2p1_WPTight_Gsf_FullRunRange.root)
+      // then this method returns the content in (1, 1)-bin.
+      // Even this phenomonon does not occur always; only some of histograms has this issue.
+      // This complicated way to get the bin indices is for avoiding this problem.
+      hData->GetBinXYZ(hData->FindBin(fBinCX, fBinCY), nIdxX, nIdxY, nIdxZ);
+
+      m_listVal[i-1].push_back(hData->GetBinContent(nIdxX, nIdxY));
+      m_listErr[i-1].push_back(hData->GetBinError(nIdxX, nIdxY));
+    }
+  }
+
+  fData->Close();
+
+  return 0;
+}
+
+double computePtEtaTable::getFactor(Float_t fPt, Float_t fEta, int direction) {
+  if (!m_bValid) return -1;
+
+  // Seeking the bin in which the value is contained
+  auto GetIdxRange = [](Float_t fX, std::vector<Float_t> listEdges) {
+    if (fX < listEdges[0]) return 0;
+    else if (fX >= listEdges[listEdges.size()-1]) return ((Int_t)listEdges.size())-2;
+    else return (int)(std::lower_bound(listEdges.begin(), listEdges.end(), fX)-listEdges.begin())-1;
+  };
+
+  Int_t nIdxPt, nIdxEta;
+
+  if (m_bUseAbsEta) fEta = std::abs(fEta);
+
+  nIdxPt  = GetIdxRange(fPt,  m_listPtBin);
+  nIdxEta = GetIdxRange(fEta, m_listEtaBin);
+
+  return m_listVal[nIdxEta][nIdxPt]+direction*m_listErr[nIdxEta][nIdxPt];
+}