Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:32:28

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 )