| 1 |
#! /usr/bin/env python
|
| 2 |
import os
|
| 3 |
import glob
|
| 4 |
import math
|
| 5 |
|
| 6 |
import ROOT
|
| 7 |
|
| 8 |
from optparse import OptionParser
|
| 9 |
|
| 10 |
parser = OptionParser()
|
| 11 |
|
| 12 |
|
| 13 |
############################################
|
| 14 |
# Job steering #
|
| 15 |
############################################
|
| 16 |
|
| 17 |
# Input inputFiles to use. This is in "glob" format, so you can use wildcards.
|
| 18 |
# If you get a "cannot find file" type of error, be sure to use "\*" instead
|
| 19 |
# of "*" to make sure you don't confuse the shell.
|
| 20 |
parser.add_option('--inputFiles', metavar='F', type='string', action='store',
|
| 21 |
default = "",
|
| 22 |
dest='inputFiles',
|
| 23 |
help='Input files')
|
| 24 |
|
| 25 |
parser.add_option('--txtfiles', metavar='F', type='string', action='store',
|
| 26 |
default = "",
|
| 27 |
dest='txtfiles',
|
| 28 |
help='Input txt files')
|
| 29 |
|
| 30 |
parser.add_option("--onDcache", action='store_true',
|
| 31 |
default=True,
|
| 32 |
dest="onDcache",
|
| 33 |
help="onDcache(1), onDcache(0)")
|
| 34 |
|
| 35 |
# Output name to use.
|
| 36 |
parser.add_option('--outputFile', metavar='F', type='string', action='store',
|
| 37 |
default='shyft_fwlite.root',
|
| 38 |
dest='outputFile',
|
| 39 |
help='output file name')
|
| 40 |
|
| 41 |
# Using MC info or not. For MC, truth information is accessed.
|
| 42 |
parser.add_option('--doMC', metavar='F', action='store_true',
|
| 43 |
default=False,
|
| 44 |
dest='doMC',
|
| 45 |
help='Check MC Information')
|
| 46 |
|
| 47 |
# Which lepton type to use
|
| 48 |
parser.add_option('--lepType', metavar='F', type='int', action='store',
|
| 49 |
default=0,
|
| 50 |
dest='lepType',
|
| 51 |
help='Lepton type. Options are 0 = muons, 1 = electrons')
|
| 52 |
|
| 53 |
# Invert MET cut?
|
| 54 |
parser.add_option('--invertMET', action='store_true',
|
| 55 |
default=False,
|
| 56 |
dest='invertMET',
|
| 57 |
help='Invert MET cut')
|
| 58 |
|
| 59 |
# Remove MET cut?
|
| 60 |
parser.add_option('--relaxMET', action='store_true',
|
| 61 |
default=False,
|
| 62 |
dest='relaxMET',
|
| 63 |
help='Relax MET cut')
|
| 64 |
|
| 65 |
# Invert PF isolation cut?
|
| 66 |
parser.add_option('--invertPFIso', action='store_true',
|
| 67 |
default=False,
|
| 68 |
dest='invertPFIso',
|
| 69 |
help='Invert PF isolation cut')
|
| 70 |
|
| 71 |
|
| 72 |
(options, args) = parser.parse_args()
|
| 73 |
|
| 74 |
argv = []
|
| 75 |
|
| 76 |
# Import everything from ROOT
|
| 77 |
import ROOT
|
| 78 |
ROOT.gROOT.Macro("rootlogon.C")
|
| 79 |
|
| 80 |
# Import stuff from FWLite
|
| 81 |
import sys
|
| 82 |
from DataFormats.FWLite import Events, Handle
|
| 83 |
|
| 84 |
#infile = open( options.inputFiles )
|
| 85 |
#infileStr = infile.read().rstrip()
|
| 86 |
|
| 87 |
#print 'Getting files from this dir: ' + infileStr
|
| 88 |
|
| 89 |
# Get the file list.
|
| 90 |
#files = glob.glob( infileStr )
|
| 91 |
|
| 92 |
# Get the file list.
|
| 93 |
if options.inputFiles:
|
| 94 |
files = glob.glob( options.files )
|
| 95 |
print 'getting files', files
|
| 96 |
elif options.txtfiles:
|
| 97 |
files = []
|
| 98 |
with open(options.txtfiles, 'r') as input_:
|
| 99 |
for line in input_:
|
| 100 |
files.append(line.strip())
|
| 101 |
else:
|
| 102 |
files = []
|
| 103 |
|
| 104 |
print 'getting files: ', files
|
| 105 |
|
| 106 |
if options.onDcache:
|
| 107 |
files = ["dcap://" + x for x in files]
|
| 108 |
#print 'new files', *files, sep='\n'
|
| 109 |
#print('new files', files[0], files[1], ..., sep='\n')
|
| 110 |
|
| 111 |
fname = options.txtfiles
|
| 112 |
fileN = fname[fname.rfind('/')+1:]
|
| 113 |
|
| 114 |
print files
|
| 115 |
|
| 116 |
|
| 117 |
# Create the output file.
|
| 118 |
f = ROOT.TFile(options.outputFile, "recreate")
|
| 119 |
f.cd()
|
| 120 |
|
| 121 |
|
| 122 |
# Make histograms
|
| 123 |
print "Creating histograms"
|
| 124 |
secvtxMassHist = ROOT.TH1F('secvtxMassHist', "Secondary Vertex Mass", 150, 0., 5.0)
|
| 125 |
secvtxMassHistB = ROOT.TH1F('secvtxMassHistB', "Secondary Vertex Mass, b jets", 150, 0., 5.0)
|
| 126 |
secvtxMassHistC = ROOT.TH1F('secvtxMassHistC', "Secondary Vertex Mass, c jets", 150, 0., 5.0)
|
| 127 |
secvtxMassHistL = ROOT.TH1F('secvtxMassHistL', "Secondary Vertex Mass, udsg jets", 150, 0., 5.0)
|
| 128 |
metVsIso = ROOT.TH2F('metVsIso', 'MET Versus PFIsolation', 15, 0., 150., 100, 0., 2.5)
|
| 129 |
|
| 130 |
jetPtHist = ROOT.TH1F('jetPtHist', 'Jet p_{T}', 150, 0., 600.)
|
| 131 |
m3Hist = ROOT.TH1F('m3Hist', 'M3 Histogram', 150, 0., 600.)
|
| 132 |
#Jet Multiplicity
|
| 133 |
nJetHist = ROOT.TH1F('nJetHist', 'Jet Multiplicity', 7, 3.5, 10.5)
|
| 134 |
#Muon pt
|
| 135 |
muPtHist = ROOT.TH1F('muPtHist', 'Muon pt', 75, 0., 300.)
|
| 136 |
#Muon eta
|
| 137 |
muEtaHist = ROOT.TH1F('muEtaHist', 'Muon eta', 40, -2.2, 2.2)
|
| 138 |
#MET
|
| 139 |
metHist = ROOT.TH1F('metHist', 'MET', 50, 0., 200.)
|
| 140 |
#HT
|
| 141 |
#this is a scalar sum of pt in the events -- we will choose whether of jets, leptons and jets, etc.
|
| 142 |
htHist = ROOT.TH1F('htHist', 'HT', 250, 0., 1000.)
|
| 143 |
#transverse mass
|
| 144 |
transMHist = ROOT.TH1F('transMHist', 'Transverse mass', 75, 0., 300)
|
| 145 |
#N btags
|
| 146 |
nBTagHist = ROOT.TH1F('nBTagHist', 'Btag jet multiplicity', 4, 0.5, 4.5)
|
| 147 |
|
| 148 |
############################################
|
| 149 |
# Physics level parameters for systematics #
|
| 150 |
############################################
|
| 151 |
|
| 152 |
# Kinematic cuts:
|
| 153 |
jetPtMin = 30.0
|
| 154 |
leadJetPtMin = 30.0
|
| 155 |
isoMax = 0.12
|
| 156 |
ssvheCut = 1.74
|
| 157 |
minJets = 4
|
| 158 |
|
| 159 |
if options.lepType == 0 :
|
| 160 |
muonPtMin = 45.0
|
| 161 |
electronPtMin = 20.0
|
| 162 |
metMin = 20.0
|
| 163 |
lepStr = 'Mu'
|
| 164 |
else:
|
| 165 |
muonPtMin = 20.0
|
| 166 |
electronPtMin = 35.0
|
| 167 |
metMin = 20.0
|
| 168 |
lepStr = 'Ele'
|
| 169 |
|
| 170 |
|
| 171 |
|
| 172 |
|
| 173 |
events = Events (files)
|
| 174 |
|
| 175 |
# Make the entirety of the handles required for the
|
| 176 |
# analysis.
|
| 177 |
postfix = ""
|
| 178 |
if options.invertPFIso :
|
| 179 |
postfix = "Loose"
|
| 180 |
|
| 181 |
puHandle = Handle( "std::vector<float>" )
|
| 182 |
puLabel = ( "PUNtupleDumper", "PUweightNominalUpDown" )
|
| 183 |
|
| 184 |
|
| 185 |
jetPtHandle = Handle( "std::vector<float>" )
|
| 186 |
jetPtLabel = ( "pfShyftTupleJets" + lepStr + postfix, "pt" )
|
| 187 |
jetEtaHandle = Handle( "std::vector<float>" )
|
| 188 |
jetEtaLabel = ( "pfShyftTupleJets" + lepStr + postfix, "eta" )
|
| 189 |
jetPhiHandle = Handle( "std::vector<float>" )
|
| 190 |
jetPhiLabel = ( "pfShyftTupleJets" + lepStr + postfix, "phi" )
|
| 191 |
jetMassHandle = Handle( "std::vector<float>" )
|
| 192 |
jetMassLabel = ( "pfShyftTupleJets" + lepStr + postfix, "mass" )
|
| 193 |
jetSecvtxMassHandle = Handle( "std::vector<float>" )
|
| 194 |
jetSecvtxMassLabel = ( "pfShyftTupleJets" + lepStr + postfix, "secvtxMass" )
|
| 195 |
jetSSVHEHandle = Handle( "std::vector<float>" )
|
| 196 |
jetSSVHELabel = ( "pfShyftTupleJets" + lepStr + postfix, "ssvhe" )
|
| 197 |
jetFlavorHandle = Handle( "std::vector<float>" )
|
| 198 |
jetFlavorLabel = ( "pfShyftTupleJets" + lepStr + postfix, "flavor" )
|
| 199 |
|
| 200 |
muonPtHandle = Handle( "std::vector<float>" )
|
| 201 |
muonPtLabel = ( "pfShyftTupleMuons"+ postfix, "pt" )
|
| 202 |
muonEtaHandle = Handle( "std::vector<float>" )
|
| 203 |
muonEtaLabel = ( "pfShyftTupleMuons"+ postfix, "eta" )
|
| 204 |
muonPhiHandle = Handle( "std::vector<float>" )
|
| 205 |
muonPhiLabel = ( "pfShyftTupleMuons"+ postfix, "phi" )
|
| 206 |
muonNhIsoHandle = Handle( "std::vector<float>" )
|
| 207 |
muonNhIsoLabel = ( "pfShyftTupleMuons"+ postfix, "nhIso" )
|
| 208 |
muonChIsoHandle = Handle( "std::vector<float>" )
|
| 209 |
muonChIsoLabel = ( "pfShyftTupleMuons"+ postfix, "chIso" )
|
| 210 |
muonPhIsoHandle = Handle( "std::vector<float>" )
|
| 211 |
muonPhIsoLabel = ( "pfShyftTupleMuons"+ postfix, "phIso" )
|
| 212 |
muonPuIsoHandle = Handle( "std::vector<float>" )
|
| 213 |
muonPuIsoLabel = ( "pfShyftTupleMuons"+ postfix, "puIso" )
|
| 214 |
|
| 215 |
electronPtHandle = Handle( "std::vector<float>" )
|
| 216 |
electronPtLabel = ( "pfShyftTupleElectrons"+ postfix, "pt" )
|
| 217 |
electronEtaHandle = Handle( "std::vector<float>" )
|
| 218 |
electronEtaLabel = ( "pfShyftTupleElectrons"+ postfix, "eta" )
|
| 219 |
electronPhiHandle = Handle( "std::vector<float>" )
|
| 220 |
electronPhiLabel = ( "pfShyftTupleElectrons"+ postfix, "phi" )
|
| 221 |
electronNhIsoHandle = Handle( "std::vector<float>" )
|
| 222 |
electronNhIsoLabel = ( "pfShyftTupleElectrons"+ postfix, "nhIso" )
|
| 223 |
electronChIsoHandle = Handle( "std::vector<float>" )
|
| 224 |
electronChIsoLabel = ( "pfShyftTupleElectrons"+ postfix, "chIso" )
|
| 225 |
electronPhIsoHandle = Handle( "std::vector<float>" )
|
| 226 |
electronPhIsoLabel = ( "pfShyftTupleElectrons"+ postfix, "phIso" )
|
| 227 |
electronPuIsoHandle = Handle( "std::vector<float>" )
|
| 228 |
electronPuIsoLabel = ( "pfShyftTupleElectrons"+ postfix, "puIso" )
|
| 229 |
|
| 230 |
|
| 231 |
metHandle = Handle( "std::vector<float>" )
|
| 232 |
metLabel = ("pfShyftTupleMET" + lepStr + postfix, "pt" )
|
| 233 |
metPhiHandle = Handle( "std::vector<float>" )
|
| 234 |
metPhiLabel = ("pfShyftTupleMET" + lepStr + postfix, "phi" )
|
| 235 |
|
| 236 |
|
| 237 |
# Keep some timing information
|
| 238 |
nEventsAnalyzed = 0
|
| 239 |
nEventsPassed4Jets = 0
|
| 240 |
nEventsPassed1Tag = 0
|
| 241 |
timer = ROOT.TStopwatch()
|
| 242 |
timer.Start()
|
| 243 |
|
| 244 |
pairs = []
|
| 245 |
|
| 246 |
# loop over events
|
| 247 |
count = 0
|
| 248 |
ntotal = events.size()
|
| 249 |
percentDone = 0.0
|
| 250 |
ipercentDone = 0
|
| 251 |
ipercentDoneLast = -1
|
| 252 |
print "Start looping"
|
| 253 |
for event in events:
|
| 254 |
nEventsAnalyzed += 1
|
| 255 |
ipercentDone = int(percentDone)
|
| 256 |
if ipercentDone != ipercentDoneLast :
|
| 257 |
ipercentDoneLast = ipercentDone
|
| 258 |
print 'Processing {0:10.0f}/{1:10.0f} : {2:5.0f}%'.format(
|
| 259 |
count, ntotal, ipercentDone )
|
| 260 |
count = count + 1
|
| 261 |
percentDone = float(count) / float(ntotal) * 100.0
|
| 262 |
|
| 263 |
|
| 264 |
################################################
|
| 265 |
# Retrieve the jet four-vector
|
| 266 |
# ------------------------------------
|
| 267 |
# The jet 4-vectors are large and hence
|
| 268 |
# take a long time to read out. If you don't
|
| 269 |
# need the other products (eta,phi,mass of jet)
|
| 270 |
# then don't read them out.
|
| 271 |
################################################
|
| 272 |
|
| 273 |
event.getByLabel( jetPtLabel, jetPtHandle )
|
| 274 |
if not jetPtHandle.isValid():
|
| 275 |
jetPts = None
|
| 276 |
else :
|
| 277 |
jetPts = jetPtHandle.product()
|
| 278 |
|
| 279 |
if jetPts is None :
|
| 280 |
continue
|
| 281 |
event.getByLabel( jetEtaLabel, jetEtaHandle )
|
| 282 |
jetEtas = jetEtaHandle.product()
|
| 283 |
event.getByLabel( jetPhiLabel, jetPhiHandle )
|
| 284 |
jetPhis = jetPhiHandle.product()
|
| 285 |
event.getByLabel( jetMassLabel, jetMassHandle )
|
| 286 |
jetMasses = jetMassHandle.product()
|
| 287 |
|
| 288 |
|
| 289 |
# Find the njet bin we're in
|
| 290 |
njets = -1
|
| 291 |
event.getByLabel( jetPtLabel, jetPtHandle )
|
| 292 |
jetPts = jetPtHandle.product()
|
| 293 |
njets = 0
|
| 294 |
for ijet in range( 0, len( jetPts ) ) :
|
| 295 |
if jetPts[ijet] > 30.0 :
|
| 296 |
njets += 1
|
| 297 |
|
| 298 |
# We're not interested in <=4 jets
|
| 299 |
if njets < minJets :
|
| 300 |
continue
|
| 301 |
nEventsPassed4Jets = nEventsPassed4Jets + 1
|
| 302 |
|
| 303 |
|
| 304 |
|
| 305 |
################################################
|
| 306 |
# Retrieve the jet vertex mass and plot the
|
| 307 |
# secondary vertex mass for tagged jets.
|
| 308 |
################################################
|
| 309 |
event.getByLabel (jetSSVHELabel, jetSSVHEHandle)
|
| 310 |
jetSSVHEs = jetSSVHEHandle.product()
|
| 311 |
|
| 312 |
# Now loop over the jets, and count tags
|
| 313 |
ntags = 0
|
| 314 |
for ijet in range(0,len(jetPts) ) :
|
| 315 |
jetPt = jetPts[ijet]
|
| 316 |
if jetPt < 30.0 :
|
| 317 |
continue
|
| 318 |
#jetEta = jetEtas[ijet]
|
| 319 |
#jetPhi = jetPhis[ijet]
|
| 320 |
#jetMass = jetMasses[ijet]
|
| 321 |
jetSSVHE = jetSSVHEs[ijet]
|
| 322 |
# plot secondary vertex mass for tagged jets
|
| 323 |
if jetSSVHE >= ssvheCut :
|
| 324 |
ntags = ntags + 1
|
| 325 |
if ntags > 0:
|
| 326 |
nEventsPassed1Tag = nEventsPassed1Tag + 1
|
| 327 |
pairs.append( [event.object().id().run(),
|
| 328 |
event.object().id().luminosityBlock(),
|
| 329 |
event.object().id().event(),
|
| 330 |
njets,
|
| 331 |
ntags] )
|
| 332 |
else :
|
| 333 |
continue
|
| 334 |
|
| 335 |
|
| 336 |
################################################
|
| 337 |
# Require exactly one lepton (e or mu)
|
| 338 |
# ------------------------------------
|
| 339 |
# Our ntuples have both muon and electron
|
| 340 |
# events, and hence we must select events
|
| 341 |
# based on one or the other type.
|
| 342 |
# To accomplish this we check the products
|
| 343 |
# for the type we're currently plotting
|
| 344 |
# (Mu or Ele), and check if the product is
|
| 345 |
# present.
|
| 346 |
################################################
|
| 347 |
muonPts = None
|
| 348 |
electronPts = None
|
| 349 |
muonEtas = None
|
| 350 |
muonPhis = None
|
| 351 |
if options.lepType == 0 :
|
| 352 |
event.getByLabel (muonPtLabel, muonPtHandle)
|
| 353 |
if not muonPtHandle.isValid():
|
| 354 |
muonPts = None
|
| 355 |
else :
|
| 356 |
muonPts = muonPtHandle.product()
|
| 357 |
event.getByLabel (muonEtaLabel, muonEtaHandle)
|
| 358 |
muonEtas = muonEtaHandle.product()
|
| 359 |
event.getByLabel( muonPhiLabel, muonPhiHandle )
|
| 360 |
muonPhis = muonPhiHandle.product()
|
| 361 |
|
| 362 |
elif options.lepType == 1 :
|
| 363 |
event.getByLabel (electronPtLabel, electronPtHandle)
|
| 364 |
if not electronPtHandle.isValid():
|
| 365 |
electronPts = None
|
| 366 |
else :
|
| 367 |
electronPts = electronPtHandle.product()
|
| 368 |
|
| 369 |
# If neither muons nor electrons are found, skip
|
| 370 |
if muonPts is None and electronPts is None :
|
| 371 |
continue
|
| 372 |
# If we are looking for muons but none are found, skip
|
| 373 |
if options.lepType == 0 and muonPts is None :
|
| 374 |
continue
|
| 375 |
# If we are looking for electrons but none are found, skip
|
| 376 |
if options.lepType == 1 and electronPts is None :
|
| 377 |
continue
|
| 378 |
|
| 379 |
# Now get the MET
|
| 380 |
event.getByLabel( metLabel, metHandle )
|
| 381 |
metRaw = metHandle.product()[0]
|
| 382 |
event.getByLabel( metPhiLabel, metPhiHandle )
|
| 383 |
metPhi = metPhiHandle.product()[0]
|
| 384 |
|
| 385 |
# Now get the PF isolation
|
| 386 |
lepIso = -1.0
|
| 387 |
if options.lepType == 0 and muonPts is not None :
|
| 388 |
events.getByLabel( muonNhIsoLabel, muonNhIsoHandle )
|
| 389 |
events.getByLabel( muonChIsoLabel, muonChIsoHandle )
|
| 390 |
events.getByLabel( muonPhIsoLabel, muonPhIsoHandle )
|
| 391 |
events.getByLabel( muonPuIsoLabel, muonPuIsoHandle )
|
| 392 |
nhIso = muonNhIsoHandle.product()[0]
|
| 393 |
chIso = muonChIsoHandle.product()[0]
|
| 394 |
phIso = muonPhIsoHandle.product()[0]
|
| 395 |
puIso = muonPuIsoHandle.product()[0]
|
| 396 |
lepIso = (chIso + max(0.0, nhIso + phIso - 0.5*puIso)) / muonPts[0]
|
| 397 |
if options.lepType == 1 and electronPts is not None :
|
| 398 |
events.getByLabel( electronNhIsoLabel, electronNhIsoHandle )
|
| 399 |
events.getByLabel( electronChIsoLabel, electronChIsoHandle )
|
| 400 |
events.getByLabel( electronPhIsoLabel, electronPhIsoHandle )
|
| 401 |
events.getByLabel( electronPuIsoLabel, electronPuIsoHandle )
|
| 402 |
nhIso = electronNhIsoHandle.product()[0]
|
| 403 |
chIso = electronChIsoHandle.product()[0]
|
| 404 |
phIso = electronPhIsoHandle.product()[0]
|
| 405 |
puIso = electronPuIsoHandle.product()[0]
|
| 406 |
lepIso = (chIso + max(0.0, nhIso + phIso - 0.5*puIso)) / electronPts[0]
|
| 407 |
|
| 408 |
|
| 409 |
# Make a plot of the MET versus ISO for normalization purposes
|
| 410 |
metVsIso.Fill( metRaw, lepIso )
|
| 411 |
|
| 412 |
# If the ISO is higher than our cut, skip, unless we want it inverted
|
| 413 |
if not options.invertPFIso :
|
| 414 |
if lepIso > isoMax :
|
| 415 |
continue
|
| 416 |
else :
|
| 417 |
if lepIso < isoMax :
|
| 418 |
continue
|
| 419 |
|
| 420 |
# Fill histograms which do not have MET cut
|
| 421 |
metHist.Fill(metRaw)
|
| 422 |
|
| 423 |
# Now compute transverse mass
|
| 424 |
muVec = ROOT.TLorentzVector()
|
| 425 |
muVec.SetPtEtaPhiM(muonPts[0],0,muonPhis[0],0)
|
| 426 |
metVec = ROOT.TLorentzVector()
|
| 427 |
metVec.SetPtEtaPhiM(metRaw,0,metPhi,0)
|
| 428 |
massVec = muVec + metVec
|
| 429 |
transMHist.Fill( massVec.M())
|
| 430 |
|
| 431 |
# If the MET is lower than our cut, skip, unless we want it inverted
|
| 432 |
if options.invertMET :
|
| 433 |
if metRaw > metMin :
|
| 434 |
continue
|
| 435 |
elif options.relaxMET :
|
| 436 |
if metRaw < 0 :
|
| 437 |
continue
|
| 438 |
else :
|
| 439 |
if metRaw < metMin :
|
| 440 |
continue
|
| 441 |
|
| 442 |
event.getByLabel (jetSecvtxMassLabel, jetSecvtxMassHandle)
|
| 443 |
jetSecvtxMasses = jetSecvtxMassHandle.product()
|
| 444 |
if options.doMC :
|
| 445 |
event.getByLabel( jetFlavorLabel, jetFlavorHandle )
|
| 446 |
jetFlavors = jetFlavorHandle.product()
|
| 447 |
# Now loop over the jets, and store the secondary vertex mass.
|
| 448 |
ntags = 0
|
| 449 |
ht_jets = 0
|
| 450 |
jets_p4 = []
|
| 451 |
for ijet in range(0,len(jetPts) ) :
|
| 452 |
jetPt = jetPts[ijet]
|
| 453 |
if jetPt < 30.0 :
|
| 454 |
continue
|
| 455 |
ijetP4 = ROOT.TLorentzVector()
|
| 456 |
ijetP4.SetPtEtaPhiM( jetPts[ijet], jetEtas[ijet], jetPhis[ijet], jetMasses[ijet] )
|
| 457 |
jets_p4.append( ijetP4 )
|
| 458 |
jetSSVHE = jetSSVHEs[ijet]
|
| 459 |
jetSecvtxMass = jetSecvtxMasses[ijet]
|
| 460 |
jetPtHist.Fill( jetPt )
|
| 461 |
ht_jets += jetPt
|
| 462 |
# plot secondary vertex mass for tagged jets
|
| 463 |
if jetSSVHE >= ssvheCut :
|
| 464 |
ntags = ntags + 1
|
| 465 |
secvtxMassHist.Fill( jetSecvtxMass )
|
| 466 |
if options.doMC :
|
| 467 |
if abs(jetFlavors[ijet]) == 5 :
|
| 468 |
secvtxMassHistB.Fill( jetSecvtxMass )
|
| 469 |
elif abs(jetFlavors[ijet]) == 4 :
|
| 470 |
secvtxMassHistC.Fill( jetSecvtxMass )
|
| 471 |
else :
|
| 472 |
secvtxMassHistL.Fill( jetSecvtxMass )
|
| 473 |
nBTagHist.Fill( ntags )
|
| 474 |
|
| 475 |
# Now compute m3
|
| 476 |
maxPt = -1.0
|
| 477 |
m3 = -1.0
|
| 478 |
for ijet in range(0, len(jets_p4) ) :
|
| 479 |
for jjet in range(ijet + 1, len(jets_p4) ) :
|
| 480 |
for kjet in range(jjet + 1, len(jets_p4) ) :
|
| 481 |
sumP4 = jets_p4[ijet] + jets_p4[jjet] + jets_p4[kjet]
|
| 482 |
if sumP4.Perp() > maxPt :
|
| 483 |
maxPt = sumP4.Perp()
|
| 484 |
m3 = sumP4.M()
|
| 485 |
if maxPt > 0.0 :
|
| 486 |
m3Hist.Fill( m3 )
|
| 487 |
|
| 488 |
# fill remaining histograms
|
| 489 |
nJetHist.Fill(njets)
|
| 490 |
|
| 491 |
htHist.Fill(ht_jets)
|
| 492 |
|
| 493 |
for imu in range( 0, len( muonPts ) ) :
|
| 494 |
muPtHist.Fill(muonPts[imu])
|
| 495 |
for imu in range( 0, len( muonEtas) ) :
|
| 496 |
muEtaHist.Fill(muonEtas[imu])
|
| 497 |
|
| 498 |
# Stop our timer
|
| 499 |
timer.Stop()
|
| 500 |
|
| 501 |
# Print out our timing information
|
| 502 |
rtime = timer.RealTime(); # Real time (or "wall time")
|
| 503 |
ctime = timer.CpuTime(); # CPU time
|
| 504 |
print("Analyzed events: {0:6d}").format(nEventsAnalyzed)
|
| 505 |
print(">=4 jet events : {0:6d}").format(nEventsPassed4Jets)
|
| 506 |
print(">=1 tag events : {0:6d}").format(nEventsPassed1Tag)
|
| 507 |
print("RealTime={0:6.2f} seconds, CpuTime={1:6.2f} seconds").format(rtime,ctime)
|
| 508 |
print("{0:4.2f} events / RealTime second .").format( nEventsAnalyzed/rtime)
|
| 509 |
print("{0:4.2f} events / CpuTime second .").format( nEventsAnalyzed/ctime)
|
| 510 |
|
| 511 |
|
| 512 |
|
| 513 |
f.cd()
|
| 514 |
f.Write()
|
| 515 |
f.Close()
|
