Merged tth-heppy-745-v2 from repository jpata

PhysicsTools/Heppy/python/analyzers/objects/JetAnalyzer.py

@@ -8,6 +8,7 @@
 
 from PhysicsTools.Heppy.physicsutils.QGLikelihoodCalculator import QGLikelihoodCalculator
 
+import copy
 def cleanNearestJetOnly(jets,leptons,deltaR):
     dr2 = deltaR**2
     good = [ True for j in jets ]
@@ -101,6 +102,12 @@ def process(self, event):
         if self.doJEC:
 #            print "\nCalibrating jets %s for lumi %d, event %d" % (self.cfg_ana.jetCol, event.lumi, event.eventId)
             self.jetReCalibrator.correctAll(allJets, rho, delta=self.shiftJEC, metShift=self.deltaMetFromJEC)
+
+        for delta, shift in [(1.0, "JECUp"), (0.0, ""), (-1.0, "JECDown")]:
+            for j1 in allJets:
+                corr = self.jetReCalibrator.getCorrection(j1, rho, delta, self.deltaMetFromJEC)
+                setattr(j1, "corr"+shift, corr)
+
         self.allJetsUsedForMET = allJets
 #        print "after. rho",self.rho,self.cfg_ana.collectionPostFix,'allJets len ',len(allJets),'pt', [j.pt() for j in allJets]
 
@@ -373,7 +380,7 @@ def matchJets(self, event, jets):
             jet.mcJet = match[jet]
 
 
-
+  
     def smearJets(self, event, jets):
         # https://twiki.cern.ch/twiki/bin/viewauth/CMS/TWikiTopRefSyst#Jet_energy_resolution
        for jet in jets:

PhysicsTools/Heppy/python/analyzers/objects/autophobj.py

@@ -165,6 +165,9 @@
     NTupleVariable("mcPt",   lambda x : x.mcJet.pt() if getattr(x,"mcJet",None) else 0., mcOnly=True, help="p_{T} of associated gen jet"),
     NTupleVariable("mcFlavour", lambda x : x.partonFlavour(), int,     mcOnly=True, help="parton flavour (physics definition, i.e. including b's from shower)"),
     NTupleVariable("mcMatchId",  lambda x : getattr(x, 'mcMatchId', -99), int, mcOnly=True, help="Match to source from hard scatter (pdgId of heaviest particle in chain, 25 for H, 6 for t, 23/24 for W/Z), zero if non-prompt or fake"),
+    NTupleVariable("corr_JECUp",  lambda x : getattr(x, 'corrJECUp', -99), float, mcOnly=True, help=""),
+    NTupleVariable("corr_JECDown",  lambda x : getattr(x, 'corrJECDown', -99), float, mcOnly=True, help=""),
+    NTupleVariable("corr",  lambda x : getattr(x, 'corr', -99), float, mcOnly=True, help=""),
 ])
 jetTypeExtra = NTupleObjectType("jetExtra",  baseObjectTypes = [ jetType ], variables = [
     NTupleVariable("area",   lambda x : x.jetArea(), help="Catchment area of jet"),

PhysicsTools/Heppy/python/physicsutils/BTagWeightCalculator.py

@@ -0,0 +1,98 @@
+import ROOT
+import numpy as np
+
+class BTagWeightCalculator:
+    def __init__(self, fn_hf, fn_lf) :
+        self.pdfs = {}
+
+        self.pt_bins_hf = np.array([20, 30, 40, 60, 100, 160, 10000])
+        self.eta_bins_hf = np.array([0, 2.41])
+
+        self.pt_bins_lf = np.array([20, 30, 40, 60, 10000])
+        self.eta_bins_lf = np.array([0, 0.8, 1.6, 2.41])
+
+        self.btag = "combinedInclusiveSecondaryVertexV2BJetTags"
+        self.init(fn_hf, fn_lf)
+
+    def getBin(self, bvec, val):
+        return int(bvec.searchsorted(val) - 1)
+
+    def init(self, fn_hf, fn_lf) :
+        print "[BTagWeightCalculator]: Initializing from files", fn_hf, fn_lf
+
+        self.pdfs["hf"] = self.getHistosFromFile(fn_hf)
+        self.pdfs["lf"] = self.getHistosFromFile(fn_lf)
+
+        return True
+
+    def getHistosFromFile(self, fn):
+        ret = {}
+        tf = ROOT.TFile(fn)
+        if not tf or tf.IsZombie():
+            raise FileError("Could not open file {0}".format(fn))
+        ROOT.gROOT.cd()
+        for k in tf.GetListOfKeys():
+            kn = k.GetName()
+            if not kn.startswith("csv_ratio"):
+                continue
+            spl = kn.split("_")
+
+            if spl[2] == "all":
+                ptbin = -1
+                etabin = -1
+                kind = "all"
+                syst = "nominal"
+            else:
+                ptbin = int(spl[2][2:])
+                etabin = int(spl[3][3:])
+                kind = spl[4]
+                if len(spl)==6:
+                    syst = spl[5]
+                else:
+                    syst = "nominal"
+            ret[(ptbin, etabin, kind, syst)] = k.ReadObj().Clone()
+        return ret
+
+    def calcJetWeight(self, jet, kind, systematic):
+        pt = jet.pt()
+        aeta = abs(jet.eta())
+        fl = abs(jet.mcFlavour)
+        csv = jet.btag(self.btag)
+
+        is_heavy = (fl == 5 or fl == 6)
+        if is_heavy:
+            ptbin = self.getBin(self.pt_bins_hf, pt)
+            etabin = self.getBin(self.eta_bins_hf, aeta)
+        else:
+            ptbin = self.getBin(self.pt_bins_lf, pt)
+            etabin = self.getBin(self.eta_bins_lf, aeta)
+
+        if ptbin < 0 or etabin < 0:
+            return 1.0
+
+        k = (ptbin, etabin, kind, systematic)
+        hdict = self.pdfs["lf"]
+        if is_heavy:
+            hdict = self.pdfs["hf"]
+        h = hdict.get(k, None)
+        if not h:
+            return 1.0
+
+        csvbin = 1
+        if csv>=0:
+            csvbin = h.FindBin(csv)
+
+        if csvbin <= 0 or csvbin > h.GetNbinsX():
+            return 1.0
+
+        w = h.GetBinContent(csvbin)
+        return w
+
+    def calcEventWeight(self, jets, kind, systematic):
+        weights = np.array(
+            [self.calcJetWeight(jet, kind, systematic)
+            for jet in jets]
+        )
+
+        wtot = np.prod(weights)
+        return wtot

PhysicsTools/Heppy/python/physicsutils/JetReCalibrator.py

@@ -25,6 +25,8 @@ def __init__(self,globalTag,jetFlavour,doResidualJECs,jecPath,upToLevel=3):
         self.JetCorrector = ROOT.FactorizedJetCorrector(self.vPar)
         if os.path.exists("%s/%s_Uncertainty_%s.txt" % (path,globalTag,jetFlavour)):
             self.JetUncertainty = ROOT.JetCorrectionUncertainty("%s/%s_Uncertainty_%s.txt" % (path,globalTag,jetFlavour));
+        elif os.path.exists("%s/Uncertainty_FAKE.txt" % path):
+            self.JetUncertainty = ROOT.JetCorrectionUncertainty("%s/Uncertainty_FAKE.txt" % path);
         else:
             print 'Missing JEC uncertainty file "%s/%s_Uncertainty_%s.txt", so jet energy uncertainties will not be available' % (path,globalTag,jetFlavour)
             self.JetUncertainty = None
@@ -38,11 +40,10 @@ def correctAll(self,jets,rho,delta=0,metShift=[0,0]):
             print "Warning: %d out of %d jets flagged bad by JEC." % (len(badJets), len(jets))
         for bj in badJets:
             jets.remove(bj)
-    def correct(self,jet,rho,delta=0,metShift=[0,0]):
-        """Corrects a jet energy (optionally shifting it also by delta times the JEC uncertainty)
-           If a two-component list is passes as 'metShift', it will be modified adding to the first and second
-           component the change to the MET along x and y due to the JEC, defined as the negative difference 
-           between the new and old jet 4-vectors, for jets with corrected pt > 10."""
+
+    def getCorrection(self,jet,rho,delta=0,metShift=[0,0]):
+        """Calculates the correction factor of a jet without modifying it
+        """
         self.JetCorrector.setJetEta(jet.eta())
         self.JetCorrector.setJetPt(jet.pt() * jet.rawFactor())
         self.JetCorrector.setJetA(jet.jetArea())
@@ -67,6 +68,14 @@ def correct(self,jet,rho,delta=0,metShift=[0,0]):
                 metShift[0] -= jet.px()*jet.rawFactor()*corr*delta*jet.jetEnergyCorrUncertainty
                 metShift[1] -= jet.py()*jet.rawFactor()*corr*delta*jet.jetEnergyCorrUncertainty
         #print "   jet with raw pt %6.2f eta %+5.3f phi %+5.3f: previous corr %.4f, my corr %.4f " % (jet.pt()*jet.rawFactor(), jet.eta(), jet.phi(), 1./jet.rawFactor(), corr)
+        return corr
+
+    def correct(self,jet,rho,delta=0,metShift=[0,0]):
+        """Corrects a jet energy (optionally shifting it also by delta times the JEC uncertainty)
+           If a two-component list is passes as 'metShift', it will be modified adding to the first and second
+           component the change to the MET along x and y due to the JEC, defined as the negative difference
+           between the new and old jet 4-vectors, for jets with corrected pt > 10."""
+        corr = self.getCorrection(jet,rho,delta,metShift)
         if corr <= 0:
             return False
         jet.setP4(jet.p4() * (corr * jet.rawFactor()))