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