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 import re
0007 
0008 from ROOT import TLorentzVector
0009 from ROOT import heppy 
0010 
0011 from PhysicsTools.Heppy.analyzers.core.Analyzer import Analyzer
0012 from PhysicsTools.HeppyCore.framework.event import Event
0013 from PhysicsTools.HeppyCore.statistics.counter import Counter, Counters
0014 from PhysicsTools.Heppy.analyzers.core.AutoHandle import AutoHandle
0015 from PhysicsTools.Heppy.physicsobjects.Photon import Photon
0016 from PhysicsTools.Heppy.physicsutils.PhotonCalibrator import Run2PhotonCalibrator
0017 
0018 from PhysicsTools.HeppyCore.utils.deltar import deltaR, deltaPhi, bestMatch, matchObjectCollection3
0019 
0020 import PhysicsTools.HeppyCore.framework.config as cfg
0021 
0022 
0023 class PhotonAnalyzer( Analyzer ):
0024 
0025     
0026     def __init__(self, cfg_ana, cfg_comp, looperName ):
0027         super(PhotonAnalyzer,self).__init__(cfg_ana,cfg_comp,looperName)
0028         self.etaCentral = self.cfg_ana.etaCentral  if hasattr(self.cfg_ana, 'etaCentral') else 9999
0029         # footprint removed isolation
0030         self.doFootprintRemovedIsolation = getattr(cfg_ana, 'doFootprintRemovedIsolation', False)
0031         if self.doFootprintRemovedIsolation:
0032             self.footprintRemovedIsolationPUCorr =  self.cfg_ana.footprintRemovedIsolationPUCorr
0033             self.IsolationComputer = heppy.IsolationComputer()
0034     #FIXME: only Embedded works
0035         if self.cfg_ana.doPhotonScaleCorrections:
0036             conf = cfg_ana.doPhotonScaleCorrections
0037             self.photonEnergyCalibrator = Run2PhotonCalibrator(
0038                 conf['data'],
0039                 cfg_comp.isMC,
0040                 conf['isSync'] if 'isSync' in conf else False,
0041             )
0042 
0043     def declareHandles(self):
0044         super(PhotonAnalyzer, self).declareHandles()
0045         self.handles['rhoPhoton'] = AutoHandle( self.cfg_ana.rhoPhoton, 'double')
0046 
0047     #----------------------------------------                                                                                                                                   
0048     # DECLARATION OF HANDLES OF PHOTONS STUFF                                                                                                                                   
0049     #----------------------------------------     
0050 
0051         self.handles['photons'] = AutoHandle( self.cfg_ana.photons,'std::vector<pat::Photon>')
0052         self.mchandles['packedGen'] = AutoHandle( 'packedGenParticles', 'std::vector<pat::PackedGenParticle>' )
0053         self.mchandles['prunedGen'] = AutoHandle( 'prunedGenParticles', 'std::vector<reco::GenParticle>' )
0054 
0055         self.handles['packedCandidates'] = AutoHandle( 'packedPFCandidates', 'std::vector<pat::PackedCandidate>')
0056         self.handles['jets'] = AutoHandle( "slimmedJets", 'std::vector<pat::Jet>' )
0057 
0058 
0059     def beginLoop(self, setup):
0060         super(PhotonAnalyzer,self).beginLoop(setup)
0061         self.counters.addCounter('events')
0062         count = self.counters.counter('events')
0063         count.register('all events')
0064         count.register('has >=1 gamma at preselection')
0065         count.register('has >=1 selected gamma')
0066 
0067     def makePhotons(self, event):
0068         event.allphotons = map( Photon, self.handles['photons'].product() )
0069         event.allphotons.sort(key = lambda l : l.pt(), reverse = True)
0070 
0071         event.selectedPhotons = []
0072         event.selectedPhotonsCentral = []
0073 
0074         if self.doFootprintRemovedIsolation:
0075             # values are taken from EGamma implementation: https://github.com/cms-sw/cmssw/blob/CMSSW_7_6_X/RecoEgamma/PhotonIdentification/plugins/PhotonIDValueMapProducer.cc#L198-L199
0076             self.IsolationComputer.setPackedCandidates(self.handles['packedCandidates'].product(), -1, 0.1, 0.2)
0077 
0078         # Photon scale calibrations
0079         if self.cfg_ana.doPhotonScaleCorrections:
0080             for gamma in event.allphotons:
0081                 self.photonEnergyCalibrator.correct(gamma, event.run)
0082 
0083         foundPhoton = False
0084         for gamma in event.allphotons:
0085             if gamma.pt() < self.cfg_ana.ptMin: continue
0086             if abs(gamma.eta()) > self.cfg_ana.etaMax: continue
0087             foundPhoton = True
0088 
0089             gamma.rho = float(self.handles['rhoPhoton'].product()[0])
0090             # https://twiki.cern.ch/twiki/bin/view/CMS/CutBasedPhotonIdentificationRun2#Selection_implementation_details
0091             if   abs(gamma.eta()) < 1.0:   gamma.EffectiveArea03 = [ 0.0234, 0.0053, 0.078  ]
0092             elif abs(gamma.eta()) < 1.479: gamma.EffectiveArea03 = [ 0.0189, 0.0103, 0.0629 ]
0093             elif abs(gamma.eta()) < 2.0:   gamma.EffectiveArea03 = [ 0.0171, 0.0057, 0.0264 ]
0094             elif abs(gamma.eta()) < 2.2:   gamma.EffectiveArea03 = [ 0.0129, 0.0070, 0.0462 ]
0095             elif abs(gamma.eta()) < 2.3:   gamma.EffectiveArea03 = [ 0.0110, 0.0152, 0.0740 ]
0096             elif abs(gamma.eta()) < 2.4:   gamma.EffectiveArea03 = [ 0.0074, 0.0232, 0.0924 ]
0097             else:                          gamma.EffectiveArea03 = [ 0.0035, 0.1709, 0.1484 ]
0098 
0099             if self.doFootprintRemovedIsolation:
0100                 self.attachFootprintRemovedIsolation(gamma)
0101 
0102             gamma.relIso = (max(gamma.chargedHadronIso()-gamma.rho*gamma.EffectiveArea03[0],0) + max(gamma.neutralHadronIso()-gamma.rho*gamma.EffectiveArea03[1],0) + max(gamma.photonIso() - gamma.rho*gamma.EffectiveArea03[2],0))/gamma.pt()
0103 
0104             def idWP(gamma,X):
0105                 """Create an integer equal to 1-2-3 for (loose,medium,tight)"""
0106 
0107                 id=0
0108                 if gamma.photonID(X%"Loose"):
0109                     id=1
0110                 #if gamma.photonID(X%"Medium"):
0111                 #    id=2 
0112                 if gamma.photonID(X%"Tight"):
0113                     id=3
0114                 return id
0115 
0116             gamma.idCutBased = idWP(gamma, "PhotonCutBasedID%s")
0117 
0118 
0119             keepThisPhoton = True
0120 
0121             if self.cfg_ana.gammaID=="PhotonCutBasedIDLoose_CSA14" or self.cfg_ana.gammaID=="PhotonCutBasedIDLoose_PHYS14" :
0122                 gamma.idCutBased = gamma.photonIDCSA14(self.cfg_ana.gammaID) 
0123                 # we're keeing sigmaietaieta sidebands:
0124                 keepThisPhoton   = gamma.photonIDCSA14(self.cfg_ana.gammaID, True) 
0125 
0126                 if gamma.hasPixelSeed():
0127                     keepThisPhoton = False
0128                     gamma.idCutBased = 0
0129             elif "NoIso" in self.cfg_ana.gammaID:
0130                 idName = re.split('_NoIso',self.cfg_ana.gammaID)
0131                 keepThisPhoton = gamma.passPhotonID(idName[0],self.cfg_ana.conversionSafe_eleVeto)
0132                 basenameID = re.split('_looseSieie',idName[0])
0133                 gamma.idCutBased = gamma.passPhotonID(basenameID[0],self.cfg_ana.conversionSafe_eleVeto)
0134             else:
0135                 # Reading from miniAOD directly
0136                 #keepThisPhoton = gamma.photonID(self.cfg_ana.gammaID)
0137 
0138                 # implement cut based ID with CMGTools
0139                 keepThisPhoton = gamma.passPhotonID(self.cfg_ana.gammaID,self.cfg_ana.conversionSafe_eleVeto) and gamma.passPhotonIso(self.cfg_ana.gammaID,self.cfg_ana.gamma_isoCorr)
0140 
0141             if keepThisPhoton:
0142                 event.selectedPhotons.append(gamma)
0143 
0144             if keepThisPhoton and abs(gamma.eta()) < self.etaCentral:
0145                 event.selectedPhotonsCentral.append(gamma)
0146 
0147         event.selectedPhotons.sort(key = lambda l : l.pt(), reverse = True)
0148         event.selectedPhotonsCentral.sort(key = lambda l : l.pt(), reverse = True)
0149 
0150         self.counters.counter('events').inc('all events')
0151         if foundPhoton: self.counters.counter('events').inc('has >=1 gamma at preselection')
0152         if len(event.selectedPhotons): self.counters.counter('events').inc('has >=1 selected gamma')
0153        
0154     def matchPhotons(self, event):
0155         event.genPhotons = [ x for x in event.genParticles if x.status() == 1 and abs(x.pdgId()) == 22 ]
0156         event.genPhotonsWithMom = [ x for x in event.genPhotons if x.numberOfMothers()>0 ]
0157         event.genPhotonsWithoutMom = [ x for x in event.genPhotons if x.numberOfMothers()==0 ]
0158         event.genPhotonsMatched = [ x for x in event.genPhotonsWithMom if abs(x.mother(0).pdgId())<23 or x.mother(0).pdgId()==2212 ]
0159         match = matchObjectCollection3(event.allphotons, event.genPhotonsMatched, deltaRMax = 0.1)
0160         matchNoMom = matchObjectCollection3(event.allphotons, event.genPhotonsWithoutMom, deltaRMax = 0.1)
0161         packedGenParts = [ p for p in self.mchandles['packedGen'].product() if abs(p.eta()) < 3.1 ]
0162         partons = [ p for p in self.mchandles['prunedGen'].product() if (p.status()==23 or p.status()==22) and abs(p.pdgId())<22 ]
0163         for gamma in event.allphotons:
0164           gen = match[gamma]
0165           gamma.mcGamma = gen
0166           if gen and gen.pt()>=0.5*gamma.pt() and gen.pt()<=2.*gamma.pt():
0167             gamma.mcMatchId = 22
0168             sumPt03 = 0.;
0169             sumPt04 = 0.;
0170             for part in packedGenParts:
0171               if abs(part.pdgId())==12: continue # exclude neutrinos
0172               if abs(part.pdgId())==14: continue
0173               if abs(part.pdgId())==16: continue
0174               if abs(part.pdgId())==18: continue
0175               deltar = deltaR(gen.eta(), gen.phi(), part.eta(), part.phi())
0176               if deltar <= 0.3:
0177                 sumPt03 += part.pt()
0178               if deltar <= 0.4:
0179                 sumPt04 += part.pt()
0180             sumPt03 -= gen.pt()
0181             sumPt04 -= gen.pt()
0182             if sumPt03<0. : sumPt03=0.
0183             if sumPt04<0. : sumPt04=0.
0184             gamma.genIso03 = sumPt03
0185             gamma.genIso04 = sumPt04
0186             # match to parton
0187             deltaRmin = 999.
0188             for p in partons:
0189               deltar = deltaR(gen.eta(), gen.phi(), p.eta(), p.phi())
0190               if deltar < deltaRmin:
0191                 deltaRmin = deltar
0192             gamma.drMinParton = deltaRmin
0193           else:
0194             genNoMom = matchNoMom[gamma]
0195             if genNoMom:
0196               gamma.mcMatchId = 7
0197               sumPt03 = 0.;
0198               sumPt04 = 0.;
0199               for part in packedGenParts:
0200                 if abs(part.pdgId())==12: continue # exclude neutrinos
0201                 if abs(part.pdgId())==14: continue
0202                 if abs(part.pdgId())==16: continue
0203                 if abs(part.pdgId())==18: continue
0204                 deltar = deltaR(genNoMom.eta(), genNoMom.phi(), part.eta(), part.phi())
0205                 if deltar <= 0.3:
0206                   sumPt03 += part.pt()
0207                 if deltar <= 0.4:
0208                   sumPt04 += part.pt()
0209               sumPt03 -= genNoMom.pt()
0210               sumPt04 -= genNoMom.pt()
0211               if sumPt03<0. : sumPt03=0.
0212               if sumPt04<0. : sumPt04=0.
0213               gamma.genIso03 = sumPt03
0214               gamma.genIso04 = sumPt04
0215               # match to parton
0216               deltaRmin = 999.
0217               for p in partons:
0218                 deltar = deltaR(genNoMom.eta(), genNoMom.phi(), p.eta(), p.phi())
0219                 if deltar < deltaRmin:
0220                   deltaRmin = deltar
0221               gamma.drMinParton = deltaRmin
0222             else:
0223               gamma.mcMatchId = 0
0224               gamma.genIso03 = -1.
0225               gamma.genIso04 = -1.
0226               gamma.drMinParton = -1.
0227 
0228 
0229 
0230 
0231 
0232     def checkMatch( self, eta, phi, particles, deltar ):
0233 
0234       for part in particles:
0235         if deltaR(eta, phi, part.eta(), part.phi()) < deltar:
0236           return True
0237 
0238       return False
0239 
0240 
0241 
0242 
0243 
0244     def computeRandomCone( self, event, eta, phi, deltarmax, charged, jets, photons ):
0245 
0246       if self.checkMatch( eta, phi, jets, 2.*deltarmax ): 
0247         return -1.
0248     
0249       if self.checkMatch( eta, phi, photons, 2.*deltarmax ): 
0250         return -1.
0251     
0252       if self.checkMatch( eta, phi, event.selectedLeptons, deltarmax ): 
0253         return -1.
0254 
0255       iso = 0.
0256 
0257       for part in charged:
0258         if deltaR(eta, phi, part.eta(), part.phi()) > deltarmax : continue
0259         #if deltaR(eta, phi, part.eta(), part.phi()) < 0.02: continue
0260         iso += part.pt()
0261 
0262       return iso
0263 
0264 
0265 
0266 
0267             
0268 
0269     def randomCone( self, event ):
0270 
0271         patcands  = self.handles['packedCandidates'].product()
0272         jets  = self.handles['jets'].product()
0273 
0274         charged   = [ p for p in patcands if ( p.charge() != 0 and abs(p.pdgId())>20 and abs(p.dz())<=0.1 and p.fromPV()>1 and p.trackHighPurity() ) ]
0275         photons10 = [ p for p in patcands if ( p.pdgId() == 22 and p.pt()>10. ) ]
0276         jets20 = [ j for j in jets if j.pt() > 20 and abs(j.eta())<2.5 ]
0277 
0278         for gamma in event.allphotons:
0279 
0280           etaPhot = gamma.eta()
0281           phiPhot = gamma.eta()
0282           pi = 3.14159
0283           phiRC = phiPhot + 0.5*pi
0284           while phiRC>pi:
0285             phiRC -= 2.*pi
0286 
0287 
0288           gamma.chHadIsoRC03 = self.computeRandomCone( event, etaPhot, phiRC, 0.3, charged, jets20, photons10 )
0289           gamma.chHadIsoRC04 = self.computeRandomCone( event, etaPhot, phiRC, 0.4, charged, jets20, photons10 )
0290           
0291           
0292           #try other side
0293           phiRC = phiPhot - 0.5*pi
0294           while phiRC<-pi:
0295             phiRC += 2.*pi
0296           
0297           if gamma.chHadIsoRC03<0. : gamma.chHadIsoRC03 = self.computeRandomCone( event, etaPhot, phiRC, 0.3, charged, jets20, photons10 )
0298           if gamma.chHadIsoRC04<0. : gamma.chHadIsoRC04 = self.computeRandomCone( event, etaPhot, phiRC, 0.4, charged, jets20, photons10 )
0299 
0300 
0301     def attachFootprintRemovedIsolation(self, gamma):
0302         # cone deltar=0.3
0303         gamma.ftprAbsIsoCharged03 = self.IsolationComputer.chargedAbsIso(gamma.physObj, 0.3, 0, 0.0);
0304         gamma.ftprAbsIsoPho03     = self.IsolationComputer.photonAbsIsoRaw( gamma.physObj, 0.3, 0, 0.0);
0305         gamma.ftprAbsIsoNHad03    = self.IsolationComputer.neutralHadAbsIsoRaw( gamma.physObj, 0.3, 0, 0.0);
0306         gamma.ftprAbsIso03 = gamma.ftprAbsIsoCharged03 + gamma.ftprAbsIsoPho03 + gamma.ftprAbsIsoNHad03
0307         if self.cfg_ana.gamma_isoCorr == "rhoArea":
0308             gamma.ftprAbsIso03 = (max(gamma.ftprAbsIsoCharged03-gamma.rho*gamma.EffectiveArea03[0],0) + max(gamma.ftprAbsIsoPho03-gamma.rho*gamma.EffectiveArea03[1],0) + max(gamma.ftprAbsIsoNHad03 - gamma.rho*gamma.EffectiveArea03[2],0))
0309         elif self.cfg_ana.gamma_isoCorr != 'raw':
0310             raise RuntimeError("Unsupported gamma_isoCorr name '" + str(self.cfg_ana.gamma_isoCorr) +  "'! For now only 'rhoArea', 'raw' are supported.")
0311         gamma.ftprRelIso03 = gamma.ftprAbsIso03/gamma.pt()
0312 
0313     def printInfo(self, event):
0314         print('----------------')
0315         if len(event.selectedPhotons)>0:
0316             print('lenght: ',len(event.selectedPhotons))
0317             print('gamma candidate pt: ',event.selectedPhotons[0].pt())
0318             print('gamma candidate eta: ',event.selectedPhotons[0].eta())
0319             print('gamma candidate phi: ',event.selectedPhotons[0].phi())
0320             print('gamma candidate mass: ',event.selectedPhotons[0].mass())
0321             print('gamma candidate HoE: ',event.selectedPhotons[0].hOVERe())
0322             print('gamma candidate r9: ',event.selectedPhotons[0].full5x5_r9())
0323             print('gamma candidate sigmaIetaIeta: ',event.selectedPhotons[0].full5x5_sigmaIetaIeta())
0324             print('gamma candidate had iso: ',event.selectedPhotons[0].chargedHadronIso())
0325             print('gamma candidate neu iso: ',event.selectedPhotons[0].neutralHadronIso())
0326             print('gamma candidate gamma iso: ',event.selectedPhotons[0].photonIso())
0327             print('gamma idCutBased',event.selectedPhotons[0].idCutBased)
0328 
0329 
0330     def process(self, event):
0331         self.readCollections( event.input )
0332         self.makePhotons(event)
0333 #        self.printInfo(event)   
0334 
0335         if self.cfg_ana.do_randomCone:
0336             self.randomCone(event)
0337 
0338         if not self.cfg_comp.isMC:
0339             return True
0340 
0341         if self.cfg_ana.do_mc_match and hasattr(event, 'genParticles'):
0342             self.matchPhotons(event)
0343 
0344 
0345         return True
0346 
0347 
0348 setattr(PhotonAnalyzer,"defaultConfig",cfg.Analyzer(
0349     class_object=PhotonAnalyzer,
0350     photons='slimmedPhotons',
0351     ptMin = 20,
0352     etaMax = 2.5,
0353     # energy scale corrections (off by default)
0354     doPhotonScaleCorrections=False, 
0355     gammaID = "PhotonCutBasedIDLoose_CSA14",
0356     rhoPhoton = 'fixedGridRhoFastjetAll',
0357     gamma_isoCorr = 'rhoArea',
0358     # Footprint-removed isolation, removing all the footprint of the photon
0359     doFootprintRemovedIsolation = False, # off by default since it requires access to all PFCandidates
0360     packedCandidates = 'packedPFCandidates',
0361     footprintRemovedIsolationPUCorr = 'rhoArea', # Allowed options: 'rhoArea', 'raw' (uncorrected)
0362     conversionSafe_eleVeto = False,
0363     do_mc_match = True,
0364     do_randomCone = False,
0365   )
0366 )
0367