| 1 |
#! /usr/bin/env python
|
| 2 |
import os
|
| 3 |
import glob
|
| 4 |
import math
|
| 5 |
import ROOT
|
| 6 |
import sys
|
| 7 |
from DataFormats.FWLite import Events, Handle
|
| 8 |
from optparse import OptionParser
|
| 9 |
|
| 10 |
parser = OptionParser()
|
| 11 |
|
| 12 |
parser.add_option('-i', '--inputDir', metavar='F', type='string', action='store',
|
| 13 |
default = 'Jet_Run2011A-May10ReReco_ttbsm_v6_ttbsmTuples_v3',
|
| 14 |
dest = 'inputDir',
|
| 15 |
help = 'Input Dir')
|
| 16 |
|
| 17 |
parser.add_option('-p', '--prefix', metavar='N', type='string', action='store',
|
| 18 |
default = './',
|
| 19 |
dest = 'prefix',
|
| 20 |
help = 'Dir prefix')
|
| 21 |
|
| 22 |
(options, args) = parser.parse_args()
|
| 23 |
|
| 24 |
files = glob.glob( options.prefix + '/' + options.inputDir + "/res/*.root" )
|
| 25 |
#print files
|
| 26 |
|
| 27 |
# We select all the events:
|
| 28 |
events = Events (files)
|
| 29 |
|
| 30 |
|
| 31 |
# Load all jets
|
| 32 |
#---------------------------------------------------------------------------------------------------------------------#
|
| 33 |
|
| 34 |
wbCandHandle = Handle ( "vector<ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<double> > >" )
|
| 35 |
wbCandLabel = ( "ttbsmAna" , "wTagP4")
|
| 36 |
|
| 37 |
wbDiscHandle = Handle( "std::vector<double>" )
|
| 38 |
wbDiscLabel = ( "ttbsmAna" , "wTagBDisc" )
|
| 39 |
|
| 40 |
#---------------------------------------------------------------------------------------------------------------------#
|
| 41 |
|
| 42 |
|
| 43 |
#Load all hemisphere 0 jets
|
| 44 |
#---------------------------------------------------------------------------------------------------------------------#
|
| 45 |
|
| 46 |
hemis0topHandle = Handle( "vector<ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<double> > >")
|
| 47 |
hemis0topLabel = ( "ttbsmAna", "topTagP4Hemis0" )
|
| 48 |
|
| 49 |
hemis0NSubJetsHandle = Handle ( "vector<double> " )
|
| 50 |
hemis0NSubJetsLabel = ( "ttbsmAna" , "topTagNSubjetsHemis0")
|
| 51 |
|
| 52 |
hemis0MinMassHandle = Handle( "std::vector<double>" )
|
| 53 |
hemis0MinMassLabel = ( "ttbsmAna", "topTagMinMassHemis0" )
|
| 54 |
|
| 55 |
hemis0TopMassHandle = Handle( "std::vector<double>" )
|
| 56 |
hemis0TopMassLabel = ( "ttbsmAna", "topTagTopMassHemis0" )
|
| 57 |
|
| 58 |
hemis0wbHandle = Handle ("vector<ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<double> > >")
|
| 59 |
hemis0wbLabel = ( "ttbsmAna", "wTagP4Hemis0" )
|
| 60 |
|
| 61 |
hemis0BDiscHandle = Handle( "std::vector<double>" )
|
| 62 |
hemis0BDiscLabel = ( "ttbsmAna" , "wTagBDiscHemis0" )
|
| 63 |
|
| 64 |
#---------------------------------------------------------------------------------------------------------------------#
|
| 65 |
|
| 66 |
|
| 67 |
# Load all Hemisphere 1 jets
|
| 68 |
#---------------------------------------------------------------------------------------------------------------------#
|
| 69 |
|
| 70 |
hemis1topHandle = Handle( "vector<ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<double> > >")
|
| 71 |
hemis1topLabel = ( "ttbsmAna", "topTagP4Hemis1" )
|
| 72 |
|
| 73 |
hemis1NSubJetsHandle = Handle ( "vector<double> " )
|
| 74 |
hemis1NSubJetsLabel = ( "ttbsmAna" , "topTagNSubjetsHemis1")
|
| 75 |
|
| 76 |
hemis1MinMassHandle = Handle( "std::vector<double>" )
|
| 77 |
hemis1MinMassLabel = ( "ttbsmAna", "topTagMinMassHemis1" )
|
| 78 |
|
| 79 |
hemis1TopMassHandle = Handle( "std::vector<double>" )
|
| 80 |
hemis1TopMassLabel = ( "ttbsmAna", "topTagTopMassHemis1" )
|
| 81 |
|
| 82 |
hemis1wbHandle = Handle( "vector<ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<double> > >")
|
| 83 |
hemis1wbLabel = ( "ttbsmAna", "wTagP4Hemis1" )
|
| 84 |
|
| 85 |
hemis1BdiscHandle = Handle( "std::vector<double>" )
|
| 86 |
hemis1BdiscLabel = ( "ttbsmAna", "wTagBDiscHemis1" )
|
| 87 |
|
| 88 |
hemis1Jet3Handle = Handle("int")
|
| 89 |
hemis1Jet3Label = ( "ttbsmAna", "jet3Hemis1" )
|
| 90 |
|
| 91 |
hemis1MuHandle = Handle( "std::vector<double>")
|
| 92 |
hemis1MuLabel = ( "ttbsmAna", "wTagMuHemis1" )
|
| 93 |
|
| 94 |
#---------------------------------------------------------------------------------------------------------------------#
|
| 95 |
|
| 96 |
f = ROOT.TFile( "TB:" + options.inputDir + ".root", "recreate" )
|
| 97 |
BBKGD = ROOT.TFile("BKGDII.root")
|
| 98 |
f.cd()
|
| 99 |
|
| 100 |
print "Creating histograms"
|
| 101 |
|
| 102 |
#Define Histograms
|
| 103 |
#---------------------------------------------------------------------------------------------------------------------#
|
| 104 |
|
| 105 |
#Raw Jet Plots
|
| 106 |
nJets = ROOT.TH1F("nJets", "Number of Jets", 10, -0.5, 9.5 )
|
| 107 |
mass2bJets = ROOT.TH1F("mass2bJets", "invariant mass of 2 b Jets", 100, 0, 2000 )
|
| 108 |
energyJets = ROOT.TH1F("energyJets", "sum energy of all jets > 30Gev", 100, 30, 3000 )
|
| 109 |
dphi = ROOT.TH1F("dphi", "delta-phi of b jets", 60, 0, 4 )
|
| 110 |
deta = ROOT.TH1F("deta", "delta-eta between b jets", 35, 0, 6 )
|
| 111 |
dR = ROOT.TH1F("dR", "delta-R of b jets", 80, 0, 6 )
|
| 112 |
|
| 113 |
#After Cuts
|
| 114 |
nJetsh1 = ROOT.TH1F("nJetsh1", "Number of Jets in h1", 10, -0.5, 9.5 )
|
| 115 |
ptbJetsh1 = ROOT.TH1F("ptbJetsh1", "pt of b-like Jets in h1", 150, 0, 800 )
|
| 116 |
ptbJetsh1wt = ROOT.TH1F("ptbJetsh1wt", "pt of b Jets in h1 with btag", 150, 0, 800 )
|
| 117 |
ptwJetsh1 = ROOT.TH1F("ptwJetsh1", "pt of w Jets in h1", 150, 0, 800 )
|
| 118 |
WJetmass = ROOT.TH1F("WJetmass", "mass of w tagged jets", 150, 0, 200 )
|
| 119 |
WJetmasswc = ROOT.TH1F("WJetmasswc", "mass of w tagged jets after mass cut", 150, 0, 200 )
|
| 120 |
topCandmass = ROOT.TH1F("topCandmass", "pairwise miss of h1 b+W", 80, 0, 400 )
|
| 121 |
topCandmasswb = ROOT.TH1F("topCandmasswb", "pairwise miss of h1 b+W with btag", 80, 0, 400 )
|
| 122 |
ptbh0 = ROOT.TH1F("ptbh0", "pt of b Jet in h0", 150, 0, 800 )
|
| 123 |
massbh0 = ROOT.TH1F("massbh0", "Mass of high pt b jet in h0", 45, 0, 200 )
|
| 124 |
|
| 125 |
masswprimet11 = ROOT.TH1F("masswprimet11", "Mtb in type b+1 channel", 40, 500, 3000 )
|
| 126 |
masswprimewb = ROOT.TH1F("masswprimewb", "Mtb in type b+2 W Btag h1", 40, 500, 3000 )
|
| 127 |
masswprimewob = ROOT.TH1F("masswprimewob", "Mtb in type b+2 WO Btag h1", 40, 500, 3000 )
|
| 128 |
mistagB11 = ROOT.TH1F("mistagB11", "Mtb type b+1 bkgd with mistag used for B", 40, 500, 3000 )
|
| 129 |
mistagB12 = ROOT.TH1F("mistagB12", "Mtb type b+2 bkgd with mistag used in hem 1", 40, 500, 3000 )
|
| 130 |
|
| 131 |
#---------------------------------------------------------------------------------------------------------------------#
|
| 132 |
|
| 133 |
# loop over events
|
| 134 |
#---------------------------------------------------------------------------------------------------------------------#
|
| 135 |
|
| 136 |
count = 0
|
| 137 |
print "Start looping"
|
| 138 |
mistag = BBKGD.Get("mistag")
|
| 139 |
for event in events:
|
| 140 |
nb = []
|
| 141 |
esum = 0.0
|
| 142 |
dRn = 0.0
|
| 143 |
dphin = 0.0
|
| 144 |
detan = 0.0
|
| 145 |
pairMass = 0.0
|
| 146 |
count = count + 1
|
| 147 |
|
| 148 |
if count % 10000 == 0 :
|
| 149 |
print '--------- Processing Event ' + str(count)
|
| 150 |
|
| 151 |
# We load up the three vectors for each event.
|
| 152 |
#---------------------------------------------------------------------------------------------------------------------#
|
| 153 |
|
| 154 |
event.getByLabel( wbCandLabel , wbCandHandle )
|
| 155 |
wbCandJets = wbCandHandle.product()
|
| 156 |
|
| 157 |
event.getByLabel ( wbDiscLabel , wbDiscHandle )
|
| 158 |
wJetBDisc = wbDiscHandle.product()
|
| 159 |
|
| 160 |
event.getByLabel (hemis1BdiscLabel, hemis1BdiscHandle)
|
| 161 |
wJetBDisch1 = hemis1BdiscHandle.product()
|
| 162 |
|
| 163 |
event.getByLabel (hemis0BDiscLabel, hemis0BDiscHandle)
|
| 164 |
wJetBDisch0 = hemis0BDiscHandle.product()
|
| 165 |
|
| 166 |
event.getByLabel (hemis1wbLabel, hemis1wbHandle)
|
| 167 |
wJetsh1 = hemis1wbHandle.product()
|
| 168 |
|
| 169 |
event.getByLabel (hemis0wbLabel, hemis0wbHandle)
|
| 170 |
wJetsh0 = hemis0wbHandle.product()
|
| 171 |
|
| 172 |
event.getByLabel (hemis1Jet3Label, hemis1Jet3Handle )
|
| 173 |
jet3 = (hemis1Jet3Handle.product())[0]
|
| 174 |
|
| 175 |
event.getByLabel (hemis1MuLabel, hemis1MuHandle)
|
| 176 |
wJetMu = hemis1MuHandle.product()
|
| 177 |
|
| 178 |
event.getByLabel (hemis1topLabel, hemis1topHandle)
|
| 179 |
topJetsh1 = hemis1topHandle.product()
|
| 180 |
|
| 181 |
event.getByLabel (hemis0topLabel, hemis0topHandle)
|
| 182 |
topJetsh0 = hemis0topHandle.product()
|
| 183 |
|
| 184 |
event.getByLabel (hemis0TopMassLabel, hemis0TopMassHandle)
|
| 185 |
topJetMassh0 = hemis0TopMassHandle.product()
|
| 186 |
|
| 187 |
event.getByLabel (hemis1TopMassLabel, hemis1TopMassHandle)
|
| 188 |
topJetMassh1 = hemis1TopMassHandle.product()
|
| 189 |
|
| 190 |
event.getByLabel (hemis0MinMassLabel, hemis0MinMassHandle)
|
| 191 |
topJetMinMassh0 = hemis0MinMassHandle.product()
|
| 192 |
|
| 193 |
event.getByLabel (hemis1MinMassLabel, hemis1MinMassHandle)
|
| 194 |
topJetMinMassh1 = hemis1MinMassHandle.product()
|
| 195 |
|
| 196 |
events.getByLabel ( hemis0NSubJetsLabel, hemis0NSubJetsHandle )
|
| 197 |
NJetsh0 = hemis0NSubJetsHandle.product()
|
| 198 |
|
| 199 |
events.getByLabel ( hemis1NSubJetsLabel , hemis1NSubJetsHandle )
|
| 200 |
NJetsh1 = hemis1NSubJetsHandle.product()
|
| 201 |
|
| 202 |
#---------------------------------------------------------------------------------------------------------------------#
|
| 203 |
|
| 204 |
#This section is used to plot overall Jet collection Information
|
| 205 |
#---------------------------------------------------------------------------------------------------------------------#
|
| 206 |
|
| 207 |
# Here we plot the number of jets, by adding hemisphere 0 + hemisphere 1.
|
| 208 |
nJets.Fill( len(wJetsh1) + len(wJetsh0))
|
| 209 |
|
| 210 |
for j in range(0,len(wbCandJets)) :
|
| 211 |
if wbCandJets[j].energy() > 30 :
|
| 212 |
esum += wbCandJets[j].energy()
|
| 213 |
if (wJetBDisc[j] > 3.3) :
|
| 214 |
nb.append(wbCandJets[j])
|
| 215 |
|
| 216 |
# If there are two b-jets, we find the various quantities.
|
| 217 |
if len(nb) == 2 :
|
| 218 |
dphin = abs(nb[0].phi()-nb[1].phi())
|
| 219 |
if dphin > math.pi :
|
| 220 |
dphin -= -2*math.pi
|
| 221 |
mass2bJets.Fill( (nb[0] + nb[1]).mass() )
|
| 222 |
dphi.Fill(dphin)
|
| 223 |
detan = abs(nb[0].eta()-nb[1].eta())
|
| 224 |
deta.Fill(detan)
|
| 225 |
dRn = math.sqrt(detan*detan + dphin*dphin)
|
| 226 |
dR.Fill(dRn)
|
| 227 |
|
| 228 |
# We plot the energy of the jet.
|
| 229 |
energyJets.Fill( esum )
|
| 230 |
|
| 231 |
#---------------------------------------------------------------------------------------------------------------------#
|
| 232 |
|
| 233 |
#This gives the conditional statements for type 1+2 and the two
|
| 234 |
#type 1+1 cases to be sure the number of jets exist
|
| 235 |
#---------------------------------------------------------------------------------------------------------------------#
|
| 236 |
|
| 237 |
njets12 = (len(wJetsh1) >= 2) and (len(wJetsh0) >= 1)
|
| 238 |
njets11b0 = ((len(topJetsh1) >= 1) and (len(wJetsh0) >= 1))
|
| 239 |
njets11b1 = ((len(topJetsh0) >= 1) and (len(wJetsh1) >= 1))
|
| 240 |
|
| 241 |
#---------------------------------------------------------------------------------------------------------------------#
|
| 242 |
|
| 243 |
#We look at type 1+1 W` first
|
| 244 |
#---------------------------------------------------------------------------------------------------------------------#
|
| 245 |
|
| 246 |
#If b is the highest pt jet
|
| 247 |
if njets11b0 :
|
| 248 |
hasBTag0 = wJetBDisch0[0] > 3.3
|
| 249 |
hasTopTag1 = topJetMassh1[0] > 140 and topJetMassh1[0] < 250 and topJetMinMassh1[0] > 50 and NJetsh1[0] > 2
|
| 250 |
highptbjh0 = wJetsh0[0].pt() > 350
|
| 251 |
highpttjh1 = topJetsh1[0].pt() > 350
|
| 252 |
if highptbjh0 and highpttjh1 and hasTopTag1:
|
| 253 |
pt1Bin = mistag.FindBin( wJetsh0[0].pt() )
|
| 254 |
mistagB11.Fill((topJetsh1[0]+wJetsh0[0]).mass(), mistag.GetBinContent(pt1Bin))
|
| 255 |
if highptbjh0 and highpttjh1 and hasTopTag1 and hasBTag0:
|
| 256 |
masswprimet11.Fill((topJetsh1[0]+wJetsh0[0]).mass())
|
| 257 |
|
| 258 |
#If t is the highest pt jet
|
| 259 |
if njets11b1 :
|
| 260 |
hasBTag1 = wJetBDisch1[0] > 3.3
|
| 261 |
hasTopTag0 = topJetMassh0[0] > 140 and topJetMassh0[0] < 250 and topJetMinMassh0[0] > 50 and NJetsh0[0] > 2
|
| 262 |
highpttjh0 = topJetsh0[0].pt() > 350
|
| 263 |
highptbjh1 = wJetsh1[0].pt() > 350
|
| 264 |
if highptbjh1 and highpttjh0 and hasTopTag0 and (not hasTopTag1):
|
| 265 |
pt1Bin = mistag.FindBin( wJetsh1[0].pt() )
|
| 266 |
mistagB11.Fill((topJetsh0[0]+wJetsh1[0]).mass(), mistag.GetBinContent(pt1Bin))
|
| 267 |
if highptbjh1 and highpttjh0 and hasTopTag0 and hasBTag1 and (not hasTopTag1):
|
| 268 |
masswprimet11.Fill((topJetsh0[0]+wJetsh1[0]).mass())
|
| 269 |
|
| 270 |
#---------------------------------------------------------------------------------------------------------------------#
|
| 271 |
|
| 272 |
#Now we look at type 1+2 W`
|
| 273 |
#---------------------------------------------------------------------------------------------------------------------#
|
| 274 |
|
| 275 |
if njets12 :
|
| 276 |
#Type2Top: Wjet pt > 200 mu < 0.4 and Bjet pt > 30
|
| 277 |
highptbjh0 = wJetsh0[0].pt() > 350
|
| 278 |
hemis1kincuts = (highptbjh0) and (wJetsh1[0].pt() > 200) and (wJetMu[0] < 0.4) and (wJetsh1[jet3].pt() > 30)
|
| 279 |
pairMass = (wJetsh1[jet3]+wJetsh1[0]).mass()
|
| 280 |
hasBTag0 = wJetBDisch0[0] > 3.3
|
| 281 |
#Toplike mass cuts
|
| 282 |
hash1Top = wJetsh1[0].mass() > 60 and wJetsh1[0].mass() < 130 and pairMass > 140 and pairMass < 250
|
| 283 |
hasBTag2 = wJetBDisch1[jet3] > 3.3
|
| 284 |
if highptbjh0 and hemis1kincuts :
|
| 285 |
nJetsh1.Fill(len(wJetsh1))
|
| 286 |
WJetmass.Fill(wJetsh1[0].mass())
|
| 287 |
if wJetsh1[0].mass() > 60 and wJetsh1[0].mass() < 130 :
|
| 288 |
WJetmasswc.Fill(wJetsh1[0].mass())
|
| 289 |
ptwJetsh1.Fill(wJetsh1[0].pt())
|
| 290 |
ptbJetsh1.Fill(wJetsh1[jet3].pt())
|
| 291 |
if hasBTag2 :
|
| 292 |
if hasBTag0:
|
| 293 |
ptbh0.Fill(wJetsh0[0].pt())
|
| 294 |
massbh0.Fill(wJetsh0[0].mass())
|
| 295 |
ptbJetsh1wt.Fill(wJetsh1[jet3].pt())
|
| 296 |
if hash1Top :
|
| 297 |
topCandmass.Fill(pairMass)
|
| 298 |
if hasBTag0:
|
| 299 |
masswprimewob.Fill((wJetsh1[jet3]+wJetsh1[0]+wJetsh0[0]).mass())
|
| 300 |
if hasBTag2 :
|
| 301 |
ptBin = mistag.FindBin( wJetsh0[0].pt() )
|
| 302 |
mistagB12.Fill((wJetsh1[jet3]+wJetsh1[0]+wJetsh0[0]).mass(), mistag.GetBinContent(ptBin))
|
| 303 |
topCandmasswb.Fill(pairMass)
|
| 304 |
if hasBTag0 :
|
| 305 |
masswprimewb.Fill((wJetsh1[jet3]+wJetsh1[0]+wJetsh0[0]).mass())
|
| 306 |
|
| 307 |
#---------------------------------------------------------------------------------------------------------------------#
|
| 308 |
|
| 309 |
f.cd()
|
| 310 |
f.Write()
|
| 311 |
f.Close()
|
