CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
LeptonAnalyzer.py
Go to the documentation of this file.
1 from PhysicsTools.Heppy.analyzers.core.Analyzer import Analyzer
2 from PhysicsTools.Heppy.analyzers.core.AutoHandle import AutoHandle
3 from PhysicsTools.Heppy.physicsobjects.Electron import Electron
4 from PhysicsTools.Heppy.physicsobjects.Muon import Muon
5 #from CMGTools.TTHAnalysis.tools.EfficiencyCorrector import EfficiencyCorrector
6 
7 from PhysicsTools.HeppyCore.utils.deltar import bestMatch
8 from PhysicsTools.Heppy.physicsutils.RochesterCorrections import rochcor
9 from PhysicsTools.Heppy.physicsutils.MuScleFitCorrector import MuScleFitCorr
10 from PhysicsTools.Heppy.physicsutils.KalmanMuonCorrector import KalmanMuonCorrector
11 from PhysicsTools.Heppy.physicsutils.ElectronCalibrator import Run2ElectronCalibrator
12 #from CMGTools.TTHAnalysis.electronCalibrator import ElectronCalibrator
13 import PhysicsTools.HeppyCore.framework.config as cfg
14 from PhysicsTools.HeppyCore.utils.deltar import *
15 from PhysicsTools.Heppy.physicsutils.genutils import *
16 
17 
18 from ROOT import heppy
19 cmgMuonCleanerBySegments = heppy.CMGMuonCleanerBySegmentsAlgo()
20 
21 class LeptonAnalyzer( Analyzer ):
22 
23 
24  def __init__(self, cfg_ana, cfg_comp, looperName ):
25  super(LeptonAnalyzer,self).__init__(cfg_ana,cfg_comp,looperName)
26  if hasattr(self.cfg_ana, 'doMuScleFitCorrections'):
27  raise RuntimeError("doMuScleFitCorrections is not supported. Please set instead doMuonScaleCorrections = ( 'MuScleFit', <name> )")
28  if hasattr(self.cfg_ana, 'doRochesterCorrections'):
29  raise RuntimeError("doRochesterCorrections is not supported. Please set instead doMuonScaleCorrections = ( 'Rochester', <name> )")
30  if self.cfg_ana.doMuonScaleCorrections:
31  algo, options = self.cfg_ana.doMuonScaleCorrections
32  if algo == "Kalman":
33  corr = options['MC' if self.cfg_comp.isMC else 'Data']
34  self.muonScaleCorrector = KalmanMuonCorrector(corr,
35  self.cfg_comp.isMC,
36  options['isSync'] if 'isSync' in options else False,
37  options['smearMode'] if 'smearMode' in options else "ebe")
38  elif algo == "Rochester":
39  print "WARNING: the Rochester correction in heppy is still from Run 1"
40  self.muonScaleCorrector = RochesterCorrections()
41  elif algo == "MuScleFit":
42  print "WARNING: the MuScleFit correction in heppy is still from Run 1 (and probably no longer functional)"
43  if options not in [ "prompt", "prompt-sync", "rereco", "rereco-sync" ]:
44  raise RuntimeError('MuScleFit correction name must be one of [ "prompt", "prompt-sync", "rereco", "rereco-sync" ] ')
45  rereco = ("prompt" not in self.cfg_ana.doMuScleFitCorrections)
46  sync = ("sync" in self.cfg_ana.doMuScleFitCorrections)
47  self.muonScaleCorrector = MuScleFitCorr(cfg_comp.isMC, rereco, sync)
48  else: raise RuntimeError("Unknown muon scale correction algorithm")
49  else:
50  self.muonScaleCorrector = None
51  #FIXME: only Embedded works
52  if self.cfg_ana.doElectronScaleCorrections:
53  conf = cfg_ana.doElectronScaleCorrections
54  self.electronEnergyCalibrator = Run2ElectronCalibrator(
55  conf['data'],
56  conf['GBRForest'],
57  cfg_comp.isMC,
58  conf['isSync'] if 'isSync' in conf else False,
59  )
60 # if hasattr(cfg_comp,'efficiency'):
61 # self.efficiency= EfficiencyCorrector(cfg_comp.efficiency)
62  # Isolation cut
63  if hasattr(cfg_ana, 'loose_electron_isoCut'):
64  self.eleIsoCut = cfg_ana.loose_electron_isoCut
65  else:
66  self.eleIsoCut = lambda ele : (
67  ele.relIso03 <= self.cfg_ana.loose_electron_relIso and
68  ele.absIso03 < getattr(self.cfg_ana,'loose_electron_absIso',9e99))
69  if hasattr(cfg_ana, 'loose_muon_isoCut'):
70  self.muIsoCut = cfg_ana.loose_muon_isoCut
71  else:
72  self.muIsoCut = lambda mu : (
73  mu.relIso03 <= self.cfg_ana.loose_muon_relIso and
74  mu.absIso03 < getattr(self.cfg_ana,'loose_muon_absIso',9e99))
75 
76 
77 
78  self.eleEffectiveArea = getattr(cfg_ana, 'ele_effectiveAreas', "Spring15_25ns_v1")
79  self.muEffectiveArea = getattr(cfg_ana, 'mu_effectiveAreas', "Spring15_25ns_v1")
80  # MiniIsolation
81  self.doMiniIsolation = getattr(cfg_ana, 'doMiniIsolation', False)
82  if self.doMiniIsolation:
83  self.miniIsolationPUCorr = self.cfg_ana.miniIsolationPUCorr
84  self.miniIsolationVetoLeptons = self.cfg_ana.miniIsolationVetoLeptons
85  if self.miniIsolationVetoLeptons not in [ None, 'any', 'inclusive' ]:
86  raise RuntimeError("miniIsolationVetoLeptons should be None, or 'any' (all reco leptons), or 'inclusive' (all inclusive leptons)")
87  if self.miniIsolationPUCorr == "weights":
88  self.IsolationComputer = heppy.IsolationComputer(0.4)
89  else:
90  self.IsolationComputer = heppy.IsolationComputer()
91 
92  self.doIsoAnnulus = getattr(cfg_ana, 'doIsoAnnulus', False)
93  if self.doIsoAnnulus:
94  if not self.doMiniIsolation:
95  self.IsolationComputer = heppy.IsolationComputer()
96 
97  self.doIsolationScan = getattr(cfg_ana, 'doIsolationScan', False)
98  if self.doIsolationScan:
99  if self.doMiniIsolation:
100  assert (self.miniIsolationPUCorr!="weights")
101  assert (self.miniIsolationVetoLeptons==None)
102  else:
103  self.IsolationComputer = heppy.IsolationComputer()
104 
105 
106  self.doMatchToPhotons = getattr(cfg_ana, 'do_mc_match_photons', False)
107 
108  #----------------------------------------
109  # DECLARATION OF HANDLES OF LEPTONS STUFF
110  #----------------------------------------
111 
112 
113  def declareHandles(self):
114  super(LeptonAnalyzer, self).declareHandles()
115 
116  #leptons
117  self.handles['muons'] = AutoHandle(self.cfg_ana.muons,"std::vector<pat::Muon>")
118  self.handles['electrons'] = AutoHandle(self.cfg_ana.electrons,"std::vector<pat::Electron>")
119 
120  #rho for muons
121  self.handles['rhoMu'] = AutoHandle( self.cfg_ana.rhoMuon, 'double')
122  #rho for electrons
123  self.handles['rhoEle'] = AutoHandle( self.cfg_ana.rhoElectron, 'double')
124 
125  if self.doMiniIsolation or self.doIsolationScan:
126  self.handles['packedCandidates'] = AutoHandle( self.cfg_ana.packedCandidates, 'std::vector<pat::PackedCandidate>')
127 
128  if self.doMatchToPhotons:
129  if self.doMatchToPhotons == "any":
130  self.mchandles['genPhotons'] = AutoHandle( 'packedGenParticles', 'std::vector<pat::PackedGenParticle>' )
131  else:
132  self.mchandles['genPhotons'] = AutoHandle( 'prunedGenParticles', 'std::vector<reco::GenParticle>' )
133 
134  def beginLoop(self, setup):
135  super(LeptonAnalyzer,self).beginLoop(setup)
136  self.counters.addCounter('events')
137  count = self.counters.counter('events')
138  count.register('all events')
139 
140  #------------------
141  # MAKE LEPTON LISTS
142  #------------------
143 
144 
145  def makeLeptons(self, event):
146  ### inclusive leptons = all leptons that could be considered somewhere in the analysis, with minimal requirements (used e.g. to match to MC)
147  event.inclusiveLeptons = []
148  ### selected leptons = subset of inclusive leptons passing some basic id definition and pt requirement
149  ### other leptons = subset of inclusive leptons failing some basic id definition and pt requirement
150  event.selectedLeptons = []
151  event.selectedMuons = []
152  event.selectedElectrons = []
153  event.otherLeptons = []
154 
155  if self.doMiniIsolation or self.doIsolationScan:
156  self.IsolationComputer.setPackedCandidates(self.handles['packedCandidates'].product())
157  if self.doMiniIsolation:
158  if self.miniIsolationVetoLeptons == "any":
159  for lep in self.handles['muons'].product():
160  self.IsolationComputer.addVetos(lep)
161  for lep in self.handles['electrons'].product():
162  self.IsolationComputer.addVetos(lep)
163 
164  #muons
165  allmuons = self.makeAllMuons(event)
166 
167  #electrons
168  allelectrons = self.makeAllElectrons(event)
169 
170  #make inclusive leptons
171  inclusiveMuons = []
172  inclusiveElectrons = []
173  for mu in allmuons:
174  if (mu.track().isNonnull() and mu.muonID(self.cfg_ana.inclusive_muon_id) and
175  mu.pt()>self.cfg_ana.inclusive_muon_pt and abs(mu.eta())<self.cfg_ana.inclusive_muon_eta and
176  abs(mu.dxy())<self.cfg_ana.inclusive_muon_dxy and abs(mu.dz())<self.cfg_ana.inclusive_muon_dz):
177  inclusiveMuons.append(mu)
178  for ele in allelectrons:
179  if ( ele.electronID(self.cfg_ana.inclusive_electron_id) and
180  ele.pt()>self.cfg_ana.inclusive_electron_pt and abs(ele.eta())<self.cfg_ana.inclusive_electron_eta and
181  abs(ele.dxy())<self.cfg_ana.inclusive_electron_dxy and abs(ele.dz())<self.cfg_ana.inclusive_electron_dz and
182  ele.lostInner()<=self.cfg_ana.inclusive_electron_lostHits ):
183  inclusiveElectrons.append(ele)
184  event.inclusiveLeptons = inclusiveMuons + inclusiveElectrons
185 
186  if self.doMiniIsolation:
187  if self.miniIsolationVetoLeptons == "inclusive":
188  for lep in event.inclusiveLeptons:
189  self.IsolationComputer.addVetos(lep.physObj)
190  for lep in event.inclusiveLeptons:
191  self.attachMiniIsolation(lep)
192 
193  if self.doIsoAnnulus:
194  for lep in event.inclusiveLeptons:
195  self.attachIsoAnnulus04(lep)
196 
197  if self.doIsolationScan:
198  for lep in event.inclusiveLeptons:
199  self.attachIsolationScan(lep)
200 
201  # make loose leptons (basic selection)
202  for mu in inclusiveMuons:
203  if (mu.muonID(self.cfg_ana.loose_muon_id) and
204  mu.pt() > self.cfg_ana.loose_muon_pt and abs(mu.eta()) < self.cfg_ana.loose_muon_eta and
205  abs(mu.dxy()) < self.cfg_ana.loose_muon_dxy and abs(mu.dz()) < self.cfg_ana.loose_muon_dz and
206  self.muIsoCut(mu)):
207  mu.looseIdSusy = True
208  event.selectedLeptons.append(mu)
209  event.selectedMuons.append(mu)
210  else:
211  mu.looseIdSusy = False
212  event.otherLeptons.append(mu)
213  looseMuons = event.selectedLeptons[:]
214  for ele in inclusiveElectrons:
215  if (ele.electronID(self.cfg_ana.loose_electron_id) and
216  ele.pt()>self.cfg_ana.loose_electron_pt and abs(ele.eta())<self.cfg_ana.loose_electron_eta and
217  abs(ele.dxy()) < self.cfg_ana.loose_electron_dxy and abs(ele.dz())<self.cfg_ana.loose_electron_dz and
218  self.eleIsoCut(ele) and
219  ele.lostInner() <= self.cfg_ana.loose_electron_lostHits and
220  ( True if getattr(self.cfg_ana,'notCleaningElectrons',False) else (bestMatch(ele, looseMuons)[1] > (self.cfg_ana.min_dr_electron_muon**2)) )):
221  event.selectedLeptons.append(ele)
222  event.selectedElectrons.append(ele)
223  ele.looseIdSusy = True
224  else:
225  event.otherLeptons.append(ele)
226  ele.looseIdSusy = False
227 
228  event.otherLeptons.sort(key = lambda l : l.pt(), reverse = True)
229  event.selectedLeptons.sort(key = lambda l : l.pt(), reverse = True)
230  event.selectedMuons.sort(key = lambda l : l.pt(), reverse = True)
231  event.selectedElectrons.sort(key = lambda l : l.pt(), reverse = True)
232  event.inclusiveLeptons.sort(key = lambda l : l.pt(), reverse = True)
233 
234  for lepton in event.selectedLeptons:
235  if hasattr(self,'efficiency'):
236  self.efficiency.attachToObject(lepton)
237 
238  def makeAllMuons(self, event):
239  """
240  make a list of all muons, and apply basic corrections to them
241  """
242  # Start from all muons
243  allmuons = map( Muon, self.handles['muons'].product() )
244 
245  # Muon scale and resolution corrections (if enabled)
246  if self.muonScaleCorrector:
247  self.muonScaleCorrector.correct_all(allmuons, event.run)
248 
249  # Clean up dulicate muons (note: has no effect unless the muon id is removed)
250  if self.cfg_ana.doSegmentBasedMuonCleaning:
251  isgood = cmgMuonCleanerBySegments.clean( self.handles['muons'].product() )
252  newmu = []
253  for i,mu in enumerate(allmuons):
254  if isgood[i]: newmu.append(mu)
255  allmuons = newmu
256 
257  # Attach EAs for isolation:
258  for mu in allmuons:
259  mu.rho = float(self.handles['rhoMu'].product()[0])
260  if self.muEffectiveArea == "Data2012":
261  if aeta < 1.0 : mu.EffectiveArea03 = 0.382;
262  elif aeta < 1.47 : mu.EffectiveArea03 = 0.317;
263  elif aeta < 2.0 : mu.EffectiveArea03 = 0.242;
264  elif aeta < 2.2 : mu.EffectiveArea03 = 0.326;
265  elif aeta < 2.3 : mu.EffectiveArea03 = 0.462;
266  else : mu.EffectiveArea03 = 0.372;
267  if aeta < 1.0 : mu.EffectiveArea04 = 0.674;
268  elif aeta < 1.47 : mu.EffectiveArea04 = 0.565;
269  elif aeta < 2.0 : mu.EffectiveArea04 = 0.442;
270  elif aeta < 2.2 : mu.EffectiveArea04 = 0.515;
271  elif aeta < 2.3 : mu.EffectiveArea04 = 0.821;
272  else : mu.EffectiveArea04 = 0.660;
273  elif self.muEffectiveArea == "Phys14_25ns_v1":
274  aeta = abs(mu.eta())
275  if aeta < 0.800: mu.EffectiveArea03 = 0.0913
276  elif aeta < 1.300: mu.EffectiveArea03 = 0.0765
277  elif aeta < 2.000: mu.EffectiveArea03 = 0.0546
278  elif aeta < 2.200: mu.EffectiveArea03 = 0.0728
279  else: mu.EffectiveArea03 = 0.1177
280  if aeta < 0.800: mu.EffectiveArea04 = 0.1564
281  elif aeta < 1.300: mu.EffectiveArea04 = 0.1325
282  elif aeta < 2.000: mu.EffectiveArea04 = 0.0913
283  elif aeta < 2.200: mu.EffectiveArea04 = 0.1212
284  else: mu.EffectiveArea04 = 0.2085
285  elif self.muEffectiveArea == "Spring15_25ns_v1":
286  aeta = abs(mu.eta())
287  if aeta < 0.800: mu.EffectiveArea03 = 0.0735
288  elif aeta < 1.300: mu.EffectiveArea03 = 0.0619
289  elif aeta < 2.000: mu.EffectiveArea03 = 0.0465
290  elif aeta < 2.200: mu.EffectiveArea03 = 0.0433
291  else: mu.EffectiveArea03 = 0.0577
292  mu.EffectiveArea04 = 0 # not computed
293  else: raise RuntimeError("Unsupported value for mu_effectiveAreas: can only use Data2012 (rho: ?) and Phys14_25ns_v1 or Spring15_25ns_v1 (rho: fixedGridRhoFastjetAll)")
294  # Attach the vertex to them, for dxy/dz calculation
295  for mu in allmuons:
296  mu.associatedVertex = event.goodVertices[0] if len(event.goodVertices)>0 else event.vertices[0]
297  mu.setTrackForDxyDz(self.cfg_ana.muon_dxydz_track)
298 
299  # Set tight id if specified
300  if hasattr(self.cfg_ana, "mu_tightId"):
301  for mu in allmuons:
302  mu.tightIdResult = mu.muonID(self.cfg_ana.mu_tightId)
303 
304  # Compute relIso in 0.3 and 0.4 cones
305  for mu in allmuons:
306  if self.cfg_ana.mu_isoCorr=="rhoArea" :
307  mu.absIso03 = (mu.pfIsolationR03().sumChargedHadronPt + max( mu.pfIsolationR03().sumNeutralHadronEt + mu.pfIsolationR03().sumPhotonEt - mu.rho * mu.EffectiveArea03,0.0))
308  mu.absIso04 = (mu.pfIsolationR04().sumChargedHadronPt + max( mu.pfIsolationR04().sumNeutralHadronEt + mu.pfIsolationR04().sumPhotonEt - mu.rho * mu.EffectiveArea04,0.0))
309  elif self.cfg_ana.mu_isoCorr=="deltaBeta" :
310  mu.absIso03 = (mu.pfIsolationR03().sumChargedHadronPt + max( mu.pfIsolationR03().sumNeutralHadronEt + mu.pfIsolationR03().sumPhotonEt - mu.pfIsolationR03().sumPUPt/2,0.0))
311  mu.absIso04 = (mu.pfIsolationR04().sumChargedHadronPt + max( mu.pfIsolationR04().sumNeutralHadronEt + mu.pfIsolationR04().sumPhotonEt - mu.pfIsolationR04().sumPUPt/2,0.0))
312  else :
313  raise RuntimeError("Unsupported mu_isoCorr name '" + str(self.cfg_ana.mu_isoCorr) + "'! For now only 'rhoArea' and 'deltaBeta' are supported.")
314  mu.relIso03 = mu.absIso03/mu.pt()
315  mu.relIso04 = mu.absIso04/mu.pt()
316  return allmuons
317 
318  def makeAllElectrons(self, event):
319  """
320  make a list of all electrons, and apply basic corrections to them
321  """
322  allelectrons = map( Electron, self.handles['electrons'].product() )
323 
324  ## Duplicate removal for fast sim (to be checked if still necessary in latest greatest 5.3.X releases)
325  allelenodup = []
326  for e in allelectrons:
327  dup = False
328  for e2 in allelenodup:
329  if abs(e.pt()-e2.pt()) < 1e-6 and abs(e.eta()-e2.eta()) < 1e-6 and abs(e.phi()-e2.phi()) < 1e-6 and e.charge() == e2.charge():
330  dup = True
331  break
332  if not dup: allelenodup.append(e)
333  allelectrons = allelenodup
334 
335  # fill EA for rho-corrected isolation
336  for ele in allelectrons:
337  ele.rho = float(self.handles['rhoEle'].product()[0])
338  if self.eleEffectiveArea == "Data2012":
339  # https://twiki.cern.ch/twiki/bin/viewauth/CMS/EgammaEARhoCorrection?rev=14
340  SCEta = abs(ele.superCluster().eta())
341  if SCEta < 1.0 : ele.EffectiveArea03 = 0.13 # 0.130;
342  elif SCEta < 1.479: ele.EffectiveArea03 = 0.14 # 0.137;
343  elif SCEta < 2.0 : ele.EffectiveArea03 = 0.07 # 0.067;
344  elif SCEta < 2.2 : ele.EffectiveArea03 = 0.09 # 0.089;
345  elif SCEta < 2.3 : ele.EffectiveArea03 = 0.11 # 0.107;
346  elif SCEta < 2.4 : ele.EffectiveArea03 = 0.11 # 0.110;
347  else : ele.EffectiveArea03 = 0.14 # 0.138;
348  if SCEta < 1.0 : ele.EffectiveArea04 = 0.208;
349  elif SCEta < 1.479: ele.EffectiveArea04 = 0.209;
350  elif SCEta < 2.0 : ele.EffectiveArea04 = 0.115;
351  elif SCEta < 2.2 : ele.EffectiveArea04 = 0.143;
352  elif SCEta < 2.3 : ele.EffectiveArea04 = 0.183;
353  elif SCEta < 2.4 : ele.EffectiveArea04 = 0.194;
354  else : ele.EffectiveArea04 = 0.261;
355  elif self.eleEffectiveArea == "Phys14_25ns_v1":
356  aeta = abs(ele.eta())
357  if aeta < 0.800: ele.EffectiveArea03 = 0.1013
358  elif aeta < 1.300: ele.EffectiveArea03 = 0.0988
359  elif aeta < 2.000: ele.EffectiveArea03 = 0.0572
360  elif aeta < 2.200: ele.EffectiveArea03 = 0.0842
361  else: ele.EffectiveArea03 = 0.1530
362  if aeta < 0.800: ele.EffectiveArea04 = 0.1830
363  elif aeta < 1.300: ele.EffectiveArea04 = 0.1734
364  elif aeta < 2.000: ele.EffectiveArea04 = 0.1077
365  elif aeta < 2.200: ele.EffectiveArea04 = 0.1565
366  else: ele.EffectiveArea04 = 0.2680
367  elif self.eleEffectiveArea == "Spring15_50ns_v1":
368  SCEta = abs(ele.superCluster().eta())
369  ## ----- https://github.com/ikrav/cmssw/blob/egm_id_747_v2/RecoEgamma/ElectronIdentification/data/Spring15/effAreaElectrons_cone03_pfNeuHadronsAndPhotons_50ns.txt
370  if SCEta < 0.800: ele.EffectiveArea03 = 0.0973
371  elif SCEta < 1.300: ele.EffectiveArea03 = 0.0954
372  elif SCEta < 2.000: ele.EffectiveArea03 = 0.0632
373  elif SCEta < 2.200: ele.EffectiveArea03 = 0.0727
374  else: ele.EffectiveArea03 = 0.1337
375  # warning: EAs not computed for cone DR=0.4 yet. Do not correct
376  ele.EffectiveArea04 = 0.0
377  elif self.eleEffectiveArea == "Spring15_25ns_v1":
378  SCEta = abs(ele.superCluster().eta())
379  ## ----- https://github.com/ikrav/cmssw/blob/egm_id_747_v2/RecoEgamma/ElectronIdentification/data/Spring15/effAreaElectrons_cone03_pfNeuHadronsAndPhotons_25ns.txt
380  if SCEta < 1.000: ele.EffectiveArea03 = 0.1752
381  elif SCEta < 1.479: ele.EffectiveArea03 = 0.1862
382  elif SCEta < 2.000: ele.EffectiveArea03 = 0.1411
383  elif SCEta < 2.200: ele.EffectiveArea03 = 0.1534
384  elif SCEta < 2.300: ele.EffectiveArea03 = 0.1903
385  elif SCEta < 2.400: ele.EffectiveArea03 = 0.2243
386  else: ele.EffectiveArea03 = 0.2687
387  # warning: EAs not computed for cone DR=0.4 yet. Do not correct
388  ele.EffectiveArea04 = 0.0
389  else: raise RuntimeError("Unsupported value for ele_effectiveAreas: can only use Data2012 (rho: ?), Phys14_v1 and Spring15_v1 (rho: fixedGridRhoFastjetAll)")
390 
391  # Electron scale calibrations
392  if self.cfg_ana.doElectronScaleCorrections:
393  for ele in allelectrons:
394  self.electronEnergyCalibrator.correct(ele, event.run)
395 
396  # Attach the vertex
397  for ele in allelectrons:
398  ele.associatedVertex = event.goodVertices[0] if len(event.goodVertices)>0 else event.vertices[0]
399 
400  # Compute relIso with R=0.3 and R=0.4 cones
401  for ele in allelectrons:
402  if self.cfg_ana.ele_isoCorr=="rhoArea" :
403  ele.absIso03 = (ele.chargedHadronIsoR(0.3) + max(ele.neutralHadronIsoR(0.3)+ele.photonIsoR(0.3)-ele.rho*ele.EffectiveArea03,0))
404  ele.absIso04 = (ele.chargedHadronIsoR(0.4) + max(ele.neutralHadronIsoR(0.4)+ele.photonIsoR(0.4)-ele.rho*ele.EffectiveArea04,0))
405  elif self.cfg_ana.ele_isoCorr=="deltaBeta" :
406  ele.absIso03 = (ele.chargedHadronIsoR(0.3) + max( ele.neutralHadronIsoR(0.3)+ele.photonIsoR(0.3) - ele.puChargedHadronIsoR(0.3)/2, 0.0))
407  ele.absIso04 = (ele.chargedHadronIsoR(0.4) + max( ele.neutralHadronIsoR(0.4)+ele.photonIsoR(0.4) - ele.puChargedHadronIsoR(0.4)/2, 0.0))
408  else :
409  raise RuntimeError("Unsupported ele_isoCorr name '" + str(self.cfg_ana.ele_isoCorr) + "'! For now only 'rhoArea' and 'deltaBeta' are supported.")
410  ele.relIso03 = ele.absIso03/ele.pt()
411  ele.relIso04 = ele.absIso04/ele.pt()
412 
413  # Set tight MVA id
414  for ele in allelectrons:
415  if self.cfg_ana.ele_tightId=="MVA" :
416  ele.tightIdResult = ele.electronID("POG_MVA_ID_Trig_full5x5")
417  elif self.cfg_ana.ele_tightId=="Cuts_2012" :
418  ele.tightIdResult = -1 + 1*ele.electronID("POG_Cuts_ID_2012_Veto_full5x5") + 1*ele.electronID("POG_Cuts_ID_2012_Loose_full5x5") + 1*ele.electronID("POG_Cuts_ID_2012_Medium_full5x5") + 1*ele.electronID("POG_Cuts_ID_2012_Tight_full5x5")
419  elif self.cfg_ana.ele_tightId=="Cuts_PHYS14_25ns_v1_ConvVetoDxyDz" :
420  ele.tightIdResult = -1 + 1*ele.electronID("POG_Cuts_ID_PHYS14_25ns_v1_ConvVetoDxyDz_Veto_full5x5") + 1*ele.electronID("POG_Cuts_ID_PHYS14_25ns_v1_ConvVetoDxyDz_Loose_full5x5") + 1*ele.electronID("POG_Cuts_ID_PHYS14_25ns_v1_ConvVetoDxyDz_Medium_full5x5") + 1*ele.electronID("POG_Cuts_ID_PHYS14_25ns_v1_ConvVetoDxyDz_Tight_full5x5")
421  elif self.cfg_ana.ele_tightId=="Cuts_SPRING15_25ns_v1_ConvVetoDxyDz" :
422  ele.tightIdResult = -1 + 1*ele.electronID("POG_Cuts_ID_SPRING15_25ns_v1_ConvVetoDxyDz_Veto_full5x5") + 1*ele.electronID("POG_Cuts_ID_SPRING15_25ns_v1_ConvVetoDxyDz_Loose_full5x5") + 1*ele.electronID("POG_Cuts_ID_SPRING15_25ns_v1_ConvVetoDxyDz_Medium_full5x5") + 1*ele.electronID("POG_Cuts_ID_SPRING15_25ns_v1_ConvVetoDxyDz_Tight_full5x5")
423 
424  else :
425  try:
426  ele.tightIdResult = ele.electronID(self.cfg_ana.ele_tightId)
427  except RuntimeError:
428  raise RuntimeError("Unsupported ele_tightId name '" + str(self.cfg_ana.ele_tightId) + "'! For now only 'MVA' and 'Cuts_2012' are supported, in addition to what provided in Electron.py.")
429 
430 
431  return allelectrons
432 
433  def attachMiniIsolation(self, mu):
434  mu.miniIsoR = 10.0/min(max(mu.pt(), 50),200)
435  # -- version with increasing cone at low pT, gives slightly better performance for tight cuts and low pt leptons
436  # mu.miniIsoR = 10.0/min(max(mu.pt(), 50),200) if mu.pt() > 20 else 4.0/min(max(mu.pt(),10),20)
437  what = "mu" if (abs(mu.pdgId()) == 13) else ("eleB" if mu.isEB() else "eleE")
438  if what == "mu":
439  mu.miniAbsIsoCharged = self.IsolationComputer.chargedAbsIso(mu.physObj, mu.miniIsoR, {"mu":0.0001,"eleB":0,"eleE":0.015}[what], 0.0);
440  else:
441  mu.miniAbsIsoCharged = self.IsolationComputer.chargedAbsIso(mu.physObj, mu.miniIsoR, {"mu":0.0001,"eleB":0,"eleE":0.015}[what], 0.0,self.IsolationComputer.selfVetoNone);
442 
443  if self.miniIsolationPUCorr == None: puCorr = self.cfg_ana.mu_isoCorr if what=="mu" else self.cfg_ana.ele_isoCorr
444  else: puCorr = self.miniIsolationPUCorr
445 
446  if puCorr == "weights":
447  if what == "mu":
448  mu.miniAbsIsoNeutral = self.IsolationComputer.neutralAbsIsoWeighted(mu.physObj, mu.miniIsoR, 0.01, 0.5);
449  else:
450  mu.miniAbsIsoNeutral = ( self.IsolationComputer.photonAbsIsoWeighted( mu.physObj, mu.miniIsoR, 0.08 if what == "eleE" else 0.0, 0.0, self.IsolationComputer.selfVetoNone) +
451  self.IsolationComputer.neutralHadAbsIsoWeighted(mu.physObj, mu.miniIsoR, 0.0, 0.0, self.IsolationComputer.selfVetoNone) )
452  else:
453  if what == "mu":
454  mu.miniAbsIsoNeutral = self.IsolationComputer.neutralAbsIsoRaw(mu.physObj, mu.miniIsoR, 0.01, 0.5);
455  else:
456  mu.miniAbsIsoPho = self.IsolationComputer.photonAbsIsoRaw( mu.physObj, mu.miniIsoR, 0.08 if what == "eleE" else 0.0, 0.0, self.IsolationComputer.selfVetoNone)
457  mu.miniAbsIsoNHad = self.IsolationComputer.neutralHadAbsIsoRaw(mu.physObj, mu.miniIsoR, 0.0, 0.0, self.IsolationComputer.selfVetoNone)
458  mu.miniAbsIsoNeutral = mu.miniAbsIsoPho + mu.miniAbsIsoNHad
459  # -- version relying on PF candidate vetos; apparently less performant, and the isolation computed at RECO level doesn't have them
460  #mu.miniAbsIsoPhoSV = self.IsolationComputer.photonAbsIsoRaw( mu.physObj, mu.miniIsoR, 0.0, 0.0)
461  #mu.miniAbsIsoNHadSV = self.IsolationComputer.neutralHadAbsIsoRaw(mu.physObj, mu.miniIsoR, 0.0, 0.0)
462  #mu.miniAbsIsoNeutral = mu.miniAbsIsoPhoSV + mu.miniAbsIsoNHadSV
463  if puCorr == "rhoArea":
464  mu.miniAbsIsoNeutral = max(0.0, mu.miniAbsIsoNeutral - mu.rho * mu.EffectiveArea03 * (mu.miniIsoR/0.3)**2)
465  elif puCorr == "deltaBeta":
466  if what == "mu":
467  mu.miniAbsIsoPU = self.IsolationComputer.puAbsIso(mu.physObj, mu.miniIsoR, 0.01, 0.5);
468  else:
469  mu.miniAbsIsoPU = self.IsolationComputer.puAbsIso(mu.physObj, mu.miniIsoR, 0.015 if what == "eleE" else 0.0, 0.0,self.IsolationComputer.selfVetoNone);
470  mu.miniAbsIsoNeutral = max(0.0, mu.miniAbsIsoNeutral - 0.5*mu.miniAbsIsoPU)
471  elif puCorr != 'raw':
472  raise RuntimeError("Unsupported miniIsolationCorr name '" + puCorr + "'! For now only 'rhoArea', 'deltaBeta', 'raw', 'weights' are supported (and 'weights' is not tested).")
473 
474  mu.miniAbsIso = mu.miniAbsIsoCharged + mu.miniAbsIsoNeutral
475  mu.miniRelIso = mu.miniAbsIso/mu.pt()
476 
477 
478  def attachIsoAnnulus04(self, mu): # annulus isolation with outer cone of 0.4 and delta beta PU correction
479  mu.miniIsoR = 10.0/min(max(mu.pt(), 50),200)
480  mu.absIsoAnCharged = self.IsolationComputer.chargedAbsIso (mu.physObj, 0.4, mu.miniIsoR, 0.0, self.IsolationComputer.selfVetoNone)
481  mu.absIsoAnPho = self.IsolationComputer.photonAbsIsoRaw (mu.physObj, 0.4, mu.miniIsoR, 0.0, self.IsolationComputer.selfVetoNone)
482  mu.absIsoAnNHad = self.IsolationComputer.neutralHadAbsIsoRaw(mu.physObj, 0.4, mu.miniIsoR, 0.0, self.IsolationComputer.selfVetoNone)
483  mu.absIsoAnPU = self.IsolationComputer.puAbsIso (mu.physObj, 0.4, mu.miniIsoR, 0.0, self.IsolationComputer.selfVetoNone)
484  mu.absIsoAnNeutral = max(0.0, mu.absIsoAnPho + mu.absIsoAnNHad - 0.5*mu.absIsoAnPU)
485 
486  mu.absIsoAn04 = mu.absIsoAnCharged + mu.absIsoAnNeutral
487  mu.relIsoAn04 = mu.absIsoAn04/mu.pt()
488 
489 
490  def attachIsolationScan(self, mu):
491 
492  what = "mu" if (abs(mu.pdgId()) == 13) else ("eleB" if mu.isEB() else "eleE")
493  vetoreg = {"mu":0.0001,"eleB":0,"eleE":0.015}[what]
494 
495  if what=="mu":
496  mu.ScanAbsIsoCharged005 = self.IsolationComputer.chargedAbsIso(mu.physObj, 0.05, vetoreg, 0.0)
497  mu.ScanAbsIsoCharged01 = self.IsolationComputer.chargedAbsIso(mu.physObj, 0.1, vetoreg, 0.0)
498  mu.ScanAbsIsoCharged02 = self.IsolationComputer.chargedAbsIso(mu.physObj, 0.2, vetoreg, 0.0)
499  mu.ScanAbsIsoCharged03 = self.IsolationComputer.chargedAbsIso(mu.physObj, 0.3, vetoreg, 0.0)
500  mu.ScanAbsIsoCharged04 = self.IsolationComputer.chargedAbsIso(mu.physObj, 0.4, vetoreg, 0.0)
501  else:
502  mu.ScanAbsIsoCharged005 = self.IsolationComputer.chargedAbsIso(mu.physObj, 0.05, vetoreg, 0.0, self.IsolationComputer.selfVetoNone)
503  mu.ScanAbsIsoCharged01 = self.IsolationComputer.chargedAbsIso(mu.physObj, 0.1, vetoreg, 0.0, self.IsolationComputer.selfVetoNone)
504  mu.ScanAbsIsoCharged02 = self.IsolationComputer.chargedAbsIso(mu.physObj, 0.2, vetoreg, 0.0, self.IsolationComputer.selfVetoNone)
505  mu.ScanAbsIsoCharged03 = self.IsolationComputer.chargedAbsIso(mu.physObj, 0.3, vetoreg, 0.0, self.IsolationComputer.selfVetoNone)
506  mu.ScanAbsIsoCharged04 = self.IsolationComputer.chargedAbsIso(mu.physObj, 0.4, vetoreg, 0.0, self.IsolationComputer.selfVetoNone)
507 
508  if what=="mu":
509  mu.ScanAbsIsoNeutral005 = self.IsolationComputer.neutralAbsIsoRaw(mu.physObj, 0.05, 0.01, 0.5)
510  mu.ScanAbsIsoNeutral01 = self.IsolationComputer.neutralAbsIsoRaw(mu.physObj, 0.1, 0.01, 0.5)
511  mu.ScanAbsIsoNeutral02 = self.IsolationComputer.neutralAbsIsoRaw(mu.physObj, 0.2, 0.01, 0.5)
512  mu.ScanAbsIsoNeutral03 = self.IsolationComputer.neutralAbsIsoRaw(mu.physObj, 0.3, 0.01, 0.5)
513  mu.ScanAbsIsoNeutral04 = self.IsolationComputer.neutralAbsIsoRaw(mu.physObj, 0.4, 0.01, 0.5)
514  else:
515  vetoreg = {"eleB":0.0,"eleE":0.08}[what]
516  mu.ScanAbsIsoNeutral005 = self.IsolationComputer.photonAbsIsoRaw(mu.physObj, 0.05, vetoreg, 0.0, self.IsolationComputer.selfVetoNone)+self.IsolationComputer.neutralHadAbsIsoRaw(mu.physObj, 0.05, 0.0, 0.0, self.IsolationComputer.selfVetoNone)
517  mu.ScanAbsIsoNeutral01 = self.IsolationComputer.photonAbsIsoRaw(mu.physObj, 0.1, vetoreg, 0.0, self.IsolationComputer.selfVetoNone)+self.IsolationComputer.neutralHadAbsIsoRaw(mu.physObj, 0.1, 0.0, 0.0, self.IsolationComputer.selfVetoNone)
518  mu.ScanAbsIsoNeutral02 = self.IsolationComputer.photonAbsIsoRaw(mu.physObj, 0.2, vetoreg, 0.0, self.IsolationComputer.selfVetoNone)+self.IsolationComputer.neutralHadAbsIsoRaw(mu.physObj, 0.2, 0.0, 0.0, self.IsolationComputer.selfVetoNone)
519  mu.ScanAbsIsoNeutral03 = self.IsolationComputer.photonAbsIsoRaw(mu.physObj, 0.3, vetoreg, 0.0, self.IsolationComputer.selfVetoNone)+self.IsolationComputer.neutralHadAbsIsoRaw(mu.physObj, 0.3, 0.0, 0.0, self.IsolationComputer.selfVetoNone)
520  mu.ScanAbsIsoNeutral04 = self.IsolationComputer.photonAbsIsoRaw(mu.physObj, 0.4, vetoreg, 0.0, self.IsolationComputer.selfVetoNone)+self.IsolationComputer.neutralHadAbsIsoRaw(mu.physObj, 0.4, 0.0, 0.0, self.IsolationComputer.selfVetoNone)
521 
522 
523  def matchLeptons(self, event):
524  def plausible(rec,gen):
525  if abs(rec.pdgId()) == 11 and abs(gen.pdgId()) != 11: return False
526  if abs(rec.pdgId()) == 13 and abs(gen.pdgId()) != 13: return False
527  dr = deltaR(rec.eta(),rec.phi(),gen.eta(),gen.phi())
528  if dr < 0.3: return True
529  if rec.pt() < 10 and abs(rec.pdgId()) == 13 and gen.pdgId() != rec.pdgId(): return False
530  if dr < 0.7: return True
531  if min(rec.pt(),gen.pt())/max(rec.pt(),gen.pt()) < 0.3: return False
532  return True
533 
534  leps = event.inclusiveLeptons if self.cfg_ana.match_inclusiveLeptons else event.selectedLeptons
535  match = matchObjectCollection3(leps,
536  event.genleps + event.gentauleps,
537  deltaRMax = 1.2, filter = plausible)
538  for lep in leps:
539  gen = match[lep]
540  lep.mcMatchId = (gen.sourceId if gen != None else 0)
541  lep.mcMatchTau = (gen in event.gentauleps if gen else -99)
542  lep.mcLep=gen
543 
544  def isFromB(self,particle,bid=5, done={}):
545  for i in xrange( particle.numberOfMothers() ):
546  mom = particle.mother(i)
547  momid = abs(mom.pdgId())
548  if momid / 1000 == bid or momid / 100 == bid or momid == bid:
549  return True
550  elif mom.status() == 2 and self.isFromB(mom, done=done, bid=bid):
551  return True
552  return False
553 
554  def matchAnyLeptons(self, event):
555  event.anyLeptons = [ x for x in event.genParticles if x.status() == 1 and abs(x.pdgId()) in [11,13] ]
556  leps = event.inclusiveLeptons if hasattr(event, 'inclusiveLeptons') else event.selectedLeptons
557  match = matchObjectCollection3(leps, event.anyLeptons, deltaRMax = 0.3, filter = lambda x,y : abs(x.pdgId()) == abs(y.pdgId()))
558  for lep in leps:
559  gen = match[lep]
560  lep.mcMatchAny_gp = gen
561  if gen:
562  if self.isFromB(gen): lep.mcMatchAny = 5 # B (inclusive of B->D)
563  elif self.isFromB(gen,bid=4): lep.mcMatchAny = 4 # Charm
564  else: lep.mcMatchAny = 1
565  if not getattr(lep, 'mcLep', None): lep.mcLep = gen
566  else:
567  if not getattr(lep, 'mcLep', None): lep.mcLep = None
568  lep.mcMatchAny = 0
569  # fix case where the matching with the only prompt leptons failed, but we still ended up with a prompt match
570  if gen != None and hasattr(lep,'mcMatchId') and lep.mcMatchId == 0:
571  if isPromptLepton(gen, False) or (gen.isPromptFinalState() or gen.isDirectPromptTauDecayProductFinalState()):
572  lep.mcMatchId = 100
573  lep.mcLep = gen
574  elif not hasattr(lep,'mcMatchId'):
575  lep.mcMatchId = 0
576  if not hasattr(lep,'mcMatchTau'): lep.mcMatchTau = 0
577 
578  def matchToPhotons(self, event):
579  event.anyPho = [ x for x in self.mchandles['genPhotons'].product() if x.status() == 1 and x.pdgId() == 22 and x.pt() > 1.0 ]
580  leps = event.inclusiveLeptons if hasattr(event, 'inclusiveLeptons') else event.selectedLeptons
581  leps = [ l for l in leps if abs(l.pdgId()) == 11 ]
582  plausible = lambda rec, gen : 0.3*gen.pt() < rec.pt() and rec.pt() < 1.5*gen.pt()
583  match = matchObjectCollection3(leps, event.anyPho, deltaRMax = 0.3, filter = plausible)
584  for lep in leps:
585  gen = match[lep]
586  lep.mcPho = gen
587  if lep.mcPho and lep.mcLep:
588  # I have both, I should keep the best one
589  def distance(rec,gen):
590  dr = deltaR(rec.eta(),rec.phi(),gen.eta(),gen.phi())
591  dptRel = abs(rec.pt()-gen.pt())/gen.pt()
592  return dr + 0.2*dptRel
593  dpho = distance(lep,lep.mcPho)
594  dlep = distance(lep,lep.mcLep)
595  if getattr(lep,'mcMatchAny_gp',None) and lep.mcMatchAny_gp != lep.mcLep:
596  dlep = min(dlep, distance(lep,lep.mcMatchAny_gp))
597  if dlep <= dpho: lep.mcPho = None
598 
599  def process(self, event):
600  self.readCollections( event.input )
601  self.counters.counter('events').inc('all events')
602 
603  #call the leptons functions
604  self.makeLeptons(event)
605 
606  if self.cfg_comp.isMC and self.cfg_ana.do_mc_match:
607  self.matchLeptons(event)
608  self.matchAnyLeptons(event)
609  if self.doMatchToPhotons:
610  self.matchToPhotons(event)
611 
612  return True
613 
614 #A default config
615 setattr(LeptonAnalyzer,"defaultConfig",cfg.Analyzer(
616  verbose=False,
617  class_object=LeptonAnalyzer,
618  # input collections
619  muons='slimmedMuons',
620  electrons='slimmedElectrons',
621  rhoMuon= 'fixedGridRhoFastjetAll',
622  rhoElectron = 'fixedGridRhoFastjetAll',
623 ## photons='slimmedPhotons',
624  # energy scale corrections and ghost muon suppression (off by default)
625  doMuonScaleCorrections=False,
626  doElectronScaleCorrections=False,
627  doSegmentBasedMuonCleaning=False,
628  # inclusive very loose muon selection
629  inclusive_muon_id = "POG_ID_Loose",
630  inclusive_muon_pt = 3,
631  inclusive_muon_eta = 2.4,
632  inclusive_muon_dxy = 0.5,
633  inclusive_muon_dz = 1.0,
634  muon_dxydz_track = "muonBestTrack",
635  # loose muon selection
636  loose_muon_id = "POG_ID_Loose",
637  loose_muon_pt = 5,
638  loose_muon_eta = 2.4,
639  loose_muon_dxy = 0.05,
640  loose_muon_dz = 0.2,
641  loose_muon_relIso = 0.4,
642  # loose_muon_isoCut = lambda muon :muon.miniRelIso < 0.2
643  # inclusive very loose electron selection
644  inclusive_electron_id = "",
645  inclusive_electron_pt = 5,
646  inclusive_electron_eta = 2.5,
647  inclusive_electron_dxy = 0.5,
648  inclusive_electron_dz = 1.0,
649  inclusive_electron_lostHits = 1.0,
650  # loose electron selection
651  loose_electron_id = "", #POG_MVA_ID_NonTrig_full5x5",
652  loose_electron_pt = 7,
653  loose_electron_eta = 2.4,
654  loose_electron_dxy = 0.05,
655  loose_electron_dz = 0.2,
656  loose_electron_relIso = 0.4,
657  # loose_electron_isoCut = lambda electron : electron.miniRelIso < 0.1
658  loose_electron_lostHits = 1.0,
659  # muon isolation correction method (can be "rhoArea" or "deltaBeta")
660  mu_isoCorr = "rhoArea" ,
661  mu_effectiveAreas = "Spring15_25ns_v1", #(can be 'Data2012' or 'Phys14_25ns_v1' or 'Spring15_25ns_v1')
662  mu_tightId = "POG_ID_Tight" ,
663  # electron isolation correction method (can be "rhoArea" or "deltaBeta")
664  ele_isoCorr = "rhoArea" ,
665  ele_effectiveAreas = "Spring15_25ns_v1" , #(can be 'Data2012' or 'Phys14_25ns_v1', or 'Spring15_50ns_v1' or 'Spring15_25ns_v1')
666  ele_tightId = "Cuts_2012" ,
667  # minimum deltaR between a loose electron and a loose muon (on overlaps, discard the electron)
668  min_dr_electron_muon = 0.02,
669  # Mini-isolation, with pT dependent cone: will fill in the miniRelIso, miniRelIsoCharged, miniRelIsoNeutral variables of the leptons (see https://indico.cern.ch/event/368826/ )
670  doMiniIsolation = False, # off by default since it requires access to all PFCandidates
671  packedCandidates = 'packedPFCandidates',
672  miniIsolationPUCorr = 'rhoArea', # Allowed options: 'rhoArea' (EAs for 03 cone scaled by R^2), 'deltaBeta', 'raw' (uncorrected), 'weights' (delta beta weights; not validated)
673  # Choose None to just use the individual object's PU correction
674  miniIsolationVetoLeptons = None, # use 'inclusive' to veto inclusive leptons and their footprint in all isolation cones
675  # Activity Annulus
676  doIsoAnnulus = False, # off by default since it requires access to all PFCandidates
677  # do MC matching
678  do_mc_match = True, # note: it will in any case try it only on MC, not on data
679  do_mc_match_photons = False, # mc match electrons to photons
680  match_inclusiveLeptons = False, # match to all inclusive leptons
681  )
682 )
deltar.bestMatch
def bestMatch
Definition: deltar.py:139
HLT_FULL_cff.distance
tuple distance
Definition: HLT_FULL_cff.py:5234
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
min
T min(T a, T b)
Definition: MathUtil.h:58
deltaR
double deltaR(double eta1, double eta2, double phi1, double phi2)
Definition: TreeUtility.cc:17
objects.LeptonAnalyzer.LeptonAnalyzer.eleEffectiveArea
eleEffectiveArea
Duplicate removal for fast sim (to be checked if still necessary in latest greatest 5...
Definition: LeptonAnalyzer.py:338
objects.LeptonAnalyzer.LeptonAnalyzer.miniIsolationVetoLeptons
miniIsolationVetoLeptons
inclusive leptons = all leptons that could be considered somewhere in the analysis, with minimal requirements (used e.g.
Definition: LeptonAnalyzer.py:158
genutils.isPromptLepton
def isPromptLepton
Definition: genutils.py:49
deltar.matchObjectCollection3
def matchObjectCollection3
Definition: deltar.py:41
Generated for CMSSW Reference Manual by doxygen 1.8.5