allow instantiating corrections with different constructors

Makefile

@@ -2,10 +2,12 @@ SOURCE_DIR=TIMBER/Framework/src/
 HEADER_DIR=TIMBER/Framework/include/
 EXT_DIR=TIMBER/Framework/ext/
 BIN_DIR=bin/libtimber/
+BOOST_DIR=bin/libboost/include/
+FMT_DIR=bin/fmt/include
 
 CC=gcc
-INCLUDE=-I/usr/include/ -I bin/ -I./ `root-config --cflags --ldflags --glibs` -I/usr/include/ -I$(EXT_DIR)
-LIBS=-lstdc++ -l boost_wserialization -l boost_filesystem -L bin/libarchive/lib/ -Wl,-rpath=bin/libarchive/lib/ -l archive
+INCLUDE=-I/usr/include/ -I bin/ -I./ `root-config --cflags --ldflags --glibs` -I/usr/include/ -I$(EXT_DIR) -I$(BOOST_DIR) -I$(FMT_DIR)
+LIBS=-lstdc++ -L bin/libarchive/lib/ -Wl,-rpath=bin/libarchive/lib/ -l archive
 CFLAGS=-g -Wno-attributes -fPIC -c
 CVMFS=/cvmfs/cms.cern.ch/$(SCRAM_ARCH)/cms/cmssw/$(CMSSW_VERSION)
 CMSSW=-I$(CVMFS)/src -I/cvmfs/cms.cern.ch/$(SCRAM_ARCH)/external/boost/1.72.0-bcolbf/include -L/cvmfs/cms.cern.ch/$(SCRAM_ARCH)/external/boost/1.72.0-bcolbf/lib -L$(CVMFS)/lib/$(SCRAM_ARCH)

TIMBER/Analyzer.py

@@ -124,7 +124,7 @@ def __init__(self,fileName,eventsTreeName="Events",runTreeName="Runs",multiSampl
             self.isData = False
         else:
             self.isData = True
- 
+
         # Count number of generated events if not data
         self.genEventSumw = 0.0
         self.genEventCount = 0
@@ -167,10 +167,11 @@ def __init__(self,fileName,eventsTreeName="Events",runTreeName="Runs",multiSampl
         if 'CMSSW_BASE' not in os.environ.keys():
             skipHeaders = ['JME_common.h','JetSmearer.h','JetRecalibrator.h','JES_weight.h','JER_weight.h','JMS_weight.h','JMR_weight.h']
 
+
         for f in glob.glob(os.environ["TIMBERPATH"]+'TIMBER/Framework/include/*.h'):
             if f.split('/')[-1] in skipHeaders: continue
             CompileCpp('#include "%s"\n'%f)
- 
+
     def _addFile(self,f):
         '''Add file to TChains being tracked.
 
@@ -1909,7 +1910,7 @@ class ModuleWorker(object):
     Writing the C++ modules requires the desired branch/column names must be specified or be used as the argument variable names
     to allow the framework to automatically determine what branch/column to use in GetCall().
     '''
-    def __init__(self,name,script,constructor=[],mainFunc='eval',columnList=None,isClone=False,cloneFuncInfo=None):
+    def __init__(self,name,script,constructor=[],mainFunc='eval',columnList=None,isClone=False,cloneFuncInfo=None,isNewConstr=False):
         '''Constructor
 
         @param name (str): Unique name to identify the instantiated worker object.
@@ -1931,8 +1932,7 @@ def __init__(self,name,script,constructor=[],mainFunc='eval',columnList=None,isC
         # Correction name
         self.name = name
         self._script = self._getScript(script)
-        if not isClone: self._funcInfo = self._getFuncInfo(mainFunc)
-        else: self._funcInfo = cloneFuncInfo
+        self._funcInfo = self._getFuncInfo(mainFunc)
         self._mainFunc = list(self._funcInfo.keys())[0]
         self._columnNames = LoadColumnNames() if columnList == None else columnList
         self._constructor = constructor 
@@ -1949,6 +1949,9 @@ def __init__(self,name,script,constructor=[],mainFunc='eval',columnList=None,isC
 
             self._instantiate(constructor)
 
+        if isNewConstr:
+            self._instantiate(constructor)            
+
     def Clone(self,name,newMainFunc=None):
         '''Makes a clone of current instance.
 
@@ -2204,7 +2207,7 @@ class Correction(ModuleWorker):
     (2) the return must be a vector ordered as {nominal, up, down} for "weight" type and 
     {up, down} for "uncert" type.
     '''
-    def __init__(self,name,script='',constructor=[],mainFunc='eval',corrtype=None,columnList=None,isClone=False,cloneFuncInfo=None):
+    def __init__(self,name,script='',constructor=[],mainFunc='eval',corrtype=None,columnList=None,isClone=False,cloneFuncInfo=None,isNewConstr=False):
         '''Constructor
 
         @param name (str): Correction name.
@@ -2231,7 +2234,7 @@ def __init__(self,name,script='',constructor=[],mainFunc='eval',corrtype=None,co
         # str
         # Name of correction 
         if script != '':
-            super(Correction,self).__init__(name,script,constructor,mainFunc,columnList,isClone,cloneFuncInfo)
+            super(Correction,self).__init__(name,script,constructor,mainFunc,columnList,isClone,cloneFuncInfo,isNewConstr)
             self.existing = False
         else:
             self.existing = True

TIMBER/Framework/Zbb_modules/AK4Btag_SF.cc

@@ -0,0 +1,130 @@
+// without CMSSW / standalone:
+#include "TIMBER/Framework/ext/BTagCalibrationStandalone.h"
+#include "TIMBER/Framework/ext/BTagCalibrationStandalone.cpp"
+#include <ROOT/RVec.hxx>
+
+using namespace ROOT::VecOps;
+
+/**
+ * @class AK4Btag_SF
+ * @brief b tagging scale factor lookup.
+ */
+class AK4Btag_SF {
+    public:
+        /**
+         * @brief Construct a new b tag scale factor lookup object
+         * 
+         * @param year 16, 17, or 18.
+         * @param tagger Ex. DeepCSV, DeepJet
+         * @param op_string "loose", "medium", "tight", "reshaping"
+         */
+        AK4Btag_SF(std::string year, std::string tagger, std::string op_string);
+        ~AK4Btag_SF(){};
+        /**
+         * @brief Per-event evaluation function
+         * 
+         * @param pt \f$p_{T}\f$ of subjet
+         * @param eta \f$\eta\f$ of subjet
+         * @return RVec<float> Nominal, up, down scale factor values.
+         */
+        RVec<float> eval(float pt, float eta,int flav, float disc);
+        RVec<float> evalCollection(int nJet, RVec<float> pt, RVec<float> eta, RVec<float> flav, RVec<float> disc,int var);
+    private:
+        std::string csv_file;
+        BTagEntry::OperatingPoint operating_point;
+        BTagCalibration calib;
+        BTagCalibrationReader reader;
+};
+
+AK4Btag_SF::AK4Btag_SF(std::string year, std::string tagger, std::string op_string){
+		std::cout<<tagger<<" "<<year<<" "<<op_string<<"\n";
+        if (op_string == "loose") {
+            operating_point = BTagEntry::OP_LOOSE;
+        } else if (op_string == "medium") {
+            operating_point = BTagEntry::OP_MEDIUM;
+        } else if (op_string == "tight") {
+            operating_point = BTagEntry::OP_TIGHT;
+        } else if (op_string == "reshaping") {
+            operating_point = BTagEntry::OP_RESHAPING;
+        }
+        else{
+            throw "Operating point type not supported!";
+        }
+
+        if (year == "2016APV") {
+            csv_file = std::string(std::getenv("TIMBERPATH"))+"TIMBER/data/OfficialSFs/FILENAME.csv";
+        } else if (year == "2016") {
+            csv_file = std::string(std::getenv("TIMBERPATH"))+"TIMBER/data/OfficialSFs/FILENAME.csv";
+        } else if (year == "2017") {
+            csv_file = std::string(std::getenv("TIMBERPATH"))+"TIMBER/data/OfficialSFs/DeepJet_106XUL17SF_V2.csv";
+        } else if (year == "2018") {
+            csv_file = std::string(std::getenv("TIMBERPATH"))+"TIMBER/data/OfficialSFs/DeepJet_106XUL18SF.csv";
+        }
+
+        std::cout<<csv_file<<"\n";
+        // setup calibration + reader
+        calib = BTagCalibration(tagger, csv_file);
+        std::cout<<"Created calib\n";
+        reader = BTagCalibrationReader(operating_point,  // operating point
+                                     "central",             // central sys type
+                                     {"up_hf", "down_hf"});      // other sys types
+        std::cout<<"Created reader\n";
+        reader.load(calib,               
+                    BTagEntry::FLAV_B,    
+                    "iterativefit");
+        reader.load(calib,
+                    BTagEntry::FLAV_C,
+                    "iterativefit");
+        reader.load(calib,
+                    BTagEntry::FLAV_UDSG,
+                    "iterativefit");
+        std::cout<<"Loaded reader\n";
+
+}
+
+
+RVec<float> AK4Btag_SF::eval(float pt, float eta, int flav, float disc) {
+    RVec<float> jet_scalefactor(3);
+
+    BTagEntry::JetFlavor fl = static_cast<BTagEntry::JetFlavor>(flav);
+
+    float nom = reader.eval_auto_bounds("central", fl, eta, pt, disc);//eta, pt, discr
+    //float up = reader.eval_auto_bounds("up_hf", fl, eta, pt, disc);
+    //float down = reader.eval_auto_bounds("down_hf", fl, eta, pt, disc);
+    float down = 0.98*nom;
+    float up = 1.02*nom;
+
+    jet_scalefactor[0] = nom;
+    jet_scalefactor[1] = down;
+    jet_scalefactor[2] = up;
+
+    return jet_scalefactor;
+};
+
+RVec<float> AK4Btag_SF::evalCollection(int nJet, RVec<float> pt, RVec<float> eta, RVec<float> flav, RVec<float> disc,int var) {
+    //int var 0,1,2 = nom,down,up
+    //evaluate at central, apply 2% unc if var>0
+    RVec<float> recalib_disc(nJet);
+    BTagEntry::JetFlavor fl;
+    for(int i=0; i<nJet;i++){
+        if(flav[i]==5){
+            fl = static_cast<BTagEntry::JetFlavor>(2);//b
+        }
+        else if(flav[i]==4){
+            fl = static_cast<BTagEntry::JetFlavor>(1);//c
+        }
+        else{
+            fl = static_cast<BTagEntry::JetFlavor>(0);//udsg
+        }
+        float sf = reader.eval_auto_bounds("central", fl, eta[i], pt[i], disc[i]);
+        if(var==1){
+            sf = sf*0.98;
+        }
+        else if(var==2){
+            sf = sf * 1.02;
+        }
+        recalib_disc[i] = sf*disc[i];
+    }
+
+    return recalib_disc;
+};

TIMBER/Framework/Zbb_modules/BranchCorrection.cc

@@ -0,0 +1,39 @@
+#include <ROOT/RVec.hxx>
+/**
+ * @class BranchCorrection
+ * @brief Trivial class to load a branch as correction in TIMBER
+ */
+using namespace ROOT::VecOps;
+class BranchCorrection {
+
+    public:
+        BranchCorrection(){};
+        ~BranchCorrection(){};
+        RVec<float> evalCorrection(float val);
+        RVec<float> evalWeight(float val,float valUp,float valDown);
+        RVec<float> evalUncert(float valUp,float valDown);
+
+};
+
+
+RVec<float> BranchCorrection::evalCorrection(float val){
+    RVec<float> correction(1);
+    correction[0]=val;
+    return correction;
+};
+
+RVec<float> BranchCorrection::evalWeight(float val,float valUp,float valDown){
+    RVec<float> weight(3);
+    weight[0]=val;
+    weight[1]=valUp;
+    weight[2]=valDown;
+    return weight;
+};
+
+RVec<float> BranchCorrection::evalUncert(float valUp,float valDown){
+    RVec<float> uncert(2);
+    uncert[0]=valUp;
+    uncert[1]=valDown;
+    return uncert;
+};
+

TIMBER/Framework/Zbb_modules/TTstitching.cc

@@ -0,0 +1,51 @@
+#include <TFile.h>
+#include <TMath.h>
+#include <stdio.h>
+#include <vector>
+#include <iostream>
+#include "ROOT/RVec.hxx"
+
+using namespace ROOT::VecOps;
+using rvec_i = ROOT::VecOps::RVec<int>;
+using rvec_f = ROOT::VecOps::RVec<float>;
+using LVector = ROOT::Math::PtEtaPhiMVector;
+
+Float_t getMTT(Int_t nGenPart, rvec_i GenPart_pdgId, rvec_f GenPart_pt, rvec_f GenPart_phi, rvec_f GenPart_eta, rvec_f GenPart_mass);
+Float_t getPartIdx(Int_t nGenPart, rvec_i GenPart_pdgId, rvec_i GenPart_statusFlags, Int_t pdgId);
+
+Float_t getPartIdx(Int_t nGenPart, rvec_i GenPart_pdgId, rvec_i GenPart_statusFlags, Int_t pdgId){
+    //returns idx of the first hard process parton with GenPart_pdgId == pdgId, -1 otherwise
+    // statusFlags bit 7 : isHardProcess, bit counting starts from zero!
+    for(Int_t i=0;i<nGenPart;i++){
+        if(GenPart_pdgId[i]==pdgId && (GenPart_statusFlags[i]&(1 << 7))){
+            return i;
+        }
+    }
+    return -1;    
+}
+
+Float_t getMTT(Int_t nGenPart, rvec_i GenPart_pdgId, rvec_f GenPart_pt, rvec_f GenPart_phi, rvec_f GenPart_eta, rvec_f GenPart_mass){
+    //Finds two first genParts with pdgId 6 or -6 and calculates their inv mass
+    //Assumes TTbar samples where there will be top and antitop one after the other in the particle chain
+    Int_t nTops = 0;
+    std::vector<LVector> topVecs;
+    Float_t invMass = 0.0;
+
+    for(Int_t i=0;i<nGenPart;i++){
+        if(nTops>1){
+            break;
+        }
+        if(GenPart_pdgId[i]==6 || GenPart_pdgId[i]==-6){
+            topVecs.push_back(LVector(GenPart_pt[i],GenPart_eta[i],GenPart_phi[i],GenPart_mass[i]));
+            nTops=nTops+1;
+        }
+    }
+
+    if(nTops==2){
+        invMass = (topVecs[0]+topVecs[1]).M();
+    }
+
+    return invMass;
+
+}
+

TIMBER/Framework/Zbb_modules/TrigEff.cc

@@ -0,0 +1,87 @@
+#include <string>
+#include "TFile.h"
+#include "TEfficiency.h"
+#include <iostream>
+/**
+ * @class TrigEff
+ * @brief Generic histogram loader with methods to return bin values.
+ */
+class TrigEff {
+    private:
+        TFile *file;
+        TEfficiency *efficiency; //!< Histogram object
+        int binx;
+        int globalbin;
+        float effval;
+        float effup;
+        float effdown;
+
+    public:
+        /**
+         * @brief Empty constructor
+         */
+        TrigEff(){};
+        /**
+         * @brief Construct a new TrigEff object
+         * 
+         * @param filename File to access.
+         * @param histname Histogram name in the file.
+         */
+        TrigEff(std::string filename, std::string histname);
+        /**
+         * @brief Evaluate by global bin number.
+         * 
+         * @param globalbin 
+         * @return std::vector<float> {nominal value, up error+nominal, down error+nominal}
+         */
+        std::vector<float> eval_byglobal(int globalbin);
+        /**
+         * @brief Evaluate by per-axis bin numbers.
+         * 
+         * @param binx 
+         * @param biny
+         * @param binz 
+         * @return std::vector<float> {nominal value, up error+nominal, down error+nominal}
+         */
+        std::vector<float> eval_bybin(int binx, int biny = 0, int binz = 0);
+        /**
+         * @brief Evaluate by axis value.
+         * 
+         * @param xval 
+         * @param yval
+         * @param zval
+         * @return std::vector<float> {nominal value, up error+nominal, down error+nominal}
+         */
+        std::vector<float> eval(float xval, float yval = 0, float zval = 0);
+
+};
+
+TrigEff::TrigEff(std::string filename, std::string effname) {
+    file = TFile::Open(filename.c_str());
+    efficiency = (TEfficiency*)file->Get(effname.c_str());
+}
+
+std::vector<float> TrigEff::eval_byglobal(int globalbin){
+    effval  = efficiency->GetEfficiency(globalbin);
+    effup   = effval + efficiency->GetEfficiencyErrorUp(globalbin)+0.01;//1% additional uncertainty to account for JetHT/SingleMuon differences
+    effdown = effval - efficiency->GetEfficiencyErrorLow(globalbin)-0.01;
+
+    if(effup>1.0){
+        effup = 1.0;
+    }
+
+    if(effdown<0.0){
+        effdown = 0.0;
+    }
+    return {effval,effup,effdown};
+}
+
+std::vector<float> TrigEff::eval_bybin(int binx, int biny, int binz){
+    globalbin    = efficiency->GetGlobalBin(binx, biny, binz);
+    return eval_byglobal(globalbin);
+}
+
+std::vector<float> TrigEff::eval(float xval, float yval, float zval){
+    globalbin    = efficiency->FindFixBin(xval, yval, zval);
+    return eval_byglobal(globalbin);
+}