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