diff --git a/analysis/data/scaleFactor/2016/2016LegacyReReco_ElectronTight.root b/analysis/data/scaleFactor/2016/2016LegacyReReco_ElectronTight.root new file mode 100644 index 00000000..bb86cf25 Binary files /dev/null and b/analysis/data/scaleFactor/2016/2016LegacyReReco_ElectronTight.root differ diff --git a/analysis/data/scaleFactor/2016/EfficienciesAndSF_Run2016.root b/analysis/data/scaleFactor/2016/EfficienciesAndSF_Run2016.root new file mode 100644 index 00000000..01855ebb Binary files /dev/null and b/analysis/data/scaleFactor/2016/EfficienciesAndSF_Run2016.root differ diff --git a/analysis/data/scaleFactor/2016/HLT_Ele32_eta2p1_WPTight_Gsf_FullRunRange.root b/analysis/data/scaleFactor/2016/HLT_Ele32_eta2p1_WPTight_Gsf_FullRunRange.root new file mode 100644 index 00000000..84087aa9 Binary files /dev/null and b/analysis/data/scaleFactor/2016/HLT_Ele32_eta2p1_WPTight_Gsf_FullRunRange.root differ diff --git a/analysis/data/scaleFactor/2016/Run2016_SF_ID.root b/analysis/data/scaleFactor/2016/Run2016_SF_ID.root new file mode 100644 index 00000000..89c46b91 Binary files /dev/null and b/analysis/data/scaleFactor/2016/Run2016_SF_ID.root differ diff --git a/analysis/interface/nanoBase.h b/analysis/interface/nanoBase.h index a251c758..9ad835bd 100644 --- a/analysis/interface/nanoBase.h +++ b/analysis/interface/nanoBase.h @@ -9,6 +9,8 @@ #include +#include + #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; }; diff --git a/analysis/src/nanoBase.cc b/analysis/src/nanoBase.cc index 5e680e2e..1f45105c 100644 --- a/analysis/src/nanoBase.cc +++ b/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() diff --git a/analysis/src/topEventSelectionSL.cc b/analysis/src/topEventSelectionSL.cc index c78cd5bf..4b26e528 100644 --- a/analysis/src/topEventSelectionSL.cc +++ b/analysis/src/topEventSelectionSL.cc @@ -145,16 +145,16 @@ 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); @@ -162,9 +162,10 @@ int topEventSelectionSL::EventSelection() 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 diff --git a/external/interface/computePtEtaTable.h b/external/interface/computePtEtaTable.h new file mode 100644 index 00000000..986e93c5 --- /dev/null +++ b/external/interface/computePtEtaTable.h @@ -0,0 +1,34 @@ +#ifndef computePtEtaTable_H +#define computePtEtaTable_H + + +#include +#include + + +class computePtEtaTable { +public: + bool m_bValid; + + std::vector> m_listVal; + std::vector> m_listErr; + + std::vector m_listPtBin; + std::vector 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 + + diff --git a/external/scripts/mergeSFroots.py b/external/scripts/mergeSFroots.py new file mode 100755 index 00000000..627dbcd3 --- /dev/null +++ b/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() + + diff --git a/external/src/computePtEtaTable.cc b/external/src/computePtEtaTable.cc new file mode 100644 index 00000000..d0c99d63 --- /dev/null +++ b/external/src/computePtEtaTable.cc @@ -0,0 +1,102 @@ +#include "nano/external/interface/computePtEtaTable.h" + +#include + + +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 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]; +}