analyzers/eventtopology/RazorAnalyzer.py

0001 import operator
0002 import itertools
0003 import copy
0004 from math import *
0005
0006 from ROOT import std
0007 from ROOT import TLorentzVector, TVector3, TVectorD
0008
0009 from PhysicsTools.Heppy.analyzers.core.Analyzer import Analyzer
0010 from PhysicsTools.HeppyCore.framework.event import Event
0011 from PhysicsTools.HeppyCore.statistics.counter import Counter, Counters
0012 from PhysicsTools.Heppy.analyzers.core.AutoHandle import AutoHandle
0013
0014 from PhysicsTools.HeppyCore.utils.deltar import deltaR
0015
0016 from ROOT.heppy import Megajet
0017 from ROOT.heppy import ReclusterJets
0018
0019 import ROOT
0020
0021 import os
0022
0023 class RazorAnalyzer( Analyzer ):
0024     def __init__(self, cfg_ana, cfg_comp, looperName ):
0025         super(RazorAnalyzer,self).__init__(cfg_ana,cfg_comp,looperName)
0026
0027     def declareHandles(self):
0028         super(RazorAnalyzer, self).declareHandles()
0029        #genJets
0030         self.handles['genJets'] = AutoHandle( 'slimmedGenJets','std::vector<reco::GenJet>')
0031
0032     def beginLoop(self, setup):
0033         super(RazorAnalyzer,self).beginLoop(setup)
0034         self.counters.addCounter('pairs')
0035         count = self.counters.counter('pairs')
0036         count.register('all events')
0037
0038     def computeMR(self, ja, jb):
0039         A = ja.P();
0040         B = jb.P();
0041         az = ja.Pz();
0042         bz = jb.Pz();
0043         mr = sqrt((A+B)*(A+B)-(az+bz)*(az+bz));
0044         return mr
0045
0046     def computeMTR(self, ja, jb, met):
0047
0048         mtr = met.Vect().Mag()*(ja.Pt()+jb.Pt()) - met.Vect().Dot(ja.Vect()+jb.Vect());
0049         mtr = sqrt(mtr/2.);
0050         return mtr;
0051
0052     def computeR(self, ja, jb, met):
0053
0054         mr = self.computeMR(ja,jb)
0055         mtr = self.computeMTR(ja,jb,met)
0056         r = 999999. if mr <= 0 else mtr/mr
0057         return r
0058
0059
0060     def makeRAZOR(self, event):
0061 #        print '==> INSIDE THE PRINT MT2'
0062 #        print 'MET=',event.met.pt()
0063
0064         import array
0065         import numpy
0066
0067
0068 ## ===> hadronic RAZOR
0069
0070         (met, metphi)  = event.met.pt(), event.met.phi()
0071         metp4 = ROOT.TLorentzVector()
0072         metp4.SetPtEtaPhiM(met,0,metphi,0)
0073
0074         objects40jc = [ j for j in event.cleanJets if j.pt() > 40 and abs(j.eta())<2.5 ]
0075
0076 #### get megajets (association method: default 1 = minimum sum of the invariant masses of the two megajets)
0077
0078         if len(objects40jc)>=2:
0079
0080             pxvec  = ROOT.std.vector(float)()
0081             pyvec  = ROOT.std.vector(float)()
0082             pzvec  = ROOT.std.vector(float)()
0083             Evec  = ROOT.std.vector(float)()
0084             grouping  = ROOT.std.vector(int)()
0085
0086             for jet in objects40jc:
0087                 pxvec.push_back(jet.px())
0088                 pyvec.push_back(jet.py())
0089                 pzvec.push_back(jet.pz())
0090                 Evec.push_back(jet.energy())
0091
0092             megajet = Megajet(pxvec, pyvec, pzvec, Evec, 1)
0093
0094             pseudoJet1px = megajet.getAxis1()[0] * megajet.getAxis1()[3]
0095             pseudoJet1py = megajet.getAxis1()[1] * megajet.getAxis1()[3]
0096             pseudoJet1pz = megajet.getAxis1()[2] * megajet.getAxis1()[3]
0097             pseudoJet1energy = megajet.getAxis1()[4]
0098
0099             pseudoJet2px = megajet.getAxis2()[0] * megajet.getAxis2()[3]
0100             pseudoJet2py = megajet.getAxis2()[1] * megajet.getAxis2()[3]
0101             pseudoJet2pz = megajet.getAxis2()[2] * megajet.getAxis2()[3]
0102             pseudoJet2energy = megajet.getAxis2()[4]
0103
0104             pseudoJet1pt2 = pseudoJet1px*pseudoJet1px + pseudoJet1py*pseudoJet1py
0105             pseudoJet2pt2 = pseudoJet2px*pseudoJet2px + pseudoJet2py*pseudoJet2py
0106
0107             if pseudoJet1pt2 >= pseudoJet2pt2:
0108                 event.pseudoJet1_had = ROOT.TLorentzVector( pseudoJet1px, pseudoJet1py, pseudoJet1pz, pseudoJet1energy)
0109                 event.pseudoJet2_had = ROOT.TLorentzVector( pseudoJet2px, pseudoJet2py, pseudoJet2pz, pseudoJet2energy)
0110             else:
0111                 event.pseudoJet2_had = ROOT.TLorentzVector( pseudoJet1px, pseudoJet1py, pseudoJet1pz, pseudoJet1energy)
0112                 event.pseudoJet1_had = ROOT.TLorentzVector( pseudoJet2px, pseudoJet2py, pseudoJet2pz, pseudoJet2energy)
0113
0114             event.mr_had = self.computeMR(event.pseudoJet1_had, event.pseudoJet2_had)
0115             event.mtr_had = self.computeMTR(event.pseudoJet1_had, event.pseudoJet2_had, metp4)
0116             event.r_had = self.computeR(event.pseudoJet1_had, event.pseudoJet2_had, metp4)
0117
0118 #### do same things for GEN
0119
0120         if self.cfg_comp.isMC:
0121             (genmet, genmetphi)  = event.met.genMET().pt(), event.met.genMET().phi()
0122             genmetp4 = ROOT.TLorentzVector()
0123             genmetp4.SetPtEtaPhiM(genmet,0,genmetphi,0)
0124
0125             allGenJets = [ x for x in self.handles['genJets'].product() ]
0126             objects40jc_Gen = [ j for j in allGenJets if j.pt() > 40 and abs(j.eta())<2.5 ]
0127
0128             if len(objects40jc_Gen)>=2:
0129
0130                 pxvec  = ROOT.std.vector(float)()
0131                 pyvec  = ROOT.std.vector(float)()
0132                 pzvec  = ROOT.std.vector(float)()
0133                 Evec  = ROOT.std.vector(float)()
0134                 grouping  = ROOT.std.vector(int)()
0135
0136                 for jet in objects40jc_Gen:
0137                     pxvec.push_back(jet.px())
0138                     pyvec.push_back(jet.py())
0139                     pzvec.push_back(jet.pz())
0140                     Evec.push_back(jet.energy())
0141
0142                 megajet = Megajet(pxvec, pyvec, pzvec, Evec, 1)
0143
0144                 pseudoJet1px = megajet.getAxis1()[0] * megajet.getAxis1()[3]
0145                 pseudoJet1py = megajet.getAxis1()[1] * megajet.getAxis1()[3]
0146                 pseudoJet1pz = megajet.getAxis1()[2] * megajet.getAxis1()[3]
0147                 pseudoJet1energy = megajet.getAxis1()[4]
0148
0149                 pseudoJet2px = megajet.getAxis2()[0] * megajet.getAxis2()[3]
0150                 pseudoJet2py = megajet.getAxis2()[1] * megajet.getAxis2()[3]
0151                 pseudoJet2pz = megajet.getAxis2()[2] * megajet.getAxis2()[3]
0152                 pseudoJet2energy = megajet.getAxis2()[4]
0153
0154                 pseudoJet1pt2 = pseudoJet1px*pseudoJet1px + pseudoJet1py*pseudoJet1py
0155                 pseudoJet2pt2 = pseudoJet2px*pseudoJet2px + pseudoJet2py*pseudoJet2py
0156
0157                 if pseudoJet1pt2 >= pseudoJet2pt2:
0158                     pseudoJet1_gen = ROOT.TLorentzVector( pseudoJet1px, pseudoJet1py, pseudoJet1pz, pseudoJet1energy)
0159                     pseudoJet2_gen = ROOT.TLorentzVector( pseudoJet2px, pseudoJet2py, pseudoJet2pz, pseudoJet2energy)
0160                 else:
0161                     pseudoJet2_gen = ROOT.TLorentzVector( pseudoJet1px, pseudoJet1py, pseudoJet1pz, pseudoJet1energy)
0162                     pseudoJet1_gen = ROOT.TLorentzVector( pseudoJet2px, pseudoJet2py, pseudoJet2pz, pseudoJet2energy)
0163
0164                 event.mr_gen = self.computeMR(pseudoJet1_gen, pseudoJet2_gen)
0165                 event.mtr_gen = self.computeMTR(pseudoJet1_gen, pseudoJet2_gen, genmetp4)
0166                 event.r_gen = self.computeR(pseudoJet1_gen, pseudoJet2_gen, genmetp4)
0167         else:
0168             event.mr_gen = -999
0169             event.mtr_gen = -999
0170             event.r_gen = -999
0171
0172
0173 ## ===> full RAZOR (jets + leptons)
0174         objects10lc = [ l for l in event.selectedLeptons if l.pt() > 10 and abs(l.eta())<2.5 ]
0175         if hasattr(event, 'selectedIsoCleanTrack'):
0176             objects10lc = [ l for l in event.selectedLeptons if l.pt() > 10 and abs(l.eta())<2.5 ] + [ t for t in event.selectedIsoCleanTrack ]
0177
0178         objects40j10lc = objects40jc + objects10lc
0179
0180         objects40j10lc.sort(key = lambda obj : obj.pt(), reverse = True)
0181
0182         if len(objects40j10lc)>=2:
0183
0184             pxvec  = ROOT.std.vector(float)()
0185             pyvec  = ROOT.std.vector(float)()
0186             pzvec  = ROOT.std.vector(float)()
0187             Evec  = ROOT.std.vector(float)()
0188             grouping  = ROOT.std.vector(int)()
0189
0190             for obj in objects40j10lc:
0191                 pxvec.push_back(obj.px())
0192                 pyvec.push_back(obj.py())
0193                 pzvec.push_back(obj.pz())
0194                 Evec.push_back(obj.energy())
0195
0196             #for obj in objects_fullmt2:
0197             #    print "pt: ", obj.pt(), ", eta: ", obj.eta(), ", phi: ", obj.phi(), ", mass: ", obj.mass()
0198
0199             #### get megajets  (association method: default 1 = minimum sum of the invariant masses of the two megajets)
0200
0201             megajet = Megajet(pxvec, pyvec, pzvec, Evec, 1)
0202
0203             pseudoJet1px = megajet.getAxis1()[0] * megajet.getAxis1()[3]
0204             pseudoJet1py = megajet.getAxis1()[1] * megajet.getAxis1()[3]
0205             pseudoJet1pz = megajet.getAxis1()[2] * megajet.getAxis1()[3]
0206             pseudoJet1energy = megajet.getAxis1()[4]
0207
0208             pseudoJet2px = megajet.getAxis2()[0] * megajet.getAxis2()[3]
0209             pseudoJet2py = megajet.getAxis2()[1] * megajet.getAxis2()[3]
0210             pseudoJet2pz = megajet.getAxis2()[2] * megajet.getAxis2()[3]
0211             pseudoJet2energy = megajet.getAxis2()[4]
0212
0213             pseudoJet1pt2 = pseudoJet1px*pseudoJet1px + pseudoJet1py*pseudoJet1py
0214             pseudoJet2pt2 = pseudoJet2px*pseudoJet2px + pseudoJet2py*pseudoJet2py
0215
0216             if pseudoJet1pt2 >= pseudoJet2pt2:
0217                 event.pseudoJet1 = ROOT.TLorentzVector( pseudoJet1px, pseudoJet1py, pseudoJet1pz, pseudoJet1energy)
0218                 event.pseudoJet2 = ROOT.TLorentzVector( pseudoJet2px, pseudoJet2py, pseudoJet2pz, pseudoJet2energy)
0219             else:
0220                 event.pseudoJet2 = ROOT.TLorentzVector( pseudoJet1px, pseudoJet1py, pseudoJet1pz, pseudoJet1energy)
0221                 event.pseudoJet1 = ROOT.TLorentzVector( pseudoJet2px, pseudoJet2py, pseudoJet2pz, pseudoJet2energy)
0222
0223             ###
0224
0225             event.mr = self.computeMR(event.pseudoJet1, event.pseudoJet2)
0226             event.mtr = self.computeMTR(event.pseudoJet1, event.pseudoJet2, metp4)
0227             event.r = self.computeR(event.pseudoJet1, event.pseudoJet2, metp4)
0228
0229
0230
0231 #### do the razor with one or two b jets (medium CSV)
0232         if len(event.bjetsMedium)>=2:
0233
0234             bJet1 = ROOT.TLorentzVector(event.bjetsMedium[0].px(), event.bjetsMedium[0].py(), event.bjetsMedium[0].pz(), event.bjetsMedium[0].energy())
0235             bJet2 = ROOT.TLorentzVector(event.bjetsMedium[1].px(), event.bjetsMedium[1].py(), event.bjetsMedium[1].pz(), event.bjetsMedium[1].energy())
0236
0237             event.mr_bb  = self.computeMR(bJet1, bJet2)
0238             event.mtr_bb = self.computeMTR(bJet1, bJet2, metp4)
0239             event.r_bb   = self.computeR(bJet1, bJet2, metp4)
0240
0241 #            print 'MR(2b)',event.mr_bb
0242         if len(event.bjetsMedium)==1:
0243
0244             objects40jcCSV = [ j for j in event.cleanJets if j.pt() > 40 and abs(j.eta())<2.5 and j.p4()!=event.bjetsMedium[0].p4() ]
0245             objects40jcCSV.sort(key = lambda l : l.btag('combinedInclusiveSecondaryVertexV2BJetTags'), reverse = True)
0246
0247             if len(objects40jcCSV)>0:
0248
0249                 bJet1 = ROOT.TLorentzVector(event.bjetsMedium[0].px(), event.bjetsMedium[0].py(), event.bjetsMedium[0].pz(), event.bjetsMedium[0].energy())
0250                 bJet2 = ROOT.TLorentzVector(objects40jcCSV[0].px(), objects40jcCSV[0].py(), objects40jcCSV[0].pz(), objects40jcCSV[0].energy())
0251
0252                 event.mr_bb = self.computeMR(bJet1, bJet2)
0253                 event.mtr_bb = self.computeMTR(bJet1, bJet2, metp4)
0254                 event.r_bb = self.computeR(bJet1, bJet2, metp4)
0255 ##                print 'MRbb(1b)',event.mr_bb
0256
0257 ## ===> leptonic MR
0258         if not self.cfg_ana.doOnlyDefault:
0259             if len(event.selectedLeptons)>=2:
0260
0261                 lep1 = ROOT.TLorentzVector(event.selectedLeptons[0].px(), event.selectedLeptons[0].py(), event.selectedLeptons[0].pz(), event.selectedLeptons[0].energy())
0262                 lep2 = ROOT.TLorentzVector(event.selectedLeptons[1].px(), event.selectedLeptons[1].py(), event.selectedLeptons[1].pz(), event.selectedLeptons[1].energy())
0263
0264                 event.mr_lept = self.computeMR(lep1, lep2)
0265                 event.mtr_lept = self.computeMTR(lep1, lep2, metp4)
0266                 event.r_lept = self.computeR(lep1, lep2, metp4)
0267
0268
0269
0270 ###
0271
0272     def process(self, event):
0273         self.readCollections( event.input )
0274
0275         event.mr_gen=-999
0276         event.mtr_gen=-999
0277         event.r_gen=-999
0278
0279         event.mr_bb=-999
0280         event.mtr_bb=-999
0281         event.r_bb=-999
0282
0283         event.mr_lept=-999
0284         event.mtr_lept=-999
0285         event.r_lept=-999
0286
0287         event.mr_had=-999
0288         event.mtr_had=-999
0289         event.r_had=-999
0290
0291         event.mr=-999
0292         event.mtr=-999
0293         event.r=-999
0294
0295         event.pseudoJet1 = ROOT.TLorentzVector( 0, 0, 0, 0 )
0296         event.pseudoJet2 = ROOT.TLorentzVector( 0, 0, 0, 0 )
0297
0298         ###
0299
0300         self.makeRAZOR(event)
0301
0302 #        print 'variables computed: MR=',event.mr_had,'R=',event.r,'MTR=',event.mtr
0303 #        print 'pseudoJet1 px=',event.pseudoJet1.px(),' py=',event.pseudoJet1.py(),' pz=',event.pseudoJet1.pz()
0304 #        print 'pseudoJet2 px=',event.pseudoJet2.px(),' py=',event.pseudoJet2.py(),' pz=',event.pseudoJet2.pz()
0305
0306         return True