ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/src/TBSideband.py
Revision: 1.1
Committed: Fri Jul 1 19:51:06 2011 UTC (13 years, 10 months ago) by knash
Content type: text/x-python
Branch: MAIN
Log Message:
test

File Contents

# User Rev Content
1 knash 1.1 #! /usr/bin/env python
2     import os
3     import glob
4     import math
5     import ROOT
6     from ROOT import *
7     import sys
8     from DataFormats.FWLite import Events, Handle
9     from optparse import OptionParser
10     gROOT.Macro("rootlogon.C")
11    
12     parser = OptionParser()
13    
14     parser.add_option('-i', '--inputDir', metavar='F', type='string', action='store',
15     default = 'Jet_Run2011A-May10ReReco_ttbsm_v6_ttbsmTuples_v3',
16     dest = 'inputDir',
17     help = 'Input Dir')
18    
19     parser.add_option('-p', '--prefix', metavar='N', type='string', action='store',
20     default = './',
21     dest = 'prefix',
22     help = 'Dir prefix')
23    
24     (options, args) = parser.parse_args()
25    
26     files = glob.glob( options.prefix + '/' + options.inputDir + "/res/*.root" )
27     #print files
28    
29     # We select all the events:
30     events = Events (files)
31    
32     #Load all hemisphere 0 jets
33     #---------------------------------------------------------------------------------------------------------------------#
34    
35     hemis0topHandle = Handle( "vector<ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<double> > >")
36     hemis0topLabel = ( "ttbsmAna", "topTagP4Hemis0" )
37    
38     hemis0NSubJetsHandle = Handle ( "vector<double> " )
39     hemis0NSubJetsLabel = ( "ttbsmAna" , "topTagNSubjetsHemis0")
40    
41     hemis0MinMassHandle = Handle( "std::vector<double>" )
42     hemis0MinMassLabel = ( "ttbsmAna", "topTagMinMassHemis0" )
43    
44     hemis0TopMassHandle = Handle( "std::vector<double>" )
45     hemis0TopMassLabel = ( "ttbsmAna", "topTagTopMassHemis0" )
46    
47     hemis0wbHandle = Handle ("vector<ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<double> > >")
48     hemis0wbLabel = ( "ttbsmAna", "wTagP4Hemis0" )
49    
50     hemis0BDiscHandle = Handle( "std::vector<double>" )
51     hemis0BDiscLabel = ( "ttbsmAna" , "wTagBDiscHemis0" )
52    
53     #---------------------------------------------------------------------------------------------------------------------#
54    
55    
56     # Load all Hemisphere 1 jets
57     #---------------------------------------------------------------------------------------------------------------------#
58    
59     hemis1topHandle = Handle( "vector<ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<double> > >")
60     hemis1topLabel = ( "ttbsmAna", "topTagP4Hemis1" )
61    
62     hemis1NSubJetsHandle = Handle ( "vector<double> " )
63     hemis1NSubJetsLabel = ( "ttbsmAna" , "topTagNSubjetsHemis1")
64    
65     hemis1MinMassHandle = Handle( "std::vector<double>" )
66     hemis1MinMassLabel = ( "ttbsmAna", "topTagMinMassHemis1" )
67    
68     hemis1TopMassHandle = Handle( "std::vector<double>" )
69     hemis1TopMassLabel = ( "ttbsmAna", "topTagTopMassHemis1" )
70    
71     hemis1wbHandle = Handle( "vector<ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<double> > >")
72     hemis1wbLabel = ( "ttbsmAna", "wTagP4Hemis1" )
73    
74     hemis1BdiscHandle = Handle( "std::vector<double>" )
75     hemis1BdiscLabel = ( "ttbsmAna", "wTagBDiscHemis1" )
76    
77     hemis1Jet3Handle = Handle("int")
78     hemis1Jet3Label = ( "ttbsmAna", "jet3Hemis1" )
79    
80     hemis1MuHandle = Handle( "std::vector<double>")
81     hemis1MuLabel = ( "ttbsmAna", "wTagMuHemis1" )
82    
83     #---------------------------------------------------------------------------------------------------------------------#
84    
85     f = ROOT.TFile( "TBSideband.root", "recreate" )
86     f.cd()
87    
88     c = TCanvas('c', 'Signal to Background Estimates', 1200, 400)
89     c1 = TCanvas('c1', 'Signal to Background Estimates', 600, 400)
90     c.cd()
91     c.Divide(3,1)
92    
93     print "Creating histograms"
94    
95     #Define Histograms
96     #---------------------------------------------------------------------------------------------------------------------#
97    
98    
99     masswprimewb = ROOT.TH1F("masswprimewb", "Mtb in type 1+2 W Btag h1", 35, 500, 3000 )
100     masswprimewob = ROOT.TH1F("masswprimewob", "Mtb in type 1+2 WO Btag h1", 35, 500, 3000 )
101    
102     masswprimet11sr = ROOT.TH1F("masswprimet11sr", "Mtb in type 1+1 channel", 35, 500, 3000 )
103     masswprimet11srnb = ROOT.TH1F("masswprimet11srnb", "Mtb in type 1+1 channel no B", 35, 500, 3000 )
104     masswprimet11sb1 = ROOT.TH1F("masswprimet11sb1", "Mtb sideband in type b+1 channel Mjet", 35, 500, 3000 )
105     masswprimet11sb2 = ROOT.TH1F("masswprimet11sb2", "Mtb sideband in type b+1 channel Minmass", 35, 500, 3000 )
106     masswprimet11sb3 = ROOT.TH1F("masswprimet11sb3", "Mtb sideband in type b+1 channel Nsj", 35, 500, 3000 )
107     masswprimet11sb1nb = ROOT.TH1F("masswprimet11sb1nb", "Mtb sideband in type b+1 channel Nsj + pm no B", 35, 500, 3000 )
108     masswprimet11sb2nb = ROOT.TH1F("masswprimet11sb2nb", "Mtb sideband in type b+1 channel pm no B", 35, 500, 3000 )
109     masswprimet11sb3nb = ROOT.TH1F("masswprimet11sb3nb", "Mtb sideband in type b+1 channel Nsj no B", 35, 500, 3000 )
110    
111     masswprimet12sb2wob = ROOT.TH1F("masswprimet12sb2wob", "Mtb sideband2 in type 1+2 channel WO b", 35, 500, 3000 )
112     wCandVsTopCandh1sb = ROOT.TH2F("wCandVsTopCandh1sb", "W Cand Vs Top Cand Mass in sideband", 50, 0, 250, 100, 0, 500 )
113     wCandVsTopCandh1sr = ROOT.TH2F("wCandVsTopCandh1sr", "W Cand Vs Top Cand Mass in signal region", 50, 0, 250, 100, 0, 500 )
114    
115    
116     #---------------------------------------------------------------------------------------------------------------------#
117    
118     # loop over events
119     #---------------------------------------------------------------------------------------------------------------------#
120     masswprimet11sr.Sumw2()
121     masswprimewob.Sumw2()
122     count = 0
123     print "Start looping"
124     for event in events:
125     #if count > 20000:
126     #break
127     pairMass = 0.0
128     count = count + 1
129     nwp = 0
130    
131     if count % 10000 == 0 :
132     print '--------- Processing Event ' + str(count)
133    
134     # We load up the three vectors for each event.
135     #---------------------------------------------------------------------------------------------------------------------#
136    
137     event.getByLabel (hemis1BdiscLabel, hemis1BdiscHandle)
138     wJetBDisch1 = hemis1BdiscHandle.product()
139    
140     event.getByLabel (hemis0BDiscLabel, hemis0BDiscHandle)
141     wJetBDisch0 = hemis0BDiscHandle.product()
142    
143     event.getByLabel (hemis1wbLabel, hemis1wbHandle)
144     wJetsh1 = hemis1wbHandle.product()
145    
146     event.getByLabel (hemis0wbLabel, hemis0wbHandle)
147     wJetsh0 = hemis0wbHandle.product()
148    
149     event.getByLabel (hemis1Jet3Label, hemis1Jet3Handle )
150     jet3 = (hemis1Jet3Handle.product())[0]
151    
152     event.getByLabel (hemis1MuLabel, hemis1MuHandle)
153     wJetMu = hemis1MuHandle.product()
154    
155     event.getByLabel (hemis1topLabel, hemis1topHandle)
156     topJetsh1 = hemis1topHandle.product()
157    
158     event.getByLabel (hemis0topLabel, hemis0topHandle)
159     topJetsh0 = hemis0topHandle.product()
160    
161     event.getByLabel (hemis0TopMassLabel, hemis0TopMassHandle)
162     topJetMassh0 = hemis0TopMassHandle.product()
163    
164     event.getByLabel (hemis1TopMassLabel, hemis1TopMassHandle)
165     topJetMassh1 = hemis1TopMassHandle.product()
166    
167     event.getByLabel (hemis0MinMassLabel, hemis0MinMassHandle)
168     topJetMinMassh0 = hemis0MinMassHandle.product()
169    
170     event.getByLabel (hemis1MinMassLabel, hemis1MinMassHandle)
171     topJetMinMassh1 = hemis1MinMassHandle.product()
172    
173     events.getByLabel ( hemis0NSubJetsLabel, hemis0NSubJetsHandle )
174     NJetsh0 = hemis0NSubJetsHandle.product()
175    
176     events.getByLabel ( hemis1NSubJetsLabel , hemis1NSubJetsHandle )
177     NJetsh1 = hemis1NSubJetsHandle.product()
178    
179     #---------------------------------------------------------------------------------------------------------------------#
180    
181     #This gives the conditional statements for type 1+2 and the two
182     #type 1+1 cases to be sure the number of jets exist
183     #---------------------------------------------------------------------------------------------------------------------#
184    
185     njets12 = (len(wJetsh1) >= 2) and (len(wJetsh0) >= 1)
186     njets11b0 = ((len(topJetsh1) >= 1) and (len(wJetsh0) >= 1))
187     njets11b1 = ((len(topJetsh0) >= 1) and (len(wJetsh1) >= 1))
188    
189     #---------------------------------------------------------------------------------------------------------------------#
190    
191     #We look at type 1+1 W` first
192     #---------------------------------------------------------------------------------------------------------------------#
193    
194     #If b is the highest pt jet
195     if njets11b0 :
196     sidebandth1t1 = (topJetMassh1[0] > 100 and topJetMassh1[0] < 140) or (topJetMassh1[0] > 250 and topJetMassh1[0] < 300)
197     sidebandth1t2 = topJetMassh1[0] > 140 and topJetMassh1[0] < 250 and topJetMinMassh1[0] <= 50 and NJetsh1[0] > 2
198     sidebandth1t3 = topJetMassh1[0] > 140 and topJetMassh1[0] < 250 and NJetsh1[0] <= 2
199     hasBTag0 = wJetBDisch0[0] > 3.3
200     hasTopTag1 = topJetMassh1[0] > 140 and topJetMassh1[0] < 250 and topJetMinMassh1[0] > 50 and NJetsh1[0] > 2
201     hasTopMassh1 = topJetMassh1[0] > 140 and topJetMassh1[0] < 250
202     hasTopLikeh1 = topJetMinMassh1[0] > 50 and NJetsh1[0] > 2
203     highptbjh0 = wJetsh0[0].pt() > 350
204     highpttjh1 = topJetsh1[0].pt() > 350
205     if highpttjh1 and highptbjh0 and sidebandth1t1 and hasTopLikeh1:
206     masswprimet11sb1nb.Fill((topJetsh1[0]+wJetsh0[0]).mass())
207     if hasBTag0 :
208     masswprimet11sb1.Fill((topJetsh1[0]+wJetsh0[0]).mass())
209     if highpttjh1 and highptbjh0 and sidebandth1t2:
210     masswprimet11sb2nb.Fill((topJetsh1[0]+wJetsh0[0]).mass())
211     if hasBTag0 :
212     masswprimet11sb2.Fill((topJetsh1[0]+wJetsh0[0]).mass())
213     if highpttjh1 and highptbjh0 and sidebandth1t3 :
214     masswprimet11sb3nb.Fill((topJetsh1[0]+wJetsh0[0]).mass())
215     if hasBTag0 :
216     masswprimet11sb3.Fill((topJetsh1[0]+wJetsh0[0]).mass())
217     if highptbjh0 and highpttjh1 and hasTopTag1 :
218     masswprimet11srnb.Fill((topJetsh1[0]+wJetsh0[0]).mass())
219     if hasBTag0 :
220     masswprimet11sr.Fill((topJetsh1[0]+wJetsh0[0]).mass())
221    
222    
223     #If t is the highest pt jet
224     if njets11b1 :
225     sidebandth0t1 = (topJetMassh0[0] > 100 and topJetMassh0[0] < 140) or (topJetMassh0[0] > 250 and topJetMassh0[0] < 300)
226     sidebandth0t2 = topJetMassh0[0] > 140 and topJetMassh0[0] < 250 and topJetMinMassh0[0] <= 50 and NJetsh0[0] > 2
227     sidebandth0t3 = topJetMassh0[0] > 140 and topJetMassh0[0] < 250 and NJetsh0[0] <= 2
228     hasBTag1 = wJetBDisch1[0] > 3.3
229     hasTopTag0 = topJetMassh0[0] > 140 and topJetMassh0[0] < 250 and topJetMinMassh0[0] > 50 and NJetsh0[0] > 2
230     hasTopMassh0 = topJetMassh0[0] > 140 and topJetMassh0[0] < 250
231     hasTopLikeh0 = topJetMinMassh0[0] > 50 and NJetsh0[0] > 2
232     highpttjh0 = topJetsh0[0].pt() > 350
233     highptbjh1 = wJetsh1[0].pt() > 350
234     if highptbjh1 and highpttjh0 and sidebandth0t1 and (not sidebandth1t1) and hasTopLikeh0:
235     masswprimet11sb1nb.Fill((topJetsh0[0]+wJetsh1[0]).mass())
236     if hasBTag1 :
237     masswprimet11sb1.Fill((topJetsh0[0]+wJetsh1[0]).mass())
238     if highptbjh1 and highpttjh0 and sidebandth0t2 and (not sidebandth1t2):
239     masswprimet11sb2nb.Fill((topJetsh0[0]+wJetsh1[0]).mass())
240     if hasBTag1 :
241     masswprimet11sb2.Fill((topJetsh0[0]+wJetsh1[0]).mass())
242     if highptbjh1 and highpttjh0 and sidebandth0t3 and (not sidebandth1t3):
243     masswprimet11sb3nb.Fill((topJetsh0[0]+wJetsh1[0]).mass())
244     if hasBTag1 :
245     masswprimet11sb3.Fill((topJetsh0[0]+wJetsh1[0]).mass())
246     if highptbjh1 and highpttjh0 and hasTopTag0 and (not hasTopTag1):
247     masswprimet11srnb.Fill((topJetsh0[0]+wJetsh1[0]).mass())
248     if hasBTag1 :
249     masswprimet11sr.Fill((topJetsh0[0]+wJetsh1[0]).mass())
250    
251     #---------------------------------------------------------------------------------------------------------------------#
252    
253    
254     #Now we look at type 1+2 W`
255     #---------------------------------------------------------------------------------------------------------------------#
256    
257     if njets12 :
258     #Type2Top: Wjet pt > 200 mu < 0.4 and Bjet pt > 30
259     highptbjh0 = wJetsh0[0].pt() > 350
260     hemis1kincuts = (highptbjh0) and (wJetsh1[0].pt() > 200) and (wJetMu[0] < 0.4) and (wJetsh1[jet3].pt() > 30)
261     pairMass = (wJetsh1[jet3]+wJetsh1[0]).mass()
262     hasBTag0 = wJetBDisch0[0] > 3.3
263     #Toplike mass cuts
264     hash1Top = wJetsh1[0].mass() > 60 and wJetsh1[0].mass() < 130 and pairMass > 140 and pairMass < 250
265     sidebandt2top = wJetsh1[0].mass() > 40 and wJetsh1[0].mass() < 130 and pairMass > 100 and pairMass < 300
266     hasBTag2 = wJetBDisch1[jet3] > 3.3
267     if highptbjh0 and hemis1kincuts :
268     if hash1Top :
269     wCandVsTopCandh1sr.Fill( wJetsh1[0].mass() , pairMass )
270     if hasBTag0:
271     masswprimewob.Fill((wJetsh1[jet3]+wJetsh1[0]+wJetsh0[0]).mass())
272     if hasBTag2 :
273     masswprimewb.Fill((wJetsh1[jet3]+wJetsh1[0]+wJetsh0[0]).mass())
274     if sidebandt2top and (not hash1Top):
275     wCandVsTopCandh1sb.Fill( wJetsh1[0].mass() , pairMass )
276     if hasBTag0:
277     masswprimet12sb2wob.Fill((wJetsh1[jet3]+wJetsh1[0]+wJetsh0[0]).mass())
278    
279     #---------------------------------------------------------------------------------------------------------------------#
280    
281     sb1norm = masswprimet11srnb.GetEntries()/masswprimet11sb1nb.GetEntries()
282     sb2norm = masswprimet11srnb.GetEntries()/masswprimet11sb2nb.GetEntries()
283     sb3norm = masswprimet11srnb.GetEntries()/masswprimet11sb3nb.GetEntries()
284     sb1t12norm = wCandVsTopCandh1sr.GetEntries()/wCandVsTopCandh1sb.GetEntries()
285    
286     masswprimet12sb2wob.Scale( sb1t12norm )
287     masswprimet11sb1.Scale( sb1norm )
288     masswprimet11sb2.Scale( sb2norm )
289     masswprimet11sb3.Scale( sb3norm )
290     print "Normalizations: SB1: " + str(sb1norm) + " SB2: " + str(sb2norm) + " SB3: " + str(sb3norm) + " Type B+2: " + str(sb1t12norm)
291     masswprimet11sb1.SetMaximum(75)
292     masswprimet11sb2.SetMaximum(75)
293     masswprimet11sb3.SetMaximum(75)
294     masswprimet12sb2wob.SetMaximum(125)
295    
296     masswprimet11sb1.SetFillColor(TROOT.kYellow)
297     masswprimet11sb2.SetFillColor(TROOT.kYellow)
298     masswprimet11sb3.SetFillColor(TROOT.kYellow)
299     masswprimet12sb2wob.SetFillColor(TROOT.kYellow)
300    
301     c.cd(1)
302     masswprimet11sb1.Draw()
303     masswprimet11sr.Draw('same')
304     c.cd(2)
305     masswprimet11sb2.Draw()
306     masswprimet11sr.Draw('same')
307     c.cd(3)
308     masswprimet11sb3.Draw()
309     masswprimet11sr.Draw('same')
310    
311     c1.cd()
312     masswprimet12sb2wob.Draw()
313     masswprimewob.Draw('same')
314    
315     #c.Print('Sidebands1.pdf', 'pdf')
316     c.Print('Sidebandstb1.root', 'root')
317     c1.Print('Sidebandstb2.root', 'root')
318    
319    
320     f.cd()
321     f.Write()
322     f.Close()