Merge pull request cms-sw#62 from arizzi/VhbbToColin

Thanks VHbb team!

PhysicsTools/Heppy/python/analyzers/core/AutoFillTreeProducer.py

@@ -125,7 +125,7 @@ def fillCoreVariables(self, tr, event, isMC):
             ## PU weights, check if a PU analyzer actually filled it
             if hasattr(event,"nPU"):
                     tr.fill("nTrueInt", event.nPU)
-                    tr.fill("puWeight", event.eventWeight)
+                    tr.fill("puWeight", event.puWeight)
             else :
                     tr.fill("nTrueInt", -1)
                     tr.fill("puWeight", 1.0)

PhysicsTools/Heppy/python/analyzers/core/PileUpAnalyzer.py

@@ -104,7 +104,7 @@ def process(self, event):
         if self.cfg_comp.isEmbed :
           return True
 
-        event.vertexWeight = 1
+        event.puWeight = 1
         event.nPU = None
         event.pileUpVertex_z = []
         event.pileUpVertex_ptHat = []
@@ -140,18 +140,18 @@ def process(self, event):
         if self.enable:
             bin = self.datahist.FindBin(event.nPU)
             if bin<1 or bin>self.datahist.GetNbinsX():
-                event.vertexWeight = 0
+                event.puWeight = 0
             else:
                 data = self.datahist.GetBinContent(bin)
                 mc = self.mchist.GetBinContent(bin)
                 #Protect 0 division!!!!
                 if mc !=0.0:
-                    event.vertexWeight = data/mc
+                    event.puWeight = data/mc
                 else:
-                    event.vertexWeight = 1
+                    event.puWeight = 1
 
-        event.eventWeight *= event.vertexWeight
-        self.averages['vertexWeight'].add( event.vertexWeight )
+        event.eventWeight *= event.puWeight
+        self.averages['puWeight'].add( event.puWeight )
         return True
 
     def write(self, setup):

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

@@ -47,6 +47,21 @@ def cleanJetsAndLeptons(jets,leptons,deltaR,arbitration):
              [ l for (il,l) in enumerate(leptons) if goodlep[il] == True ] )
 
 
+
+def shiftJERfactor(JERShift, aeta):
+        factor = 1.079 + JERShift*0.026
+        if   aeta > 3.2: factor = 1.056 + JERShift * 0.191
+        elif aeta > 2.8: factor = 1.395 + JERShift * 0.063
+        elif aeta > 2.3: factor = 1.254 + JERShift * 0.062
+        elif aeta > 1.7: factor = 1.208 + JERShift * 0.046
+        elif aeta > 1.1: factor = 1.121 + JERShift * 0.029
+        elif aeta > 0.5: factor = 1.099 + JERShift * 0.028
+        return factor 
+
+
+
+
+
 class JetAnalyzer( Analyzer ):
     """Taken from RootTools.JetAnalyzer, simplified, modified, added corrections    """
     def __init__(self, cfg_ana, cfg_comp, looperName):
@@ -55,7 +70,7 @@ def __init__(self, cfg_ana, cfg_comp, looperName):
         dataGT = cfg_ana.dataGT if hasattr(cfg_ana,'dataGT') else "GR_70_V2_AN1"
         self.shiftJEC = self.cfg_ana.shiftJEC if hasattr(self.cfg_ana, 'shiftJEC') else 0
         self.recalibrateJets = self.cfg_ana.recalibrateJets
-        self.addJECShifts = self.cfg_ana.addJECShifts
+        self.addJECShifts = self.cfg_ana.addJECShifts if hasattr(self.cfg_ana, 'addJECShifts') else 0
         if   self.recalibrateJets == "MC"  : self.recalibrateJets =     self.cfg_comp.isMC
         elif self.recalibrateJets == "Data": self.recalibrateJets = not self.cfg_comp.isMC
         elif self.recalibrateJets not in [True,False]: raise RuntimeError, "recalibrateJets must be any of { True, False, 'MC', 'Data' }, while it is %r " % self.recalibrateJets
@@ -86,6 +101,7 @@ def declareHandles(self):
         self.handles['jets']   = AutoHandle( self.cfg_ana.jetCol, 'std::vector<pat::Jet>' )
         self.handles['genJet'] = AutoHandle( self.cfg_ana.genJetCol, 'vector<reco::GenJet>' )
         self.shiftJER = self.cfg_ana.shiftJER if hasattr(self.cfg_ana, 'shiftJER') else 0
+        self.addJERShifts = self.cfg_ana.addJERShifts if hasattr(self.cfg_ana, 'addJERShifts') else 0
         self.handles['rho'] = AutoHandle( self.cfg_ana.rho, 'double' )
 
     def beginLoop(self, setup):
@@ -122,6 +138,9 @@ def process(self, event):
                 self.matchJets(event, allJets)
             if getattr(self.cfg_ana, 'smearJets', False):
                 self.smearJets(event, allJets)
+
+
+
 
 	##Sort Jets by pT 
         allJets.sort(key = lambda j : j.pt(), reverse = True)
@@ -365,7 +384,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:
@@ -374,22 +393,25 @@ def smearJets(self, event, jets):
                genpt, jetpt, aeta = gen.pt(), jet.pt(), abs(jet.eta())
                # from https://twiki.cern.ch/twiki/bin/view/CMS/JetResolution
                #8 TeV tables
-
-               factor = 1.079 + self.shiftJER*0.026
-               if   aeta > 3.2: factor = 1.056 + self.shiftJER * 0.191
-               elif aeta > 2.8: factor = 1.395 + self.shiftJER * 0.063
-               elif aeta > 2.3: factor = 1.254 + self.shiftJER * 0.062
-               elif aeta > 1.7: factor = 1.208 + self.shiftJER * 0.046
-               elif aeta > 1.1: factor = 1.121 + self.shiftJER * 0.029
-               elif aeta > 0.5: factor = 1.099 + self.shiftJER * 0.028
-               ptscale = max(0.0, (jetpt + (factor-1)*(jetpt-genpt))/jetpt)
+               factor = shiftJERfactor(self.shiftJER, aeta)
+               ptscale = max(0.0, (jetpt + (factor-1)*(jetpt-genpt))/jetpt)             
                #print "get with pt %.1f (gen pt %.1f, ptscale = %.3f)" % (jetpt,genpt,ptscale)
                jet.deltaMetFromJetSmearing = [ -(ptscale-1)*jet.rawFactor()*jet.px(), -(ptscale-1)*jet.rawFactor()*jet.py() ]
                if ptscale != 0:
                   jet.setP4(jet.p4()*ptscale)
                # leave the uncorrected unchanged for sync
                jet._rawFactorMultiplier *= (1.0/ptscale) if ptscale != 0 else 1
             #else: print "jet with pt %.1d, eta %.2f is unmatched" % (jet.pt(), jet.eta())
+               if (self.shiftJER==0) and (self.addJERShifts):
+                   setattr(jet, "corrJER", ptscale )
+                   factorJERUp= shiftJERfactor(1, aeta)
+                   ptscaleJERUp = max(0.0, (jetpt + (factorJERUp-1)*(jetpt-genpt))/jetpt)
+                   setattr(jet, "corrJERUp", ptscaleJERUp)
+                   factorJERDown= shiftJERfactor(-1, aeta)
+                   ptscaleJERDown = max(0.0, (jetpt + (factorJERDown-1)*(jetpt-genpt))/jetpt)
+                   setattr(jet, "corrJERDown", ptscaleJERDown)
+
+
 
 
 
@@ -417,6 +439,7 @@ def smearJets(self, event, jets):
     addJECShifts = False, # if true, add  "corr", "corrJECUp", and "corrJECDown" for each jet (requires uncertainties to be available!)
     smearJets = True,
     shiftJER = 0, # set to +1 or -1 to get +/-1 sigma shifts    
+    addJERShifts = 0, # add +/-1 sigma shifts to jets, intended to be used with shiftJER=0
     cleanJetsFromFirstPhoton = False,
     cleanJetsFromTaus = False,
     cleanJetsFromIsoTracks = False,

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

@@ -176,9 +176,12 @@
     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=""),
+    NTupleVariable("corr_JECUp",  lambda x : getattr(x, 'corrJECUp', -99), float,  help=""),
+    NTupleVariable("corr_JECDown",  lambda x : getattr(x, 'corrJECDown', -99), float,help=""),
+    NTupleVariable("corr",  lambda x : getattr(x, 'corr', -99), float, help=""),
+    NTupleVariable("corr_JERUp",  lambda x : getattr(x, 'corrJERUp', -99), float, mcOnly=True, help=""),
+    NTupleVariable("corr_JERDown", lambda x : getattr(x, 'corrJERDown', -99), float, mcOnly=True, help=""),
+    NTupleVariable("corr_JER", lambda x : getattr(x, 'corrJER', -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

@@ -33,34 +33,51 @@ def getHistosFromFile(self, fn):
         ROOT.gROOT.cd()
         for k in tf.GetListOfKeys():
             kn = k.GetName()
-            if not kn.startswith("csv_ratio"):
+            if not (kn.startswith("csv_ratio") or kn.startswith("c_csv_ratio") ):
                 continue
             spl = kn.split("_")
+            is_c = 1 if kn.startswith("c_csv_ratio") else 0
 
-            if spl[2] == "all":
+            if spl[2+is_c] == "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]
+                ptbin = int(spl[2+is_c][2:])
+                etabin = int(spl[3+is_c][3:])
+                kind = spl[4+is_c]
+                if len(spl)==(6+is_c):
+                    syst = spl[5+is_c]
                 else:
                     syst = "nominal"
             ret[(ptbin, etabin, kind, syst)] = k.ReadObj().Clone()
         return ret
 
     def calcJetWeight(self, jet, kind, systematic):
-        pt = jet.pt()
+        pt   = jet.pt()
         aeta = abs(jet.eta())
-        fl = abs(jet.mcFlavour)
-        csv = jet.btag(self.btag)
+        fl   = abs(jet.hadronFlavour())
+        csv  = jet.btag(self.btag)
 
-        is_heavy = (fl == 5 or fl == 6)
-        if is_heavy:
+        is_b = (fl == 5)
+        is_c = (fl == 4)
+        is_l = (fl < 4)
+
+        if is_b and not (systematic in ["JESUp", "JESDown", "LFUp", "LFDown", 
+                                        "Stats1Up", "Stats1Down", "Stats2Up", "Stats2Down", 
+                                        "nominal"]):
+            return 1.0
+        if is_c and not (systematic in ["cErr1Up", "cErr1Down", "cErr2Up", "cErr2Down", 
+                                        "nominal"]):
+            return 1.0
+        if is_l and not (systematic in ["JESUp", "JESDown", "HFUp", "HFDown", 
+                                        "Stats1Up", "Stats1Down", "Stats2Up", "Stats2Down", 
+                                        "nominal"]):
+            return 1.0
+
+
+        if is_b or is_c:
             ptbin = self.getBin(self.pt_bins_hf, pt)
             etabin = self.getBin(self.eta_bins_hf, aeta)
         else:
@@ -72,7 +89,7 @@ def calcJetWeight(self, jet, kind, systematic):
 
         k = (ptbin, etabin, kind, systematic)
         hdict = self.pdfs["lf"]
-        if is_heavy:
+        if is_b or is_c:
             hdict = self.pdfs["hf"]
         h = hdict.get(k, None)
         if not h: