File indexing completed on 2024-04-06 12:23:27
0001 from __future__ import print_function
0002 from builtins import range
0003 from PhysicsTools.Heppy.analyzers.core.Analyzer import Analyzer
0004 from PhysicsTools.Heppy.analyzers.core.AutoHandle import AutoHandle
0005 from PhysicsTools.Heppy.physicsobjects.Electron import Electron
0006 from PhysicsTools.Heppy.physicsobjects.Muon import Muon
0007
0008
0009 from PhysicsTools.HeppyCore.utils.deltar import bestMatch
0010 from PhysicsTools.Heppy.physicsutils.RochesterCorrections import rochcor
0011 from PhysicsTools.Heppy.physicsutils.MuScleFitCorrector import MuScleFitCorr
0012 from PhysicsTools.Heppy.physicsutils.KalmanMuonCorrector import KalmanMuonCorrector
0013 from PhysicsTools.Heppy.physicsutils.ElectronCalibrator import Run2ElectronCalibrator
0014
0015 import PhysicsTools.HeppyCore.framework.config as cfg
0016 from PhysicsTools.HeppyCore.utils.deltar import *
0017 from PhysicsTools.Heppy.physicsutils.genutils import *
0018
0019
0020 from ROOT import heppy
0021 cmgMuonCleanerBySegments = heppy.CMGMuonCleanerBySegmentsAlgo()
0022
0023 class LeptonAnalyzer( Analyzer ):
0024
0025
0026 def __init__(self, cfg_ana, cfg_comp, looperName ):
0027 super(LeptonAnalyzer,self).__init__(cfg_ana,cfg_comp,looperName)
0028 if hasattr(self.cfg_ana, 'doMuScleFitCorrections'):
0029 raise RuntimeError("doMuScleFitCorrections is not supported. Please set instead doMuonScaleCorrections = ( 'MuScleFit', <name> )")
0030 if hasattr(self.cfg_ana, 'doRochesterCorrections'):
0031 raise RuntimeError("doRochesterCorrections is not supported. Please set instead doMuonScaleCorrections = ( 'Rochester', <name> )")
0032 if self.cfg_ana.doMuonScaleCorrections:
0033 algo, options = self.cfg_ana.doMuonScaleCorrections
0034 if algo == "Kalman":
0035 corr = options['MC' if self.cfg_comp.isMC else 'Data']
0036 self.muonScaleCorrector = KalmanMuonCorrector(corr,
0037 self.cfg_comp.isMC,
0038 options['isSync'] if 'isSync' in options else False,
0039 options['smearMode'] if 'smearMode' in options else "ebe")
0040 elif algo == "Rochester":
0041 print("WARNING: the Rochester correction in heppy is still from Run 1")
0042 self.muonScaleCorrector = RochesterCorrections()
0043 elif algo == "MuScleFit":
0044 print("WARNING: the MuScleFit correction in heppy is still from Run 1 (and probably no longer functional)")
0045 if options not in [ "prompt", "prompt-sync", "rereco", "rereco-sync" ]:
0046 raise RuntimeError('MuScleFit correction name must be one of [ "prompt", "prompt-sync", "rereco", "rereco-sync" ] ')
0047 rereco = ("prompt" not in self.cfg_ana.doMuScleFitCorrections)
0048 sync = ("sync" in self.cfg_ana.doMuScleFitCorrections)
0049 self.muonScaleCorrector = MuScleFitCorr(cfg_comp.isMC, rereco, sync)
0050 else: raise RuntimeError("Unknown muon scale correction algorithm")
0051 else:
0052 self.muonScaleCorrector = None
0053
0054 if self.cfg_ana.doElectronScaleCorrections:
0055 conf = cfg_ana.doElectronScaleCorrections
0056 self.electronEnergyCalibrator = Run2ElectronCalibrator(
0057 conf['data'],
0058 conf['GBRForest'],
0059 cfg_comp.isMC,
0060 conf['isSync'] if 'isSync' in conf else False,
0061 )
0062
0063
0064
0065 if hasattr(cfg_ana, 'loose_electron_isoCut'):
0066 self.eleIsoCut = cfg_ana.loose_electron_isoCut
0067 else:
0068 self.eleIsoCut = lambda ele : (
0069 ele.relIso03 <= self.cfg_ana.loose_electron_relIso and
0070 ele.absIso03 < getattr(self.cfg_ana,'loose_electron_absIso',9e99))
0071 if hasattr(cfg_ana, 'loose_muon_isoCut'):
0072 self.muIsoCut = cfg_ana.loose_muon_isoCut
0073 else:
0074 self.muIsoCut = lambda mu : (
0075 mu.relIso03 <= self.cfg_ana.loose_muon_relIso and
0076 mu.absIso03 < getattr(self.cfg_ana,'loose_muon_absIso',9e99))
0077
0078
0079
0080 self.eleEffectiveArea = getattr(cfg_ana, 'ele_effectiveAreas', "Spring15_25ns_v1")
0081 self.muEffectiveArea = getattr(cfg_ana, 'mu_effectiveAreas', "Spring15_25ns_v1")
0082
0083 self.doMiniIsolation = getattr(cfg_ana, 'doMiniIsolation', False)
0084 if self.doMiniIsolation:
0085 self.miniIsolationPUCorr = self.cfg_ana.miniIsolationPUCorr
0086 self.miniIsolationVetoLeptons = self.cfg_ana.miniIsolationVetoLeptons
0087 if self.miniIsolationVetoLeptons not in [ None, 'any', 'inclusive' ]:
0088 raise RuntimeError("miniIsolationVetoLeptons should be None, or 'any' (all reco leptons), or 'inclusive' (all inclusive leptons)")
0089 if self.miniIsolationPUCorr == "weights":
0090 self.IsolationComputer = heppy.IsolationComputer(0.4)
0091 else:
0092 self.IsolationComputer = heppy.IsolationComputer()
0093
0094 self.doIsoAnnulus = getattr(cfg_ana, 'doIsoAnnulus', False)
0095 if self.doIsoAnnulus:
0096 if not self.doMiniIsolation:
0097 self.IsolationComputer = heppy.IsolationComputer()
0098
0099 self.doIsolationScan = getattr(cfg_ana, 'doIsolationScan', False)
0100 if self.doIsolationScan:
0101 if self.doMiniIsolation:
0102 assert (self.miniIsolationPUCorr!="weights")
0103 assert (self.miniIsolationVetoLeptons==None)
0104 else:
0105 self.IsolationComputer = heppy.IsolationComputer()
0106
0107
0108 self.doMatchToPhotons = getattr(cfg_ana, 'do_mc_match_photons', False)
0109
0110
0111
0112
0113
0114
0115 def declareHandles(self):
0116 super(LeptonAnalyzer, self).declareHandles()
0117
0118
0119 self.handles['muons'] = AutoHandle(self.cfg_ana.muons,"std::vector<pat::Muon>")
0120 self.handles['electrons'] = AutoHandle(self.cfg_ana.electrons,"std::vector<pat::Electron>")
0121
0122
0123 self.handles['rhoMu'] = AutoHandle( self.cfg_ana.rhoMuon, 'double')
0124
0125 self.handles['rhoEle'] = AutoHandle( self.cfg_ana.rhoElectron, 'double')
0126
0127 if self.doMiniIsolation or self.doIsolationScan:
0128 self.handles['packedCandidates'] = AutoHandle( self.cfg_ana.packedCandidates, 'std::vector<pat::PackedCandidate>')
0129
0130 if self.doMatchToPhotons:
0131 if self.doMatchToPhotons == "any":
0132 self.mchandles['genPhotons'] = AutoHandle( 'packedGenParticles', 'std::vector<pat::PackedGenParticle>' )
0133 else:
0134 self.mchandles['genPhotons'] = AutoHandle( 'prunedGenParticles', 'std::vector<reco::GenParticle>' )
0135
0136 def beginLoop(self, setup):
0137 super(LeptonAnalyzer,self).beginLoop(setup)
0138 self.counters.addCounter('events')
0139 count = self.counters.counter('events')
0140 count.register('all events')
0141
0142
0143
0144
0145
0146
0147 def makeLeptons(self, event):
0148
0149 event.inclusiveLeptons = []
0150
0151
0152 event.selectedLeptons = []
0153 event.selectedMuons = []
0154 event.selectedElectrons = []
0155 event.otherLeptons = []
0156
0157 if self.doMiniIsolation or self.doIsolationScan:
0158 self.IsolationComputer.setPackedCandidates(self.handles['packedCandidates'].product())
0159 if self.doMiniIsolation:
0160 if self.miniIsolationVetoLeptons == "any":
0161 for lep in self.handles['muons'].product():
0162 self.IsolationComputer.addVetos(lep)
0163 for lep in self.handles['electrons'].product():
0164 self.IsolationComputer.addVetos(lep)
0165
0166
0167 allmuons = self.makeAllMuons(event)
0168
0169
0170 allelectrons = self.makeAllElectrons(event)
0171
0172
0173 inclusiveMuons = []
0174 inclusiveElectrons = []
0175 for mu in allmuons:
0176 if (mu.track().isNonnull() and mu.muonID(self.cfg_ana.inclusive_muon_id) and
0177 mu.pt()>self.cfg_ana.inclusive_muon_pt and abs(mu.eta())<self.cfg_ana.inclusive_muon_eta and
0178 abs(mu.dxy())<self.cfg_ana.inclusive_muon_dxy and abs(mu.dz())<self.cfg_ana.inclusive_muon_dz):
0179 inclusiveMuons.append(mu)
0180 for ele in allelectrons:
0181 if ( ele.electronID(self.cfg_ana.inclusive_electron_id) and
0182 ele.pt()>self.cfg_ana.inclusive_electron_pt and abs(ele.eta())<self.cfg_ana.inclusive_electron_eta and
0183 abs(ele.dxy())<self.cfg_ana.inclusive_electron_dxy and abs(ele.dz())<self.cfg_ana.inclusive_electron_dz and
0184 ele.lostInner()<=self.cfg_ana.inclusive_electron_lostHits ):
0185 inclusiveElectrons.append(ele)
0186 event.inclusiveLeptons = inclusiveMuons + inclusiveElectrons
0187
0188 if self.doMiniIsolation:
0189 if self.miniIsolationVetoLeptons == "inclusive":
0190 for lep in event.inclusiveLeptons:
0191 self.IsolationComputer.addVetos(lep.physObj)
0192 for lep in event.inclusiveLeptons:
0193 self.attachMiniIsolation(lep)
0194
0195 if self.doIsoAnnulus:
0196 for lep in event.inclusiveLeptons:
0197 self.attachIsoAnnulus04(lep)
0198
0199 if self.doIsolationScan:
0200 for lep in event.inclusiveLeptons:
0201 self.attachIsolationScan(lep)
0202
0203
0204 for mu in inclusiveMuons:
0205 if (mu.muonID(self.cfg_ana.loose_muon_id) and
0206 mu.pt() > self.cfg_ana.loose_muon_pt and abs(mu.eta()) < self.cfg_ana.loose_muon_eta and
0207 abs(mu.dxy()) < self.cfg_ana.loose_muon_dxy and abs(mu.dz()) < self.cfg_ana.loose_muon_dz and
0208 self.muIsoCut(mu)):
0209 mu.looseIdSusy = True
0210 event.selectedLeptons.append(mu)
0211 event.selectedMuons.append(mu)
0212 else:
0213 mu.looseIdSusy = False
0214 event.otherLeptons.append(mu)
0215 looseMuons = event.selectedLeptons[:]
0216 for ele in inclusiveElectrons:
0217 if (ele.electronID(self.cfg_ana.loose_electron_id) and
0218 ele.pt()>self.cfg_ana.loose_electron_pt and abs(ele.eta())<self.cfg_ana.loose_electron_eta and
0219 abs(ele.dxy()) < self.cfg_ana.loose_electron_dxy and abs(ele.dz())<self.cfg_ana.loose_electron_dz and
0220 self.eleIsoCut(ele) and
0221 ele.lostInner() <= self.cfg_ana.loose_electron_lostHits and
0222 ( True if getattr(self.cfg_ana,'notCleaningElectrons',False) else (bestMatch(ele, looseMuons)[1] > (self.cfg_ana.min_dr_electron_muon**2)) )):
0223 event.selectedLeptons.append(ele)
0224 event.selectedElectrons.append(ele)
0225 ele.looseIdSusy = True
0226 else:
0227 event.otherLeptons.append(ele)
0228 ele.looseIdSusy = False
0229
0230 event.otherLeptons.sort(key = lambda l : l.pt(), reverse = True)
0231 event.selectedLeptons.sort(key = lambda l : l.pt(), reverse = True)
0232 event.selectedMuons.sort(key = lambda l : l.pt(), reverse = True)
0233 event.selectedElectrons.sort(key = lambda l : l.pt(), reverse = True)
0234 event.inclusiveLeptons.sort(key = lambda l : l.pt(), reverse = True)
0235
0236 for lepton in event.selectedLeptons:
0237 if hasattr(self,'efficiency'):
0238 self.efficiency.attachToObject(lepton)
0239
0240 def makeAllMuons(self, event):
0241 """
0242 make a list of all muons, and apply basic corrections to them
0243 """
0244
0245 allmuons = map( Muon, self.handles['muons'].product() )
0246
0247
0248 if self.muonScaleCorrector:
0249 self.muonScaleCorrector.correct_all(allmuons, event.run)
0250
0251
0252 if self.cfg_ana.doSegmentBasedMuonCleaning:
0253 isgood = cmgMuonCleanerBySegments.clean( self.handles['muons'].product() )
0254 newmu = []
0255 for i,mu in enumerate(allmuons):
0256 if isgood[i]: newmu.append(mu)
0257 allmuons = newmu
0258
0259
0260 for mu in allmuons:
0261 mu.rho = float(self.handles['rhoMu'].product()[0])
0262 if self.muEffectiveArea == "Data2012":
0263 if aeta < 1.0 : mu.EffectiveArea03 = 0.382;
0264 elif aeta < 1.47 : mu.EffectiveArea03 = 0.317;
0265 elif aeta < 2.0 : mu.EffectiveArea03 = 0.242;
0266 elif aeta < 2.2 : mu.EffectiveArea03 = 0.326;
0267 elif aeta < 2.3 : mu.EffectiveArea03 = 0.462;
0268 else : mu.EffectiveArea03 = 0.372;
0269 if aeta < 1.0 : mu.EffectiveArea04 = 0.674;
0270 elif aeta < 1.47 : mu.EffectiveArea04 = 0.565;
0271 elif aeta < 2.0 : mu.EffectiveArea04 = 0.442;
0272 elif aeta < 2.2 : mu.EffectiveArea04 = 0.515;
0273 elif aeta < 2.3 : mu.EffectiveArea04 = 0.821;
0274 else : mu.EffectiveArea04 = 0.660;
0275 elif self.muEffectiveArea == "Phys14_25ns_v1":
0276 aeta = abs(mu.eta())
0277 if aeta < 0.800: mu.EffectiveArea03 = 0.0913
0278 elif aeta < 1.300: mu.EffectiveArea03 = 0.0765
0279 elif aeta < 2.000: mu.EffectiveArea03 = 0.0546
0280 elif aeta < 2.200: mu.EffectiveArea03 = 0.0728
0281 else: mu.EffectiveArea03 = 0.1177
0282 if aeta < 0.800: mu.EffectiveArea04 = 0.1564
0283 elif aeta < 1.300: mu.EffectiveArea04 = 0.1325
0284 elif aeta < 2.000: mu.EffectiveArea04 = 0.0913
0285 elif aeta < 2.200: mu.EffectiveArea04 = 0.1212
0286 else: mu.EffectiveArea04 = 0.2085
0287 elif self.muEffectiveArea == "Spring15_25ns_v1":
0288 aeta = abs(mu.eta())
0289 if aeta < 0.800: mu.EffectiveArea03 = 0.0735
0290 elif aeta < 1.300: mu.EffectiveArea03 = 0.0619
0291 elif aeta < 2.000: mu.EffectiveArea03 = 0.0465
0292 elif aeta < 2.200: mu.EffectiveArea03 = 0.0433
0293 else: mu.EffectiveArea03 = 0.0577
0294 mu.EffectiveArea04 = 0
0295 else: raise RuntimeError("Unsupported value for mu_effectiveAreas: can only use Data2012 (rho: ?) and Phys14_25ns_v1 or Spring15_25ns_v1 (rho: fixedGridRhoFastjetAll)")
0296
0297 for mu in allmuons:
0298 mu.associatedVertex = event.goodVertices[0] if len(event.goodVertices)>0 else event.vertices[0]
0299 mu.setTrackForDxyDz(self.cfg_ana.muon_dxydz_track)
0300
0301
0302 if hasattr(self.cfg_ana, "mu_tightId"):
0303 for mu in allmuons:
0304 mu.tightIdResult = mu.muonID(self.cfg_ana.mu_tightId)
0305
0306
0307 for mu in allmuons:
0308 if self.cfg_ana.mu_isoCorr=="rhoArea" :
0309 mu.absIso03 = (mu.pfIsolationR03().sumChargedHadronPt + max( mu.pfIsolationR03().sumNeutralHadronEt + mu.pfIsolationR03().sumPhotonEt - mu.rho * mu.EffectiveArea03,0.0))
0310 mu.absIso04 = (mu.pfIsolationR04().sumChargedHadronPt + max( mu.pfIsolationR04().sumNeutralHadronEt + mu.pfIsolationR04().sumPhotonEt - mu.rho * mu.EffectiveArea04,0.0))
0311 elif self.cfg_ana.mu_isoCorr=="deltaBeta" :
0312 mu.absIso03 = (mu.pfIsolationR03().sumChargedHadronPt + max( mu.pfIsolationR03().sumNeutralHadronEt + mu.pfIsolationR03().sumPhotonEt - mu.pfIsolationR03().sumPUPt/2,0.0))
0313 mu.absIso04 = (mu.pfIsolationR04().sumChargedHadronPt + max( mu.pfIsolationR04().sumNeutralHadronEt + mu.pfIsolationR04().sumPhotonEt - mu.pfIsolationR04().sumPUPt/2,0.0))
0314 else :
0315 raise RuntimeError("Unsupported mu_isoCorr name '" + str(self.cfg_ana.mu_isoCorr) + "'! For now only 'rhoArea' and 'deltaBeta' are supported.")
0316 mu.relIso03 = mu.absIso03/mu.pt()
0317 mu.relIso04 = mu.absIso04/mu.pt()
0318 return allmuons
0319
0320 def makeAllElectrons(self, event):
0321 """
0322 make a list of all electrons, and apply basic corrections to them
0323 """
0324 allelectrons = map( Electron, self.handles['electrons'].product() )
0325
0326
0327 allelenodup = []
0328 for e in allelectrons:
0329 dup = False
0330 for e2 in allelenodup:
0331 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():
0332 dup = True
0333 break
0334 if not dup: allelenodup.append(e)
0335 allelectrons = allelenodup
0336
0337
0338 for ele in allelectrons:
0339 ele.rho = float(self.handles['rhoEle'].product()[0])
0340 if self.eleEffectiveArea == "Data2012":
0341
0342 SCEta = abs(ele.superCluster().eta())
0343 if SCEta < 1.0 : ele.EffectiveArea03 = 0.13
0344 elif SCEta < 1.479: ele.EffectiveArea03 = 0.14
0345 elif SCEta < 2.0 : ele.EffectiveArea03 = 0.07
0346 elif SCEta < 2.2 : ele.EffectiveArea03 = 0.09
0347 elif SCEta < 2.3 : ele.EffectiveArea03 = 0.11
0348 elif SCEta < 2.4 : ele.EffectiveArea03 = 0.11
0349 else : ele.EffectiveArea03 = 0.14
0350 if SCEta < 1.0 : ele.EffectiveArea04 = 0.208;
0351 elif SCEta < 1.479: ele.EffectiveArea04 = 0.209;
0352 elif SCEta < 2.0 : ele.EffectiveArea04 = 0.115;
0353 elif SCEta < 2.2 : ele.EffectiveArea04 = 0.143;
0354 elif SCEta < 2.3 : ele.EffectiveArea04 = 0.183;
0355 elif SCEta < 2.4 : ele.EffectiveArea04 = 0.194;
0356 else : ele.EffectiveArea04 = 0.261;
0357 elif self.eleEffectiveArea == "Phys14_25ns_v1":
0358 aeta = abs(ele.eta())
0359 if aeta < 0.800: ele.EffectiveArea03 = 0.1013
0360 elif aeta < 1.300: ele.EffectiveArea03 = 0.0988
0361 elif aeta < 2.000: ele.EffectiveArea03 = 0.0572
0362 elif aeta < 2.200: ele.EffectiveArea03 = 0.0842
0363 else: ele.EffectiveArea03 = 0.1530
0364 if aeta < 0.800: ele.EffectiveArea04 = 0.1830
0365 elif aeta < 1.300: ele.EffectiveArea04 = 0.1734
0366 elif aeta < 2.000: ele.EffectiveArea04 = 0.1077
0367 elif aeta < 2.200: ele.EffectiveArea04 = 0.1565
0368 else: ele.EffectiveArea04 = 0.2680
0369 elif self.eleEffectiveArea == "Spring15_50ns_v1":
0370 SCEta = abs(ele.superCluster().eta())
0371
0372 if SCEta < 0.800: ele.EffectiveArea03 = 0.0973
0373 elif SCEta < 1.300: ele.EffectiveArea03 = 0.0954
0374 elif SCEta < 2.000: ele.EffectiveArea03 = 0.0632
0375 elif SCEta < 2.200: ele.EffectiveArea03 = 0.0727
0376 else: ele.EffectiveArea03 = 0.1337
0377
0378 ele.EffectiveArea04 = 0.0
0379 elif self.eleEffectiveArea == "Spring15_25ns_v1":
0380 SCEta = abs(ele.superCluster().eta())
0381
0382 if SCEta < 1.000: ele.EffectiveArea03 = 0.1752
0383 elif SCEta < 1.479: ele.EffectiveArea03 = 0.1862
0384 elif SCEta < 2.000: ele.EffectiveArea03 = 0.1411
0385 elif SCEta < 2.200: ele.EffectiveArea03 = 0.1534
0386 elif SCEta < 2.300: ele.EffectiveArea03 = 0.1903
0387 elif SCEta < 2.400: ele.EffectiveArea03 = 0.2243
0388 else: ele.EffectiveArea03 = 0.2687
0389
0390 ele.EffectiveArea04 = 0.0
0391 else: raise RuntimeError("Unsupported value for ele_effectiveAreas: can only use Data2012 (rho: ?), Phys14_v1 and Spring15_v1 (rho: fixedGridRhoFastjetAll)")
0392
0393
0394 if self.cfg_ana.doElectronScaleCorrections:
0395 for ele in allelectrons:
0396 self.electronEnergyCalibrator.correct(ele, event.run)
0397
0398
0399 for ele in allelectrons:
0400 ele.associatedVertex = event.goodVertices[0] if len(event.goodVertices)>0 else event.vertices[0]
0401
0402
0403 for ele in allelectrons:
0404 if self.cfg_ana.ele_isoCorr=="rhoArea" :
0405 ele.absIso03 = (ele.chargedHadronIsoR(0.3) + max(ele.neutralHadronIsoR(0.3)+ele.photonIsoR(0.3)-ele.rho*ele.EffectiveArea03,0))
0406 ele.absIso04 = (ele.chargedHadronIsoR(0.4) + max(ele.neutralHadronIsoR(0.4)+ele.photonIsoR(0.4)-ele.rho*ele.EffectiveArea04,0))
0407 elif self.cfg_ana.ele_isoCorr=="deltaBeta" :
0408 ele.absIso03 = (ele.chargedHadronIsoR(0.3) + max( ele.neutralHadronIsoR(0.3)+ele.photonIsoR(0.3) - ele.puChargedHadronIsoR(0.3)/2, 0.0))
0409 ele.absIso04 = (ele.chargedHadronIsoR(0.4) + max( ele.neutralHadronIsoR(0.4)+ele.photonIsoR(0.4) - ele.puChargedHadronIsoR(0.4)/2, 0.0))
0410 else :
0411 raise RuntimeError("Unsupported ele_isoCorr name '" + str(self.cfg_ana.ele_isoCorr) + "'! For now only 'rhoArea' and 'deltaBeta' are supported.")
0412 ele.relIso03 = ele.absIso03/ele.pt()
0413 ele.relIso04 = ele.absIso04/ele.pt()
0414
0415
0416 for ele in allelectrons:
0417 if self.cfg_ana.ele_tightId=="MVA" :
0418 ele.tightIdResult = ele.electronID("POG_MVA_ID_Trig_full5x5")
0419 elif self.cfg_ana.ele_tightId=="Cuts_2012" :
0420 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")
0421 elif self.cfg_ana.ele_tightId=="Cuts_PHYS14_25ns_v1_ConvVetoDxyDz" :
0422 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")
0423 elif self.cfg_ana.ele_tightId=="Cuts_SPRING15_25ns_v1_ConvVetoDxyDz" :
0424 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")
0425
0426 else :
0427 try:
0428 ele.tightIdResult = ele.electronID(self.cfg_ana.ele_tightId)
0429 except RuntimeError:
0430 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.")
0431
0432
0433 return allelectrons
0434
0435 def attachMiniIsolation(self, mu):
0436 mu.miniIsoR = 10.0/min(max(mu.pt(), 50),200)
0437
0438
0439 what = "mu" if (abs(mu.pdgId()) == 13) else ("eleB" if mu.isEB() else "eleE")
0440 if what == "mu":
0441 mu.miniAbsIsoCharged = self.IsolationComputer.chargedAbsIso(mu.physObj, mu.miniIsoR, {"mu":0.0001,"eleB":0,"eleE":0.015}[what], 0.0);
0442 else:
0443 mu.miniAbsIsoCharged = self.IsolationComputer.chargedAbsIso(mu.physObj, mu.miniIsoR, {"mu":0.0001,"eleB":0,"eleE":0.015}[what], 0.0,self.IsolationComputer.selfVetoNone);
0444
0445 if self.miniIsolationPUCorr == None: puCorr = self.cfg_ana.mu_isoCorr if what=="mu" else self.cfg_ana.ele_isoCorr
0446 else: puCorr = self.miniIsolationPUCorr
0447
0448 if puCorr == "weights":
0449 if what == "mu":
0450 mu.miniAbsIsoNeutral = self.IsolationComputer.neutralAbsIsoWeighted(mu.physObj, mu.miniIsoR, 0.01, 0.5);
0451 else:
0452 mu.miniAbsIsoNeutral = ( self.IsolationComputer.photonAbsIsoWeighted( mu.physObj, mu.miniIsoR, 0.08 if what == "eleE" else 0.0, 0.0, self.IsolationComputer.selfVetoNone) +
0453 self.IsolationComputer.neutralHadAbsIsoWeighted(mu.physObj, mu.miniIsoR, 0.0, 0.0, self.IsolationComputer.selfVetoNone) )
0454 else:
0455 if what == "mu":
0456 mu.miniAbsIsoNeutral = self.IsolationComputer.neutralAbsIsoRaw(mu.physObj, mu.miniIsoR, 0.01, 0.5);
0457 else:
0458 mu.miniAbsIsoPho = self.IsolationComputer.photonAbsIsoRaw( mu.physObj, mu.miniIsoR, 0.08 if what == "eleE" else 0.0, 0.0, self.IsolationComputer.selfVetoNone)
0459 mu.miniAbsIsoNHad = self.IsolationComputer.neutralHadAbsIsoRaw(mu.physObj, mu.miniIsoR, 0.0, 0.0, self.IsolationComputer.selfVetoNone)
0460 mu.miniAbsIsoNeutral = mu.miniAbsIsoPho + mu.miniAbsIsoNHad
0461
0462
0463
0464
0465 if puCorr == "rhoArea":
0466 mu.miniAbsIsoNeutral = max(0.0, mu.miniAbsIsoNeutral - mu.rho * mu.EffectiveArea03 * (mu.miniIsoR/0.3)**2)
0467 elif puCorr == "deltaBeta":
0468 if what == "mu":
0469 mu.miniAbsIsoPU = self.IsolationComputer.puAbsIso(mu.physObj, mu.miniIsoR, 0.01, 0.5);
0470 else:
0471 mu.miniAbsIsoPU = self.IsolationComputer.puAbsIso(mu.physObj, mu.miniIsoR, 0.015 if what == "eleE" else 0.0, 0.0,self.IsolationComputer.selfVetoNone);
0472 mu.miniAbsIsoNeutral = max(0.0, mu.miniAbsIsoNeutral - 0.5*mu.miniAbsIsoPU)
0473 elif puCorr != 'raw':
0474 raise RuntimeError("Unsupported miniIsolationCorr name '" + puCorr + "'! For now only 'rhoArea', 'deltaBeta', 'raw', 'weights' are supported (and 'weights' is not tested).")
0475
0476 mu.miniAbsIso = mu.miniAbsIsoCharged + mu.miniAbsIsoNeutral
0477 mu.miniRelIso = mu.miniAbsIso/mu.pt()
0478
0479
0480 def attachIsoAnnulus04(self, mu):
0481 mu.miniIsoR = 10.0/min(max(mu.pt(), 50),200)
0482 mu.absIsoAnCharged = self.IsolationComputer.chargedAbsIso (mu.physObj, 0.4, mu.miniIsoR, 0.0, self.IsolationComputer.selfVetoNone)
0483 mu.absIsoAnPho = self.IsolationComputer.photonAbsIsoRaw (mu.physObj, 0.4, mu.miniIsoR, 0.0, self.IsolationComputer.selfVetoNone)
0484 mu.absIsoAnNHad = self.IsolationComputer.neutralHadAbsIsoRaw(mu.physObj, 0.4, mu.miniIsoR, 0.0, self.IsolationComputer.selfVetoNone)
0485 mu.absIsoAnPU = self.IsolationComputer.puAbsIso (mu.physObj, 0.4, mu.miniIsoR, 0.0, self.IsolationComputer.selfVetoNone)
0486 mu.absIsoAnNeutral = max(0.0, mu.absIsoAnPho + mu.absIsoAnNHad - 0.5*mu.absIsoAnPU)
0487
0488 mu.absIsoAn04 = mu.absIsoAnCharged + mu.absIsoAnNeutral
0489 mu.relIsoAn04 = mu.absIsoAn04/mu.pt()
0490
0491
0492 def attachIsolationScan(self, mu):
0493
0494 what = "mu" if (abs(mu.pdgId()) == 13) else ("eleB" if mu.isEB() else "eleE")
0495 vetoreg = {"mu":0.0001,"eleB":0,"eleE":0.015}[what]
0496
0497 if what=="mu":
0498 mu.ScanAbsIsoCharged005 = self.IsolationComputer.chargedAbsIso(mu.physObj, 0.05, vetoreg, 0.0)
0499 mu.ScanAbsIsoCharged01 = self.IsolationComputer.chargedAbsIso(mu.physObj, 0.1, vetoreg, 0.0)
0500 mu.ScanAbsIsoCharged02 = self.IsolationComputer.chargedAbsIso(mu.physObj, 0.2, vetoreg, 0.0)
0501 mu.ScanAbsIsoCharged03 = self.IsolationComputer.chargedAbsIso(mu.physObj, 0.3, vetoreg, 0.0)
0502 mu.ScanAbsIsoCharged04 = self.IsolationComputer.chargedAbsIso(mu.physObj, 0.4, vetoreg, 0.0)
0503 else:
0504 mu.ScanAbsIsoCharged005 = self.IsolationComputer.chargedAbsIso(mu.physObj, 0.05, vetoreg, 0.0, self.IsolationComputer.selfVetoNone)
0505 mu.ScanAbsIsoCharged01 = self.IsolationComputer.chargedAbsIso(mu.physObj, 0.1, vetoreg, 0.0, self.IsolationComputer.selfVetoNone)
0506 mu.ScanAbsIsoCharged02 = self.IsolationComputer.chargedAbsIso(mu.physObj, 0.2, vetoreg, 0.0, self.IsolationComputer.selfVetoNone)
0507 mu.ScanAbsIsoCharged03 = self.IsolationComputer.chargedAbsIso(mu.physObj, 0.3, vetoreg, 0.0, self.IsolationComputer.selfVetoNone)
0508 mu.ScanAbsIsoCharged04 = self.IsolationComputer.chargedAbsIso(mu.physObj, 0.4, vetoreg, 0.0, self.IsolationComputer.selfVetoNone)
0509
0510 if what=="mu":
0511 mu.ScanAbsIsoNeutral005 = self.IsolationComputer.neutralAbsIsoRaw(mu.physObj, 0.05, 0.01, 0.5)
0512 mu.ScanAbsIsoNeutral01 = self.IsolationComputer.neutralAbsIsoRaw(mu.physObj, 0.1, 0.01, 0.5)
0513 mu.ScanAbsIsoNeutral02 = self.IsolationComputer.neutralAbsIsoRaw(mu.physObj, 0.2, 0.01, 0.5)
0514 mu.ScanAbsIsoNeutral03 = self.IsolationComputer.neutralAbsIsoRaw(mu.physObj, 0.3, 0.01, 0.5)
0515 mu.ScanAbsIsoNeutral04 = self.IsolationComputer.neutralAbsIsoRaw(mu.physObj, 0.4, 0.01, 0.5)
0516 else:
0517 vetoreg = {"eleB":0.0,"eleE":0.08}[what]
0518 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)
0519 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)
0520 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)
0521 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)
0522 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)
0523
0524
0525 def matchLeptons(self, event):
0526 def plausible(rec,gen):
0527 if abs(rec.pdgId()) == 11 and abs(gen.pdgId()) != 11: return False
0528 if abs(rec.pdgId()) == 13 and abs(gen.pdgId()) != 13: return False
0529 dr = deltaR(rec.eta(),rec.phi(),gen.eta(),gen.phi())
0530 if dr < 0.3: return True
0531 if rec.pt() < 10 and abs(rec.pdgId()) == 13 and gen.pdgId() != rec.pdgId(): return False
0532 if dr < 0.7: return True
0533 if min(rec.pt(),gen.pt())/max(rec.pt(),gen.pt()) < 0.3: return False
0534 return True
0535
0536 leps = event.inclusiveLeptons if self.cfg_ana.match_inclusiveLeptons else event.selectedLeptons
0537 match = matchObjectCollection3(leps,
0538 event.genleps + event.gentauleps,
0539 deltaRMax = 1.2, filter = plausible)
0540 for lep in leps:
0541 gen = match[lep]
0542 lep.mcMatchId = (gen.sourceId if gen != None else 0)
0543 lep.mcMatchTau = (gen in event.gentauleps if gen else -99)
0544 lep.mcLep=gen
0545
0546 def isFromB(self,particle,bid=5, done={}):
0547 for i in range( particle.numberOfMothers() ):
0548 mom = particle.mother(i)
0549 momid = abs(mom.pdgId())
0550 if momid / 1000 == bid or momid / 100 == bid or momid == bid:
0551 return True
0552 elif mom.status() == 2 and self.isFromB(mom, done=done, bid=bid):
0553 return True
0554 return False
0555
0556 def matchAnyLeptons(self, event):
0557 event.anyLeptons = [ x for x in event.genParticles if x.status() == 1 and abs(x.pdgId()) in [11,13] ]
0558 leps = event.inclusiveLeptons if hasattr(event, 'inclusiveLeptons') else event.selectedLeptons
0559 match = matchObjectCollection3(leps, event.anyLeptons, deltaRMax = 0.3, filter = lambda x,y : abs(x.pdgId()) == abs(y.pdgId()))
0560 for lep in leps:
0561 gen = match[lep]
0562 lep.mcMatchAny_gp = gen
0563 if gen:
0564 if self.isFromB(gen): lep.mcMatchAny = 5
0565 elif self.isFromB(gen,bid=4): lep.mcMatchAny = 4
0566 else: lep.mcMatchAny = 1
0567 if not getattr(lep, 'mcLep', None): lep.mcLep = gen
0568 else:
0569 if not getattr(lep, 'mcLep', None): lep.mcLep = None
0570 lep.mcMatchAny = 0
0571
0572 if gen != None and hasattr(lep,'mcMatchId') and lep.mcMatchId == 0:
0573 if isPromptLepton(gen, False) or (gen.isPromptFinalState() or gen.isDirectPromptTauDecayProductFinalState()):
0574 lep.mcMatchId = 100
0575 lep.mcLep = gen
0576 elif not hasattr(lep,'mcMatchId'):
0577 lep.mcMatchId = 0
0578 if not hasattr(lep,'mcMatchTau'): lep.mcMatchTau = 0
0579
0580 def matchToPhotons(self, event):
0581 event.anyPho = [ x for x in self.mchandles['genPhotons'].product() if x.status() == 1 and x.pdgId() == 22 and x.pt() > 1.0 ]
0582 leps = event.inclusiveLeptons if hasattr(event, 'inclusiveLeptons') else event.selectedLeptons
0583 leps = [ l for l in leps if abs(l.pdgId()) == 11 ]
0584 plausible = lambda rec, gen : 0.3*gen.pt() < rec.pt() and rec.pt() < 1.5*gen.pt()
0585 match = matchObjectCollection3(leps, event.anyPho, deltaRMax = 0.3, filter = plausible)
0586 for lep in leps:
0587 gen = match[lep]
0588 lep.mcPho = gen
0589 if lep.mcPho and lep.mcLep:
0590
0591 def distance(rec,gen):
0592 dr = deltaR(rec.eta(),rec.phi(),gen.eta(),gen.phi())
0593 dptRel = abs(rec.pt()-gen.pt())/gen.pt()
0594 return dr + 0.2*dptRel
0595 dpho = distance(lep,lep.mcPho)
0596 dlep = distance(lep,lep.mcLep)
0597 if getattr(lep,'mcMatchAny_gp',None) and lep.mcMatchAny_gp != lep.mcLep:
0598 dlep = min(dlep, distance(lep,lep.mcMatchAny_gp))
0599 if dlep <= dpho: lep.mcPho = None
0600
0601 def process(self, event):
0602 self.readCollections( event.input )
0603 self.counters.counter('events').inc('all events')
0604
0605
0606 self.makeLeptons(event)
0607
0608 if self.cfg_comp.isMC and self.cfg_ana.do_mc_match:
0609 self.matchLeptons(event)
0610 self.matchAnyLeptons(event)
0611 if self.doMatchToPhotons:
0612 self.matchToPhotons(event)
0613
0614 return True
0615
0616
0617 setattr(LeptonAnalyzer,"defaultConfig",cfg.Analyzer(
0618 verbose=False,
0619 class_object=LeptonAnalyzer,
0620
0621 muons='slimmedMuons',
0622 electrons='slimmedElectrons',
0623 rhoMuon= 'fixedGridRhoFastjetAll',
0624 rhoElectron = 'fixedGridRhoFastjetAll',
0625
0626
0627 doMuonScaleCorrections=False,
0628 doElectronScaleCorrections=False,
0629 doSegmentBasedMuonCleaning=False,
0630
0631 inclusive_muon_id = "POG_ID_Loose",
0632 inclusive_muon_pt = 3,
0633 inclusive_muon_eta = 2.4,
0634 inclusive_muon_dxy = 0.5,
0635 inclusive_muon_dz = 1.0,
0636 muon_dxydz_track = "muonBestTrack",
0637
0638 loose_muon_id = "POG_ID_Loose",
0639 loose_muon_pt = 5,
0640 loose_muon_eta = 2.4,
0641 loose_muon_dxy = 0.05,
0642 loose_muon_dz = 0.2,
0643 loose_muon_relIso = 0.4,
0644
0645
0646 inclusive_electron_id = "",
0647 inclusive_electron_pt = 5,
0648 inclusive_electron_eta = 2.5,
0649 inclusive_electron_dxy = 0.5,
0650 inclusive_electron_dz = 1.0,
0651 inclusive_electron_lostHits = 1.0,
0652
0653 loose_electron_id = "",
0654 loose_electron_pt = 7,
0655 loose_electron_eta = 2.4,
0656 loose_electron_dxy = 0.05,
0657 loose_electron_dz = 0.2,
0658 loose_electron_relIso = 0.4,
0659
0660 loose_electron_lostHits = 1.0,
0661
0662 mu_isoCorr = "rhoArea" ,
0663 mu_effectiveAreas = "Spring15_25ns_v1",
0664 mu_tightId = "POG_ID_Tight" ,
0665
0666 ele_isoCorr = "rhoArea" ,
0667 ele_effectiveAreas = "Spring15_25ns_v1" ,
0668 ele_tightId = "Cuts_2012" ,
0669
0670 min_dr_electron_muon = 0.02,
0671
0672 doMiniIsolation = False,
0673 packedCandidates = 'packedPFCandidates',
0674 miniIsolationPUCorr = 'rhoArea',
0675
0676 miniIsolationVetoLeptons = None,
0677
0678 doIsoAnnulus = False,
0679
0680 do_mc_match = True,
0681 do_mc_match_photons = False,
0682 match_inclusiveLeptons = False,
0683 )
0684 )