Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-11-25 02:29:49

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