Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Jet Energy Scale systematics can be split by sources and grouped #750

Open
wants to merge 3 commits into
base: heppy_94X_dev
Choose a base branch
from
Open

Jet Energy Scale systematics can be split by sources and grouped #750

Changes from 2 commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
465e94e
jet calibration can add specific shifts by uncertainty source
GaelTouquet Mar 25, 2019
1347266
harmonizing naming
GaelTouquet Apr 5, 2019
a22f2eb
fix source grouping when only 1 source in JEC
GaelTouquet Apr 16, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
49 changes: 37 additions & 12 deletions PhysicsTools/Heppy/python/physicsutils/JetReCalibrator.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,8 @@
class JetReCalibrator:
def __init__(self,globalTag,jetFlavour,doResidualJECs,jecPath,upToLevel=3,
calculateSeparateCorrections=False,
calculateType1METCorrection=False, type1METParams={'jetPtThreshold':15., 'skipEMfractionThreshold':0.9, 'skipMuons':True} ):
calculateType1METCorrection=False, type1METParams={'jetPtThreshold':15., 'skipEMfractionThreshold':0.9, 'skipMuons':True},
groupForUncertaintySources = {}):
"""Create a corrector object that reads the payloads from the text dumps of a global tag under
CMGTools/RootTools/data/jec (see the getJec.py there to make the dumps).
It will apply the L1,L2,L3 and possibly the residual corrections to the jets.
Expand All @@ -18,6 +19,7 @@ def __init__(self,globalTag,jetFlavour,doResidualJECs,jecPath,upToLevel=3,
self.upToLevel = upToLevel
self.calculateType1METCorr = calculateType1METCorrection
self.type1METParams = type1METParams
self.groupForUncertaintySources = groupForUncertaintySources
# Make base corrections
path = os.path.expandvars(jecPath) #"%s/src/CMGTools/RootTools/data/jec" % os.environ['CMSSW_BASE'];
self.L1JetPar = ROOT.JetCorrectorParameters("%s/%s_L1FastJet_%s.txt" % (path,globalTag,jetFlavour),"");
Expand All @@ -35,6 +37,12 @@ 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));
if os.path.exists("%s/%s_Uncertainty%s_%s.txt" % (path,globalTag,'Sources',jetFlavour)) and self.groupForUncertaintySources:
self.JetUncertaintyBySources = {}
for group, sources in groupForUncertaintySources.iteritems():
for source in sources:
params = ROOT.JetCorrectorParameters("%s/%s_UncertaintySources_%s.txt" % (path,globalTag,jetFlavour),source)
self.JetUncertaintyBySources[source] = ROOT.JetCorrectionUncertainty(params)
elif os.path.exists("%s/Uncertainty_FAKE.txt" % path):
self.JetUncertainty = ROOT.JetCorrectionUncertainty("%s/Uncertainty_FAKE.txt" % path);
else:
Expand All @@ -58,7 +66,7 @@ def __init__(self,globalTag,jetFlavour,doResidualJECs,jecPath,upToLevel=3,
for i in [self.L1JetPar,self.L2JetPar,self.L3JetPar,self.ResJetPar]: self.vParL3Res.push_back(i)
self.separateJetCorrectors["L1L2L3Res"] = ROOT.FactorizedJetCorrector(self.vParL3Res)

def getCorrection(self,jet,rho,delta=0,corrector=None):
def getCorrection(self,jet,rho,delta=0,corrector=None, sources=[]):
if not corrector: corrector = self.JetCorrector
if corrector != self.JetCorrector and delta!=0: raise RuntimeError('Configuration not supported')
corrector.setJetEta(jet.eta())
Expand All @@ -67,16 +75,27 @@ def getCorrection(self,jet,rho,delta=0,corrector=None):
corrector.setRho(rho)
corr = corrector.getCorrection()
if delta != 0:
if not self.JetUncertainty: raise RuntimeError("Jet energy scale uncertainty shifts requested, but not available")
self.JetUncertainty.setJetEta(jet.eta())
self.JetUncertainty.setJetPt(corr * jet.pt() * jet.rawFactor())
try:
jet.jetEnergyCorrUncertainty = self.JetUncertainty.getUncertainty(True)
except RuntimeError as r:
print "Caught %s when getting uncertainty for jet of pt %.1f, eta %.2f\n" % (r,corr * jet.pt() * jet.rawFactor(),jet.eta())
jet.jetEnergyCorrUncertainty = 0.5
#print " jet with corr pt %6.2f has an uncertainty %.2f " % (jet.pt()*jet.rawFactor()*corr, jet.jetEnergyCorrUncertainty)
corr *= max(0, 1+delta*jet.jetEnergyCorrUncertainty)
if sources:
group_shift = 0.
for source in sources:
if source not in self.JetUncertaintyBySources:
raise RuntimeError("Jet energy scale uncertainty source requested, but not available")
self.JetUncertaintyBySources[source].setJetEta(jet.eta())
self.JetUncertaintyBySources[source].setJetPt(corr * jet.pt() * jet.rawFactor())
unc = self.JetUncertaintyBySources[source].getUncertainty(True if delta>0. else False)
group_shift += unc*unc
GaelTouquet marked this conversation as resolved.
Show resolved Hide resolved
group_shift = sqrt(group_shift)
corr *= max(0, 1+delta*group_shift)
else:
if not self.JetUncertainty: raise RuntimeError("Jet energy scale uncertainty shifts requested, but not available")
self.JetUncertainty.setJetEta(jet.eta())
self.JetUncertainty.setJetPt(corr * jet.pt() * jet.rawFactor())
try:
jet.jetEnergyCorrUncertainty = self.JetUncertainty.getUncertainty(True if delta>0. else False)
except RuntimeError as r:
print "Caught %s when getting uncertainty for jet of pt %.1f, eta %.2f\n" % (r,corr * jet.pt() * jet.rawFactor(),jet.eta())
jet.jetEnergyCorrUncertainty = 0.5
corr *= max(0, 1+delta*jet.jetEnergyCorrUncertainty)
return corr

def rawP4forType1MET_(self, jet):
Expand Down Expand Up @@ -125,6 +144,12 @@ def correct(self,jet,rho,delta=0,addCorr=False,addShifts=False, metShift=None, t
for cdelta,shift in [(1.0, "JECUp"), (-1.0, "JECDown")]:
cshift = self.getCorrection(jet,rho,delta+cdelta)
setattr(jet, "corr"+shift, cshift)
if self.groupForUncertaintySources:
for group, sources in self.groupForUncertaintySources.iteritems():
shift_up = self.getCorrection(jet,rho,delta+1.0,sources=sources)
shift_down = self.getCorrection(jet,rho,delta-1.0,sources=sources)
setattr(jet, "corr_"+group+"_JEC_up", shift_up)
setattr(jet, "corr_"+group+"_JEC_down", shift_down)
if corr <= 0:
return False
newpt = jet.pt()*raw*corr
Expand Down