1 |
#!/usr/bin/env python
|
2 |
import sys
|
3 |
import os
|
4 |
import re
|
5 |
import time
|
6 |
from math import *
|
7 |
from array import *
|
8 |
from decimal import *
|
9 |
from optparse import OptionParser
|
10 |
from OSUT3Analysis.Configuration.configurationOptions import *
|
11 |
from OSUT3Analysis.Configuration.processingUtilities import *
|
12 |
from OSUT3Analysis.Configuration.formattingUtilities import *
|
13 |
|
14 |
|
15 |
|
16 |
### parse the command-line options
|
17 |
|
18 |
parser = OptionParser()
|
19 |
parser = set_commandline_arguments(parser)
|
20 |
|
21 |
parser.add_option("-f", "--fancy", action="store_true", dest="makeFancy", default=False,
|
22 |
help="removes the title and replaces it with the official CMS plot heading")
|
23 |
parser.add_option("--ylog", action="store_true", dest="setLogY", default=False,
|
24 |
help="Set logarithmic scale on vertical axis on all plots")
|
25 |
parser.add_option("--ymin", dest="setYMin",
|
26 |
help="Minimum of y axis")
|
27 |
parser.add_option("--ymax", dest="setYMax",
|
28 |
help="Maximum of y axis")
|
29 |
parser.add_option("-E", "--ratioRelErrMax", dest="ratioRelErrMax",
|
30 |
help="maximum error used in rebinning the ratio histogram")
|
31 |
|
32 |
(arguments, args) = parser.parse_args()
|
33 |
|
34 |
|
35 |
if arguments.localConfig:
|
36 |
sys.path.append(os.getcwd())
|
37 |
exec("from " + arguments.localConfig.rstrip('.py') + " import *")
|
38 |
|
39 |
#### deal with conflicting arguments
|
40 |
if arguments.normalizeToData and arguments.normalizeToUnitArea:
|
41 |
print "Conflicting normalizations requsted, will normalize to unit area"
|
42 |
arguments.normalizeToData = False
|
43 |
if arguments.normalizeToData and arguments.noStack:
|
44 |
print "You have asked to scale non-stacked backgrounds to data. This is a very strange request. Will normalize to unit area instead"
|
45 |
arguments.normalizeToData = False
|
46 |
arguments.normalizeToUnitArea = True
|
47 |
if arguments.makeRatioPlots and arguments.makeDiffPlots:
|
48 |
print "You have requested both ratio and difference plots. Will make just ratio plots instead"
|
49 |
arguments.makeRatioPlots = False
|
50 |
if arguments.makeRatioPlots and arguments.noStack:
|
51 |
print "You have asked to make a ratio plot and to not stack the backgrounds. This is a very strange request. Will skip making the ratio plot."
|
52 |
arguments.makeRatioPlots = False
|
53 |
if arguments.makeDiffPlots and arguments.noStack:
|
54 |
print "You have asked to make a difference plot and to not stack the backgrounds. This is a very strange request. Will skip making the difference plot."
|
55 |
arguments.makeDiffPlots = False
|
56 |
|
57 |
|
58 |
from ROOT import TFile, gROOT, gStyle, gDirectory, TStyle, THStack, TH1F, TCanvas, TString, TLegend, TLegendEntry, THStack, TIter, TKey, TPaveLabel, gPad
|
59 |
|
60 |
|
61 |
### setting ROOT options so our plots will look awesome and everyone will love us
|
62 |
|
63 |
gROOT.SetBatch()
|
64 |
gStyle.SetOptStat(0)
|
65 |
gStyle.SetCanvasBorderMode(0)
|
66 |
gStyle.SetPadBorderMode(0)
|
67 |
gStyle.SetPadColor(0)
|
68 |
gStyle.SetCanvasColor(0)
|
69 |
gStyle.SetTextFont(42)
|
70 |
gStyle.SetCanvasDefH(600)
|
71 |
gStyle.SetCanvasDefW(600)
|
72 |
gStyle.SetCanvasDefX(0)
|
73 |
gStyle.SetCanvasDefY(0)
|
74 |
gStyle.SetPadTopMargin(0.07)
|
75 |
gStyle.SetPadBottomMargin(0.13)
|
76 |
gStyle.SetPadLeftMargin(0.15)
|
77 |
gStyle.SetPadRightMargin(0.05)
|
78 |
gStyle.SetTitleColor(1, "XYZ")
|
79 |
gStyle.SetTitleFont(42, "XYZ")
|
80 |
gStyle.SetTitleSize(0.04, "XYZ")
|
81 |
gStyle.SetTitleXOffset(1.1)
|
82 |
gStyle.SetTitleYOffset(2)
|
83 |
gStyle.SetTextAlign(12)
|
84 |
gStyle.SetLabelColor(1, "XYZ")
|
85 |
gStyle.SetLabelFont(42, "XYZ")
|
86 |
gStyle.SetLabelOffset(0.007, "XYZ")
|
87 |
gStyle.SetLabelSize(0.04, "XYZ")
|
88 |
gStyle.SetAxisColor(1, "XYZ")
|
89 |
gStyle.SetStripDecimals(True)
|
90 |
gStyle.SetTickLength(0.03, "XYZ")
|
91 |
gStyle.SetNdivisions(510, "XYZ")
|
92 |
gStyle.SetPadTickX(1)
|
93 |
gStyle.SetPadTickY(1)
|
94 |
gROOT.ForceStyle()
|
95 |
|
96 |
|
97 |
#set the text for the luminosity label
|
98 |
if(intLumi < 1000.):
|
99 |
LumiInPb = intLumi
|
100 |
LumiText = "L_{int} = " + str(intLumi) + " pb^{-1}"
|
101 |
LumiText = "L_{int} = " + str.format('{0:.1f}', LumiInPb) + " pb^{-1}"
|
102 |
else:
|
103 |
LumiInFb = intLumi/1000.
|
104 |
LumiText = "L_{int} = " + str.format('{0:.1f}', LumiInFb) + " fb^{-1}"
|
105 |
|
106 |
#bestest place for lumi. label, in top left corner
|
107 |
topLeft_x_left = 0.1375839
|
108 |
topLeft_x_right = 0.4580537
|
109 |
topLeft_y_bottom = 0.8479021
|
110 |
topLeft_y_top = 0.9475524
|
111 |
topLeft_y_offset = 0.035
|
112 |
|
113 |
#set the text for the fancy heading
|
114 |
HeaderText = "CMS Preliminary: " + LumiText + " at #sqrt{s} = 8 TeV"
|
115 |
|
116 |
#position for header
|
117 |
header_x_left = 0.2181208
|
118 |
header_x_right = 0.9562937
|
119 |
header_y_bottom = 0.9479866
|
120 |
header_y_top = 0.9947552
|
121 |
|
122 |
|
123 |
|
124 |
##########################################################################################################################################
|
125 |
##########################################################################################################################################
|
126 |
##########################################################################################################################################
|
127 |
|
128 |
# some fancy-ass code from Andrzej Zuranski to merge bins in the ratio plot until the error goes below some threshold
|
129 |
def ratioHistogram( dataHist, mcHist, relErrMax=0.1):
|
130 |
|
131 |
if not dataHist:
|
132 |
print "Error: trying to run ratioHistogram but dataHist is invalid"
|
133 |
return
|
134 |
|
135 |
if not mcHist:
|
136 |
print "Error: trying to run ratioHistogram but mcHist is invalid"
|
137 |
return
|
138 |
|
139 |
def groupR(group):
|
140 |
Data,MC = [float(sum(hist.GetBinContent(i) for i in group)) for hist in [dataHist,mcHist]]
|
141 |
return (Data-MC)/MC if MC else 0
|
142 |
|
143 |
def groupErr(group):
|
144 |
Data,MC = [float(sum(hist.GetBinContent(i) for i in group)) for hist in [dataHist,mcHist]]
|
145 |
dataErr2,mcErr2 = [sum(hist.GetBinError(i)**2 for i in group) for hist in [dataHist,mcHist]]
|
146 |
if Data > 0 and MC > 0 and Data != MC:
|
147 |
return abs(math.sqrt( (dataErr2+mcErr2)/(Data-MC)**2 + mcErr2/MC**2 ) * (Data-MC)/MC)
|
148 |
else:
|
149 |
return 0
|
150 |
|
151 |
def regroup(groups):
|
152 |
err,iG = max( (groupErr(g),groups.index(g)) for g in groups )
|
153 |
if err < relErrMax or len(groups)<3 : return groups
|
154 |
iH = max( [iG-1,iG+1], key = lambda i: groupErr(groups[i]) if 0<=i<len(groups) else -1 )
|
155 |
iLo,iHi = sorted([iG,iH])
|
156 |
return regroup(groups[:iLo] + [groups[iLo]+groups[iHi]] + groups[iHi+1:])
|
157 |
|
158 |
#don't rebin the histograms of the number of a given object (except for the pileup ones)
|
159 |
if ((dataHist.GetName().find("num") is not -1 and dataHist.GetName().find("Primaryvertexs") is -1) or
|
160 |
dataHist.GetName().find("CutFlow") is not -1 or
|
161 |
dataHist.GetName().find("GenMatch") is not -1):
|
162 |
ratio = dataHist.Clone()
|
163 |
ratio.Add(mcHist,-1)
|
164 |
ratio.Divide(mcHist)
|
165 |
ratio.SetTitle("")
|
166 |
else:
|
167 |
groups = regroup( [(i,) for i in range(1,1+dataHist.GetNbinsX())] )
|
168 |
ratio = TH1F("ratio","",len(groups), array('d', [dataHist.GetBinLowEdge(min(g)) for g in groups ] + [dataHist.GetXaxis().GetBinUpEdge(dataHist.GetNbinsX())]) )
|
169 |
for i,g in enumerate(groups) :
|
170 |
ratio.SetBinContent(i+1,groupR(g))
|
171 |
ratio.SetBinError(i+1,groupErr(g))
|
172 |
|
173 |
ratio.GetYaxis().SetTitle("#frac{Data-MC}{MC}")
|
174 |
ratio.SetLineColor(1)
|
175 |
ratio.SetLineWidth(2)
|
176 |
return ratio
|
177 |
|
178 |
##########################################################################################################################################
|
179 |
##########################################################################################################################################
|
180 |
##########################################################################################################################################
|
181 |
|
182 |
|
183 |
def MakeOneDHist(pathToDir,histogramName):
|
184 |
|
185 |
blindData = False
|
186 |
# To blind histograms, define a list of histsToBlind in the localOptions.py file
|
187 |
try: # Use "try" in case histsToBlind does not exist
|
188 |
if histsToBlind:
|
189 |
for histToBlind in histsToBlind:
|
190 |
if histToBlind in histogramName:
|
191 |
print "Blinding data for histogram " + histogramName
|
192 |
blindData = True
|
193 |
except NameError:
|
194 |
time.sleep(0.000001) # Do nothing if histsToBlind does not exist
|
195 |
|
196 |
numBgMCSamples = 0
|
197 |
numDataSamples = 0
|
198 |
numSignalSamples = 0
|
199 |
|
200 |
Stack = THStack("stack",histogramName)
|
201 |
|
202 |
HeaderLabel = TPaveLabel(header_x_left,header_y_bottom,header_x_right,header_y_top,HeaderText,"NDC")
|
203 |
HeaderLabel.SetTextAlign(32)
|
204 |
HeaderLabel.SetBorderSize(0)
|
205 |
HeaderLabel.SetFillColor(0)
|
206 |
HeaderLabel.SetFillStyle(0)
|
207 |
|
208 |
LumiLabel = TPaveLabel(topLeft_x_left,topLeft_y_bottom,topLeft_x_right,topLeft_y_top,LumiText,"NDC")
|
209 |
LumiLabel.SetBorderSize(0)
|
210 |
LumiLabel.SetFillColor(0)
|
211 |
LumiLabel.SetFillStyle(0)
|
212 |
|
213 |
NormLabel = TPaveLabel()
|
214 |
NormLabel.SetDrawOption("NDC")
|
215 |
NormLabel.SetX1NDC(topLeft_x_left)
|
216 |
NormLabel.SetX2NDC(topLeft_x_right)
|
217 |
|
218 |
NormLabel.SetBorderSize(0)
|
219 |
NormLabel.SetFillColor(0)
|
220 |
NormLabel.SetFillStyle(0)
|
221 |
|
222 |
NormText = ""
|
223 |
if arguments.normalizeToUnitArea:
|
224 |
NormText = "Scaled to unit area"
|
225 |
elif arguments.normalizeToData:
|
226 |
NormText = "MC scaled to data"
|
227 |
NormLabel.SetLabel(NormText)
|
228 |
|
229 |
|
230 |
BgMCLegend = TLegend()
|
231 |
BgTitle = BgMCLegend.AddEntry(0, "Data & Bkgd. MC", "H")
|
232 |
BgTitle.SetTextAlign(22)
|
233 |
BgTitle.SetTextFont(62)
|
234 |
BgMCLegend.SetBorderSize(0)
|
235 |
BgMCLegend.SetFillColor(0)
|
236 |
BgMCLegend.SetFillStyle(0)
|
237 |
SignalMCLegend = TLegend()
|
238 |
SignalTitle = SignalMCLegend.AddEntry(0, "Signal MC", "H")
|
239 |
SignalTitle.SetTextAlign(22)
|
240 |
SignalTitle.SetTextFont(62)
|
241 |
SignalMCLegend.SetBorderSize(0)
|
242 |
SignalMCLegend.SetFillColor(0)
|
243 |
SignalMCLegend.SetFillStyle(0)
|
244 |
|
245 |
outputFile.cd(pathToDir)
|
246 |
Canvas = TCanvas(histogramName)
|
247 |
BgMCHistograms = []
|
248 |
BgMCLegendEntries = []
|
249 |
SignalMCHistograms = []
|
250 |
SignalMCLegendEntries = []
|
251 |
DataHistograms = []
|
252 |
DataLegendEntries = []
|
253 |
|
254 |
|
255 |
backgroundIntegral = 0
|
256 |
dataIntegral = 0
|
257 |
scaleFactor = 1
|
258 |
|
259 |
for sample in processed_datasets: # loop over different samples as listed in configurationOptions.py
|
260 |
dataset_file = "%s/%s.root" % (condor_dir,sample)
|
261 |
inputFile = TFile(dataset_file)
|
262 |
HistogramObj = inputFile.Get(pathToDir+"/"+histogramName)
|
263 |
if not HistogramObj:
|
264 |
print "WARNING: Could not find histogram " + pathToDir + "/" + histogramName + " in file " + dataset_file + ". Will skip it and continue."
|
265 |
continue
|
266 |
Histogram = HistogramObj.Clone()
|
267 |
Histogram.SetDirectory(0)
|
268 |
inputFile.Close()
|
269 |
if arguments.rebinFactor:
|
270 |
RebinFactor = int(arguments.rebinFactor)
|
271 |
#don't rebin histograms which will have less than 5 bins or any gen-matching histograms
|
272 |
if Histogram.GetNbinsX() >= RebinFactor*5 and Histogram.GetName().find("GenMatch") is -1:
|
273 |
Histogram.Rebin(RebinFactor)
|
274 |
|
275 |
xAxisLabel = Histogram.GetXaxis().GetTitle()
|
276 |
unitBeginIndex = xAxisLabel.find("[")
|
277 |
unitEndIndex = xAxisLabel.find("]")
|
278 |
|
279 |
if unitBeginIndex is not -1 and unitEndIndex is not -1: #x axis has a unit
|
280 |
yAxisLabel = "Entries / " + str(Histogram.GetXaxis().GetBinWidth(1)) + " " + xAxisLabel[unitBeginIndex+1:unitEndIndex]
|
281 |
else:
|
282 |
yAxisLabel = "Entries per bin (" + str(Histogram.GetXaxis().GetBinWidth(1)) + " width)"
|
283 |
|
284 |
if not arguments.makeFancy:
|
285 |
histoTitle = Histogram.GetTitle()
|
286 |
else:
|
287 |
histoTitle = ""
|
288 |
|
289 |
|
290 |
legLabel = labels[sample]
|
291 |
if (arguments.printYields):
|
292 |
yieldHist = Histogram.Integral()
|
293 |
legLabel = legLabel + " (%.1f)" % yieldHist
|
294 |
|
295 |
if( types[sample] == "bgMC"):
|
296 |
|
297 |
numBgMCSamples += 1
|
298 |
backgroundIntegral += Histogram.Integral()
|
299 |
|
300 |
Histogram.SetLineStyle(1)
|
301 |
if(arguments.noStack):
|
302 |
Histogram.SetFillStyle(0)
|
303 |
Histogram.SetLineColor(colors[sample])
|
304 |
Histogram.SetLineWidth(2)
|
305 |
else:
|
306 |
Histogram.SetFillStyle(1001)
|
307 |
Histogram.SetFillColor(colors[sample])
|
308 |
Histogram.SetLineColor(1)
|
309 |
Histogram.SetLineWidth(1)
|
310 |
|
311 |
BgMCLegendEntries.append(legLabel)
|
312 |
BgMCHistograms.append(Histogram)
|
313 |
|
314 |
|
315 |
elif( types[sample] == "signalMC"):
|
316 |
|
317 |
numSignalSamples += 1
|
318 |
|
319 |
Histogram.SetFillStyle(0)
|
320 |
Histogram.SetLineColor(colors[sample])
|
321 |
Histogram.SetLineStyle(1)
|
322 |
Histogram.SetLineWidth(2)
|
323 |
if(arguments.normalizeToUnitArea and Histogram.Integral() > 0):
|
324 |
Histogram.Scale(1./Histogram.Integral())
|
325 |
|
326 |
SignalMCLegendEntries.append(legLabel)
|
327 |
SignalMCHistograms.append(Histogram)
|
328 |
|
329 |
elif( types[sample] == "data" and not blindData):
|
330 |
|
331 |
numDataSamples += 1
|
332 |
dataIntegral += Histogram.Integral()
|
333 |
|
334 |
Histogram.SetMarkerStyle(20)
|
335 |
Histogram.SetMarkerSize(0.8)
|
336 |
Histogram.SetFillStyle(0)
|
337 |
Histogram.SetLineColor(colors[sample])
|
338 |
Histogram.SetLineStyle(1)
|
339 |
Histogram.SetLineWidth(2)
|
340 |
if(arguments.normalizeToUnitArea and Histogram.Integral() > 0):
|
341 |
Histogram.Scale(1./Histogram.Integral())
|
342 |
|
343 |
DataLegendEntries.append(legLabel)
|
344 |
DataHistograms.append(Histogram)
|
345 |
|
346 |
#scaling histograms as per user's specifications
|
347 |
if dataIntegral > 0 and backgroundIntegral > 0:
|
348 |
scaleFactor = dataIntegral/backgroundIntegral
|
349 |
for bgMCHist in BgMCHistograms:
|
350 |
if arguments.normalizeToData:
|
351 |
bgMCHist.Scale(scaleFactor)
|
352 |
|
353 |
if arguments.normalizeToUnitArea and not arguments.noStack and backgroundIntegral > 0:
|
354 |
bgMCHist.Scale(1./backgroundIntegral)
|
355 |
elif arguments.normalizeToUnitArea and arguments.noStack and bgMCHist.Integral() > 0:
|
356 |
bgMCHist.Scale(1./bgMCHist.Integral())
|
357 |
|
358 |
if not arguments.noStack:
|
359 |
Stack.Add(bgMCHist)
|
360 |
|
361 |
|
362 |
|
363 |
### formatting data histograms and adding to legend
|
364 |
legendIndex = 0
|
365 |
for Histogram in DataHistograms:
|
366 |
BgMCLegend.AddEntry(Histogram,DataLegendEntries[legendIndex],"LEP")
|
367 |
legendIndex = legendIndex+1
|
368 |
|
369 |
|
370 |
### creating the histogram to represent the statistical errors on the stack
|
371 |
if numBgMCSamples is not 0 and not arguments.noStack:
|
372 |
ErrorHisto = BgMCHistograms[0].Clone("errors")
|
373 |
ErrorHisto.SetFillStyle(3001)
|
374 |
ErrorHisto.SetFillColor(13)
|
375 |
ErrorHisto.SetLineWidth(0)
|
376 |
BgMCLegend.AddEntry(ErrorHisto,"Stat. Errors","F")
|
377 |
for Histogram in BgMCHistograms:
|
378 |
if Histogram is not BgMCHistograms[0]:
|
379 |
ErrorHisto.Add(Histogram)
|
380 |
|
381 |
|
382 |
### formatting bgMC histograms and adding to legend
|
383 |
legendIndex = numBgMCSamples-1
|
384 |
for Histogram in reversed(BgMCHistograms):
|
385 |
if(arguments.noStack):
|
386 |
BgMCLegend.AddEntry(Histogram,BgMCLegendEntries[legendIndex],"L")
|
387 |
else:
|
388 |
BgMCLegend.AddEntry(Histogram,BgMCLegendEntries[legendIndex],"F")
|
389 |
legendIndex = legendIndex-1
|
390 |
|
391 |
|
392 |
### formatting signalMC histograms and adding to legend
|
393 |
legendIndex = 0
|
394 |
for Histogram in SignalMCHistograms:
|
395 |
SignalMCLegend.AddEntry(Histogram,SignalMCLegendEntries[legendIndex],"L")
|
396 |
legendIndex = legendIndex+1
|
397 |
|
398 |
|
399 |
### finding the maximum value of anything going on the canvas, so we know how to set the y-axis
|
400 |
finalMax = 0
|
401 |
if numBgMCSamples is not 0 and not arguments.noStack:
|
402 |
finalMax = ErrorHisto.GetMaximum() + ErrorHisto.GetBinError(ErrorHisto.GetMaximumBin())
|
403 |
else:
|
404 |
for bgMCHist in BgMCHistograms:
|
405 |
if(bgMCHist.GetMaximum() > finalMax):
|
406 |
finalMax = bgMCHist.GetMaximum()
|
407 |
for signalMCHist in SignalMCHistograms:
|
408 |
if(signalMCHist.GetMaximum() > finalMax):
|
409 |
finalMax = signalMCHist.GetMaximum()
|
410 |
for dataHist in DataHistograms:
|
411 |
if(dataHist.GetMaximum() > finalMax):
|
412 |
finalMax = dataHist.GetMaximum() + dataHist.GetBinError(dataHist.GetMaximumBin())
|
413 |
finalMax = 1.15*finalMax
|
414 |
if arguments.setYMax:
|
415 |
finalMax = float(arguments.setYMax)
|
416 |
|
417 |
|
418 |
### Drawing histograms to canvas
|
419 |
|
420 |
outputFile.cd(pathToDir)
|
421 |
|
422 |
makeRatioPlots = arguments.makeRatioPlots
|
423 |
makeDiffPlots = arguments.makeDiffPlots
|
424 |
|
425 |
yAxisMin = 0.0001
|
426 |
if arguments.setYMin:
|
427 |
yAxisMin = float(arguments.setYMin)
|
428 |
|
429 |
if numBgMCSamples is 0 or numDataSamples is not 1:
|
430 |
makeRatioPlots = False
|
431 |
makeDiffPlots = False
|
432 |
if makeRatioPlots or makeDiffPlots:
|
433 |
Canvas.SetFillStyle(0)
|
434 |
Canvas.Divide(1,2)
|
435 |
Canvas.cd(1)
|
436 |
gPad.SetPad(0,0.25,1,1)
|
437 |
gPad.SetMargin(0.15,0.05,0.01,0.07)
|
438 |
gPad.SetFillStyle(0)
|
439 |
gPad.Update()
|
440 |
gPad.Draw()
|
441 |
if arguments.setLogY:
|
442 |
gPad.SetLogy()
|
443 |
Canvas.cd(2)
|
444 |
gPad.SetPad(0,0,1,0.25)
|
445 |
#format: gPad.SetMargin(l,r,b,t)
|
446 |
gPad.SetMargin(0.15,0.05,0.4,0.01)
|
447 |
gPad.SetFillStyle(0)
|
448 |
gPad.SetGridy(1)
|
449 |
gPad.Update()
|
450 |
gPad.Draw()
|
451 |
|
452 |
Canvas.cd(1)
|
453 |
|
454 |
if numBgMCSamples is not 0: # the first thing to draw to the canvas is a bgMC sample
|
455 |
|
456 |
if not arguments.noStack: # draw stacked background samples
|
457 |
Stack.SetTitle(histoTitle)
|
458 |
Stack.Draw("HIST")
|
459 |
Stack.GetXaxis().SetTitle(xAxisLabel)
|
460 |
Stack.GetYaxis().SetTitle(yAxisLabel)
|
461 |
Stack.SetMaximum(finalMax)
|
462 |
Stack.SetMinimum(yAxisMin)
|
463 |
if makeRatioPlots or makeDiffPlots:
|
464 |
Stack.GetHistogram().GetXaxis().SetLabelSize(0)
|
465 |
#draw shaded error bands
|
466 |
ErrorHisto.Draw("A E2 SAME")
|
467 |
|
468 |
else: #draw the unstacked backgrounds
|
469 |
BgMCHistograms[0].SetTitle(histoTitle)
|
470 |
BgMCHistograms[0].Draw("HIST")
|
471 |
BgMCHistograms[0].GetXaxis().SetTitle(xAxisLabel)
|
472 |
BgMCHistograms[0].GetYaxis().SetTitle(yAxisLabel)
|
473 |
BgMCHistograms[0].SetMaximum(finalMax)
|
474 |
BgMCHistograms[0].SetMinimum(yAxisMin)
|
475 |
for bgMCHist in BgMCHistograms:
|
476 |
bgMCHist.Draw("A HIST SAME")
|
477 |
|
478 |
for signalMCHist in SignalMCHistograms:
|
479 |
signalMCHist.Draw("A HIST SAME")
|
480 |
for dataHist in DataHistograms:
|
481 |
dataHist.Draw("A E X0 SAME")
|
482 |
|
483 |
|
484 |
elif numSignalSamples is not 0: # the first thing to draw to the canvas is a signalMC sample
|
485 |
SignalMCHistograms[0].SetTitle(histoTitle)
|
486 |
SignalMCHistograms[0].Draw("HIST")
|
487 |
SignalMCHistograms[0].GetXaxis().SetTitle(xAxisLabel)
|
488 |
SignalMCHistograms[0].GetYaxis().SetTitle(yAxisLabel)
|
489 |
SignalMCHistograms[0].SetMaximum(finalMax)
|
490 |
SignalMCHistograms[0].SetMinimum(yAxisMin)
|
491 |
|
492 |
for signalMCHist in SignalMCHistograms:
|
493 |
if(signalMCHist is not SignalMCHistograms[0]):
|
494 |
signalMCHist.Draw("A HIST SAME")
|
495 |
for dataHist in DataHistograms:
|
496 |
dataHist.Draw("A E X0 SAME")
|
497 |
|
498 |
|
499 |
elif(numDataSamples is not 0): # the first thing to draw to the canvas is a data sample
|
500 |
DataHistograms[0].SetTitle(histoTitle)
|
501 |
DataHistograms[0].Draw("E")
|
502 |
DataHistograms[0].GetXaxis().SetTitle(xAxisLabel)
|
503 |
DataHistograms[0].GetYaxis().SetTitle(yAxisLabel)
|
504 |
DataHistograms[0].SetMaximum(finalMax)
|
505 |
DataHistograms[0].SetMinimum(yAxisMin)
|
506 |
for dataHist in DataHistograms:
|
507 |
if(dataHist is not DataHistograms[0]):
|
508 |
dataHist.Draw("A E X0 SAME")
|
509 |
|
510 |
|
511 |
|
512 |
#legend coordinates, empirically determined :-)
|
513 |
x_left = 0.6761745
|
514 |
x_right = 0.9328859
|
515 |
x_width = x_right - x_left
|
516 |
y_max = 0.9
|
517 |
entry_height = 0.05
|
518 |
|
519 |
if(numBgMCSamples is not 0 or numDataSamples is not 0): #then draw the data & bgMC legend
|
520 |
|
521 |
numExtraEntries = 1 # count the legend title
|
522 |
BgMCLegend.SetX1NDC(x_left)
|
523 |
if numBgMCSamples > 0:
|
524 |
numExtraEntries = numExtraEntries + 1 # count the stat. errors entry
|
525 |
|
526 |
BgMCLegend.SetY1NDC(y_max-entry_height*(numExtraEntries+numBgMCSamples+numDataSamples))
|
527 |
BgMCLegend.SetX2NDC(x_right)
|
528 |
BgMCLegend.SetY2NDC(y_max)
|
529 |
BgMCLegend.Draw()
|
530 |
|
531 |
if(numSignalSamples is not 0): #then draw the signalMC legend to the left of the other one
|
532 |
SignalMCLegend.SetX1NDC(x_left-x_width)
|
533 |
SignalMCLegend.SetY1NDC(y_max-entry_height*(1+numSignalSamples)) # add one for the title
|
534 |
SignalMCLegend.SetX2NDC(x_left)
|
535 |
SignalMCLegend.SetY2NDC(y_max)
|
536 |
SignalMCLegend.Draw()
|
537 |
|
538 |
elif numSignalSamples is not 0: #draw the signalMC legend in the upper right corner
|
539 |
SignalMCLegend.SetX1NDC(x_left)
|
540 |
SignalMCLegend.SetY1NDC(y_max-entry_height*(1+numSignalSamples)) # add one for the title
|
541 |
SignalMCLegend.SetX2NDC(x_right)
|
542 |
SignalMCLegend.SetY2NDC(y_max)
|
543 |
SignalMCLegend.Draw()
|
544 |
|
545 |
|
546 |
# Deciding which text labels to draw and drawing them
|
547 |
drawLumiLabel = False
|
548 |
drawNormLabel = False
|
549 |
offsetNormLabel = False
|
550 |
drawHeaderLabel = False
|
551 |
|
552 |
if not arguments.normalizeToUnitArea or numDataSamples > 0: #don't draw the lumi label if there's no data and it's scaled to unit area
|
553 |
drawLumiLabel = True
|
554 |
#move the normalization label down before drawing if we drew the lumi. label
|
555 |
offsetNormLabel = True
|
556 |
if arguments.normalizeToUnitArea or arguments.normalizeToData:
|
557 |
drawNormLabel = True
|
558 |
if arguments.makeFancy:
|
559 |
drawHeaderLabel = True
|
560 |
drawLumiLabel = False
|
561 |
|
562 |
#now that flags are set, draw the appropriate labels
|
563 |
|
564 |
if drawLumiLabel:
|
565 |
LumiLabel.Draw()
|
566 |
|
567 |
if drawNormLabel:
|
568 |
if offsetNormLabel:
|
569 |
NormLabel.SetY1NDC(topLeft_y_bottom-topLeft_y_offset)
|
570 |
NormLabel.SetY2NDC(topLeft_y_top-topLeft_y_offset)
|
571 |
else:
|
572 |
NormLabel.SetY1NDC(topLeft_y_bottom)
|
573 |
NormLabel.SetY2NDC(topLeft_y_top)
|
574 |
NormLabel.Draw()
|
575 |
|
576 |
if drawHeaderLabel:
|
577 |
HeaderLabel.Draw()
|
578 |
|
579 |
|
580 |
|
581 |
|
582 |
#drawing the ratio or difference plot if requested
|
583 |
|
584 |
if (makeRatioPlots or makeDiffPlots):
|
585 |
Canvas.cd(2)
|
586 |
BgSum = Stack.GetStack().Last()
|
587 |
if makeRatioPlots:
|
588 |
if arguments.ratioRelErrMax:
|
589 |
Comparison = ratioHistogram(DataHistograms[0],BgSum,arguments.ratioRelErrMax)
|
590 |
else:
|
591 |
Comparison = ratioHistogram(DataHistograms[0],BgSum)
|
592 |
elif makeDiffPlots:
|
593 |
Comparison = DataHistograms[0].Clone("diff")
|
594 |
Comparison.Add(BgSum,-1)
|
595 |
Comparison.SetTitle("")
|
596 |
Comparison.GetYaxis().SetTitle("Data-MC")
|
597 |
Comparison.GetXaxis().SetTitle(xAxisLabel)
|
598 |
Comparison.GetYaxis().CenterTitle()
|
599 |
Comparison.GetYaxis().SetTitleSize(0.1)
|
600 |
Comparison.GetYaxis().SetTitleOffset(0.5)
|
601 |
Comparison.GetXaxis().SetTitleSize(0.15)
|
602 |
Comparison.GetYaxis().SetLabelSize(0.1)
|
603 |
Comparison.GetXaxis().SetLabelSize(0.15)
|
604 |
if makeRatioPlots:
|
605 |
RatioYRange = 1.15
|
606 |
if arguments.ratioYRange:
|
607 |
RatioYRange = float(arguments.ratioYRange)
|
608 |
Comparison.GetYaxis().SetRangeUser(-1*RatioYRange, RatioYRange)
|
609 |
elif makeDiffPlots:
|
610 |
YMax = Comparison.GetMaximum()
|
611 |
YMin = Comparison.GetMinimum()
|
612 |
if YMax <= 0 and YMin <= 0:
|
613 |
Comparison.GetYaxis().SetRangeUser(-1.2*YMin,0)
|
614 |
elif YMax >= 0 and YMin >= 0:
|
615 |
Comparison.GetYaxis().SetRangeUser(0,1.2*YMax)
|
616 |
else: #axis crosses y=0
|
617 |
if abs(YMax) > abs(YMin):
|
618 |
Comparison.GetYaxis().SetRangeUser(-1.2*YMax,1.2*YMax)
|
619 |
else:
|
620 |
Comparison.GetYaxis().SetRangeUser(-1.2*YMin,1.2*YMin)
|
621 |
|
622 |
Comparison.GetYaxis().SetNdivisions(205)
|
623 |
Comparison.Draw()
|
624 |
|
625 |
Canvas.Write()
|
626 |
if arguments.savePDFs:
|
627 |
pathToDirString = plainTextString(pathToDir)
|
628 |
Canvas.SaveAs(condor_dir+"/stacked_histograms_pdfs/"+pathToDirString+"/"+histogramName+".pdf")
|
629 |
|
630 |
|
631 |
##########################################################################################################################################
|
632 |
##########################################################################################################################################
|
633 |
##########################################################################################################################################
|
634 |
|
635 |
def MakeTwoDHist(pathToDir,histogramName):
|
636 |
blindData = False
|
637 |
# To blind histograms, define a list of histsToBlind in the localOptions.py file
|
638 |
try: # Use "try" in case histsToBlind does not exist
|
639 |
if histsToBlind:
|
640 |
for histToBlind in histsToBlind:
|
641 |
if histToBlind in histogramName:
|
642 |
print "Blinding data for histogram " + histogramName
|
643 |
blindData = True
|
644 |
except NameError:
|
645 |
time.sleep(0.000001) # Do nothing if histsToBlind does not exist
|
646 |
|
647 |
numBgMCSamples = 0
|
648 |
numDataSamples = 0
|
649 |
numSignalSamples = 0
|
650 |
|
651 |
|
652 |
HeaderLabel = TPaveLabel(header_x_left,header_y_bottom,header_x_right,header_y_top,HeaderText,"NDC")
|
653 |
HeaderLabel.SetTextAlign(32)
|
654 |
HeaderLabel.SetBorderSize(0)
|
655 |
HeaderLabel.SetFillColor(0)
|
656 |
HeaderLabel.SetFillStyle(0)
|
657 |
|
658 |
LumiLabel = TPaveLabel(topLeft_x_left,topLeft_y_bottom,topLeft_x_right,topLeft_y_top,LumiText,"NDC")
|
659 |
LumiLabel.SetBorderSize(0)
|
660 |
LumiLabel.SetFillColor(0)
|
661 |
LumiLabel.SetFillStyle(0)
|
662 |
|
663 |
NormLabel = TPaveLabel()
|
664 |
NormLabel.SetDrawOption("NDC")
|
665 |
NormLabel.SetX1NDC(topLeft_x_left)
|
666 |
NormLabel.SetX2NDC(topLeft_x_right)
|
667 |
|
668 |
NormLabel.SetBorderSize(0)
|
669 |
NormLabel.SetFillColor(0)
|
670 |
NormLabel.SetFillStyle(0)
|
671 |
|
672 |
NormText = ""
|
673 |
if arguments.normalizeToUnitArea:
|
674 |
NormText = "Scaled to unit area"
|
675 |
elif arguments.normalizeToData:
|
676 |
NormText = "MC scaled to data"
|
677 |
NormLabel.SetLabel(NormText)
|
678 |
|
679 |
BgMCLegend = TLegend(0.76,0.65,0.99,0.9)
|
680 |
BgMCLegend.AddEntry (0, "Data & Bkgd. MC", "H").SetTextFont (62)
|
681 |
BgMCLegend.SetBorderSize(0)
|
682 |
BgMCLegend.SetFillColor(0)
|
683 |
BgMCLegend.SetFillStyle(0)
|
684 |
SignalMCLegend = TLegend(0.76,0.135,0.99,0.377)
|
685 |
SignalMCLegend.AddEntry (0, "Signal MC", "H").SetTextFont (62)
|
686 |
SignalMCLegend.SetBorderSize(0)
|
687 |
SignalMCLegend.SetFillColor(0)
|
688 |
SignalMCLegend.SetFillStyle(0)
|
689 |
|
690 |
outputFile.cd(pathToDir)
|
691 |
Canvas = TCanvas(histogramName)
|
692 |
Canvas.SetRightMargin(0.2413793);
|
693 |
BgMCHistograms = []
|
694 |
SignalMCHistograms = []
|
695 |
DataHistograms = []
|
696 |
|
697 |
for sample in processed_datasets: # loop over different samples as listed in configurationOptions.py
|
698 |
dataset_file = "%s/%s.root" % (condor_dir,sample)
|
699 |
inputFile = TFile(dataset_file)
|
700 |
HistogramObj = inputFile.Get(pathToDir+"/"+histogramName)
|
701 |
if not HistogramObj:
|
702 |
print "WARNING: Could not find histogram " + pathToDir + "/" + histogramName + " in file " + dataset_file + ". Will skip it and continue."
|
703 |
continue
|
704 |
Histogram = HistogramObj.Clone()
|
705 |
Histogram.SetDirectory(0)
|
706 |
inputFile.Close()
|
707 |
if arguments.rebinFactor:
|
708 |
RebinFactor = int(arguments.rebinFactor)
|
709 |
#don't rebin histograms which will have less than 5 bins or any gen-matching histograms
|
710 |
if Histogram.GetNbinsX() >= RebinFactor*5 and Histogram.GetName().find("GenMatch") is -1:
|
711 |
Histogram.Rebin(RebinFactor)
|
712 |
xAxisLabel = Histogram.GetXaxis().GetTitle()
|
713 |
yAxisLabel = Histogram.GetYaxis().GetTitle()
|
714 |
if not arguments.makeFancy:
|
715 |
histoTitle = Histogram.GetTitle()
|
716 |
else:
|
717 |
histoTitle = ""
|
718 |
|
719 |
if( types[sample] == "bgMC"):
|
720 |
|
721 |
numBgMCSamples += 1
|
722 |
Histogram.SetMarkerColor(colors[sample])
|
723 |
Histogram.SetMarkerStyle(24)
|
724 |
Histogram.SetMarkerSize(1.2)
|
725 |
Histogram.SetFillColor(colors[sample])
|
726 |
BgMCLegend.AddEntry(Histogram,labels[sample],"P").SetTextFont (42)
|
727 |
BgMCHistograms.append(Histogram)
|
728 |
|
729 |
elif( types[sample] == "signalMC"):
|
730 |
|
731 |
numSignalSamples += 1
|
732 |
Histogram.SetMarkerColor(colors[sample])
|
733 |
Histogram.SetMarkerStyle(20)
|
734 |
Histogram.SetMarkerSize(1.2)
|
735 |
Histogram.SetFillColor(colors[sample])
|
736 |
BgMCLegend.AddEntry(Histogram,labels[sample],"P").SetTextFont (42)
|
737 |
# SignalMCLegend.AddEntry(Histogram,labels[sample],"P").SetTextFont (42)
|
738 |
SignalMCHistograms.append(Histogram)
|
739 |
|
740 |
elif( types[sample] == "data" and not blindData):
|
741 |
|
742 |
numDataSamples += 1
|
743 |
Histogram.SetMarkerColor(colors[sample])
|
744 |
Histogram.SetMarkerStyle(34)
|
745 |
Histogram.SetMarkerSize(1.2)
|
746 |
Histogram.SetFillColor(colors[sample])
|
747 |
BgMCLegend.AddEntry(Histogram,labels[sample],"P").SetTextFont (42)
|
748 |
DataHistograms.append(Histogram)
|
749 |
|
750 |
|
751 |
outputFile.cd(pathToDir)
|
752 |
|
753 |
if(numBgMCSamples is not 0):
|
754 |
BgMCHistograms[0].SetTitle(histoTitle)
|
755 |
BgMCHistograms[0].GetXaxis().SetTitle(xAxisLabel)
|
756 |
BgMCHistograms[0].GetYaxis().SetTitle(yAxisLabel)
|
757 |
BgMCHistograms[0].Draw()
|
758 |
for signalMCHist in SignalMCHistograms:
|
759 |
signalMCHist.Draw("SAME")
|
760 |
for dataHist in DataHistograms:
|
761 |
dataHist.Draw("SAME")
|
762 |
|
763 |
elif(numSignalSamples is not 0):
|
764 |
SignalMCHistograms[0].SetTitle(histoTitle)
|
765 |
SignalMCHistograms[0].Draw()
|
766 |
SignalMCHistograms[0].GetXaxis().SetTitle(xAxisLabel)
|
767 |
SignalMCHistograms[0].GetYaxis().SetTitle(yAxisLabel)
|
768 |
for signalMCHist in SignalMCHistograms:
|
769 |
if(signalMCHist is not SignalMCHistograms[0]):
|
770 |
signalMCHist.Draw("SAME")
|
771 |
for dataHist in DataHistograms:
|
772 |
dataHist.Draw("SAME")
|
773 |
|
774 |
elif(numDataSamples is not 0):
|
775 |
DataHistograms[0].SetTitle(histoTitle)
|
776 |
DataHistograms[0].GetXaxis().SetTitle(xAxisLabel)
|
777 |
DataHistograms[0].GetYaxis().SetTitle(yAxisLabel)
|
778 |
DataHistograms[0].Draw()
|
779 |
for dataHist in DataHistograms:
|
780 |
if(dataHist is not DataHistograms[0]):
|
781 |
dataHist.Draw("SAME")
|
782 |
|
783 |
|
784 |
|
785 |
# Deciding which text labels to draw and drawing them
|
786 |
drawLumiLabel = False
|
787 |
drawNormLabel = False
|
788 |
offsetNormLabel = False
|
789 |
drawHeaderLabel = False
|
790 |
|
791 |
if not arguments.normalizeToUnitArea or numDataSamples > 0: #don't draw the lumi label if there's no data and it's scaled to unit area
|
792 |
drawLumiLabel = True
|
793 |
#move the normalization label down before drawing if we drew the lumi. label
|
794 |
offsetNormLabel = True
|
795 |
if arguments.normalizeToUnitArea or arguments.normalizeToData:
|
796 |
drawNormLabel = True
|
797 |
if arguments.makeFancy:
|
798 |
drawHeaderLabel = True
|
799 |
drawLumiLabel = False
|
800 |
|
801 |
#now that flags are set, draw the appropriate labels
|
802 |
|
803 |
if drawLumiLabel:
|
804 |
LumiLabel.Draw()
|
805 |
|
806 |
if drawNormLabel:
|
807 |
if offsetNormLabel:
|
808 |
NormLabel.SetY1NDC(topLeft_y_bottom-topLeft_y_offset)
|
809 |
NormLabel.SetY2NDC(topLeft_y_top-topLeft_y_offset)
|
810 |
else:
|
811 |
NormLabel.SetY1NDC(topLeft_y_bottom)
|
812 |
NormLabel.SetY2NDC(topLeft_y_top)
|
813 |
NormLabel.Draw()
|
814 |
|
815 |
if drawHeaderLabel:
|
816 |
HeaderLabel.Draw()
|
817 |
|
818 |
|
819 |
|
820 |
|
821 |
if(numBgMCSamples is not 0 or numDataSamples is not 0):
|
822 |
BgMCLegend.Draw()
|
823 |
if(numSignalSamples is not 0):
|
824 |
SignalMCLegend.Draw()
|
825 |
|
826 |
Canvas.Write()
|
827 |
|
828 |
|
829 |
|
830 |
|
831 |
##########################################################################################################################################
|
832 |
##########################################################################################################################################
|
833 |
##########################################################################################################################################
|
834 |
|
835 |
processed_datasets = []
|
836 |
|
837 |
condor_dir = set_condor_output_dir(arguments)
|
838 |
|
839 |
#### check which input datasets have valid output files
|
840 |
for sample in datasets:
|
841 |
fileName = condor_dir + "/" + sample + ".root"
|
842 |
if not os.path.exists(fileName):
|
843 |
continue
|
844 |
testFile = TFile(fileName)
|
845 |
if testFile.IsZombie() or not testFile.GetNkeys():
|
846 |
continue
|
847 |
processed_datasets.append(sample)
|
848 |
|
849 |
if len(processed_datasets) is 0:
|
850 |
sys.exit("No datasets have been processed")
|
851 |
|
852 |
|
853 |
#### make output file
|
854 |
outputFileName = "stacked_histograms.root"
|
855 |
if arguments.outputFileName:
|
856 |
outputFileName = arguments.outputFileName
|
857 |
|
858 |
outputFile = TFile(condor_dir + "/" + outputFileName, "RECREATE")
|
859 |
|
860 |
|
861 |
|
862 |
#### use the first input file as a template and make stacked versions of all its histograms
|
863 |
inputFile = TFile(condor_dir + "/" + processed_datasets[0] + ".root")
|
864 |
inputFile.cd()
|
865 |
outputFile.cd()
|
866 |
|
867 |
if arguments.savePDFs:
|
868 |
os.system("rm -rf %s/stacked_histograms_pdfs" % (condor_dir))
|
869 |
os.system("mkdir %s/stacked_histograms_pdfs" % (condor_dir))
|
870 |
|
871 |
|
872 |
#get root directory in the first layer, generally "OSUAnalysis"
|
873 |
for key in inputFile.GetListOfKeys():
|
874 |
if (key.GetClassName() != "TDirectoryFile"):
|
875 |
continue
|
876 |
rootDirectory = key.GetName()
|
877 |
outputFile.mkdir(rootDirectory)
|
878 |
if arguments.savePDFs:
|
879 |
os.system("mkdir %s/stacked_histograms_pdfs/%s" % (condor_dir,plainTextString(rootDirectory)))
|
880 |
|
881 |
#cd to root directory and look for histograms
|
882 |
inputFile.cd(rootDirectory)
|
883 |
for key2 in gDirectory.GetListOfKeys():
|
884 |
|
885 |
if re.match ('TH1', key2.GetClassName()): # found a 1-D histogram
|
886 |
MakeOneDHist(rootDirectory,key2.GetName())
|
887 |
elif re.match ('TH2', key2.GetClassName()) and arguments.draw2DPlots: # found a 2-D histogram
|
888 |
MakeTwoDHist(rootDirectory,key2.GetName())
|
889 |
|
890 |
elif (key2.GetClassName() == "TDirectoryFile"): # found a directory, cd there and look for histograms
|
891 |
level2Directory = rootDirectory+"/"+key2.GetName()
|
892 |
|
893 |
#make a corresponding directory in the output file
|
894 |
outputFile.cd(rootDirectory)
|
895 |
gDirectory.mkdir(key2.GetName())
|
896 |
if arguments.savePDFs:
|
897 |
os.system("mkdir %s/stacked_histograms_pdfs/%s" % (condor_dir,plainTextString(level2Directory)))
|
898 |
|
899 |
#####################################################
|
900 |
### This layer is typically the "channels" layer ###
|
901 |
#####################################################
|
902 |
|
903 |
inputFile.cd(level2Directory)
|
904 |
for key3 in gDirectory.GetListOfKeys():
|
905 |
if re.match ('TH1', key3.GetClassName()): # found a 1-D histogram
|
906 |
MakeOneDHist(level2Directory,key3.GetName())
|
907 |
elif re.match ('TH2', key3.GetClassName()) and arguments.draw2DPlots: # found a 2-D histogram
|
908 |
MakeTwoDHist(level2Directory,key3.GetName())
|
909 |
|
910 |
elif (key3.GetClassName() == "TDirectoryFile"): # found a directory, cd there and look for histograms
|
911 |
level3Directory = level2Directory+"/"+key3.GetName()
|
912 |
|
913 |
#make a corresponding directory in the output file
|
914 |
outputFile.cd(level2Directory)
|
915 |
gDirectory.mkdir(key3.GetName())
|
916 |
if arguments.savePDFs:
|
917 |
os.system("mkdir %s/stacked_histograms_pdfs/%s" % (condor_dir,plainTextString(level3Directory)))
|
918 |
|
919 |
#################################################
|
920 |
### This layer is typically the "cuts" layer ###
|
921 |
#################################################
|
922 |
|
923 |
inputFile.cd(level3Directory)
|
924 |
for key3 in gDirectory.GetListOfKeys():
|
925 |
if re.match ('TH1', key3.GetClassName()): # found a 1-D histogram
|
926 |
MakeOneDHist(level3Directory,key3.GetName())
|
927 |
elif re.match ('TH2', key3.GetClassName()) and arguments.draw2DPlots: # found a 2-D histogram
|
928 |
MakeTwoDHist(level3Directory,key3.GetName())
|
929 |
|
930 |
|
931 |
outputFile.Close()
|