| 1 |
#! /usr/bin/env python2
|
| 2 |
# -*- coding: utf-8 -*-
|
| 3 |
|
| 4 |
import ROOT
|
| 5 |
import argparse
|
| 6 |
import ConfigParser
|
| 7 |
import os.path
|
| 8 |
|
| 9 |
# to use user defined help message, sys.arv has to be sent to python and not
|
| 10 |
# to TApplication
|
| 11 |
ROOT.PyConfig.IgnoreCommandLineOptions = True
|
| 12 |
|
| 13 |
# use tdr style defined in own macro
|
| 14 |
import Styles
|
| 15 |
Styles.tdrStyle()
|
| 16 |
|
| 17 |
# this is the configuration file for plots including label, etc.
|
| 18 |
axisConf = ConfigParser.SafeConfigParser()
|
| 19 |
axisConf.read("axis.cfg")
|
| 20 |
|
| 21 |
# global variables:
|
| 22 |
integratedLuminosity = 19.3 #fb
|
| 23 |
|
| 24 |
def readTree( filename, treename = "susyTree" ):
|
| 25 |
"""
|
| 26 |
filename: name of file containing the tree
|
| 27 |
treename: name of the tree
|
| 28 |
returns: TTree Object
|
| 29 |
"""
|
| 30 |
if not os.path.isfile(filename):
|
| 31 |
print( "File %s does not exist"%filename)
|
| 32 |
tree = ROOT.TChain( treename )
|
| 33 |
tree.AddFile( filename )
|
| 34 |
return tree
|
| 35 |
|
| 36 |
def createHistoFromTree2D(tree, variable, weight="", nBin=50):
|
| 37 |
"""
|
| 38 |
tree: tree to create histo from
|
| 39 |
variable: variable to plot (must be a branch of the tree)
|
| 40 |
weight: weights to apply (e.g. "var1*(var2 > 15)" will use weights from var1 and cut on var2 > 15
|
| 41 |
nBins, firstBin, lastBin: number of bins, first bin and last bin (same as in TH1F constructor)
|
| 42 |
nEvents: number of events to process (-1 = all)
|
| 43 |
returns: histogram
|
| 44 |
"""
|
| 45 |
from random import randint
|
| 46 |
from sys import maxint
|
| 47 |
#make a random name you could give something meaningfull here,
|
| 48 |
#but that would make this less readable
|
| 49 |
name = "%x"%(randint(0, maxint))
|
| 50 |
|
| 51 |
# automatic binning if lastbin < firstbin
|
| 52 |
result = ROOT.TH2F(name, variable, nBin, 0, -1, nBin, 0, -1)
|
| 53 |
result.Sumw2()
|
| 54 |
tree.Draw("%s>>%s"%(variable, name), weight, "goff,colz")
|
| 55 |
return result
|
| 56 |
|
| 57 |
def createHistoFromTree(tree, variable, weight="", nBins=100, firstBin=None, lastBin=None, nEvents=-1):
|
| 58 |
"""
|
| 59 |
tree: tree to create histo from
|
| 60 |
variable: variable to plot (must be a branch of the tree)
|
| 61 |
weight: weights to apply (e.g. "var1*(var2 > 15)" will use weights from var1 and cut on var2 > 15
|
| 62 |
nBins, firstBin, lastBin: number of bins, first bin and last bin (same as in TH1F constructor)
|
| 63 |
nEvents: number of events to process (-1 = all)
|
| 64 |
returns: histogram
|
| 65 |
"""
|
| 66 |
if ":" in variable:
|
| 67 |
return createHistoFromTree2D(tree, variable, weight="")
|
| 68 |
from ROOT import TH1F
|
| 69 |
from random import randint
|
| 70 |
from sys import maxint
|
| 71 |
if nEvents < 0:
|
| 72 |
nEvents = maxint
|
| 73 |
if firstBin == None:
|
| 74 |
firstBin = tree.GetMinimum( variable )
|
| 75 |
if lastBin == None:
|
| 76 |
lastBin = tree.GetMaximum( variable )
|
| 77 |
#make a random name you could give something meaningfull here,
|
| 78 |
#but that would make this less readable
|
| 79 |
name = "%x"%(randint(0, maxint))
|
| 80 |
result = TH1F(name, variable, nBins, firstBin, lastBin)
|
| 81 |
result.Sumw2()
|
| 82 |
tree.Draw("%s>>%s"%(variable, name), weight, "goff", nEvents)
|
| 83 |
return result
|
| 84 |
|
| 85 |
def plot2D( tree, plot, cut, save=False ):
|
| 86 |
# read axis config file
|
| 87 |
var = plot.split(":")
|
| 88 |
try:
|
| 89 |
xlabel = axisConf.get( var[0], "label" )
|
| 90 |
except:
|
| 91 |
xlabel = ""
|
| 92 |
try:
|
| 93 |
ylabel = axisConf.get( var[1], "label" )
|
| 94 |
except:
|
| 95 |
ylabel = ""
|
| 96 |
try:
|
| 97 |
xunit = axisConf.get( var[0], "unit" )
|
| 98 |
except:
|
| 99 |
xunit = ""
|
| 100 |
try:
|
| 101 |
yunit = axisConf.get( var[1], "unit" )
|
| 102 |
except:
|
| 103 |
yunit = ""
|
| 104 |
if xunit != "":
|
| 105 |
xunit = " ["+xunit+"]"
|
| 106 |
if yunit != "":
|
| 107 |
yunit = " ["+yunit+"]"
|
| 108 |
|
| 109 |
# modify histo
|
| 110 |
histo = createHistoFromTree2D( tree, plot, cut )
|
| 111 |
histo.SetTitle(";%s%s;%s%s"%(ylabel,yunit,xlabel,xunit))
|
| 112 |
|
| 113 |
# draw histo
|
| 114 |
canvas = ROOT.TCanvas()
|
| 115 |
canvas.cd()
|
| 116 |
histo.Draw("colz")
|
| 117 |
|
| 118 |
# draw information
|
| 119 |
cutText = ROOT.TPaveText(.5,.8,.9,.9,"ndc")
|
| 120 |
cutText.SetBorderSize(0)
|
| 121 |
cutText.SetFillColor(0)
|
| 122 |
for line in cut.split("&&"):
|
| 123 |
cutText.AddText( line )
|
| 124 |
cutText.Draw()
|
| 125 |
|
| 126 |
topInfo = ROOT.TPaveText(.1,.95,1,1,"ndc")
|
| 127 |
topInfo.SetBorderSize(0)
|
| 128 |
topInfo.SetFillColor(0)
|
| 129 |
topInfo.AddText("Work in progress, %sfb^{-1}, #sqrt{s}=8TeV"%"?")
|
| 130 |
topInfo.Draw()
|
| 131 |
|
| 132 |
if save:
|
| 133 |
canvas.SaveAs("%.pdf"%variable)
|
| 134 |
else:
|
| 135 |
raw_input()
|
| 136 |
|
| 137 |
def fillWeights( tree, jetCut, photonCut ):
|
| 138 |
|
| 139 |
# this is the binning with witch the fo reweighting will be done and has to be studied
|
| 140 |
nBins = 50
|
| 141 |
xMin = 80
|
| 142 |
xMax = 300
|
| 143 |
h_photonPt = createHistoFromTree( tree, "photon[0].pt", photonCut, nBins, xMin, xMax )
|
| 144 |
h_jetPt = createHistoFromTree( tree, "photon[0].pt", jetCut, nBins, xMin, xMax )
|
| 145 |
|
| 146 |
# h_jetPt is now histogram with w^{-1} = h_jetPt / h_photonPt
|
| 147 |
h_jetPt.Divide( h_photonPt )
|
| 148 |
import array
|
| 149 |
weight = array.array( "f", [0] )
|
| 150 |
tree.Branch( "weight2", weight, "weight2/F" )
|
| 151 |
|
| 152 |
for event in tree:
|
| 153 |
pt = tree.photon.at(0).pt
|
| 154 |
#print "p_T = %i GeV"%pt
|
| 155 |
weight = h_jetPt.GetBinContent( h_jetPt.FindBin( pt ) )
|
| 156 |
print weight
|
| 157 |
#print tree.weight2
|
| 158 |
#tree.Fill()
|
| 159 |
|
| 160 |
|
| 161 |
|
| 162 |
def plot( tree, plot, cut, save=False ):
|
| 163 |
# read axis config file
|
| 164 |
try:
|
| 165 |
label = axisConf.get( plot, "label" )
|
| 166 |
unit = axisConf.get( plot, "unit" )
|
| 167 |
if unit != "":
|
| 168 |
unit = " ["+unit+"]"
|
| 169 |
except:
|
| 170 |
print "please specify label and unit in axis configuration file"
|
| 171 |
label = plot
|
| 172 |
unit = ""
|
| 173 |
|
| 174 |
# modify histo
|
| 175 |
histo = createHistoFromTree( tree, plot, cut )
|
| 176 |
histo.SetTitle(";%s%s;Entries"%(label,unit))
|
| 177 |
histo.SetLineColor(1)
|
| 178 |
|
| 179 |
# draw histo
|
| 180 |
canvas = ROOT.TCanvas()
|
| 181 |
canvas.cd()
|
| 182 |
canvas.SetLogy()
|
| 183 |
histo.Draw("hist")
|
| 184 |
|
| 185 |
# draw information
|
| 186 |
cutText = ROOT.TPaveText(.5,.8,.9,.9,"ndc")
|
| 187 |
cutText.SetBorderSize(0)
|
| 188 |
cutText.SetFillColor(0)
|
| 189 |
for line in cut.split("&&"):
|
| 190 |
cutText.AddText( line )
|
| 191 |
cutText.Draw()
|
| 192 |
|
| 193 |
topInfo = ROOT.TPaveText(.1,.95,1,1,"ndc")
|
| 194 |
topInfo.SetBorderSize(0)
|
| 195 |
topInfo.SetFillColor(0)
|
| 196 |
topInfo.AddText("Work in progress, %sfb^{-1}, #sqrt{s}=8TeV"%"?")
|
| 197 |
topInfo.Draw()
|
| 198 |
|
| 199 |
if save:
|
| 200 |
canvas.SaveAs("%s.pdf"%plot)
|
| 201 |
else:
|
| 202 |
raw_input()
|
| 203 |
|
| 204 |
|
| 205 |
if __name__ == "__main__":
|
| 206 |
# include knowledge about objects saved in the tree
|
| 207 |
ROOT.gSystem.Load("../treeWriter/libTreeObjects.so")
|
| 208 |
|
| 209 |
arguments = argparse.ArgumentParser( description="Calculate weighting "
|
| 210 |
+"factors for QCD background estimation." )
|
| 211 |
arguments.add_argument("-f", "--file", default="../treeWriter/myTree.root",
|
| 212 |
help="ROOT file containing a TTree produced by 'treeWriter'.")
|
| 213 |
arguments.add_argument("-c", "--cut", default="", help="TCut string.")
|
| 214 |
arguments.add_argument("-d", "--distribution", default =['met'], nargs="+",
|
| 215 |
help="Distributions to plot.")
|
| 216 |
arguments.add_argument("--save", action="store_true", help="Save canvas as pdf.")
|
| 217 |
opts = arguments.parse_args()
|
| 218 |
|
| 219 |
tree = readTree( opts.file )
|
| 220 |
|
| 221 |
# see https://twiki.cern.ch/twiki/bin/viewauth/CMS/CutBasedPhotonID2012
|
| 222 |
# for more information on 2012 Photon ID
|
| 223 |
# TODO: Single tower H/E
|
| 224 |
# TODO: correct isolation with ρ
|
| 225 |
photonCut2012ID = "photon.sigmaIetaIeta<0.012 \
|
| 226 |
&& photon.chargedIso < 2.6 \
|
| 227 |
&& photon.neutralIso < 3.5 + 0.04*photon.pt \
|
| 228 |
&& photon.photonIso < 1.3 + 0.005*photon.pt"
|
| 229 |
|
| 230 |
jetPhotonCut = "photon.sigmaIetaIeta > 0.012 && photon.sigmaIetaIeta<0.02 \
|
| 231 |
&& photon.chargedIso < 2.6 \
|
| 232 |
&& photon.neutralIso < 3.5 + 0.04*photon.pt \
|
| 233 |
&& photon.photonIso < 1.3 + 0.005*photon.pt"
|
| 234 |
|
| 235 |
fillWeights( tree, jetPhotonCut, photonCut2012ID )
|
| 236 |
|
| 237 |
if opts.cut == "photon":
|
| 238 |
opts.cut = photonCut2012ID
|
| 239 |
elif opts.cut == "jetphoton":
|
| 240 |
opts.cut = jetPhotonCut
|
| 241 |
|
| 242 |
if "all" in opts.distribution:
|
| 243 |
opts.distribution = axisConf.sections()
|
| 244 |
|
| 245 |
for plots in opts.distribution:
|
| 246 |
if ":" in plots:
|
| 247 |
plot2D( tree, plots, opts.cut, opts.save )
|
| 248 |
else:
|
| 249 |
plot( tree, plots, opts.cut, opts.save )
|
| 250 |
|
