test
CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
JetAnalyzer.py
Go to the documentation of this file.
1 import math, os
2 from PhysicsTools.Heppy.analyzers.core.Analyzer import Analyzer
3 from PhysicsTools.Heppy.analyzers.core.AutoHandle import AutoHandle
4 from PhysicsTools.Heppy.physicsobjects.PhysicsObjects import Jet
5 from PhysicsTools.HeppyCore.utils.deltar import deltaR2, deltaPhi, matchObjectCollection, matchObjectCollection2, bestMatch,matchObjectCollection3
6 from PhysicsTools.Heppy.physicsutils.JetReCalibrator import JetReCalibrator
7 import PhysicsTools.HeppyCore.framework.config as cfg
8 
9 from PhysicsTools.Heppy.physicsutils.QGLikelihoodCalculator import QGLikelihoodCalculator
10 
11 import copy
12 def cleanNearestJetOnly(jets,leptons,deltaR):
13  dr2 = deltaR**2
14  good = [ True for j in jets ]
15  for l in leptons:
16  ibest, d2m = -1, dr2
17  for i,j in enumerate(jets):
18  d2i = deltaR2(l.eta(),l.phi(), j.eta(),j.phi())
19  if d2i < d2m:
20  ibest, d2m = i, d2i
21  if ibest != -1: good[ibest] = False
22  return [ j for (i,j) in enumerate(jets) if good[i] == True ]
23 
24 def cleanJetsAndLeptons(jets,leptons,deltaR,arbitration):
25  dr2 = deltaR**2
26  goodjet = [ True for j in jets ]
27  goodlep = [ True for l in leptons ]
28  for il, l in enumerate(leptons):
29  ibest, d2m = -1, dr2
30  for i,j in enumerate(jets):
31  d2i = deltaR2(l.eta(),l.phi(), j.eta(),j.phi())
32  if d2i < dr2:
33  choice = arbitration(j,l)
34  if choice == j:
35  # if the two match, and we prefer the jet, then drop the lepton and be done
36  goodlep[il] = False
37  break
38  elif choice == (j,l) or choice == (l,j):
39  # asked to keep both, so we don't consider this match
40  continue
41  if d2i < d2m:
42  ibest, d2m = i, d2i
43  # this lepton has been killed by a jet, then we clean the jet that best matches it
44  if not goodlep[il]: continue
45  if ibest != -1: goodjet[ibest] = False
46  return ( [ j for (i ,j) in enumerate(jets) if goodjet[i ] == True ],
47  [ l for (il,l) in enumerate(leptons) if goodlep[il] == True ] )
48 
49 
50 
51 def shiftJERfactor(JERShift, aeta):
52  factor = 1.079 + JERShift*0.026
53  if aeta > 3.2: factor = 1.056 + JERShift * 0.191
54  elif aeta > 2.8: factor = 1.395 + JERShift * 0.063
55  elif aeta > 2.3: factor = 1.254 + JERShift * 0.062
56  elif aeta > 1.7: factor = 1.208 + JERShift * 0.046
57  elif aeta > 1.1: factor = 1.121 + JERShift * 0.029
58  elif aeta > 0.5: factor = 1.099 + JERShift * 0.028
59  return factor
60 
61 
62 
63 
64 
65 class JetAnalyzer( Analyzer ):
66  """Taken from RootTools.JetAnalyzer, simplified, modified, added corrections """
67  def __init__(self, cfg_ana, cfg_comp, looperName):
68  super(JetAnalyzer,self).__init__(cfg_ana, cfg_comp, looperName)
69  mcGT = cfg_ana.mcGT if hasattr(cfg_ana,'mcGT') else "PHYS14_25_V2"
70  dataGT = cfg_ana.dataGT if hasattr(cfg_ana,'dataGT') else "GR_70_V2_AN1"
71  self.shiftJEC = self.cfg_ana.shiftJEC if hasattr(self.cfg_ana, 'shiftJEC') else 0
72  self.recalibrateJets = self.cfg_ana.recalibrateJets
73  self.addJECShifts = self.cfg_ana.addJECShifts if hasattr(self.cfg_ana, 'addJECShifts') else 0
74  if self.recalibrateJets == "MC" : self.recalibrateJets = self.cfg_comp.isMC
75  elif self.recalibrateJets == "Data": self.recalibrateJets = not self.cfg_comp.isMC
76  elif self.recalibrateJets not in [True,False]: raise RuntimeError, "recalibrateJets must be any of { True, False, 'MC', 'Data' }, while it is %r " % self.recalibrateJets
77 
78  calculateSeparateCorrections = getattr(cfg_ana,"calculateSeparateCorrections", False);
79  calculateType1METCorrection = getattr(cfg_ana,"calculateType1METCorrection", False);
80  self.doJEC = self.recalibrateJets or (self.shiftJEC != 0) or self.addJECShifts or calculateSeparateCorrections or calculateType1METCorrection
81  if self.doJEC:
82  doResidual = getattr(cfg_ana, 'applyL2L3Residual', 'Data')
83  if doResidual == "MC": doResidual = self.cfg_comp.isMC
84  elif doResidual == "Data": doResidual = not self.cfg_comp.isMC
85  elif doResidual not in [True,False]: raise RuntimeError, "If specified, applyL2L3Residual must be any of { True, False, 'MC', 'Data'(default)}"
86  GT = getattr(cfg_comp, 'jecGT', mcGT if self.cfg_comp.isMC else dataGT)
87  # Now take care of the optional arguments
88  kwargs = { 'calculateSeparateCorrections':calculateSeparateCorrections,
89  'calculateType1METCorrection' :calculateType1METCorrection, }
90  if kwargs['calculateType1METCorrection']: kwargs['type1METParams'] = cfg_ana.type1METParams
91  # instantiate the jet re-calibrator
92  self.jetReCalibrator = JetReCalibrator(GT, cfg_ana.recalibrationType, doResidual, cfg_ana.jecPath, **kwargs)
93  self.doPuId = getattr(self.cfg_ana, 'doPuId', True)
94  self.jetLepDR = getattr(self.cfg_ana, 'jetLepDR', 0.4)
95  self.jetLepArbitration = getattr(self.cfg_ana, 'jetLepArbitration', lambda jet,lepton: lepton)
96  self.lepPtMin = getattr(self.cfg_ana, 'minLepPt', -1)
97  self.lepSelCut = getattr(self.cfg_ana, 'lepSelCut', lambda lep : True)
98  self.jetGammaDR = getattr(self.cfg_ana, 'jetGammaDR', 0.4)
99  self.jetGammaLepDR = getattr(self.cfg_ana, 'jetGammaLepDR', 0.4)
100  self.cleanFromLepAndGammaSimultaneously = getattr(self.cfg_ana, 'cleanFromLepAndGammaSimultaneously', False)
101  if self.cleanFromLepAndGammaSimultaneously:
102  if hasattr(self.cfg_ana, 'jetGammaLepDR'):
103  self.jetGammaLepDR = self.jetGammaLepDR
104  elif (self.jetGammaDR == self.jetLepDR):
105  self.jetGammaLepDR = self.jetGammaDR
106  else:
107  raise RuntimeError, "DR for simultaneous cleaning of jets from leptons and photons is not defined, and dR(gamma, jet)!=dR(lep, jet)"
108  if(self.cfg_ana.doQG):
109  qgdefname="{CMSSW_BASE}/src/PhysicsTools/Heppy/data/pdfQG_AK4chs_13TeV_v2b.root"
110  self.qglcalc = QGLikelihoodCalculator(getattr(self.cfg_ana,"QGpath",qgdefname).format(CMSSW_BASE= os.environ['CMSSW_BASE']))
111  if not hasattr(self.cfg_ana ,"collectionPostFix"):self.cfg_ana.collectionPostFix=""
112 
113  def declareHandles(self):
114  super(JetAnalyzer, self).declareHandles()
115  self.handles['jets'] = AutoHandle( self.cfg_ana.jetCol, 'std::vector<pat::Jet>' )
116  self.handles['genJet'] = AutoHandle( self.cfg_ana.genJetCol, 'vector<reco::GenJet>' )
117  self.shiftJER = self.cfg_ana.shiftJER if hasattr(self.cfg_ana, 'shiftJER') else 0
118  self.addJERShifts = self.cfg_ana.addJERShifts if hasattr(self.cfg_ana, 'addJERShifts') else 0
119  self.handles['rho'] = AutoHandle( self.cfg_ana.rho, 'double' )
120 
121  def beginLoop(self, setup):
122  super(JetAnalyzer,self).beginLoop(setup)
123 
124  def process(self, event):
125  self.readCollections( event.input )
126  rho = float(self.handles['rho'].product()[0])
127  self.rho = rho
128 
129  ## Read jets, if necessary recalibrate and shift MET
130  if self.cfg_ana.copyJetsByValue:
131  import ROOT
132  #from ROOT.heppy import JetUtils
133  allJets = map(lambda j:Jet(ROOT.heppy.JetUtils.copyJet(j)), self.handles['jets'].product()) #copy-by-value is safe if JetAnalyzer is ran more than once
134  else:
135  allJets = map(Jet, self.handles['jets'].product())
136 
137  #set dummy MC flavour for all jets in case we want to ntuplize discarded jets later
138  for jet in allJets:
139  jet.mcFlavour = 0
140 
141  self.deltaMetFromJEC = [0.,0.]
142  self.type1METCorr = [0.,0.,0.]
143 # print "before. rho",self.rho,self.cfg_ana.collectionPostFix,'allJets len ',len(allJets),'pt', [j.pt() for j in allJets]
144  if self.doJEC:
145  if not self.recalibrateJets: # check point that things won't change
146  jetsBefore = [ (j.pt(),j.eta(),j.phi(),j.rawFactor()) for j in allJets ]
147  self.jetReCalibrator.correctAll(allJets, rho, delta=self.shiftJEC,
148  addCorr=True, addShifts=self.addJECShifts,
149  metShift=self.deltaMetFromJEC, type1METCorr=self.type1METCorr )
150  if not self.recalibrateJets:
151  jetsAfter = [ (j.pt(),j.eta(),j.phi(),j.rawFactor()) for j in allJets ]
152  if len(jetsBefore) != len(jetsAfter):
153  print "ERROR: I had to recompute jet corrections, and they rejected some of the jets:\nold = %s\n new = %s\n" % (jetsBefore,jetsAfter)
154  else:
155  for (told, tnew) in zip(jetsBefore,jetsAfter):
156  if (deltaR2(told[1],told[2],tnew[1],tnew[2])) > 0.0001:
157  print "ERROR: I had to recompute jet corrections, and one jet direction moved: old = %s, new = %s\n" % (told, tnew)
158  elif abs(told[0]-tnew[0])/(told[0]+tnew[0]) > 0.5e-3 or abs(told[3]-tnew[3]) > 0.5e-3:
159  print "ERROR: I had to recompute jet corrections, and one jet pt or corr changed: old = %s, new = %s\n" % (told, tnew)
160  self.allJetsUsedForMET = allJets
161 # print "after. rho",self.rho,self.cfg_ana.collectionPostFix,'allJets len ',len(allJets),'pt', [j.pt() for j in allJets]
162 
163  if self.cfg_comp.isMC:
164  self.genJets = [ x for x in self.handles['genJet'].product() ]
165  if self.cfg_ana.do_mc_match:
166  for igj, gj in enumerate(self.genJets):
167  gj.index = igj
168 # self.matchJets(event, allJets)
169  self.matchJets(event, [ j for j in allJets if j.pt()>self.cfg_ana.jetPt ]) # To match only jets above chosen threshold
170  if getattr(self.cfg_ana, 'smearJets', False):
171  self.smearJets(event, allJets)
172 
173 
174 
175 
176  ##Sort Jets by pT
177  allJets.sort(key = lambda j : j.pt(), reverse = True)
178 
179  leptons = []
180  if hasattr(event, 'selectedLeptons'):
181  leptons = [ l for l in event.selectedLeptons if l.pt() > self.lepPtMin and self.lepSelCut(l) ]
182  if self.cfg_ana.cleanJetsFromTaus and hasattr(event, 'selectedTaus'):
183  leptons = leptons[:] + event.selectedTaus
184  if self.cfg_ana.cleanJetsFromIsoTracks and hasattr(event, 'selectedIsoCleanTrack'):
185  leptons = leptons[:] + event.selectedIsoCleanTrack
186 
187  ## Apply jet selection
188  self.jets = []
189  self.jetsFailId = []
190  self.jetsAllNoID = []
191  self.jetsIdOnly = []
192  for jet in allJets:
193  #Check if lepton and jet have overlapping PF candidates
194  leps_with_overlaps = []
195  if getattr(self.cfg_ana, 'checkLeptonPFOverlap', True):
196  for i in range(jet.numberOfSourceCandidatePtrs()):
197  p1 = jet.sourceCandidatePtr(i) #Ptr<Candidate> p1
198  for lep in leptons:
199  for j in range(lep.numberOfSourceCandidatePtrs()):
200  p2 = lep.sourceCandidatePtr(j)
201  has_overlaps = p1.key() == p2.key() and p1.refCore().id().productIndex() == p2.refCore().id().productIndex() and p1.refCore().id().processIndex() == p2.refCore().id().processIndex()
202  if has_overlaps:
203  leps_with_overlaps += [lep]
204  if len(leps_with_overlaps)>0:
205  for lep in leps_with_overlaps:
206  lep.jetOverlap = jet
207  if self.testJetNoID( jet ):
208  self.jetsAllNoID.append(jet)
209  if(self.cfg_ana.doQG):
210  jet.qgl_calc = self.qglcalc.computeQGLikelihood
211  jet.qgl_rho = rho
212  if self.testJetID( jet ):
213  self.jets.append(jet)
214  self.jetsIdOnly.append(jet)
215  else:
216  self.jetsFailId.append(jet)
217  elif self.testJetID (jet ):
218  self.jetsIdOnly.append(jet)
219 
220  jetsEtaCut = [j for j in self.jets if abs(j.eta()) < self.cfg_ana.jetEta ]
221  self.cleanJetsAll, cleanLeptons = cleanJetsAndLeptons(jetsEtaCut, leptons, self.jetLepDR, self.jetLepArbitration)
222 
223  self.cleanJets = [j for j in self.cleanJetsAll if abs(j.eta()) < self.cfg_ana.jetEtaCentral ]
224  self.cleanJetsFwd = [j for j in self.cleanJetsAll if abs(j.eta()) >= self.cfg_ana.jetEtaCentral ]
225  self.discardedJets = [j for j in self.jets if j not in self.cleanJetsAll]
226  if hasattr(event, 'selectedLeptons') and self.cfg_ana.cleanSelectedLeptons:
227  event.discardedLeptons = [ l for l in leptons if l not in cleanLeptons ]
228  event.selectedLeptons = [ l for l in event.selectedLeptons if l not in event.discardedLeptons ]
229  for lep in leptons:
230  if hasattr(lep, "jetOverlap"):
231  if lep.jetOverlap in self.cleanJetsAll:
232  #print "overlap reco", lep.p4().pt(), lep.p4().eta(), lep.p4().phi(), lep.jetOverlap.p4().pt(), lep.jetOverlap.p4().eta(), lep.jetOverlap.p4().phi()
233  lep.jetOverlapIdx = self.cleanJetsAll.index(lep.jetOverlap)
234  elif lep.jetOverlap in self.discardedJets:
235  #print "overlap discarded", lep.p4().pt(), lep.p4().eta(), lep.p4().phi(), lep.jetOverlap.p4().pt(), lep.jetOverlap.p4().eta(), lep.jetOverlap.p4().phi()
236  lep.jetOverlapIdx = 1000 + self.discardedJets.index(lep.jetOverlap)
237 
238  ## First cleaning, then Jet Id
239  self.noIdCleanJetsAll, cleanLeptons = cleanJetsAndLeptons(self.jetsAllNoID, leptons, self.jetLepDR, self.jetLepArbitration)
240  self.noIdCleanJets = [j for j in self.noIdCleanJetsAll if abs(j.eta()) < self.cfg_ana.jetEtaCentral ]
241  self.noIdCleanJetsFwd = [j for j in self.noIdCleanJetsAll if abs(j.eta()) >= self.cfg_ana.jetEtaCentral ]
242  self.noIdDiscardedJets = [j for j in self.jetsAllNoID if j not in self.noIdCleanJetsAll]
243 
244  ## Clean Jets from photons (first cleaning, then Jet Id)
245  photons = []
246  if hasattr(event, 'selectedPhotons'):
247  if self.cfg_ana.cleanJetsFromFirstPhoton:
248  photons = event.selectedPhotons[:1]
249  else:
250  photons = [ g for g in event.selectedPhotons ]
251 
252  self.gamma_cleanJetaAll = []
253  self.gamma_noIdCleanJetsAll = []
254 
255  if self.cleanFromLepAndGammaSimultaneously:
256  self.gamma_cleanJetsAll = cleanNearestJetOnly(jetsEtaCut, photons+leptons, self.jetGammaLepDR)
257  self.gamma_noIdCleanJetsAll = cleanNearestJetOnly(self.jetsAllNoID, photons+leptons, self.jetGammaLepDR)
258  else:
259  self.gamma_cleanJetsAll = cleanNearestJetOnly(self.cleanJetsAll, photons, self.jetGammaDR)
260  self.gamma_noIdCleanJetsAll = cleanNearestJetOnly(self.noIdCleanJetsAll, photons, self.jetGammaDR)
261 
262  self.gamma_cleanJets = [j for j in self.gamma_cleanJetsAll if abs(j.eta()) < self.cfg_ana.jetEtaCentral ]
263  self.gamma_cleanJetsFwd = [j for j in self.gamma_cleanJetsAll if abs(j.eta()) >= self.cfg_ana.jetEtaCentral ]
264 
265  self.gamma_noIdCleanJets = [j for j in self.gamma_noIdCleanJetsAll if abs(j.eta()) < self.cfg_ana.jetEtaCentral ]
266  self.gamma_noIdCleanJetsFwd = [j for j in self.gamma_noIdCleanJetsAll if abs(j.eta()) >= self.cfg_ana.jetEtaCentral ]
267  ###
268 
269  if self.cfg_ana.alwaysCleanPhotons:
270  self.cleanJets = self.gamma_cleanJets
271  self.cleanJetsAll = self.gamma_cleanJetsAll
272  self.cleanJetsFwd = self.gamma_cleanJetsFwd
273  #
274  self.noIdCleanJets = self.gamma_noIdCleanJets
275  self.noIdCleanJetsAll = self.gamma_noIdCleanJetsAll
276  self.noIdCleanJetsFwd = self.gamma_noIdCleanJetsFwd
277 
278  ## Jet Id, after jet/lepton cleaning
279  self.cleanJetsFailIdAll = []
280  for jet in self.noIdCleanJetsAll:
281  if not self.testJetID( jet ):
282  self.cleanJetsFailIdAll.append(jet)
283 
284  self.cleanJetsFailId = [j for j in self.cleanJetsFailIdAll if abs(j.eta()) < self.cfg_ana.jetEtaCentral ]
285 
286  ## Jet Id, after jet/photon cleaning
287  self.gamma_cleanJetsFailIdAll = []
288  for jet in self.gamma_noIdCleanJetsAll:
289  if not self.testJetID( jet ):
290  self.gamma_cleanJetsFailIdAll.append(jet)
291 
292  self.gamma_cleanJetsFailId = [j for j in self.gamma_cleanJetsFailIdAll if abs(j.eta()) < self.cfg_ana.jetEtaCentral ]
293 
294  ## Associate jets to leptons
295  incleptons = event.inclusiveLeptons if hasattr(event, 'inclusiveLeptons') else event.selectedLeptons
296  jlpairs = matchObjectCollection(incleptons, allJets, self.jetLepDR**2)
297 
298  for jet in allJets:
299  jet.leptons = [l for l in jlpairs if jlpairs[l] == jet ]
300  for lep in incleptons:
301  jet = jlpairs[lep]
302  if jet is None:
303  setattr(lep,"jet"+self.cfg_ana.collectionPostFix,lep)
304  else:
305  setattr(lep,"jet"+self.cfg_ana.collectionPostFix,jet)
306  ## Associate jets to taus
307  taus = getattr(event,'selectedTaus',[])
308  jtaupairs = matchObjectCollection( taus, allJets, self.jetLepDR**2)
309 
310  for jet in allJets:
311  jet.taus = [l for l in jtaupairs if jtaupairs[l] == jet ]
312  for tau in taus:
313  setattr(tau,"jet"+self.cfg_ana.collectionPostFix,jtaupairs[tau])
314 
315  #MC stuff
316  if self.cfg_comp.isMC:
317  self.deltaMetFromJetSmearing = [0, 0]
318  for j in self.cleanJetsAll:
319  if hasattr(j, 'deltaMetFromJetSmearing'):
320  self.deltaMetFromJetSmearing[0] += j.deltaMetFromJetSmearing[0]
321  self.deltaMetFromJetSmearing[1] += j.deltaMetFromJetSmearing[1]
322 
323  self.cleanGenJets = cleanNearestJetOnly(self.genJets, leptons, self.jetLepDR)
324 
325  if self.cfg_ana.cleanGenJetsFromPhoton:
326  self.cleanGenJets = cleanNearestJetOnly(self.cleanGenJets, photons, self.jetLepDR)
327 
328  if getattr(self.cfg_ana, 'attachNeutrinos', True) and hasattr(self.cfg_ana,"genNuSelection") :
329  jetNus=[x for x in event.genParticles if abs(x.pdgId()) in [12,14,16] and self.cfg_ana.genNuSelection(x) ]
330  pairs= matchObjectCollection (jetNus, self.genJets, 0.4**2)
331 
332  for (nu,genJet) in pairs.iteritems() :
333  if genJet is not None :
334  if not hasattr(genJet,"nu") :
335  genJet.nu=nu.p4()
336  else :
337  genJet.nu+=nu.p4()
338 
339 
340  if self.cfg_ana.do_mc_match:
341  self.jetFlavour(event)
342 
343  if hasattr(event,"jets"+self.cfg_ana.collectionPostFix): raise RuntimeError("Event already contains a jet collection with the following postfix: "+self.cfg_ana.collectionPostFix)
344  setattr(event,"rho" +self.cfg_ana.collectionPostFix, self.rho )
345  setattr(event,"deltaMetFromJEC" +self.cfg_ana.collectionPostFix, self.deltaMetFromJEC )
346  setattr(event,"type1METCorr" +self.cfg_ana.collectionPostFix, self.type1METCorr )
347  setattr(event,"allJetsUsedForMET" +self.cfg_ana.collectionPostFix, self.allJetsUsedForMET )
348  setattr(event,"jets" +self.cfg_ana.collectionPostFix, self.jets )
349  setattr(event,"jetsFailId" +self.cfg_ana.collectionPostFix, self.jetsFailId )
350  setattr(event,"jetsAllNoID" +self.cfg_ana.collectionPostFix, self.jetsAllNoID )
351  setattr(event,"jetsIdOnly" +self.cfg_ana.collectionPostFix, self.jetsIdOnly )
352  setattr(event,"cleanJetsAll" +self.cfg_ana.collectionPostFix, self.cleanJetsAll )
353  setattr(event,"cleanJets" +self.cfg_ana.collectionPostFix, self.cleanJets )
354  setattr(event,"cleanJetsFwd" +self.cfg_ana.collectionPostFix, self.cleanJetsFwd )
355  setattr(event,"cleanJetsFailIdAll" +self.cfg_ana.collectionPostFix, self.cleanJetsFailIdAll )
356  setattr(event,"cleanJetsFailId" +self.cfg_ana.collectionPostFix, self.cleanJetsFailId )
357  setattr(event,"discardedJets" +self.cfg_ana.collectionPostFix, self.discardedJets )
358  setattr(event,"gamma_cleanJetsAll" +self.cfg_ana.collectionPostFix, self.gamma_cleanJetsAll )
359  setattr(event,"gamma_cleanJets" +self.cfg_ana.collectionPostFix, self.gamma_cleanJets )
360  setattr(event,"gamma_cleanJetsFwd" +self.cfg_ana.collectionPostFix, self.gamma_cleanJetsFwd )
361  setattr(event,"gamma_cleanJetsFailIdAll" +self.cfg_ana.collectionPostFix, self.gamma_cleanJetsFailIdAll )
362  setattr(event,"gamma_cleanJetsFailId" +self.cfg_ana.collectionPostFix, self.gamma_cleanJetsFailId )
363 
364 
365  if self.cfg_comp.isMC:
366  setattr(event,"deltaMetFromJetSmearing"+self.cfg_ana.collectionPostFix, self.deltaMetFromJetSmearing)
367  setattr(event,"cleanGenJets" +self.cfg_ana.collectionPostFix, self.cleanGenJets )
368  setattr(event,"genJets" +self.cfg_ana.collectionPostFix, self.genJets )
369  if self.cfg_ana.do_mc_match:
370  setattr(event,"bqObjects" +self.cfg_ana.collectionPostFix, self.bqObjects )
371  setattr(event,"cqObjects" +self.cfg_ana.collectionPostFix, self.cqObjects )
372  setattr(event,"partons" +self.cfg_ana.collectionPostFix, self.partons )
373  setattr(event,"heaviestQCDFlavour" +self.cfg_ana.collectionPostFix, self.heaviestQCDFlavour )
374 
375 
376  return True
377 
378 
379 
380  def testJetID(self, jet):
381  jet.puJetIdPassed = jet.puJetId()
382  jet.pfJetIdPassed = jet.jetID('POG_PFID_Loose')
383  if self.cfg_ana.relaxJetId:
384  return True
385  else:
386  return jet.pfJetIdPassed and (jet.puJetIdPassed or not(self.doPuId))
387 
388  def testJetNoID( self, jet ):
389  # 2 is loose pile-up jet id
390  return jet.pt() > self.cfg_ana.jetPt and \
391  abs( jet.eta() ) < self.cfg_ana.jetEta;
392 
393  def jetFlavour(self,event):
394  def isFlavour(x,f):
395  id = abs(x.pdgId())
396  if id > 999: return (id/1000)%10 == f
397  if id > 99: return (id/100)%10 == f
398  return id % 100 == f
399 
400 
401 
402  self.bqObjects = [ p for p in event.genParticles if (p.status() == 2 and isFlavour(p,5)) ]
403  self.cqObjects = [ p for p in event.genParticles if (p.status() == 2 and isFlavour(p,4)) ]
404 
405  self.partons = [ p for p in event.genParticles if ((p.status() == 23 or p.status() == 3) and abs(p.pdgId())>0 and (abs(p.pdgId()) in [1,2,3,4,5,21]) ) ]
406  match = matchObjectCollection2(self.cleanJetsAll,
407  self.partons,
408  deltaRMax = 0.3)
409 
410  for jet in self.cleanJetsAll:
411  parton = match[jet]
412  jet.partonId = (parton.pdgId() if parton != None else 0)
413  jet.partonMotherId = (parton.mother(0).pdgId() if parton != None and parton.numberOfMothers()>0 else 0)
414 
415  for jet in self.jets:
416  (bmatch, dr) = bestMatch(jet, self.bqObjects)
417  if dr < 0.4:
418  jet.mcFlavour = 5
419  else:
420  (cmatch, dr) = bestMatch(jet, self.cqObjects)
421  if dr < 0.4:
422  jet.mcFlavour = 4
423  else:
424  jet.mcFlavour = 0
425 
426  self.heaviestQCDFlavour = 5 if len(self.bqObjects) else (4 if len(self.cqObjects) else 1);
427 
428  def matchJets(self, event, jets):
429  match = matchObjectCollection2(jets,
430  event.genbquarks + event.genwzquarks,
431  deltaRMax = 0.3)
432  for jet in jets:
433  gen = match[jet]
434  jet.mcParton = gen
435  jet.mcMatchId = (gen.sourceId if gen != None else 0)
436  jet.mcMatchFlav = (abs(gen.pdgId()) if gen != None else 0)
437 
438  match = matchObjectCollection2(jets,
439  self.genJets,
440  deltaRMax = 0.3)
441  for jet in jets:
442  jet.mcJet = match[jet]
443 
444 
445 
446  def smearJets(self, event, jets):
447  # https://twiki.cern.ch/twiki/bin/viewauth/CMS/TWikiTopRefSyst#Jet_energy_resolution
448  for jet in jets:
449  gen = jet.mcJet
450  if gen != None:
451  genpt, jetpt, aeta = gen.pt(), jet.pt(), abs(jet.eta())
452  # from https://twiki.cern.ch/twiki/bin/view/CMS/JetResolution
453  #8 TeV tables
454  factor = shiftJERfactor(self.shiftJER, aeta)
455  ptscale = max(0.0, (jetpt + (factor-1)*(jetpt-genpt))/jetpt)
456  #print "get with pt %.1f (gen pt %.1f, ptscale = %.3f)" % (jetpt,genpt,ptscale)
457  jet.deltaMetFromJetSmearing = [ -(ptscale-1)*jet.rawFactor()*jet.px(), -(ptscale-1)*jet.rawFactor()*jet.py() ]
458  if ptscale != 0:
459  jet.setP4(jet.p4()*ptscale)
460  # leave the uncorrected unchanged for sync
461  jet.setRawFactor(jet.rawFactor()/ptscale)
462  #else: print "jet with pt %.1d, eta %.2f is unmatched" % (jet.pt(), jet.eta())
463  if (self.shiftJER==0) and (self.addJERShifts):
464  setattr(jet, "corrJER", ptscale )
465  factorJERUp= shiftJERfactor(1, aeta)
466  ptscaleJERUp = max(0.0, (jetpt + (factorJERUp-1)*(jetpt-genpt))/jetpt)
467  setattr(jet, "corrJERUp", ptscaleJERUp)
468  factorJERDown= shiftJERfactor(-1, aeta)
469  ptscaleJERDown = max(0.0, (jetpt + (factorJERDown-1)*(jetpt-genpt))/jetpt)
470  setattr(jet, "corrJERDown", ptscaleJERDown)
471 
472 
473 
474 
475 
476 setattr(JetAnalyzer,"defaultConfig", cfg.Analyzer(
477  class_object = JetAnalyzer,
478  jetCol = 'slimmedJets',
479  copyJetsByValue = False, #Whether or not to copy the input jets or to work with references (should be 'True' if JetAnalyzer is run more than once)
480  genJetCol = 'slimmedGenJets',
481  rho = ('fixedGridRhoFastjetAll','',''),
482  jetPt = 25.,
483  jetEta = 4.7,
484  jetEtaCentral = 2.4,
485  jetLepDR = 0.4,
486  jetLepArbitration = (lambda jet,lepton : lepton), # you can decide which to keep in case of overlaps; e.g. if the jet is b-tagged you might want to keep the jet
487  cleanSelectedLeptons = True, #Whether to clean 'selectedLeptons' after disambiguation. Treat with care (= 'False') if running Jetanalyzer more than once
488  minLepPt = 10,
489  lepSelCut = lambda lep : True,
490  relaxJetId = False,
491  doPuId = False, # Not commissioned in 7.0.X
492  doQG = False,
493  checkLeptonPFOverlap = True,
494  recalibrateJets = False,
495  applyL2L3Residual = 'Data', # if recalibrateJets, apply L2L3Residual to Data only
496  recalibrationType = "AK4PFchs",
497  shiftJEC = 0, # set to +1 or -1 to apply +/-1 sigma shift to the nominal jet energies
498  addJECShifts = False, # if true, add "corr", "corrJECUp", and "corrJECDown" for each jet (requires uncertainties to be available!)
499  smearJets = True,
500  shiftJER = 0, # set to +1 or -1 to get +/-1 sigma shifts
501  jecPath = "",
502  calculateSeparateCorrections = False,
503  calculateType1METCorrection = False,
504  type1METParams = { 'jetPtThreshold':15., 'skipEMfractionThreshold':0.9, 'skipMuons':True },
505  addJERShifts = 0, # add +/-1 sigma shifts to jets, intended to be used with shiftJER=0
506  cleanJetsFromFirstPhoton = False,
507  cleanJetsFromTaus = False,
508  cleanJetsFromIsoTracks = False,
509  alwaysCleanPhotons = False,
510  do_mc_match=True,
511  cleanGenJetsFromPhoton = False,
512  jetGammaDR=0.4,
513  cleanFromLepAndGammaSimultaneously = False,
514  jetGammaLepDR=0.4,
515  attachNeutrinos = True,
516  genNuSelection = lambda nu : True, #FIXME: add here check for ispromptfinalstate
517  collectionPostFix = ""
518  )
519 )
deltar.bestMatch
def bestMatch
Definition: deltar.py:139
objects.JetAnalyzer.JetAnalyzer.deltaMetFromJEC
deltaMetFromJEC
Read jets, if necessary recalibrate and shift MET.
Definition: JetAnalyzer.py:141
SingleMuPt40Fwdv2_cfi_GEN_SIM.pdgId
tuple pdgId
Definition: SingleMuPt40Fwdv2_cfi_GEN_SIM.py:85
deltar.matchObjectCollection
def matchObjectCollection
Definition: deltar.py:152
Jet
Definition: Jet.py:1
deltar.matchObjectCollection2
def matchObjectCollection2
Definition: deltar.py:167
ComparisonHelper::zip
OutputIterator zip(InputIterator1 first1, InputIterator1 last1, InputIterator2 first2, InputIterator2 last2, OutputIterator result, Compare comp)
Definition: L1TStage2CaloLayer1.h:38
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
reco::deltaR2
T1 deltaR2(T1 eta1, T2 phi1, T3 eta2, T4 phi2)
Definition: deltaR.h:36
reco::if
if(dp >Float(M_PI)) dp-
Generated for CMSSW Reference Manual by doxygen 1.8.5