CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
RazorAnalyzer.py
Go to the documentation of this file.
1 import operator
2 import itertools
3 import copy
4 from math import *
5 
6 from ROOT import std
7 from ROOT import TLorentzVector, TVector3, TVectorD
8 
9 from PhysicsTools.Heppy.analyzers.core.Analyzer import Analyzer
10 from PhysicsTools.HeppyCore.framework.event import Event
11 from PhysicsTools.HeppyCore.statistics.counter import Counter, Counters
12 from PhysicsTools.Heppy.analyzers.core.AutoHandle import AutoHandle
13 
14 from PhysicsTools.HeppyCore.utils.deltar import deltaR
15 
16 from ROOT.heppy import Megajet
17 from ROOT.heppy import ReclusterJets
18 
19 import ROOT
20 
21 import os
22 
23 class RazorAnalyzer( Analyzer ):
24  def __init__(self, cfg_ana, cfg_comp, looperName ):
25  super(RazorAnalyzer,self).__init__(cfg_ana,cfg_comp,looperName)
26 
27  def declareHandles(self):
28  super(RazorAnalyzer, self).declareHandles()
29  #genJets
30  self.handles['genJets'] = AutoHandle( 'slimmedGenJets','std::vector<reco::GenJet>')
31 
32  def beginLoop(self, setup):
33  super(RazorAnalyzer,self).beginLoop(setup)
34  self.counters.addCounter('pairs')
35  count = self.counters.counter('pairs')
36  count.register('all events')
37 
38  def computeMR(self, ja, jb):
39  A = ja.P();
40  B = jb.P();
41  az = ja.Pz();
42  bz = jb.Pz();
43  mr = sqrt((A+B)*(A+B)-(az+bz)*(az+bz));
44  return mr
45 
46  def computeMTR(self, ja, jb, met):
47 
48  mtr = met.Vect().Mag()*(ja.Pt()+jb.Pt()) - met.Vect().Dot(ja.Vect()+jb.Vect());
49  mtr = sqrt(mtr/2.);
50  return mtr;
51 
52  def computeR(self, ja, jb, met):
53 
54  mr = self.computeMR(ja,jb)
55  mtr = self.computeMTR(ja,jb,met)
56  r = 999999. if mr <= 0 else mtr/mr
57  return r
58 
59 
60  def makeRAZOR(self, event):
61 # print '==> INSIDE THE PRINT MT2'
62 # print 'MET=',event.met.pt()
63 
64  import array
65  import numpy
66 
67 
68 ## ===> hadronic RAZOR
69 
70  (met, metphi) = event.met.pt(), event.met.phi()
71  metp4 = ROOT.TLorentzVector()
72  metp4.SetPtEtaPhiM(met,0,metphi,0)
73 
74  objects40jc = [ j for j in event.cleanJets if j.pt() > 40 and abs(j.eta())<2.5 ]
75 
76 #### get megajets (association method: default 1 = minimum sum of the invariant masses of the two megajets)
77 
78  if len(objects40jc)>=2:
79 
80  pxvec = ROOT.std.vector(float)()
81  pyvec = ROOT.std.vector(float)()
82  pzvec = ROOT.std.vector(float)()
83  Evec = ROOT.std.vector(float)()
84  grouping = ROOT.std.vector(int)()
85 
86  for jet in objects40jc:
87  pxvec.push_back(jet.px())
88  pyvec.push_back(jet.py())
89  pzvec.push_back(jet.pz())
90  Evec.push_back(jet.energy())
91 
92  megajet = Megajet(pxvec, pyvec, pzvec, Evec, 1)
93 
94  pseudoJet1px = megajet.getAxis1()[0] * megajet.getAxis1()[3]
95  pseudoJet1py = megajet.getAxis1()[1] * megajet.getAxis1()[3]
96  pseudoJet1pz = megajet.getAxis1()[2] * megajet.getAxis1()[3]
97  pseudoJet1energy = megajet.getAxis1()[4]
98 
99  pseudoJet2px = megajet.getAxis2()[0] * megajet.getAxis2()[3]
100  pseudoJet2py = megajet.getAxis2()[1] * megajet.getAxis2()[3]
101  pseudoJet2pz = megajet.getAxis2()[2] * megajet.getAxis2()[3]
102  pseudoJet2energy = megajet.getAxis2()[4]
103 
104  pseudoJet1pt2 = pseudoJet1px*pseudoJet1px + pseudoJet1py*pseudoJet1py
105  pseudoJet2pt2 = pseudoJet2px*pseudoJet2px + pseudoJet2py*pseudoJet2py
106 
107  if pseudoJet1pt2 >= pseudoJet2pt2:
108  event.pseudoJet1_had = ROOT.TLorentzVector( pseudoJet1px, pseudoJet1py, pseudoJet1pz, pseudoJet1energy)
109  event.pseudoJet2_had = ROOT.TLorentzVector( pseudoJet2px, pseudoJet2py, pseudoJet2pz, pseudoJet2energy)
110  else:
111  event.pseudoJet2_had = ROOT.TLorentzVector( pseudoJet1px, pseudoJet1py, pseudoJet1pz, pseudoJet1energy)
112  event.pseudoJet1_had = ROOT.TLorentzVector( pseudoJet2px, pseudoJet2py, pseudoJet2pz, pseudoJet2energy)
113 
114  event.mr_had = self.computeMR(event.pseudoJet1_had, event.pseudoJet2_had)
115  event.mtr_had = self.computeMTR(event.pseudoJet1_had, event.pseudoJet2_had, metp4)
116  event.r_had = self.computeR(event.pseudoJet1_had, event.pseudoJet2_had, metp4)
117 
118 #### do same things for GEN
119 
120  if self.cfg_comp.isMC:
121  (genmet, genmetphi) = event.met.genMET().pt(), event.met.genMET().phi()
122  genmetp4 = ROOT.TLorentzVector()
123  genmetp4.SetPtEtaPhiM(genmet,0,genmetphi,0)
124 
125  allGenJets = [ x for x in self.handles['genJets'].product() ]
126  objects40jc_Gen = [ j for j in allGenJets if j.pt() > 40 and abs(j.eta())<2.5 ]
127 
128  if len(objects40jc_Gen)>=2:
129 
130  pxvec = ROOT.std.vector(float)()
131  pyvec = ROOT.std.vector(float)()
132  pzvec = ROOT.std.vector(float)()
133  Evec = ROOT.std.vector(float)()
134  grouping = ROOT.std.vector(int)()
135 
136  for jet in objects40jc_Gen:
137  pxvec.push_back(jet.px())
138  pyvec.push_back(jet.py())
139  pzvec.push_back(jet.pz())
140  Evec.push_back(jet.energy())
141 
142  megajet = Megajet(pxvec, pyvec, pzvec, Evec, 1)
143 
144  pseudoJet1px = megajet.getAxis1()[0] * megajet.getAxis1()[3]
145  pseudoJet1py = megajet.getAxis1()[1] * megajet.getAxis1()[3]
146  pseudoJet1pz = megajet.getAxis1()[2] * megajet.getAxis1()[3]
147  pseudoJet1energy = megajet.getAxis1()[4]
148 
149  pseudoJet2px = megajet.getAxis2()[0] * megajet.getAxis2()[3]
150  pseudoJet2py = megajet.getAxis2()[1] * megajet.getAxis2()[3]
151  pseudoJet2pz = megajet.getAxis2()[2] * megajet.getAxis2()[3]
152  pseudoJet2energy = megajet.getAxis2()[4]
153 
154  pseudoJet1pt2 = pseudoJet1px*pseudoJet1px + pseudoJet1py*pseudoJet1py
155  pseudoJet2pt2 = pseudoJet2px*pseudoJet2px + pseudoJet2py*pseudoJet2py
156 
157  if pseudoJet1pt2 >= pseudoJet2pt2:
158  pseudoJet1_gen = ROOT.TLorentzVector( pseudoJet1px, pseudoJet1py, pseudoJet1pz, pseudoJet1energy)
159  pseudoJet2_gen = ROOT.TLorentzVector( pseudoJet2px, pseudoJet2py, pseudoJet2pz, pseudoJet2energy)
160  else:
161  pseudoJet2_gen = ROOT.TLorentzVector( pseudoJet1px, pseudoJet1py, pseudoJet1pz, pseudoJet1energy)
162  pseudoJet1_gen = ROOT.TLorentzVector( pseudoJet2px, pseudoJet2py, pseudoJet2pz, pseudoJet2energy)
163 
164  event.mr_gen = self.computeMR(pseudoJet1_gen, pseudoJet2_gen)
165  event.mtr_gen = self.computeMTR(pseudoJet1_gen, pseudoJet2_gen, genmetp4)
166  event.r_gen = self.computeR(pseudoJet1_gen, pseudoJet2_gen, genmetp4)
167  else:
168  event.mr_gen = -999
169  event.mtr_gen = -999
170  event.r_gen = -999
171 
172 
173 ## ===> full RAZOR (jets + leptons)
174  objects10lc = [ l for l in event.selectedLeptons if l.pt() > 10 and abs(l.eta())<2.5 ]
175  if hasattr(event, 'selectedIsoCleanTrack'):
176  objects10lc = [ l for l in event.selectedLeptons if l.pt() > 10 and abs(l.eta())<2.5 ] + [ t for t in event.selectedIsoCleanTrack ]
177 
178  objects40j10lc = objects40jc + objects10lc
179 
180  objects40j10lc.sort(key = lambda obj : obj.pt(), reverse = True)
181 
182  if len(objects40j10lc)>=2:
183 
184  pxvec = ROOT.std.vector(float)()
185  pyvec = ROOT.std.vector(float)()
186  pzvec = ROOT.std.vector(float)()
187  Evec = ROOT.std.vector(float)()
188  grouping = ROOT.std.vector(int)()
189 
190  for obj in objects40j10lc:
191  pxvec.push_back(obj.px())
192  pyvec.push_back(obj.py())
193  pzvec.push_back(obj.pz())
194  Evec.push_back(obj.energy())
195 
196  #for obj in objects_fullmt2:
197  # print "pt: ", obj.pt(), ", eta: ", obj.eta(), ", phi: ", obj.phi(), ", mass: ", obj.mass()
198 
199  #### get megajets (association method: default 1 = minimum sum of the invariant masses of the two megajets)
200 
201  megajet = Megajet(pxvec, pyvec, pzvec, Evec, 1)
202 
203  pseudoJet1px = megajet.getAxis1()[0] * megajet.getAxis1()[3]
204  pseudoJet1py = megajet.getAxis1()[1] * megajet.getAxis1()[3]
205  pseudoJet1pz = megajet.getAxis1()[2] * megajet.getAxis1()[3]
206  pseudoJet1energy = megajet.getAxis1()[4]
207 
208  pseudoJet2px = megajet.getAxis2()[0] * megajet.getAxis2()[3]
209  pseudoJet2py = megajet.getAxis2()[1] * megajet.getAxis2()[3]
210  pseudoJet2pz = megajet.getAxis2()[2] * megajet.getAxis2()[3]
211  pseudoJet2energy = megajet.getAxis2()[4]
212 
213  pseudoJet1pt2 = pseudoJet1px*pseudoJet1px + pseudoJet1py*pseudoJet1py
214  pseudoJet2pt2 = pseudoJet2px*pseudoJet2px + pseudoJet2py*pseudoJet2py
215 
216  if pseudoJet1pt2 >= pseudoJet2pt2:
217  event.pseudoJet1 = ROOT.TLorentzVector( pseudoJet1px, pseudoJet1py, pseudoJet1pz, pseudoJet1energy)
218  event.pseudoJet2 = ROOT.TLorentzVector( pseudoJet2px, pseudoJet2py, pseudoJet2pz, pseudoJet2energy)
219  else:
220  event.pseudoJet2 = ROOT.TLorentzVector( pseudoJet1px, pseudoJet1py, pseudoJet1pz, pseudoJet1energy)
221  event.pseudoJet1 = ROOT.TLorentzVector( pseudoJet2px, pseudoJet2py, pseudoJet2pz, pseudoJet2energy)
222 
223  ###
224 
225  event.mr = self.computeMR(event.pseudoJet1, event.pseudoJet2)
226  event.mtr = self.computeMTR(event.pseudoJet1, event.pseudoJet2, metp4)
227  event.r = self.computeR(event.pseudoJet1, event.pseudoJet2, metp4)
228 
229 
230 
231 #### do the razor with one or two b jets (medium CSV)
232  if len(event.bjetsMedium)>=2:
233 
234  bJet1 = ROOT.TLorentzVector(event.bjetsMedium[0].px(), event.bjetsMedium[0].py(), event.bjetsMedium[0].pz(), event.bjetsMedium[0].energy())
235  bJet2 = ROOT.TLorentzVector(event.bjetsMedium[1].px(), event.bjetsMedium[1].py(), event.bjetsMedium[1].pz(), event.bjetsMedium[1].energy())
236 
237  event.mr_bb = self.computeMR(bJet1, bJet2)
238  event.mtr_bb = self.computeMTR(bJet1, bJet2, metp4)
239  event.r_bb = self.computeR(bJet1, bJet2, metp4)
240 
241 # print 'MR(2b)',event.mr_bb
242  if len(event.bjetsMedium)==1:
243 
244  objects40jcCSV = [ j for j in event.cleanJets if j.pt() > 40 and abs(j.eta())<2.5 and j.p4()!=event.bjetsMedium[0].p4() ]
245  objects40jcCSV.sort(key = lambda l : l.btag('combinedInclusiveSecondaryVertexV2BJetTags'), reverse = True)
246 
247  if len(objects40jcCSV)>0:
248 
249  bJet1 = ROOT.TLorentzVector(event.bjetsMedium[0].px(), event.bjetsMedium[0].py(), event.bjetsMedium[0].pz(), event.bjetsMedium[0].energy())
250  bJet2 = ROOT.TLorentzVector(objects40jcCSV[0].px(), objects40jcCSV[0].py(), objects40jcCSV[0].pz(), objects40jcCSV[0].energy())
251 
252  event.mr_bb = self.computeMR(bJet1, bJet2)
253  event.mtr_bb = self.computeMTR(bJet1, bJet2, metp4)
254  event.r_bb = self.computeR(bJet1, bJet2, metp4)
255 ## print 'MRbb(1b)',event.mr_bb
256 
257 ## ===> leptonic MR
258  if not self.cfg_ana.doOnlyDefault:
259  if len(event.selectedLeptons)>=2:
260 
261  lep1 = ROOT.TLorentzVector(event.selectedLeptons[0].px(), event.selectedLeptons[0].py(), event.selectedLeptons[0].pz(), event.selectedLeptons[0].energy())
262  lep2 = ROOT.TLorentzVector(event.selectedLeptons[1].px(), event.selectedLeptons[1].py(), event.selectedLeptons[1].pz(), event.selectedLeptons[1].energy())
263 
264  event.mr_lept = self.computeMR(lep1, lep2)
265  event.mtr_lept = self.computeMTR(lep1, lep2, metp4)
266  event.r_lept = self.computeR(lep1, lep2, metp4)
267 
268 
269 
270 ###
271 
272  def process(self, event):
273  self.readCollections( event.input )
274 
275  event.mr_gen=-999
276  event.mtr_gen=-999
277  event.r_gen=-999
278 
279  event.mr_bb=-999
280  event.mtr_bb=-999
281  event.r_bb=-999
282 
283  event.mr_lept=-999
284  event.mtr_lept=-999
285  event.r_lept=-999
286 
287  event.mr_had=-999
288  event.mtr_had=-999
289  event.r_had=-999
290 
291  event.mr=-999
292  event.mtr=-999
293  event.r=-999
294 
295  event.pseudoJet1 = ROOT.TLorentzVector( 0, 0, 0, 0 )
296  event.pseudoJet2 = ROOT.TLorentzVector( 0, 0, 0, 0 )
297 
298  ###
299 
300  self.makeRAZOR(event)
301 
302 # print 'variables computed: MR=',event.mr_had,'R=',event.r,'MTR=',event.mtr
303 # print 'pseudoJet1 px=',event.pseudoJet1.px(),' py=',event.pseudoJet1.py(),' pz=',event.pseudoJet1.pz()
304 # print 'pseudoJet2 px=',event.pseudoJet2.px(),' py=',event.pseudoJet2.py(),' pz=',event.pseudoJet2.pz()
305 
306  return True
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:48
p4
double p4[4]
Definition: TauolaWrapper.h:92
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
relval_parameters_module.energy
string energy
Definition: relval_parameters_module.py:29
LaserDQM_cfg.process
tuple process
Definition: LaserDQM_cfg.py:3
Generated for CMSSW Reference Manual by doxygen 1.8.5