Run3 update

TIMBER/Analyzer.py

@@ -8,6 +8,8 @@
 from TIMBER.Tools.Common import GenerateHash, GetHistBinningTuple, CompileCpp, ConcatCols, GetStandardFlags, ExecuteCmd, LoadColumnNames, ProgressBar
 from clang import cindex
 from collections import OrderedDict
+import array
+import numpy
 
 import ROOT
 import pprint, copy, os, subprocess, textwrap, re, glob
@@ -875,7 +877,7 @@ def _checkCorrections(self,node,correctionNames,dropList):
 
         return correctionsToApply
 
-    def MakeWeightCols(self,name='',node=None,correctionNames=None,dropList=[],correlations=[],extraNominal=''):
+    def MakeWeightCols(self,name='',node=None,correctionNames=None,dropList=[],correlations=[], uncerts_to_corr={},extraNominal=''):
         '''Makes columns/variables to store total weights based on the Corrections that have been added.
 
         This function automates the calculation of the columns that store the nominal weight and the 
@@ -957,6 +959,11 @@ def MakeWeightCols(self,name='',node=None,correctionNames=None,dropList=[],corre
                         weights[corrname+'_up'] += ' * '+correctionName+'__up'
                         weights[corrname+'_down'] += ' * '+correctionName+'__down'
 
+                    for correction in uncerts_to_corr: #From Michael Hesford
+                        if corrname in uncerts_to_corr[correction]: #remove nominal correction from up/down uncert columns
+                            weights[corrname+'_up'] = weights[corrname+'_up'].replace(correction+'__nom * ','')
+                            weights[corrname+'_down'] = weights[corrname+'_down'].replace(correction+'__nom * ','')
+
                 elif corr.GetType() == 'corr':
                     continue
 
@@ -1023,6 +1030,8 @@ def MakeTemplateHistos(self,templateHist,variables,node=None,lazy=True):
             histtitle = '%s__%s'%(baseTitle,cname.replace('weight__','').replace('__nominal',''))
 
             # Build the tuple to give as argument for template
+            if len(binningTuple) == 4:
+                binningTuple = (binningTuple[0], array.array("d", numpy.array(binningTuple[1])), binningTuple[2], array.array("d", numpy.array(binningTuple[3])))
             template_attr = (histname,histtitle) + binningTuple
 
             if dimension == 1: 

TIMBER/Framework/ecalBadCalibFilterRecipe.cc

@@ -0,0 +1,14 @@
+
+#include "../include/ecalBadCalibFilterRecipe.h"
+bool ecalBadCalibFilterRecipe(int run,  float PuppiMET_pt, float PuppiMET_phi,  int nJet, RVec<float> Jet_pt, RVec<float> Jet_eta, RVec<float> Jet_phi, RVec<float> Jet_neEmEF, RVec<float> Jet_chEmEF) {
+    if(run >= 362433 && run <= 367144){
+        if (PuppiMET_pt > 100){
+            for (int i = 0; i < nJet; i++){
+                float deltaPhi = TMath::Abs(Jet_phi.at(i) - PuppiMET_phi) < TMath::Pi() ? TMath::Abs(Jet_phi.at(i) - PuppiMET_phi) : 2*TMath::Pi() - TMath::Abs(Jet_phi.at(i) - PuppiMET_phi);
+                if(Jet_pt.at(i) > 50 && Jet_eta.at(i) > -0.5 && Jet_eta.at(i) < -0.1 && && Jet_phi.at(i) > -2.1 && Jet_phi.at(i) < -1.8 && (Jet_neEmEF.at(i) > 0.9 || Jet_chEmEF.at(i) > 0.9) && deltaPhi > 2.9)
+                    return 0;
+            }
+        }
+    }
+    return 1;
+}

TIMBER/Framework/include/calBadCalibFilterRecipe.h

@@ -0,0 +1,6 @@
+#include <TMath.h>
+#include <ROOT/RVec.hxx>
+using ROOT::VecOps::RVec;
+
+bool calBadCalibFilterRecipe(int run,  float PuppiMET_pt, float PuppiMET_phi,  int nJet, RVec<float> Jet_pt, RVec<float> Jet_eta, RVec<float> Jet_phi, RVec<float> Jet_neEmEF, RVec<float> Jet_chEmEF);
+

TIMBER/Framework/include/common.h

@@ -378,4 +378,4 @@ class TempDir {
         std::string Hash();
 
 };
-#endif
+#endif

TIMBER/Framework/include/ecalBadCalibFilterRecipe.h

@@ -0,0 +1,6 @@
+#include <TMath.h>
+#include <ROOT/RVec.hxx>
+using ROOT::VecOps::RVec;
+
+bool ecalBadCalibFilterRecipe(int run,  float PuppiMET_pt, float PuppiMET_phi,  int nJet, RVec<float> Jet_pt, RVec<float> Jet_eta, RVec<float> Jet_phi, RVec<float> Jet_neEmEF, RVec<float> Jet_chEmEF);
+

TIMBER/Framework/src/JERC_JetVeto.cc

@@ -42,23 +42,18 @@ class JERC_JetVeto {
          *
          * @param jets TIMBER-created structure for the AK4 jet collection
          */
-        template <class T>
-        int eval(std::vector<T> jets) {
-            int veto = 0;   // Start out assuming we don't veto the event
+        int eval(int nJet, RVec<float> Jet_pt, RVec<float> Jet_eta, RVec<float> Jet_phi, RVec<float> Jet_id, RVec<float> Jet_chEmEF, RVec<float> Jet_neEmEF) {
             _total++;       // Increment the number of events studied
-            for (size_t ijet = 0; ijet < jets.size(); ijet++) {
+            for (size_t ijet = 0; ijet < nJet; ijet++) {
                 // First check whether we the jet passes the nominal "loose selection"
-                bool jet_passes = (jets[ijet].pt > 15) && (jets[ijet].jetId == 6) && ((jets[ijet].chEmEF + jets[ijet].neEmEF) < 0.9);
+                bool jet_passes = (Jet_pt[ijet] > 15) && (Jet_id[ijet] >= 6) && ((Jet_chEmEF[ijet] + Jet_neEmEF[ijet]) < 0.9) && abs(Jet_eta[ijet]) < 5.2;
                 if (jet_passes) {   // Consult jet veto map
-                    // First impose checks on valid eta and phi. The eta/phi variables have a certain allowed range in correctionlib and will fail if passed a value outside that range.
-                    if (abs(jets[ijet].eta) > 5.191) {continue;}
-                    if (abs(jets[ijet].phi) > 3.1415926536) {continue;}
 
                     float veto;
                     std::map<std::string, correction::Variable::Type> map {
                         {"type", "jetvetomap"},  // name of the type of veto map. The recommended map for analyses is 'jetvetomap'. Other possible values: jetvetomap, jetvetomap_all, jetvetomap_hbp2m1, jetvetomap_hem1516, jetvetomap_hot
-                        {"eta", jets[ijet].eta},
-                        {"phi", jets[ijet].phi}
+                        {"eta", Jet_eta[ijet]},
+                        {"phi", Jet_phi[ijet]}
                     };
                     correction::Correction::Ref ref = _cset->at(_key);
                     std::vector<correction::Variable::Type> inputs;
@@ -72,7 +67,7 @@ class JERC_JetVeto {
                     }
                 }
             }
-            return veto;    // Return 0 if no vetoed jet found
+            return 0;    // Return 0 if no vetoed jet found
         };
 };
 
@@ -88,4 +83,4 @@ JERC_JetVeto::~JERC_JetVeto() {
     std::cout << "[JERC_JetVeto] Fraction of events vetoed using jet veto maps and nominal jet selection = " << fraction << "%" << std::endl;
     std::cout << "\tTotal number of events analyzed: " << _total << std::endl;
     std::cout << "\tTotal number of events vetoed:   " << _vetoed << std::endl;
-}
+}

TIMBER/Framework/src/JER_correctionlib_weight.cc

@@ -1,5 +1,6 @@
 // Requires CMSSW
-
+#ifndef JER_CORR
+#define JER_CORR
 #include <correction.h>
 #include <ROOT/RVec.hxx>
 #include <Math/GenVector/LorentzVector.h>
@@ -15,19 +16,19 @@ class JER_correctionlib_weight {
     private:
         std::unique_ptr<correction::CorrectionSet> _cset;
         // JSON keys for each of the correction objects
-        std::string _key_jes;
         std::string _key_res;
-        std::string _key_sf;        
+        std::string _key_sf;
+        std::string _fname;        
         // Variables for JER correction algorithm
         float _dRMax, _dPtMaxFactor;
         std::mt19937 _rnd;
-        static constexpr const double MIN_JET_ENERGY = 1e-2;
+        //static constexpr const double MIN_JET_ENERGY = 1e-2;
         // Variables to keep track of GEN matching info
         float _total = 0.;
         float _matched = 0.; 
 
     public:
-        JER_correctionlib_weight(std::string fname, std::string key_jes, std::string key_res, std::string key_sf, float dRMax, float dPtMaxFactor = 3.);
+        JER_correctionlib_weight(std::string fname, std::string key_res, std::string key_sf, float dRMax, float dPtMaxFactor = 3.);
         ~JER_correctionlib_weight();
 
         /**
@@ -38,106 +39,12 @@ class JER_correctionlib_weight {
          * @param resolution pT resolution to consider.
          * @return LorentzV The gen jet that matches the reco jet.
          */
-        LorentzV match(LorentzV& jet, RVec<LorentzV> genJets, float resolution);
-
+        float match( float pt, float eta, float phi, int nGenJet, RVec<float> genJet_pt, RVec<float> genJet_eta, RVec<float> genJet_phi, float resolution);
 
-        /**
-         * @brief Main function called by JER correction class.
-         *
-         * @param jets TIMBER-created structure for the jet collection (AK4/AK8). Can be a custom collection name
-         * @param genJets TIMBER-created structure for the generator jet collection. By default, this is automatically created in TIMBER as "GenJet"
-         * @param fixedGridRhoFastjetAll NanoAOD column of the same name storing the density rho
-         */
-        template <class Tjet, class TgenJet>
-        RVec<RVec<float>> eval(std::vector<Tjet> jets, std::vector<TgenJet> genJets, float fixedGridRhoFastjetAll){
-            RVec< RVec<float> > out (jets.size());
-            for (size_t ijet = 0; ijet < jets.size(); ijet++) {
-                // Book {nom, up, down} SFs for the jet at index "ijet"
-                RVec<float> ijet_out {1.0, 1.0, 1.0};
-
-                // Get the resolution for this jet given its eta, pt, and event rho
-                float res;
-                std::map<std::string, correction::Variable::Type> resMap {
-                    {"JetPt", jets[ijet].pt},
-                    {"JetEta", jets[ijet].eta},
-                    {"Rho", fixedGridRhoFastjetAll}
-                };
-                correction::Correction::Ref ref_res = _cset->at(_key_res);
-                std::vector<correction::Variable::Type> inputs_res;
-                for (const correction::Variable& input : ref_res->inputs()) {
-                    inputs_res.push_back(resMap.at(input.name()));
-                }
-                res = ref_res->evaluate(inputs_res);
-
-                // Now attempt to match the RECO jet to a GEN jet
-                LorentzV reco_jet = hardware::TLvector(jets[ijet]);
-                LorentzV genJet = match(reco_jet, hardware::TLvector(genJets), jets[ijet].pt * res);
-
-                // Keep track of the number of total and matched jets to print upon object destruction
-                _total++;
-                if (genJet.Pt() > -1) {_matched++;}
-
-                // Now begin a loop over the variations {0: "nom", 1: "up", 2: "down"}
-                float smearFactor, dpT, sigma, jet_sf;
-                for (size_t i=0; i<3; i++) {
-                    // Determine the variation
-                    std::string variation;
-                    if (i == 0) {
-                        variation = "nom";
-                    }
-                    else if (i == 1) {
-                        variation = "up";
-                    }
-                    else {
-                        variation = "down";
-                    }
-                    // Get the JER scale factor
-                    float jet_sf;
-                    std::map<std::string, correction::Variable::Type> sfMap {
-                        {"JetEta", jets[ijet].eta},
-                         {"JetPt", jets[ijet].pt},
-                        {"systematic", variation}
-                    };
-                    correction::Correction::Ref ref_sf = _cset->at(_key_sf);
-                    std::vector<correction::Variable::Type> inputs_sf;
-                    for (const correction::Variable& input : ref_sf->inputs()) {
-                        inputs_sf.push_back(sfMap.at(input.name()));
-                    }
-                    jet_sf = ref_sf->evaluate(inputs_sf);
-
-                    // Now determine how to handle the various cases.
-                    //      Case 1: we have a "good" gen jet matched to the reco jet.
-                    //      Case 2: we don't have a gen jet. Smear jet pt using a random gaussian variation.
-                    if (genJet.Pt() > -1) { // Case 1
-                        dpT = reco_jet.Pt() - genJet.Pt();
-                        smearFactor = 1. + (jet_sf -1.) * dpT / reco_jet.Pt();
-                    }
-                    else if (jet_sf > 1) {  // Case 2
-                        std::normal_distribution<> d(0, res);
-                        smearFactor = 1. + d(_rnd) * std::sqrt(jet_sf * jet_sf - 1.);
-                    }
-                    else {
-                        smearFactor = 1.;
-                    }
-
-                    if (reco_jet.E() * smearFactor < MIN_JET_ENERGY) {
-                        // Negative or too small smearFactor. We would change direction of the jet
-                        // and this is not what we want.
-                        // Recompute the smearing factor in order to have jet.energy() == MIN_JET_ENERGY
-                        smearFactor = MIN_JET_ENERGY / reco_jet.E();
-                    }
-
-                    ijet_out[i] = smearFactor;
-
-                } // end loop over variations
-                // Return {nom,up,down} values for the jet at index "ijet"
-                out[ijet] = ijet_out;
-            }
-            return out;
-        }
+        RVec<RVec<float>> eval(int nJet, RVec<float> jet_pt, RVec<float> jet_eta, RVec<float> jet_phi, int nGenJet, RVec<float> genJet_pt, RVec<float> genJet_eta, RVec<float> genJet_phi, float fixedGridRhoFastjetAll);
 };
 
-JER_correctionlib_weight::JER_correctionlib_weight(std::string fname, std::string key_jes, std::string key_res, std::string key_sf, float dRMax, float dPtMaxFactor) : _dRMax(dRMax), _dPtMaxFactor(dPtMaxFactor), _key_jes(key_jes), _key_res(key_res), _key_sf(key_sf) {
+JER_correctionlib_weight::JER_correctionlib_weight(std::string fname, std::string key_res, std::string key_sf, float dRMax, float dPtMaxFactor) : _dRMax(dRMax), _dPtMaxFactor(dPtMaxFactor), _key_res(key_res), _key_sf(key_sf), _fname(fname) {
     _cset = correction::CorrectionSet::from_file(fname.c_str());
 }
 
@@ -147,27 +54,115 @@ JER_correctionlib_weight::~JER_correctionlib_weight() {
     //       not be triggered. If you see -nan% and did not trigger the RDF event loop, don't be alarmed!!
     float fraction = _matched/_total * 100.;
     std::cout << "[JER_correctionlib_weight] Gen matching efficiency = " << fraction << "%" << std::endl;
+    std::cout << "\tTEST"<<_fname << std::endl;
     std::cout << "\tTotal number of RECO jets matched to a GEN jet: " << _matched << std::endl;
     std::cout << "\tTotal number of RECO jets analyzed:             " << _total << std::endl;
 }
 
-LorentzV JER_correctionlib_weight::match(LorentzV& jet, RVec<LorentzV> genJets, float resolution) {
-    // Match if dR < _dRMax and dPt < dPtMaxFactor
-    double min_dR = std::numeric_limits<double>::infinity();
-    LorentzV out (-1,0,0,0);
 
-    for (const LorentzV & genJet : genJets) {
-        float dR = hardware::DeltaR(genJet, jet);
+RVec<RVec<float>> JER_correctionlib_weight::eval( int nJet, RVec<float> jet_pt, RVec<float> jet_eta, RVec<float> jet_phi, int nGenJet, RVec<float> genJet_pt, RVec<float> genJet_eta, RVec<float> genJet_phi, float fixedGridRhoFastjetAll){
+    RVec< RVec<float> > out (nJet);
+    for (size_t ijet = 0; ijet < nJet; ijet++) {
+        // Book {nom, up, down} SFs for the jet at index "ijet"
+        RVec<float> ijet_out {1.0, 1.0, 1.0};
+
+        // Get the resolution for this jet given its eta, pt, and event rho
+        float res;
+        std::map<std::string, correction::Variable::Type> resMap {
+            {"JetPt", jet_pt[ijet]},
+            {"JetEta", jet_eta[ijet]},
+            {"Rho", fixedGridRhoFastjetAll}
+        };
+        correction::Correction::Ref ref_res = _cset->at(_key_res);
+        std::vector<correction::Variable::Type> inputs_res;
+        for (const correction::Variable& input : ref_res->inputs()) {
+            inputs_res.push_back(resMap.at(input.name()));
+        }
+        res = ref_res->evaluate(inputs_res);
+
+        // Now attempt to match the RECO jet to a GEN jet
+        float genPt = match(jet_pt[ijet], jet_eta[ijet], jet_phi[ijet], nGenJet, genJet_pt, genJet_eta, genJet_phi, jet_pt[ijet] * res);
+
+        // Keep track of the number of total and matched jets to print upon object destruction
+        _total++;
+        if (genPt >= 0) {_matched++;}
+        // Now begin a loop over the variations {0: "nom", 1: "up", 2: "down"}
+        float smearFactor, dpT, sigma, jet_sf;
+        for (size_t i=0; i<3; i++) {
+            // Determine the variation
+            std::string variation;
+            if (i == 0) {
+                variation = "nom";
+            }
+            else if (i == 1) {
+                variation = "up";
+            }
+            else {
+                variation = "down";
+            }
+            // Get the JER scale factor
+            float jet_sf;
+            std::map<std::string, correction::Variable::Type> sfMap {
+                {"JetEta", jet_eta[ijet]},
+                 {"JetPt", jet_pt[ijet]},
+                {"systematic", variation}
+            };
+            correction::Correction::Ref ref_sf = _cset->at(_key_sf);
+            std::vector<correction::Variable::Type> inputs_sf;
+            for (const correction::Variable& input : ref_sf->inputs()) {
+                inputs_sf.push_back(sfMap.at(input.name()));
+            }
+            jet_sf = ref_sf->evaluate(inputs_sf);
+
+            // Now determine how to handle the various cases.
+            //      Case 1: we have a "good" gen jet matched to the reco jet.
+            //      Case 2: we don't have a gen jet. Smear jet pt using a random gaussian variation.
+            if (genPt >= 0) { // Case 1
+                dpT = jet_pt.at(ijet) - genPt;
+                smearFactor = 1. + (jet_sf -1.) * dpT / jet_pt.at(ijet);
+            }
+            else if (jet_sf > 1) {  // Case 2
+                std::normal_distribution<> d(0, res);
+                smearFactor = 1. + d(_rnd) * std::sqrt(jet_sf * jet_sf - 1.);
+            }
+            else {
+                smearFactor = 1.;
+            }
+            smearFactor = smearFactor < 0 ? 0 : smearFactor;
+
+            ijet_out[i] = smearFactor;
+
+        } // end loop over variations
+        // Return {nom,up,down} values for the jet at index "ijet"
+        out[ijet] = ijet_out;
+    }
+    return out;
+}
+
+float JER_correctionlib_weight::match( float pt, float eta, float phi, int nGenJet, RVec<float> genJet_pt, RVec<float> genJet_eta, RVec<float> genJet_phi, float resolution) {
+    double min_dR = std::numeric_limits<double>::infinity();
+    // Match if dR < _dRMax and dPt < dPtMaxFactor
+    float out_pt = -1.;
+    for (int i = 0; i < nGenJet; i ++){
+        float deltaEta = std::abs(eta - genJet_eta.at(i));
+        float deltaPhi = std::abs(phi - genJet_phi.at(i)) < M_PI ? std::abs(phi - genJet_phi.at(i)) : 2*M_PI - std::abs(phi - genJet_phi.at(i));
+        float dR = sqrt(deltaEta * deltaEta + deltaPhi * deltaPhi);
         if (dR > min_dR) {
             continue;
         }
         if (dR < _dRMax) {
-            double dPt = std::abs(genJet.pt() - jet.pt());
+            double dPt = std::abs(pt - genJet_pt.at(i));
             if ((resolution == -1) || (dPt <= _dPtMaxFactor * resolution)) {
                 min_dR = dR;
-                out = genJet;
+                out_pt = genJet_pt.at(i);
             }
         }
     }
-    return out;
+    return out_pt;
 }
+
+
+
+
+
+#endif