Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:32:27

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