CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
JetReCalibrator.py
Go to the documentation of this file.
1 import ROOT
2 import os
3 from math import *
4 from PhysicsTools.HeppyCore.utils.deltar import *
5 
6 class JetReCalibrator:
7  def __init__(self,globalTag,jetFlavour,doResidualJECs,jecPath,upToLevel=3):
8  """Create a corrector object that reads the payloads from the text dumps of a global tag under
9  CMGTools/RootTools/data/jec (see the getJec.py there to make the dumps).
10  It will apply the L1,L2,L3 and possibly the residual corrections to the jets."""
11  # Make base corrections
12  path = jecPath #"%s/src/CMGTools/RootTools/data/jec" % os.environ['CMSSW_BASE'];
13  self.L1JetPar = ROOT.JetCorrectorParameters("%s/%s_L1FastJet_%s.txt" % (path,globalTag,jetFlavour),"");
14  self.L2JetPar = ROOT.JetCorrectorParameters("%s/%s_L2Relative_%s.txt" % (path,globalTag,jetFlavour),"");
15  self.L3JetPar = ROOT.JetCorrectorParameters("%s/%s_L3Absolute_%s.txt" % (path,globalTag,jetFlavour),"");
16  self.vPar = ROOT.vector(ROOT.JetCorrectorParameters)()
17  self.vPar.push_back(self.L1JetPar);
18  if upToLevel >= 2: self.vPar.push_back(self.L2JetPar);
19  if upToLevel >= 3: self.vPar.push_back(self.L3JetPar);
20  # Add residuals if needed
21  if doResidualJECs :
22  self.ResJetPar = ROOT.JetCorrectorParameters("%s/%s_L2L3Residual_%s.txt" % (path,globalTag,jetFlavour))
23  self.vPar.push_back(self.ResJetPar);
24  #Step3 (Construct a FactorizedJetCorrector object)
25  self.JetCorrector = ROOT.FactorizedJetCorrector(self.vPar)
26  if os.path.exists("%s/%s_Uncertainty_%s.txt" % (path,globalTag,jetFlavour)):
27  self.JetUncertainty = ROOT.JetCorrectionUncertainty("%s/%s_Uncertainty_%s.txt" % (path,globalTag,jetFlavour));
28  else:
29  print 'Missing JEC uncertainty file "%s/%s_Uncertainty_%s.txt", so jet energy uncertainties will not be available' % (path,globalTag,jetFlavour)
30  self.JetUncertainty = None
31  def correctAll(self,jets,rho,delta=0,metShift=[0,0]):
32  """Applies 'correct' to all the jets, discard the ones that have bad corrections (corrected pt <= 0)"""
33  badJets = []
34  for j in jets:
35  ok = self.correct(j,rho,delta,metShift)
36  if not ok: badJets.append(j)
37  if len(badJets) > 0:
38  print "Warning: %d out of %d jets flagged bad by JEC." % (len(badJets), len(jets))
39  for bj in badJets:
40  jets.remove(bj)
41  def correct(self,jet,rho,delta=0,metShift=[0,0]):
42  """Corrects a jet energy (optionally shifting it also by delta times the JEC uncertainty)
43  If a two-component list is passes as 'metShift', it will be modified adding to the first and second
44  component the change to the MET along x and y due to the JEC, defined as the negative difference
45  between the new and old jet 4-vectors, for jets with corrected pt > 10."""
46  self.JetCorrector.setJetEta(jet.eta())
47  self.JetCorrector.setJetPt(jet.pt() * jet.rawFactor())
48  self.JetCorrector.setJetA(jet.jetArea())
49  self.JetCorrector.setRho(rho)
50  corr = self.JetCorrector.getCorrection()
51  if delta != 0:
52  if not self.JetUncertainty: raise RuntimeError, "Jet energy scale uncertainty shifts requested, but not available"
53  self.JetUncertainty.setJetEta(jet.eta())
54  self.JetUncertainty.setJetPt(corr * jet.pt() * jet.rawFactor())
55  try:
56  jet.jetEnergyCorrUncertainty = self.JetUncertainty.getUncertainty(True)
57  except RuntimeError, r:
58  print "Caught %s when getting uncertainty for jet of pt %.1f, eta %.2f\n" % (r,corr * jet.pt() * jet.rawFactor(),jet.eta())
59  jet.jetEnergyCorrUncertainty = 0.5
60  if jet.photonEnergyFraction() < 0.9 and jet.pt()*corr*jet.rawFactor() > 10:
61  metShift[0] -= jet.px()*(corr*jet.rawFactor() - 1)*(1-jet.muonEnergyFraction())
62  metShift[1] -= jet.py()*(corr*jet.rawFactor() - 1)*(1-jet.muonEnergyFraction())
63  if delta != 0:
64  #print " jet with corr pt %6.2f has an uncertainty %.2f " % (jet.pt()*jet.rawFactor()*corr, jet.jetEnergyCorrUncertainty)
65  corr *= max(0, 1+delta*jet.jetEnergyCorrUncertainty)
66  if jet.pt()*jet.rawFactor()*corr > 10:
67  metShift[0] -= jet.px()*jet.rawFactor()*corr*delta*jet.jetEnergyCorrUncertainty
68  metShift[1] -= jet.py()*jet.rawFactor()*corr*delta*jet.jetEnergyCorrUncertainty
69  #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)
70  if corr <= 0:
71  return False
72  jet.setP4(jet.p4() * (corr * jet.rawFactor()))
73  jet.setRawFactor(1.0/corr)
74  return True
75 
76 class Type1METCorrection:
77  def __init__(self,globalTag,jetFlavour,doResidualJECs):
78  """Create a corrector object that reads the payloads from the text dumps of a global tag under
79  CMGTools/RootTools/data/jec (see the getJec.py there to make the dumps).
80  It will make the Type1 MET."""
81  # Make base corrections
82  path = "%s/src/CMGTools/RootTools/data/jec" % os.environ['CMSSW_BASE'];
83  self.L1JetPar = ROOT.JetCorrectorParameters("%s/%s_L1FastJet_%s.txt" % (path,globalTag,jetFlavour));
84  self.L2JetPar = ROOT.JetCorrectorParameters("%s/%s_L2Relative_%s.txt" % (path,globalTag,jetFlavour));
85  self.L3JetPar = ROOT.JetCorrectorParameters("%s/%s_L3Absolute_%s.txt" % (path,globalTag,jetFlavour));
86  self.vPar = ROOT.vector(ROOT.JetCorrectorParameters)()
87  self.vPar.push_back(self.L1JetPar);
88  self.vPar.push_back(self.L2JetPar);
89  self.vPar.push_back(self.L3JetPar);
90  # Add residuals if needed
91  if doResidualJECs :
92  self.ResJetPar = ROOT.JetCorrectorParameters("%s/%s_L2L3Residual_%s.txt" % (path,globalTag,jetFlavour))
93  self.vPar.push_back(self.ResJetPar);
94  #Step3 (Construct a FactorizedJetCorrector object)
95  self.JetCorrector = ROOT.FactorizedJetCorrector(self.vPar)
96  self.JetUncertainty = ROOT.JetCorrectionUncertainty("%s/%s_Uncertainty_%s.txt" % (path,globalTag,jetFlavour));
97  self.vPar1 = ROOT.vector(ROOT.JetCorrectorParameters)()
98  self.vPar1.push_back(self.L1JetPar);
99  self.JetCorrectorL1 = ROOT.FactorizedJetCorrector(self.vPar1)
100  def getMETCorrection(self,jets,rho,muons):
101  px0, py0 = 0, 0
102  for j in jets:
103  if j.component(4).fraction() > 0.9: continue
104  jp4 = j.p4().__class__(j.p4());
105  jp4 *= j.rawFactor()
106  mup4 = j.p4().__class__(0,0,0,0)
107  for mu in muons:
108  if mu.sourcePtr().track().isNonnull() and (mu.sourcePtr().isGlobalMuon() or mu.sourcePtr().isStandAloneMuon()) and deltaR(mu.eta(),mu.phi(),j.eta(),j.phi()) < 0.5:
109  mup4 += mu.p4()
110  jp4 -= mup4
111  self.JetCorrector.setJetEta(j.eta())
112  self.JetCorrector.setJetPt(j.pt()*j.rawFactor()) # NOTE: we don't subtract the muon here!
113  self.JetCorrector.setJetA(j.jetArea())
114  self.JetCorrector.setRho(rho)
115  corr3 = self.JetCorrector.getCorrection()
116  self.JetCorrectorL1.setJetEta(j.eta())
117  self.JetCorrectorL1.setJetPt(j.pt()*j.rawFactor()) # NOTE: we don't subtract the muon here!
118  self.JetCorrectorL1.setJetA(j.jetArea())
119  self.JetCorrectorL1.setRho(rho)
120  corr1 = self.JetCorrectorL1.getCorrection()
121  if jp4.pt() * corr3 < 10:
122  #print " no correction from jet of pt %8.2f / %8.2f , phi %+5.3f / %+5.3f (corr3 %5.3f, corr1 %5.3f): %+7.3f %+7.3f" % (j.pt(),jp4.pt(),j.phi(), jp4.phi(), corr3, corr1, -jp4.X()*(corr3 -corr1),-jp4.Y()*(corr3-corr1))
123  continue
124  #print " correction from jet of pt %8.2f, phi %+5.3f / %+5.3f (corr3 %5.3f, corr1 %5.3f): %+7.3f %+7.3f" % (j.pt(),j.phi(), jp4.phi(), corr3, corr1, -jp4.X()*(corr3-corr1), -jp4.Y()*(corr3-
125  px0 -= jp4.X()*(corr3-corr1)
126  py0 -= jp4.Y()*(corr3-corr1)
127  return [px0, py0]
__class__
Definition: c++11.h:13
deltaR
double deltaR(double eta1, double eta2, double phi1, double phi2)
Definition: TreeUtility.cc:17
Generated for CMSSW Reference Manual by doxygen 1.8.5