analyzers/objects/PhotonAnalyzer.py

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