| 1 |
kiesel |
1.1 |
#! /usr/bin/env python2
|
| 2 |
|
|
# -*- coding: utf-8 -*-
|
| 3 |
|
|
|
| 4 |
|
|
import ROOT
|
| 5 |
|
|
import argparse
|
| 6 |
|
|
import ConfigParser
|
| 7 |
|
|
import os.path
|
| 8 |
|
|
import Styles
|
| 9 |
|
|
Styles.tdrStyle()
|
| 10 |
|
|
|
| 11 |
|
|
axisConf = ConfigParser.SafeConfigParser()
|
| 12 |
|
|
axisConf.read("axis.cfg")
|
| 13 |
|
|
|
| 14 |
|
|
integratedLuminosity = 19.3 #fb
|
| 15 |
|
|
|
| 16 |
kiesel |
1.2 |
def readTree( filename, treename = "susyTree" ):
|
| 17 |
kiesel |
1.1 |
"""
|
| 18 |
|
|
filename: name of file containing the tree
|
| 19 |
|
|
treename: name of the tree
|
| 20 |
|
|
returns: TTree Object
|
| 21 |
|
|
"""
|
| 22 |
|
|
if not os.path.isfile(filename):
|
| 23 |
|
|
print( "File %s does not exist"%filename)
|
| 24 |
|
|
tree = ROOT.TChain( treename )
|
| 25 |
|
|
tree.AddFile( filename )
|
| 26 |
|
|
return tree
|
| 27 |
|
|
|
| 28 |
|
|
def createHistoFromTree2D(tree, variable, weight="", nBin=50):
|
| 29 |
|
|
"""
|
| 30 |
|
|
tree: tree to create histo from
|
| 31 |
|
|
variable: variable to plot (must be a branch of the tree)
|
| 32 |
|
|
weight: weights to apply (e.g. "var1*(var2 > 15)" will use weights from var1 and cut on var2 > 15
|
| 33 |
|
|
nBins, firstBin, lastBin: number of bins, first bin and last bin (same as in TH1F constructor)
|
| 34 |
|
|
nEvents: number of events to process (-1 = all)
|
| 35 |
|
|
returns: histogram
|
| 36 |
|
|
"""
|
| 37 |
|
|
from random import randint
|
| 38 |
|
|
from sys import maxint
|
| 39 |
|
|
#make a random name you could give something meaningfull here,
|
| 40 |
|
|
#but that would make this less readable
|
| 41 |
|
|
name = "%x"%(randint(0, maxint))
|
| 42 |
|
|
|
| 43 |
|
|
# automatic binning if lastbin < firstbin
|
| 44 |
|
|
result = ROOT.TH2F(name, variable, nBin, 0, -1, nBin, 0, -1)
|
| 45 |
|
|
result.Sumw2()
|
| 46 |
|
|
tree.Draw("%s>>%s"%(variable, name), weight, "goff,colz")
|
| 47 |
|
|
return result
|
| 48 |
|
|
|
| 49 |
|
|
def createHistoFromTree(tree, variable, weight="", nBins=100, firstBin=None, lastBin=None, nEvents=-1):
|
| 50 |
|
|
"""
|
| 51 |
|
|
tree: tree to create histo from
|
| 52 |
|
|
variable: variable to plot (must be a branch of the tree)
|
| 53 |
|
|
weight: weights to apply (e.g. "var1*(var2 > 15)" will use weights from var1 and cut on var2 > 15
|
| 54 |
|
|
nBins, firstBin, lastBin: number of bins, first bin and last bin (same as in TH1F constructor)
|
| 55 |
|
|
nEvents: number of events to process (-1 = all)
|
| 56 |
|
|
returns: histogram
|
| 57 |
|
|
"""
|
| 58 |
|
|
if ":" in variable:
|
| 59 |
|
|
return createHistoFromTree2D(tree, variable, weight="")
|
| 60 |
|
|
from ROOT import TH1F
|
| 61 |
|
|
from random import randint
|
| 62 |
|
|
from sys import maxint
|
| 63 |
kiesel |
1.2 |
import re
|
| 64 |
kiesel |
1.1 |
if nEvents < 0:
|
| 65 |
|
|
nEvents = maxint
|
| 66 |
kiesel |
1.2 |
# Get Minimum and Maximum for x-axis. The replacement is needed for e.g.
|
| 67 |
|
|
# photon[0].pt. TODO: find better min/max for this case.
|
| 68 |
kiesel |
1.1 |
if firstBin == None:
|
| 69 |
kiesel |
1.2 |
firstBin = tree.GetMinimum( re.sub("\[\d\]", "", variable ) )
|
| 70 |
kiesel |
1.1 |
if lastBin == None:
|
| 71 |
kiesel |
1.2 |
lastBin = tree.GetMaximum( re.sub("\[\d\]", "", variable ) )
|
| 72 |
kiesel |
1.1 |
#make a random name you could give something meaningfull here,
|
| 73 |
|
|
#but that would make this less readable
|
| 74 |
|
|
name = "%x"%(randint(0, maxint))
|
| 75 |
|
|
result = TH1F(name, variable, nBins, firstBin, lastBin)
|
| 76 |
|
|
result.Sumw2()
|
| 77 |
|
|
tree.Draw("%s>>%s"%(variable, name), weight, "goff", nEvents)
|
| 78 |
|
|
return result
|
| 79 |
|
|
|
| 80 |
|
|
def plot2D( tree, plot, cut, save=False ):
|
| 81 |
|
|
# read axis config file
|
| 82 |
|
|
var = plot.split(":")
|
| 83 |
|
|
try:
|
| 84 |
|
|
xlabel = axisConf.get( var[0], "label" )
|
| 85 |
|
|
except:
|
| 86 |
|
|
xlabel = ""
|
| 87 |
|
|
try:
|
| 88 |
|
|
ylabel = axisConf.get( var[1], "label" )
|
| 89 |
|
|
except:
|
| 90 |
|
|
ylabel = ""
|
| 91 |
|
|
try:
|
| 92 |
|
|
xunit = axisConf.get( var[0], "unit" )
|
| 93 |
|
|
except:
|
| 94 |
|
|
xunit = ""
|
| 95 |
|
|
try:
|
| 96 |
|
|
yunit = axisConf.get( var[1], "unit" )
|
| 97 |
|
|
except:
|
| 98 |
|
|
yunit = ""
|
| 99 |
|
|
if xunit != "":
|
| 100 |
|
|
xunit = " ["+xunit+"]"
|
| 101 |
|
|
if yunit != "":
|
| 102 |
|
|
yunit = " ["+yunit+"]"
|
| 103 |
|
|
|
| 104 |
|
|
# modify histo
|
| 105 |
|
|
histo = createHistoFromTree2D( tree, plot, cut )
|
| 106 |
|
|
histo.SetTitle(";%s%s;%s%s"%(ylabel,yunit,xlabel,xunit))
|
| 107 |
|
|
|
| 108 |
|
|
# draw histo
|
| 109 |
|
|
canvas = ROOT.TCanvas()
|
| 110 |
|
|
canvas.cd()
|
| 111 |
|
|
histo.Draw("colz")
|
| 112 |
|
|
|
| 113 |
|
|
# draw information
|
| 114 |
|
|
cutText = ROOT.TPaveText(.5,.8,.9,.9,"ndc")
|
| 115 |
|
|
cutText.SetBorderSize(0)
|
| 116 |
|
|
cutText.SetFillColor(0)
|
| 117 |
|
|
for line in cut.split("&&"):
|
| 118 |
|
|
cutText.AddText( line )
|
| 119 |
|
|
cutText.Draw()
|
| 120 |
|
|
|
| 121 |
|
|
topInfo = ROOT.TPaveText(.1,.95,1,1,"ndc")
|
| 122 |
|
|
topInfo.SetBorderSize(0)
|
| 123 |
|
|
topInfo.SetFillColor(0)
|
| 124 |
|
|
topInfo.AddText("Work in progress, %sfb^{-1}, #sqrt{s}=8TeV"%"?")
|
| 125 |
|
|
topInfo.Draw()
|
| 126 |
|
|
|
| 127 |
|
|
if save:
|
| 128 |
|
|
canvas.SaveAs("%.pdf"%variable)
|
| 129 |
|
|
else:
|
| 130 |
|
|
raw_input()
|
| 131 |
|
|
|
| 132 |
|
|
def plot( tree, plot, cut, save=False ):
|
| 133 |
|
|
# read axis config file
|
| 134 |
|
|
try:
|
| 135 |
|
|
label = axisConf.get( plot, "label" )
|
| 136 |
|
|
unit = axisConf.get( plot, "unit" )
|
| 137 |
|
|
if unit != "":
|
| 138 |
|
|
unit = " ["+unit+"]"
|
| 139 |
|
|
except:
|
| 140 |
|
|
print "please specify label and unit in axis configuration file"
|
| 141 |
|
|
label = plot
|
| 142 |
|
|
unit = ""
|
| 143 |
|
|
|
| 144 |
|
|
# modify histo
|
| 145 |
|
|
histo = createHistoFromTree( tree, plot, cut )
|
| 146 |
|
|
histo.SetTitle(";%s%s;Entries"%(label,unit))
|
| 147 |
|
|
histo.SetLineColor(1)
|
| 148 |
|
|
|
| 149 |
|
|
# draw histo
|
| 150 |
|
|
canvas = ROOT.TCanvas()
|
| 151 |
|
|
canvas.cd()
|
| 152 |
|
|
canvas.SetLogy()
|
| 153 |
|
|
histo.Draw("hist")
|
| 154 |
|
|
|
| 155 |
|
|
# draw information
|
| 156 |
|
|
cutText = ROOT.TPaveText(.5,.8,.9,.9,"ndc")
|
| 157 |
|
|
cutText.SetBorderSize(0)
|
| 158 |
|
|
cutText.SetFillColor(0)
|
| 159 |
|
|
for line in cut.split("&&"):
|
| 160 |
|
|
cutText.AddText( line )
|
| 161 |
|
|
cutText.Draw()
|
| 162 |
|
|
|
| 163 |
|
|
topInfo = ROOT.TPaveText(.1,.95,1,1,"ndc")
|
| 164 |
|
|
topInfo.SetBorderSize(0)
|
| 165 |
|
|
topInfo.SetFillColor(0)
|
| 166 |
|
|
topInfo.AddText("Work in progress, %sfb^{-1}, #sqrt{s}=8TeV"%"?")
|
| 167 |
|
|
topInfo.Draw()
|
| 168 |
|
|
|
| 169 |
|
|
if save:
|
| 170 |
|
|
canvas.SaveAs("%.pdf"%variable)
|
| 171 |
|
|
else:
|
| 172 |
|
|
raw_input()
|
| 173 |
|
|
|
| 174 |
|
|
|
| 175 |
|
|
if __name__ == "__main__":
|
| 176 |
|
|
arguments = argparse.ArgumentParser()
|
| 177 |
|
|
arguments.add_argument("-f", "--file", default="myTree.root", help="ROOT file containing tree")
|
| 178 |
|
|
arguments.add_argument("-c", "--cut", default="", help="Cut string")
|
| 179 |
|
|
arguments.add_argument("-d", "--distribution", default =['met'], nargs="+",
|
| 180 |
|
|
help="Distribution which shall be plotted.")
|
| 181 |
|
|
arguments.add_argument("--save", action="store_true", help="Save canvas as pdf.")
|
| 182 |
|
|
|
| 183 |
|
|
opts = arguments.parse_args()
|
| 184 |
|
|
|
| 185 |
|
|
tree = readTree( opts.file )
|
| 186 |
|
|
|
| 187 |
|
|
# see https://twiki.cern.ch/twiki/bin/viewauth/CMS/CutBasedPhotonID2012
|
| 188 |
|
|
# for more information on 2012 Photon ID
|
| 189 |
|
|
# TODO: Conversion safe electron veto
|
| 190 |
|
|
# TODO: Single tower H/E
|
| 191 |
|
|
# TODO: correct isolation with ρ
|
| 192 |
|
|
photonCut2012ID = "photon.sigmaIetaIeta<0.012 \
|
| 193 |
|
|
&& photon.chargedIso < 2.6 \
|
| 194 |
|
|
&& photon.neutralIso < 3.5 + 0.04*photon.pt \
|
| 195 |
|
|
&& photon.photonIso < 1.3 + 0.005*photon.pt"
|
| 196 |
|
|
|
| 197 |
|
|
jetPhotonCut = "photon.sigmaIetaIeta > 0.012 && photon.sigmaIetaIeta<0.02 \
|
| 198 |
|
|
&& photon.chargedIso < 2.6 \
|
| 199 |
|
|
&& photon.neutralIso < 3.5 + 0.04*photon.pt \
|
| 200 |
|
|
&& photon.photonIso < 1.3 + 0.005*photon.pt"
|
| 201 |
|
|
|
| 202 |
|
|
if opts.cut == "photon":
|
| 203 |
|
|
opts.cut = photonCut2012ID
|
| 204 |
|
|
elif opts.cut == "jetphoton":
|
| 205 |
|
|
opts.cut = jetPhotonCut
|
| 206 |
|
|
|
| 207 |
|
|
if "all" in opts.distribution:
|
| 208 |
|
|
opts.distribution = axisConf.sections()
|
| 209 |
|
|
|
| 210 |
|
|
for plots in opts.distribution:
|
| 211 |
|
|
if ":" in plots:
|
| 212 |
|
|
plot2D( tree, plots, opts.cut, opts.save )
|
| 213 |
|
|
else:
|
| 214 |
|
|
plot( tree, plots, opts.cut, opts.save )
|
| 215 |
|
|
|
