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