| 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()
|
