analyzers/objects/IsoTrackAnalyzer.py

0001 from __future__ import print_function
0002 import operator
0003 import itertools
0004 import copy
0005 import types
0006
0007 from ROOT import TLorentzVector
0008
0009 from PhysicsTools.Heppy.analyzers.core.Analyzer import Analyzer
0010 from PhysicsTools.HeppyCore.framework.event import Event
0011 from PhysicsTools.HeppyCore.statistics.counter import Counter, Counters
0012 from PhysicsTools.Heppy.analyzers.core.AutoHandle import AutoHandle
0013 from PhysicsTools.Heppy.physicsobjects.Lepton import Lepton
0014 from PhysicsTools.Heppy.physicsobjects.Tau import Tau
0015 from PhysicsTools.Heppy.physicsobjects.IsoTrack import IsoTrack
0016
0017 from PhysicsTools.HeppyCore.utils.deltar import deltaR, deltaPhi, bestMatch , matchObjectCollection3
0018
0019 import PhysicsTools.HeppyCore.framework.config as cfg
0020
0021 from ROOT import heppy
0022
0023
0024
0025
0026 def mtw(x1,x2):
0027     import math
0028     return math.sqrt(2*x1.pt()*x2.pt()*(1-math.cos(x1.phi()-x2.phi())))
0029
0030 def makeNearestLeptons(leptons,track, event):
0031
0032     minDeltaR = 99999
0033
0034     nearestLepton = []
0035     ibest=-1
0036     for i,lepton in enumerate(leptons):
0037         minDeltaRtemp=deltaR(lepton.eta(),lepton.phi(),track.eta(),track.phi())
0038         if minDeltaRtemp < minDeltaR:
0039             minDeltaR = minDeltaRtemp
0040             ibest=i
0041
0042     if len(leptons) > 0 and ibest!=-1:
0043         nearestLepton.append(leptons[ibest])
0044
0045     return nearestLepton
0046
0047 class IsoTrackAnalyzer( Analyzer ):
0048
0049     def __init__(self, cfg_ana, cfg_comp, looperName ):
0050         super(IsoTrackAnalyzer,self).__init__(cfg_ana,cfg_comp,looperName)
0051         self.IsoTrackIsolationComputer = heppy.IsolationComputer()
0052
0053         self.doIsoAnnulus = getattr(cfg_ana, 'doIsoAnnulus', False)
0054
0055
0056     #----------------------------------------
0057     # DECLARATION OF HANDLES OF LEPTONS STUFF
0058     #----------------------------------------
0059     def declareHandles(self):
0060         super(IsoTrackAnalyzer, self).declareHandles()
0061         self.handles['met'] = AutoHandle( 'slimmedMETs', 'std::vector<pat::MET>' )
0062         self.handles['packedCandidates'] = AutoHandle( 'packedPFCandidates', 'std::vector<pat::PackedCandidate>')
0063
0064     def beginLoop(self, setup):
0065         super(IsoTrackAnalyzer,self).beginLoop(setup)
0066         self.counters.addCounter('events')
0067         count = self.counters.counter('events')
0068         count.register('all events')
0069         count.register('has >=1 selected Track')
0070         count.register('has >=1 selected Iso Track')
0071
0072     #------------------
0073     # MAKE LIST
0074     #------------------
0075     def makeIsoTrack(self, event):
0076
0077         event.selectedIsoTrack = []
0078         event.selectedIsoCleanTrack = []
0079         #event.preIsoTrack = []
0080
0081         patcands = self.handles['packedCandidates'].product()
0082
0083         charged = [ p for p in patcands if ( p.charge() != 0 and p.fromPV() > 1 ) ]
0084
0085         self.IsoTrackIsolationComputer.setPackedCandidates(patcands, 1, 9999, 9999.)
0086
0087
0088         alltrack = map( IsoTrack, charged )
0089
0090
0091         for track in alltrack:
0092
0093             if ( abs(track.dz()) > self.cfg_ana.dzMax ): continue
0094             if ( (abs(track.pdgId())!=11) and (abs(track.pdgId())!=13) and (track.pt() < self.cfg_ana.ptMin) ): continue
0095             if ( track.pt() < self.cfg_ana.ptMinEMU ): continue
0096
0097             foundNonIsoTrack = False
0098
0099 ## ===> require is not the leading lepton and opposite to the leading lepton
0100             if( (self.cfg_ana.doSecondVeto) and len(event.selectedLeptons)>0) :
0101                if( deltaR(event.selectedLeptons[0].eta(), event.selectedLeptons[0].phi(), track.eta(), track.phi()) <0.01) : continue
0102                if ( (abs(track.pdgId())!=11) and (abs(track.pdgId())!=13) and (track.charge()*event.selectedLeptons[0].charge()) ): continue
0103
0104
0105 ## ===> Redundant:: require the Track Candidate with a  minimum dz
0106             track.associatedVertex = event.goodVertices[0] if len(event.goodVertices)>0 else event.vertices[0]
0107
0108 ## ===> compute the isolation and find the most isolated track
0109
0110             isoSum = self.IsoTrackIsolationComputer.chargedAbsIso(track.physObj, self.cfg_ana.isoDR, 0., self.cfg_ana.ptPartMin)
0111             if( abs(track.pdgId())==211 ): isoSum = isoSum - track.pt() #BM: this is an ugly hack and it is error prone. It needs to be fixed using the addVeto method properly
0112
0113             if self.cfg_ana.doRelIsolation:
0114                 relIso = (isoSum)/track.pt()
0115                 if ( (abs(track.pdgId())!=11) and (abs(track.pdgId())!=13) and (relIso > self.cfg_ana.MaxIsoSum) ): continue
0116                 elif((relIso > self.cfg_ana.MaxIsoSumEMU)): continue
0117             else:
0118                 if(isoSum > (self.cfg_ana.maxAbsIso)): continue
0119
0120             if self.doIsoAnnulus:
0121                 self.attachIsoAnnulus04(track)
0122
0123             track.absIso = isoSum
0124
0125             #### store a preIso track
0126             #event.preIsoTrack.append(track)
0127
0128 #            if (isoSum < minIsoSum ) :
0129             if self.cfg_ana.doRelIsolation or (track.absIso < min(0.2*track.pt(), self.cfg_ana.maxAbsIso)):
0130                 event.selectedIsoTrack.append(track)
0131
0132                 if self.cfg_ana.doPrune:
0133                     myMet = self.handles['met'].product()[0]
0134                     mtwIsoTrack = mtw(track, myMet)
0135                     if mtwIsoTrack < 100:
0136                         if abs(track.pdgId()) == 11 or abs(track.pdgId()) == 13:
0137                             if track.pt()>5 and track.absIso/track.pt()<0.2:
0138
0139                                 myLeptons = [ l for l in event.selectedLeptons if l.pt() > 10 ]
0140                                 nearestSelectedLeptons = makeNearestLeptons(myLeptons,track, event)
0141                                 if len(nearestSelectedLeptons) > 0:
0142                                     for lep in nearestSelectedLeptons:
0143                                         if deltaR(lep.eta(), lep.phi(), track.eta(), track.phi()) > 0.01:
0144                                             event.selectedIsoCleanTrack.append(track)
0145                                 else:
0146                                     event.selectedIsoCleanTrack.append(track)
0147
0148
0149
0150 ##        alltrack = map( IsoTrack, charged )
0151
0152 ##        for track in alltrack:
0153 ##
0154 ##            foundNonIsoTrack = False
0155 ##
0156 #### ===> require Track Candidate above some pt and charged
0157 ##            if ( (abs(track.pdgId())!=11) and (abs(track.pdgId())!=13) and (track.pt() < self.cfg_ana.ptMin) ): continue
0158 ##            if ( track.pt() < self.cfg_ana.ptMinEMU ): continue
0159 ##
0160 ##
0161 #### ===> require is not the leading lepton and opposite to the leading lepton
0162 ##            if( (self.cfg_ana.doSecondVeto) and len(event.selectedLeptons)>0) :
0163 ##               if( deltaR(event.selectedLeptons[0].eta(), event.selectedLeptons[0].phi(), track.eta(), track.phi()) <0.01) : continue
0164 ##               if ( (abs(track.pdgId())!=11) and (abs(track.pdgId())!=13) and (track.charge()*event.selectedLeptons[0].charge()) ): continue
0165 ##
0166 #### ===> Redundant:: require the Track Candidate with a  minimum dz
0167 ##            track.associatedVertex = event.goodVertices[0]
0168 ##
0169 #### ===> compute the isolation and find the most isolated track
0170 ##
0171 ##            othertracks = [ p for p in charged if( deltaR(p.eta(), p.phi(), track.eta(), track.phi()) < self.cfg_ana.isoDR and p.pt()>self.cfg_ana.ptPartMin ) ]
0172 ##            #othertracks = alltrack
0173 ##
0174 ##            isoSum=0
0175 ##            for part in othertracks:
0176 ##                #### ===> skip pfcands with a pt min (this should be 0)
0177 ##                #if part.pt()<self.cfg_ana.ptPartMin : continue
0178 ##                #### ===> skip pfcands outside the cone (this should be 0.3)
0179 ##                #if deltaR(part.eta(), part.phi(), track.eta(), track.phi()) > self.cfg_ana.isoDR : continue
0180 ##                isoSum += part.pt()
0181 ##                ### break the loop to save time
0182 ##                if(isoSum > (self.cfg_ana.maxAbsIso + track.pt())):
0183 ##                    foundNonIsoTrack = True
0184 ##                    break
0185 ##
0186 ##            if foundNonIsoTrack: continue
0187 ##
0188 ##               ## reset
0189 ##               #isoSum=0
0190 ##               #for part in othertracks :
0191 ##               #### ===> skip pfcands with a pt min (this should be 0)
0192 ##               #    if part.pt()<self.cfg_ana.ptPartMin : continue
0193 ##               #### ===> skip pfcands outside the cone (this should be 0.3)
0194 ##               #    if deltaR(part.eta(), part.phi(), track.eta(), track.phi()) > self.cfg_ana.isoDR : continue
0195 ##               #    isoSum += part.pt()
0196 ##
0197 ##            #    ###            isoSum = isoSum/track.pt()  ## <--- this is for relIso
0198 ##
0199 ##            ### ===> the sum should not contain the track candidate
0200 ##
0201 ##            track.absIso = isoSum - track.pt()
0202 ##
0203 ##            #### store a preIso track
0204 ##            #event.preIsoTrack.append(track)
0205 ##
0206 ###            if (isoSum < minIsoSum ) :
0207 ##            if(track.absIso < min(0.2*track.pt(), self.cfg_ana.maxAbsIso)):
0208 ##                event.selectedIsoTrack.append(track)
0209 ##
0210 ##                if self.cfg_ana.doPrune:
0211 ##                    myMet = self.handles['met'].product()[0]
0212 ##                    mtwIsoTrack = mtw(track, myMet)
0213 ##                    if mtwIsoTrack < 100:
0214 ##                        if abs(track.pdgId()) == 11 or abs(track.pdgId()) == 13:
0215 ##                            if track.pt()>5 and track.absIso/track.pt()<0.2:
0216 ##
0217 ##                                myLeptons = [ l for l in event.selectedLeptons if l.pt() > 10 ]
0218 ##                                nearestSelectedLeptons = makeNearestLeptons(myLeptons,track, event)
0219 ##                                if len(nearestSelectedLeptons) > 0:
0220 ##                                    for lep in nearestSelectedLeptons:
0221 ##                                        if deltaR(lep.eta(), lep.phi(), track.eta(), track.phi()) > 0.1:
0222 ##                                            event.selectedIsoCleanTrack.append(track)
0223 ##                                else:
0224 ##                                    event.selectedIsoCleanTrack.append(track)
0225
0226         event.selectedIsoTrack.sort(key = lambda l : l.pt(), reverse = True)
0227         event.selectedIsoCleanTrack.sort(key = lambda l : l.pt(), reverse = True)
0228
0229         self.counters.counter('events').inc('all events')
0230         #if(len(event.preIsoTrack)): self.counters.counter('events').inc('has >=1 selected Track')
0231         if(len(event.selectedIsoTrack)): self.counters.counter('events').inc('has >=1 selected Iso Track')
0232
0233
0234     def attachIsoAnnulus04(self, mu):
0235         mu.absIsoAnCharged = self.IsoTrackIsolationComputer.chargedAbsIso      (mu.physObj, 0.4, self.cfg_ana.isoDR, 0.0,self.IsoTrackIsolationComputer.selfVetoNone)
0236         mu.absIsoAnPho     = self.IsoTrackIsolationComputer.photonAbsIsoRaw    (mu.physObj, 0.4, self.cfg_ana.isoDR, 0.0,self.IsoTrackIsolationComputer.selfVetoNone)
0237         mu.absIsoAnNHad    = self.IsoTrackIsolationComputer.neutralHadAbsIsoRaw(mu.physObj, 0.4, self.cfg_ana.isoDR, 0.0,self.IsoTrackIsolationComputer.selfVetoNone)
0238         mu.absIsoAnPU      = self.IsoTrackIsolationComputer.puAbsIso           (mu.physObj, 0.4, self.cfg_ana.isoDR, 0.0,self.IsoTrackIsolationComputer.selfVetoNone)
0239         mu.absIsoAnNeutral = max(0.0, mu.absIsoAnPho + mu.absIsoAnNHad - 0.5*mu.absIsoAnPU)
0240
0241         mu.absIsoAn04 = mu.absIsoAnCharged + mu.absIsoAnNeutral
0242         mu.relIsoAn04 = mu.absIsoAn04/mu.pt()
0243
0244
0245     def matchIsoTrack(self, event):
0246         matchTau = matchObjectCollection3(event.selectedIsoTrack, event.gentaus + event.gentauleps + event.genleps, deltaRMax = 0.5)
0247         for lep in event.selectedIsoTrack:
0248             gen = matchTau[lep]
0249             lep.mcMatchId = 1 if gen else 0
0250
0251
0252     def printInfo(self, event):
0253         print('event to Veto')
0254         print('----------------')
0255
0256         if len(event.selectedIsoTrack)>0:
0257             print('lenght: ',len(event.selectedIsoTrack))
0258             print('track candidate pt: ',event.selectedIsoTrack[0].pt())
0259             print('track candidate eta: ',event.selectedIsoTrack[0].eta())
0260             print('track candidate phi: ',event.selectedIsoTrack[0].phi())
0261             print('track candidate mass: ',event.selectedIsoTrack[0].mass())
0262             print('pdgId candidate : ',event.selectedIsoTrack[0].pdgId())
0263             print('dz: ',event.selectedIsoTrack[0].dz())
0264             print('iso: ',event.selectedIsoTrack[0].absIso)
0265             print('matchId: ',event.selectedIsoTrack[0].mcMatchId)
0266
0267 #        for lepton in event.selectedLeptons:
0268 #            print 'good lepton type: ',lepton.pdgId()
0269 #            print 'pt: ',lepton.pt()
0270
0271 #        for tau in event.selectedTaus:
0272 #            print 'good lepton type: ',tau.pdgId()
0273 #            print 'pt: ',tau.pt()
0274
0275         print('----------------')
0276
0277
0278     def process(self, event):
0279
0280         if self.cfg_ana.setOff:
0281             return True
0282
0283         self.readCollections( event.input )
0284         self.makeIsoTrack(event)
0285
0286         if len(event.selectedIsoTrack)==0 : return True
0287
0288 ##        event.pdgIdIsoTrack.append(event.selectedIsoTrack[0].pdgId())
0289 ##        event.isoIsoTrack.append(minIsoSum)
0290 ##        event.dzIsoTrack.append(abs(dz(event.selectedIsoTrack[0])))
0291
0292 ### ===> do matching
0293
0294         if not self.cfg_comp.isMC:
0295             return True
0296
0297         if hasattr(event, 'gentaus') and hasattr(event, 'gentauleps') and hasattr(event, 'genleps') and self.cfg_ana.do_mc_match :
0298             self.matchIsoTrack(event)
0299
0300 ###        self.printInfo(event)
0301
0302 ### ===> do veto if needed
0303
0304 #        if (self.cfg_ana.doSecondVeto and (event.selectedIsoTrack[0].pdgId()!=11) and (event.selectedIsoTrack[0].pdgId()!=12) and event.isoIsoTrack < self.cfg_ana.MaxIsoSum ) :
0305 ###            self.printInfo(event)
0306 #            return False
0307
0308 #        if ((self.cfg_ana.doSecondVeto and event.selectedIsoTrack[0].pdgId()==11 or event.selectedIsoTrack[0].pdgId()==12) and event.isoIsoTrack < self.cfg_ana.MaxIsoSumEMU ) :
0309 ##            self.printInfo(event)
0310 #            return False
0311
0312
0313         return True
0314
0315
0316 setattr(IsoTrackAnalyzer,"defaultConfig",cfg.Analyzer(
0317     class_object=IsoTrackAnalyzer,
0318     setOff=True,
0319     #####
0320     candidates='packedPFCandidates',
0321     candidatesTypes='std::vector<pat::PackedCandidate>',
0322     ptMin = 5, # for pion
0323     ptMinEMU = 5, # for EMU
0324     dzMax = 0.1,
0325     #####
0326     isoDR = 0.3,
0327     ptPartMin = 0,
0328     dzPartMax = 0.1,
0329     maxAbsIso = 8,
0330     #####
0331     doRelIsolation = False,
0332     MaxIsoSum = 0.1, ### unused
0333     MaxIsoSumEMU = 0.2, ### unused
0334     doSecondVeto = False,
0335     #####
0336     doIsoAnnulus= False,
0337     ###
0338     doPrune = True,
0339     do_mc_match = True, # note: it will in any case try it only on MC, not on data
0340   )
0341 )