File indexing completed on 2024-04-06 12:23:26
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
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
0062
0063
0064 import array
0065 import numpy
0066
0067
0068
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
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
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
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
0197
0198
0199
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
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
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
0256
0257
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
0303
0304
0305
0306 return True