analyzers/objects/PhotonAnalyzer.py

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