analyzers/objects/METAnalyzer.py

0001 from PhysicsTools.Heppy.analyzers.core.Analyzer import Analyzer
0002 from PhysicsTools.Heppy.analyzers.core.AutoHandle import AutoHandle
0003 from PhysicsTools.Heppy.physicsobjects.PhysicsObjects import Jet
0004 from PhysicsTools.HeppyCore.utils.deltar import *
0005 from PhysicsTools.HeppyCore.statistics.counter import Counter, Counters
0006 from PhysicsTools.Heppy.physicsutils.JetReCalibrator import Type1METCorrector, setFakeRawMETOnOldMiniAODs
0007 import PhysicsTools.HeppyCore.framework.config as cfg
0008
0009 import copy
0010 import ROOT
0011 from math import hypot
0012
0013 from copy import deepcopy
0014
0015 class METAnalyzer( Analyzer ):
0016     def __init__(self, cfg_ana, cfg_comp, looperName ):
0017         super(METAnalyzer,self).__init__(cfg_ana,cfg_comp,looperName)
0018         self.recalibrateMET   = cfg_ana.recalibrate
0019         self.applyJetSmearing = cfg_comp.isMC and cfg_ana.applyJetSmearing
0020         self.old74XMiniAODs         = cfg_ana.old74XMiniAODs
0021         self.jetAnalyzerPostFix = getattr(cfg_ana, 'jetAnalyzerPostFix', '')
0022         if self.recalibrateMET in [ "type1", True ]:
0023             if self.recalibrateMET == "type1":
0024                 self.type1METCorrector = Type1METCorrector(self.old74XMiniAODs)
0025         elif self.recalibrateMET != False:
0026             raise RuntimeError("Unsupported value %r for option 'recalibrate': allowed are True, False, 'type1'" % self.recalibrateMET)
0027
0028     def declareHandles(self):
0029         super(METAnalyzer, self).declareHandles()
0030         self.handles['met'] = AutoHandle( self.cfg_ana.metCollection, 'std::vector<pat::MET>' )
0031         if self.cfg_ana.doMetNoPU:
0032             self.handles['nopumet'] = AutoHandle( self.cfg_ana.noPUMetCollection, 'std::vector<pat::MET>' )
0033         if self.cfg_ana.doTkMet:
0034             self.handles['cmgCand'] = AutoHandle( self.cfg_ana.candidates, self.cfg_ana.candidatesTypes )
0035             #self.handles['vertices'] =  AutoHandle( "offlineSlimmedPrimaryVertices", 'std::vector<reco::Vertex>', fallbackLabel="offlinePrimaryVertices" )
0036             self.mchandles['packedGen'] = AutoHandle( 'packedGenParticles', 'std::vector<pat::PackedGenParticle>' )
0037
0038     def beginLoop(self, setup):
0039         super(METAnalyzer,self).beginLoop(setup)
0040         self.counters.addCounter('events')
0041         count = self.counters.counter('events')
0042         count.register('all events')
0043
0044     def applyDeltaMet(self, met, deltaMet):
0045         px,py = self.met.px()+deltaMet[0], self.met.py()+deltaMet[1]
0046         met.setP4(ROOT.reco.Particle.LorentzVector(px,py, 0, hypot(px,py)))
0047
0048     def adduParaPerp(self, met, boson, postfix):
0049
0050         upara = 0
0051         uperp = 0
0052         uX = - met.px() - boson.px()
0053         uY = - met.py() - boson.py()
0054         u1 = (uX*boson.px() + uY*boson.py())/boson.pt()
0055         u2 = (uX*boson.py() - uY*boson.px())/boson.pt()
0056
0057         setattr(met, "upara"+postfix, u1)
0058         setattr(met, "uperp"+postfix, u2)
0059
0060     def makeTkMETs(self, event):
0061         charged = []
0062         chargedchs = []
0063         chargedPVLoose = []
0064         chargedPVTight = []
0065         dochs=getattr(self.cfg_ana,"includeTkMetCHS",True)
0066         dotight=getattr(self.cfg_ana,"includeTkMetPVTight",True)
0067         doloose=getattr(self.cfg_ana,"includeTkMetPVLoose",True)
0068         pfcands = self.handles['cmgCand'].product()
0069
0070         for pfcand in pfcands:
0071
0072 ## ===> require the Track Candidate charge and with a  minimum dz
0073             if (pfcand.charge()!=0):
0074
0075                 pvflag = pfcand.fromPV()
0076                 pxy = pfcand.px(), pfcand.py()
0077
0078                 if abs(pfcand.dz())<=self.cfg_ana.dzMax:
0079                     charged.append(pxy)
0080
0081                 if dochs and  pvflag>0:
0082                     chargedchs.append(pxy)
0083
0084                 if doloose and pvflag>1:
0085                     chargedPVLoose.append(pxy)
0086
0087                 if dotight and pvflag>2:
0088                     chargedPVTight.append(pxy)
0089
0090         def sumXY(pxys):
0091             px, py = sum(x[0] for x in pxys), sum(x[1] for x in pxys)
0092             return ROOT.reco.Particle.LorentzVector(-px, -py, 0, hypot(px,py))
0093         setattr(event, "tkMet"+self.cfg_ana.collectionPostFix, sumXY(charged))
0094         setattr(event, "tkMetPVchs"+self.cfg_ana.collectionPostFix, sumXY(chargedchs))
0095         setattr(event, "tkMetPVLoose"+self.cfg_ana.collectionPostFix, sumXY(chargedPVLoose))
0096         setattr(event, "tkMetPVTight"+self.cfg_ana.collectionPostFix, sumXY(chargedPVTight))
0097         getattr(event,"tkMet"+self.cfg_ana.collectionPostFix).sumEt = sum([hypot(x[0],x[1]) for x in charged])
0098         getattr(event,"tkMetPVchs"+self.cfg_ana.collectionPostFix).sumEt = sum([hypot(x[0],x[1]) for x in chargedchs])
0099         getattr(event,"tkMetPVLoose"+self.cfg_ana.collectionPostFix).sumEt = sum([hypot(x[0],x[1]) for x in chargedPVLoose])
0100         getattr(event,"tkMetPVTight"+self.cfg_ana.collectionPostFix).sumEt = sum([hypot(x[0],x[1]) for x in chargedPVTight])
0101
0102         if  hasattr(event,'zll_p4'):
0103             self.adduParaPerp(getattr(event,"tkMet"+self.cfg_ana.collectionPostFix), event.zll_p4,"_zll")
0104             self.adduParaPerp(getattr(event,"tkMetPVchs"+self.cfg_ana.collectionPostFix), event.zll_p4,"_zll")
0105             self.adduParaPerp(getattr(event,"tkMetPVLoose"+self.cfg_ana.collectionPostFix), event.zll_p4,"_zll")
0106             self.adduParaPerp(getattr(event,"tkMetPVTight"+self.cfg_ana.collectionPostFix), event.zll_p4,"_zll")
0107
0108
0109     def makeGenTkMet(self, event):
0110         genCharged = [ (x.px(),x.py()) for x in self.mchandles['packedGen'].product() if x.charge() != 0 and abs(x.eta()) < 2.4 ]
0111         px, py = sum(x[0] for x in genCharged), sum(x[1] for x in genCharged)
0112         setattr(event,"tkGenMet"+self.cfg_ana.collectionPostFix, ROOT.reco.Particle.LorentzVector(-px , -py, 0, hypot(px,py)))
0113
0114     def makeMETNoMu(self, event):
0115         self.metNoMu = copy.deepcopy(self.met)
0116         if self.cfg_ana.doMetNoPU: self.metNoMuNoPU = copy.deepcopy(self.metNoPU)
0117
0118         mupx = 0
0119         mupy = 0
0120         #sum muon momentum
0121         for mu in event.selectedMuons:
0122             mupx += mu.px()
0123             mupy += mu.py()
0124
0125         #subtract muon momentum and construct met
0126         px,py = self.metNoMu.px()+mupx, self.metNoMu.py()+mupy
0127         self.metNoMu.setP4(ROOT.reco.Particle.LorentzVector(px,py, 0, hypot(px,py)))
0128         px,py = self.metNoMuNoPU.px()+mupx, self.metNoMuNoPU.py()+mupy
0129         self.metNoMuNoPU.setP4(ROOT.reco.Particle.LorentzVector(px,py, 0, hypot(px,py)))
0130         setattr(event, "metNoMu"+self.cfg_ana.collectionPostFix, self.metNoMu)
0131         if self.cfg_ana.doMetNoPU: setattr(event, "metNoMuNoPU"+self.cfg_ana.collectionPostFix, self.metNoMuNoPU)
0132
0133
0134     def makeMETNoEle(self, event):
0135         self.metNoEle = copy.deepcopy(self.met)
0136         if self.cfg_ana.doMetNoPU: self.metNoEleNoPU = copy.deepcopy(self.metNoPU)
0137
0138         elepx = 0
0139         elepy = 0
0140         #sum electron momentum
0141         for ele in event.selectedElectrons:
0142             elepx += ele.px()
0143             elepy += ele.py()
0144
0145         #subtract electron momentum and construct met
0146         px,py = self.metNoEle.px()+elepx, self.metNoEle.py()+elepy
0147         self.metNoEle.setP4(ROOT.reco.Particle.LorentzVector(px,py, 0, hypot(px,py)))
0148
0149         px,py = self.metNoEleNoPU.px()+elepx, self.metNoEleNoPU.py()+elepy
0150         self.metNoEleNoPU.setP4(ROOT.reco.Particle.LorentzVector(px,py, 0, hypot(px,py)))
0151         setattr(event, "metNoEle"+self.cfg_ana.collectionPostFix, self.metNoEle)
0152         if self.cfg_ana.doMetNoPU: setattr(event, "metNoEleNoPU"+self.cfg_ana.collectionPostFix, self.metNoEleNoPU)
0153
0154     def makeMETNoPhoton(self, event):
0155         self.metNoPhoton = copy.deepcopy(self.met)
0156
0157         phopx = 0
0158         phopy = 0
0159         #sum photon momentum
0160         for pho in event.selectedPhotons:
0161             phopx += pho.px()
0162             phopy += pho.py()
0163
0164         #subtract photon momentum and construct met
0165         px,py = self.metNoPhoton.px()+phopx, self.metNoPhoton.py()+phopy
0166         self.metNoPhoton.setP4(ROOT.reco.Particle.LorentzVector(px,py, 0, hypot(px,py)))
0167         setattr(event, "metNoPhoton"+self.cfg_ana.collectionPostFix, self.metNoPhoton)
0168         if self.cfg_ana.doMetNoPU:
0169           self.metNoPhotonNoPU = copy.deepcopy(self.metNoPU)
0170           px,py = self.metNoPhotonNoPU.px()+phopx, self.metNoPhotonNoPU.py()+phopy
0171           self.metNoPhotonNoPU.setP4(ROOT.reco.Particle.LorentzVector(px,py, 0, hypot(px,py)))
0172           setattr(event, "metNoPhotonNoPU"+self.cfg_ana.collectionPostFix, self.metNoPhotonNoPU)
0173
0174
0175     def makeMETs(self, event):
0176         import ROOT
0177         if self.cfg_ana.copyMETsByValue:
0178           self.met = ROOT.pat.MET(self.handles['met'].product()[0])
0179           if self.cfg_ana.doMetNoPU: self.metNoPU = ROOT.pat.MET(self.handles['nopumet'].product()[0])
0180         else:
0181           self.met = self.handles['met'].product()[0]
0182           if self.cfg_ana.doMetNoPU: self.metNoPU = self.handles['nopumet'].product()[0]
0183
0184         if self.recalibrateMET == "type1":
0185           type1METCorr = getattr(event, 'type1METCorr'+self.jetAnalyzerPostFix)
0186           self.type1METCorrector.correct(self.met, type1METCorr)
0187         elif self.recalibrateMET == True:
0188           deltaMetJEC = getattr(event, 'deltaMetFromJEC'+self.jetAnalyzerPostFix)
0189           self.applyDeltaMet(self.met, deltaMetJEC)
0190         if self.applyJetSmearing:
0191           deltaMetSmear = getattr(event, 'deltaMetFromJetSmearing'+self.jetAnalyzerPostFix)
0192           self.applyDeltaMet(self.met, deltaMetSmear)
0193
0194
0195         if (not self.cfg_ana.copyMETsByValue) and getattr(self.cfg_ana, 'makeShiftedMETs', True):
0196           shifts = []
0197           for obj in 'JetEn', 'JetRes', 'MuonEn', 'ElectronEn', 'PhotonEn', 'TauEn', 'UnclusteredEn':
0198             for sh in 'Up', 'Down':
0199                 shifts.append( (obj+sh, getattr(self.met,obj+sh)) )
0200           shifts.append( ('NoShift', self.met.NoShift) )
0201           for name,i in shifts:
0202                key = i
0203                m = ROOT.pat.MET(self.met)
0204                if self.old74XMiniAODs:
0205                     if key > 12:   key = 12
0206                     elif key <= 3: key = { 'JetEnUp':0, 'JetEnDown':1, 'JetResUp':2, 'JetResDown':3 }[name]
0207                m.setP4(self.met.shiftedP4(key))
0208                setattr(event, "met{0}_shifted_{1}".format(self.cfg_ana.collectionPostFix, i),m)
0209                setattr(event, "met{0}_shifted_{1}".format(self.cfg_ana.collectionPostFix, name),m)
0210
0211         self.met_sig = self.met.significance()
0212         self.met_sumet = self.met.sumEt()
0213
0214         if self.old74XMiniAODs and self.recalibrateMET != "type1":
0215            oldraw = self.met.shiftedP2_74x(12,0);
0216            setFakeRawMETOnOldMiniAODs( self.met, oldraw.px, oldraw.py, self.met.shiftedSumEt_74x(12,0) )
0217            px, py = oldraw.px, oldraw.py
0218         else:
0219            px, py = self.met.uncorPx(), self.met.uncorPy()
0220         self.met_raw = ROOT.reco.Particle.LorentzVector(px,py,0,hypot(px,py))
0221
0222         if hasattr(event,'zll_p4'):
0223             self.adduParaPerp(self.met,event.zll_p4,"_zll")
0224             self.adduParaPerp(self.met_raw, event.zll_p4,"_zll")
0225             setattr(event,"met_raw"+self.cfg_ana.collectionPostFix, self.met_raw)
0226             setattr(event,"met_raw.upara_zll"+self.cfg_ana.collectionPostFix, self.met_raw.upara_zll)
0227             setattr(event,"met_raw.uperp_zll"+self.cfg_ana.collectionPostFix, self.met_raw.uperp_zll)
0228
0229         if hasattr(event,'gamma_p4'):
0230             self.adduParaPerp(self.met,event.gamma_p4,"_gamma")
0231             self.adduParaPerp(self.met_raw, event.gamma_p4,"_gamma")
0232             setattr(event,"met_raw"+self.cfg_ana.collectionPostFix, self.met_raw)
0233             setattr(event,"met_raw.upara_gamma"+self.cfg_ana.collectionPostFix, self.met_raw.upara_gamma)
0234             setattr(event,"met_raw.uperp_gamma"+self.cfg_ana.collectionPostFix, self.met_raw.uperp_gamma)
0235
0236         if hasattr(event,"met"+self.cfg_ana.collectionPostFix): raise RuntimeError("Event already contains met with the following postfix: "+self.cfg_ana.collectionPostFix)
0237         setattr(event, "met"+self.cfg_ana.collectionPostFix, self.met)
0238         if self.cfg_ana.doMetNoPU: setattr(event, "metNoPU"+self.cfg_ana.collectionPostFix, self.metNoPU)
0239         setattr(event, "met_sig"+self.cfg_ana.collectionPostFix, self.met_sig)
0240         setattr(event, "met_sumet"+self.cfg_ana.collectionPostFix, self.met_sumet)
0241
0242         genMET = self.met.genMET()
0243         if genMET:
0244           setattr(event, "met_genPt"+self.cfg_ana.collectionPostFix, genMET.pt())
0245           setattr(event, "met_genPhi"+self.cfg_ana.collectionPostFix, genMET.phi())
0246         else:
0247           setattr(event, "met_genPt"+self.cfg_ana.collectionPostFix, float('nan'))
0248           setattr(event, "met_genPhi"+self.cfg_ana.collectionPostFix, float('nan'))
0249
0250         if self.cfg_ana.doMetNoMu and hasattr(event, 'selectedMuons'):
0251             self.makeMETNoMu(event)
0252
0253         if self.cfg_ana.doMetNoEle and hasattr(event, 'selectedElectrons'):
0254             self.makeMETNoEle(event)
0255
0256         if self.cfg_ana.doMetNoPhoton and hasattr(event, 'selectedPhotons'):
0257             self.makeMETNoPhoton(event)
0258
0259     def process(self, event):
0260         self.readCollections( event.input)
0261         self.counters.counter('events').inc('all events')
0262
0263         self.makeMETs(event)
0264
0265         if self.cfg_ana.doTkMet:
0266             self.makeTkMETs(event);
0267
0268         if getattr(self.cfg_ana,"doTkGenMet",self.cfg_ana.doTkMet) and self.cfg_comp.isMC and hasattr(event, 'genParticles'):
0269             self.makeGenTkMet(event)
0270
0271         return True
0272
0273
0274 setattr(METAnalyzer,"defaultConfig", cfg.Analyzer(
0275     class_object = METAnalyzer,
0276     metCollection     = "slimmedMETs",
0277     noPUMetCollection = "slimmedMETs",
0278     copyMETsByValue = False,
0279     recalibrate = True,
0280     applyJetSmearing = True,
0281     jetAnalyzerPostFix = "",
0282     old74XMiniAODs = False, # need to set to True to get proper Raw MET on plain 74X MC produced with CMSSW <= 7_4_12
0283     makeShiftedMETs = True,
0284     doTkMet = False,
0285     includeTkMetCHS = True,
0286     includeTkMetPVLoose = True,
0287     includeTkMetPVTight = True,
0288     doMetNoPU = False, # Not existing in MiniAOD at the moment
0289     doMetNoMu = False,
0290     doMetNoEle = False,
0291     doMetNoPhoton = False,
0292     candidates='packedPFCandidates',
0293     candidatesTypes='std::vector<pat::PackedCandidate>',
0294     dzMax = 0.1,
0295     collectionPostFix = "",
0296     )
0297 )