| 1 |
dittmer |
1.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 |
dittmer |
1.2 |
# Remove MET cut?
|
| 60 |
|
|
parser.add_option('--relaxMET', action='store_true',
|
| 61 |
|
|
default=False,
|
| 62 |
|
|
dest='relaxMET',
|
| 63 |
|
|
help='Relax MET cut')
|
| 64 |
|
|
|
| 65 |
dittmer |
1.1 |
# 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 |
dittmer |
1.6 |
muEtaHist = ROOT.TH1F('muEtaHist', 'Muon eta', 40, -2.2, 2.2)
|
| 138 |
dittmer |
1.1 |
#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 |
dittmer |
1.5 |
isoMax = 0.12
|
| 156 |
dittmer |
1.1 |
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 |
dittmer |
1.3 |
# 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 |
dittmer |
1.1 |
# If the MET is lower than our cut, skip, unless we want it inverted
|
| 432 |
dittmer |
1.2 |
if options.invertMET :
|
| 433 |
|
|
if metRaw > metMin :
|
| 434 |
|
|
continue
|
| 435 |
|
|
elif options.relaxMET :
|
| 436 |
|
|
if metRaw < 0 :
|
| 437 |
dittmer |
1.1 |
continue
|
| 438 |
|
|
else :
|
| 439 |
dittmer |
1.4 |
if metRaw < metMin :
|
| 440 |
dittmer |
1.1 |
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 |
dittmer |
1.3 |
|
| 491 |
dittmer |
1.1 |
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()
|
