1 |
#!/usr/bin/env python
|
2 |
import sys
|
3 |
import os
|
4 |
import re
|
5 |
from math import *
|
6 |
from array import *
|
7 |
from decimal import *
|
8 |
from optparse import OptionParser
|
9 |
from OSUT3Analysis.Configuration.configurationOptions import *
|
10 |
from OSUT3Analysis.Configuration.processingUtilities import *
|
11 |
from OSUT3Analysis.Configuration.formattingUtilities import *
|
12 |
|
13 |
|
14 |
|
15 |
### parse the command-line options
|
16 |
|
17 |
parser = OptionParser()
|
18 |
parser = set_commandline_arguments(parser)
|
19 |
|
20 |
parser.remove_option("-c")
|
21 |
parser.remove_option("-e")
|
22 |
parser.remove_option("-n")
|
23 |
parser.remove_option("--2D")
|
24 |
parser.remove_option("-y")
|
25 |
|
26 |
parser.add_option("-f", "--fancy", action="store_true", dest="makeFancy", default=False,
|
27 |
help="removes the title and replaces it with the official CMS plot heading")
|
28 |
parser.add_option("--dontRebinRatio", action="store_true", dest="dontRebinRatio", default=False,
|
29 |
help="don't do the rebinning of ratio plots")
|
30 |
|
31 |
parser.add_option("--ylog", action="store_true", dest="setLogY", default=False,
|
32 |
help="Set logarithmic scale on vertical axis on all plots")
|
33 |
parser.add_option("--ymin", dest="setYMin",
|
34 |
help="Minimum of y axis")
|
35 |
parser.add_option("--ymax", dest="setYMax",
|
36 |
help="Maximum of y axis")
|
37 |
parser.add_option("--hist", action="store_true", dest="plot_hist", default=False,
|
38 |
help="plot as hollow histograms instead of error bar crosses")
|
39 |
parser.add_option("--line-width", dest="line_width",
|
40 |
help="set line width (default is 2)")
|
41 |
|
42 |
(arguments, args) = parser.parse_args()
|
43 |
|
44 |
|
45 |
if arguments.localConfig:
|
46 |
sys.path.append(os.getcwd())
|
47 |
exec("from " + arguments.localConfig.rstrip('.py') + " import *")
|
48 |
|
49 |
|
50 |
|
51 |
#### deal with conflicting arguments
|
52 |
|
53 |
if arguments.makeRatioPlots or arguments.makeDiffPlots:
|
54 |
if len(input_sources) is not 2:
|
55 |
print "You need to have exactly two input sources to produce ratio or difference plots, turning them off"
|
56 |
arguments.makeRatioPlots = False
|
57 |
arguments.makeDiffPlots = False
|
58 |
|
59 |
if arguments.makeRatioPlots and arguments.makeDiffPlots:
|
60 |
print "You have requested both ratio and difference plots. Will make just ratio plots instead"
|
61 |
arguments.makeRatioPlots = False
|
62 |
|
63 |
|
64 |
|
65 |
from ROOT import TFile, gROOT, gStyle, gDirectory, TStyle, THStack, TH1F, TCanvas, TString, TLegend, TLegendEntry, THStack, TIter, TKey, TPaveLabel, gPad
|
66 |
|
67 |
|
68 |
### setting ROOT options so our plots will look awesome and everyone will love us
|
69 |
|
70 |
gROOT.SetBatch()
|
71 |
gStyle.SetOptStat(0)
|
72 |
gStyle.SetCanvasBorderMode(0)
|
73 |
gStyle.SetPadBorderMode(0)
|
74 |
gStyle.SetPadColor(0)
|
75 |
gStyle.SetCanvasColor(0)
|
76 |
gStyle.SetTextFont(42)
|
77 |
gStyle.SetCanvasDefH(600)
|
78 |
gStyle.SetCanvasDefW(600)
|
79 |
gStyle.SetCanvasDefX(0)
|
80 |
gStyle.SetCanvasDefY(0)
|
81 |
gStyle.SetPadTopMargin(0.07)
|
82 |
gStyle.SetPadBottomMargin(0.13)
|
83 |
gStyle.SetPadLeftMargin(0.15)
|
84 |
gStyle.SetPadRightMargin(0.05)
|
85 |
gStyle.SetTitleColor(1, "XYZ")
|
86 |
gStyle.SetTitleFont(42, "XYZ")
|
87 |
gStyle.SetTitleSize(0.04, "XYZ")
|
88 |
gStyle.SetTitleXOffset(1.1)
|
89 |
gStyle.SetTitleYOffset(2)
|
90 |
gStyle.SetTextAlign(12)
|
91 |
gStyle.SetLabelColor(1, "XYZ")
|
92 |
gStyle.SetLabelFont(42, "XYZ")
|
93 |
gStyle.SetLabelOffset(0.007, "XYZ")
|
94 |
gStyle.SetLabelSize(0.04, "XYZ")
|
95 |
gStyle.SetAxisColor(1, "XYZ")
|
96 |
gStyle.SetStripDecimals(True)
|
97 |
gStyle.SetTickLength(0.03, "XYZ")
|
98 |
gStyle.SetNdivisions(510, "XYZ")
|
99 |
gStyle.SetPadTickX(1)
|
100 |
gStyle.SetPadTickY(1)
|
101 |
gROOT.ForceStyle()
|
102 |
|
103 |
line_width = 2
|
104 |
if arguments.line_width:
|
105 |
line_width = arguments.line_width
|
106 |
|
107 |
#set the text for the luminosity label
|
108 |
if(intLumi < 1000.):
|
109 |
LumiText = "L_{int} = " + str(intLumi) + " pb^{-1}"
|
110 |
LumiText = "L_{int} = " + str.format('{0:.1f}', LumiInPb) + " pb^{-1}"
|
111 |
else:
|
112 |
LumiInFb = intLumi/1000.
|
113 |
LumiText = "L_{int} = " + str.format('{0:.1f}', LumiInFb) + " fb^{-1}"
|
114 |
|
115 |
#bestest place for lumi. label, in top left corner
|
116 |
topLeft_x_left = 0.1375839
|
117 |
topLeft_x_right = 0.4580537
|
118 |
topLeft_y_bottom = 0.8479021
|
119 |
topLeft_y_top = 0.9475524
|
120 |
topLeft_y_offset = 0.035
|
121 |
|
122 |
#set the text for the fancy heading
|
123 |
HeaderText = "CMS Preliminary: " + LumiText + " at #sqrt{s} = 8 TeV"
|
124 |
|
125 |
#position for header
|
126 |
header_x_left = 0.2181208
|
127 |
header_x_right = 0.9562937
|
128 |
header_y_bottom = 0.9479866
|
129 |
header_y_top = 0.9947552
|
130 |
|
131 |
|
132 |
colors = {
|
133 |
'black' : 1,
|
134 |
'red' : 623,
|
135 |
'green' : 407,
|
136 |
'purple' : 871,
|
137 |
'blue' : 591,
|
138 |
'yellow' : 393,
|
139 |
}
|
140 |
|
141 |
colorList = [
|
142 |
'black',
|
143 |
'red',
|
144 |
'green',
|
145 |
'purple',
|
146 |
'blue',
|
147 |
'yellow',
|
148 |
]
|
149 |
|
150 |
|
151 |
markers = {
|
152 |
'circle' : 20,
|
153 |
'square' : 21,
|
154 |
'triangle' : 22,
|
155 |
}
|
156 |
|
157 |
fills = {
|
158 |
'solid' : 0,
|
159 |
'hollow' : 4,
|
160 |
}
|
161 |
|
162 |
##########################################################################################################################################
|
163 |
##########################################################################################################################################
|
164 |
##########################################################################################################################################
|
165 |
|
166 |
# some fancy-ass code from Andrzej Zuranski to merge bins in the ratio plot until the error goes below some threshold
|
167 |
def ratioHistogram( dataHist, mcHist, relErrMax=0.10):
|
168 |
|
169 |
def groupR(group):
|
170 |
Data,MC = [float(sum(hist.GetBinContent(i) for i in group)) for hist in [dataHist,mcHist]]
|
171 |
return (Data-MC)/MC if MC else 0
|
172 |
|
173 |
def groupErr(group):
|
174 |
Data,MC = [float(sum(hist.GetBinContent(i) for i in group)) for hist in [dataHist,mcHist]]
|
175 |
dataErr2,mcErr2 = [sum(hist.GetBinError(i)**2 for i in group) for hist in [dataHist,mcHist]]
|
176 |
if MC > 0 and Data > 0 and Data != MC:
|
177 |
return abs(math.sqrt( (dataErr2+mcErr2)/(Data-MC)**2 + mcErr2/MC**2 ) * (Data-MC)/MC)
|
178 |
else:
|
179 |
return 0
|
180 |
|
181 |
|
182 |
def regroup(groups):
|
183 |
err,iG = max( (groupErr(g),groups.index(g)) for g in groups )
|
184 |
if err < relErrMax or len(groups)<3 : return groups
|
185 |
iH = max( [iG-1,iG+1], key = lambda i: groupErr(groups[i]) if 0<=i<len(groups) else -1 )
|
186 |
iLo,iHi = sorted([iG,iH])
|
187 |
return regroup(groups[:iLo] + [groups[iLo]+groups[iHi]] + groups[iHi+1:])
|
188 |
|
189 |
#don't rebin the histograms of the number of a given object (except for the pileup ones)
|
190 |
if ((dataHist.GetName().find("num") is not -1 and dataHist.GetName().find("Primaryvertexs") is -1) or
|
191 |
dataHist.GetName().find("CutFlow") is not -1 or
|
192 |
dataHist.GetName().find("GenMatch") is not -1 or
|
193 |
arguments.dontRebinRatio):
|
194 |
ratio = dataHist.Clone()
|
195 |
ratio.Add(mcHist,-1)
|
196 |
ratio.Divide(mcHist)
|
197 |
ratio.SetTitle("")
|
198 |
else:
|
199 |
groups = regroup( [(i,) for i in range(1,1+dataHist.GetNbinsX())] )
|
200 |
ratio = TH1F("ratio","",len(groups), array('d', [dataHist.GetBinLowEdge(min(g)) for g in groups ] + [dataHist.GetXaxis().GetBinUpEdge(dataHist.GetNbinsX())]) )
|
201 |
for i,g in enumerate(groups) :
|
202 |
ratio.SetBinContent(i+1,groupR(g))
|
203 |
ratio.SetBinError(i+1,groupErr(g))
|
204 |
|
205 |
ratio.GetYaxis().SetTitle("#frac{hist1-hist2}{hist2}")
|
206 |
ratio.GetXaxis().SetLabelOffset(0.03)
|
207 |
ratio.SetLineColor(1)
|
208 |
ratio.SetMarkerColor(1)
|
209 |
ratio.SetLineWidth(2)
|
210 |
return ratio
|
211 |
|
212 |
##########################################################################################################################################
|
213 |
##########################################################################################################################################
|
214 |
##########################################################################################################################################
|
215 |
|
216 |
|
217 |
def MakeOneDHist(histogramName):
|
218 |
|
219 |
HeaderLabel = TPaveLabel(header_x_left,header_y_bottom,header_x_right,header_y_top,HeaderText,"NDC")
|
220 |
HeaderLabel.SetTextAlign(32)
|
221 |
HeaderLabel.SetBorderSize(0)
|
222 |
HeaderLabel.SetFillColor(0)
|
223 |
HeaderLabel.SetFillStyle(0)
|
224 |
|
225 |
LumiLabel = TPaveLabel(topLeft_x_left,topLeft_y_bottom,topLeft_x_right,topLeft_y_top,LumiText,"NDC")
|
226 |
LumiLabel.SetBorderSize(0)
|
227 |
LumiLabel.SetFillColor(0)
|
228 |
LumiLabel.SetFillStyle(0)
|
229 |
|
230 |
NormLabel = TPaveLabel()
|
231 |
NormLabel.SetDrawOption("NDC")
|
232 |
NormLabel.SetX1NDC(topLeft_x_left)
|
233 |
NormLabel.SetX2NDC(topLeft_x_right)
|
234 |
|
235 |
NormLabel.SetBorderSize(0)
|
236 |
NormLabel.SetFillColor(0)
|
237 |
NormLabel.SetFillStyle(0)
|
238 |
|
239 |
NormText = ""
|
240 |
if arguments.normalizeToUnitArea:
|
241 |
NormText = "Scaled to unit area"
|
242 |
|
243 |
Legend = TLegend()
|
244 |
Legend.SetBorderSize(0)
|
245 |
Legend.SetFillColor(0)
|
246 |
Legend.SetFillStyle(0)
|
247 |
|
248 |
Canvas = TCanvas(histogramName)
|
249 |
Histograms = []
|
250 |
LegendEntries = []
|
251 |
|
252 |
colorIndex = 0
|
253 |
|
254 |
for source in input_sources: # loop over different input sources in config file
|
255 |
dataset_file = "condor/%s/%s.root" % (source['condor_dir'],source['dataset'])
|
256 |
inputFile = TFile(dataset_file)
|
257 |
HistogramObj = inputFile.Get("OSUAnalysis/" + source['channel'] + "/" +histogramName)
|
258 |
if not HistogramObj:
|
259 |
print "WARNING: Could not find histogram " + source['channel'] + "/" + histogramName + " in file " + dataset_file + ". Will skip it and continue."
|
260 |
return
|
261 |
Histogram = HistogramObj.Clone()
|
262 |
Histogram.SetDirectory(0)
|
263 |
inputFile.Close()
|
264 |
if arguments.rebinFactor:
|
265 |
RebinFactor = int(arguments.rebinFactor)
|
266 |
#don't rebin histograms which will have less than 5 bins or any gen-matching histograms
|
267 |
if Histogram.GetNbinsX() >= RebinFactor*5 and Histogram.GetTitle().find("GenMatch") is -1:
|
268 |
Histogram.Rebin(RebinFactor)
|
269 |
|
270 |
xAxisLabel = Histogram.GetXaxis().GetTitle()
|
271 |
unitBeginIndex = xAxisLabel.find("[")
|
272 |
unitEndIndex = xAxisLabel.find("]")
|
273 |
|
274 |
if unitBeginIndex is not -1 and unitEndIndex is not -1: #x axis has a unit
|
275 |
yAxisLabel = "Entries / " + str(Histogram.GetXaxis().GetBinWidth(1)) + " " + xAxisLabel[unitBeginIndex+1:unitEndIndex]
|
276 |
else:
|
277 |
yAxisLabel = "Entries per bin (" + str(Histogram.GetXaxis().GetBinWidth(1)) + " width)"
|
278 |
if arguments.normalizeToUnitArea:
|
279 |
yAxisLabel = yAxisLabel + " (Unit Area Norm.)"
|
280 |
|
281 |
|
282 |
if not arguments.makeFancy:
|
283 |
fullTitle = Histogram.GetTitle()
|
284 |
splitTitle = fullTitle.split(":")
|
285 |
# print splitTitle
|
286 |
histoTitle = splitTitle[1].lstrip(" ")
|
287 |
else:
|
288 |
histoTitle = ""
|
289 |
|
290 |
if 'color' in source:
|
291 |
Histogram.SetMarkerColor(colors[source['color']])
|
292 |
Histogram.SetLineColor(colors[source['color']])
|
293 |
else:
|
294 |
Histogram.SetMarkerColor(colors[colorList[colorIndex]])
|
295 |
Histogram.SetLineColor(colors[colorList[colorIndex]])
|
296 |
colorIndex = colorIndex + 1
|
297 |
if colorIndex is len(colorList):
|
298 |
colorIndex = 0
|
299 |
|
300 |
markerStyle = 20
|
301 |
if 'marker' in source:
|
302 |
markerStyle = markers[source['marker']]
|
303 |
if 'fill' in source:
|
304 |
markerStyle = markerStyle + fills[source['fill']]
|
305 |
|
306 |
Histogram.SetMarkerStyle(markerStyle)
|
307 |
|
308 |
Histogram.SetLineWidth(line_width)
|
309 |
Histogram.SetFillStyle(0)
|
310 |
|
311 |
LegendEntries.append(source['legend_entry'])
|
312 |
Histograms.append(Histogram)
|
313 |
|
314 |
### scaling histograms as per user's specifications
|
315 |
for histogram in Histograms:
|
316 |
if arguments.normalizeToUnitArea and histogram.Integral() > 0:
|
317 |
histogram.Scale(1./histogram.Integral())
|
318 |
|
319 |
|
320 |
### formatting histograms and adding to legend
|
321 |
legendIndex = 0
|
322 |
for histogram in Histograms:
|
323 |
Legend.AddEntry(histogram,LegendEntries[legendIndex],"LEP")
|
324 |
legendIndex = legendIndex+1
|
325 |
|
326 |
### finding the maximum value of anything going on the canvas, so we know how to set the y-axis
|
327 |
finalMax = 0
|
328 |
for histogram in Histograms:
|
329 |
currentMax = histogram.GetMaximum() + histogram.GetBinError(histogram.GetMaximumBin())
|
330 |
if(currentMax > finalMax):
|
331 |
finalMax = currentMax
|
332 |
finalMax = 1.5*finalMax
|
333 |
if arguments.setYMax:
|
334 |
finalMax = float(arguments.setYMax)
|
335 |
|
336 |
### Drawing histograms to canvas
|
337 |
|
338 |
makeRatioPlots = arguments.makeRatioPlots
|
339 |
makeDiffPlots = arguments.makeDiffPlots
|
340 |
|
341 |
yAxisMin = 0.0001
|
342 |
if arguments.setYMin:
|
343 |
yAxisMin = float(arguments.setYMin)
|
344 |
|
345 |
if makeRatioPlots or makeDiffPlots:
|
346 |
Canvas.SetFillStyle(0)
|
347 |
Canvas.Divide(1,2)
|
348 |
Canvas.cd(1)
|
349 |
gPad.SetPad(0,0.25,1,1)
|
350 |
gPad.SetMargin(0.15,0.05,0.01,0.07)
|
351 |
gPad.SetFillStyle(0)
|
352 |
gPad.Update()
|
353 |
gPad.Draw()
|
354 |
if arguments.setLogY:
|
355 |
gPad.SetLogy()
|
356 |
Canvas.cd(2)
|
357 |
gPad.SetPad(0,0,1,0.25)
|
358 |
#format: gPad.SetMargin(l,r,b,t)
|
359 |
gPad.SetMargin(0.15,0.05,0.4,0.01)
|
360 |
gPad.SetFillStyle(0)
|
361 |
gPad.SetGridy(1)
|
362 |
gPad.Update()
|
363 |
gPad.Draw()
|
364 |
|
365 |
Canvas.cd(1)
|
366 |
|
367 |
histCounter = 0
|
368 |
plotting_options = ""
|
369 |
if arguments.plot_hist:
|
370 |
plotting_options = "HIST"
|
371 |
|
372 |
for histogram in Histograms:
|
373 |
histogram.SetTitle(histoTitle)
|
374 |
histogram.Draw(plotting_options)
|
375 |
histogram.GetXaxis().SetTitle(xAxisLabel)
|
376 |
histogram.GetYaxis().SetTitle(yAxisLabel)
|
377 |
histogram.SetMaximum(finalMax)
|
378 |
histogram.SetMinimum(yAxisMin)
|
379 |
if makeRatioPlots or makeDiffPlots:
|
380 |
histogram.GetXaxis().SetLabelSize(0)
|
381 |
if histCounter is 0:
|
382 |
plotting_options = plotting_options + " SAME"
|
383 |
histCounter = histCounter + 1
|
384 |
|
385 |
#legend coordinates, empirically determined :-)
|
386 |
|
387 |
x_left = 0.1677852
|
388 |
x_right = 0.9647651
|
389 |
y_min = 0.6765734
|
390 |
y_max = 0.9
|
391 |
|
392 |
Legend.SetX1NDC(x_left)
|
393 |
Legend.SetY1NDC(y_min)
|
394 |
Legend.SetX2NDC(x_right)
|
395 |
Legend.SetY2NDC(y_max)
|
396 |
Legend.Draw()
|
397 |
|
398 |
|
399 |
# Deciding which text labels to draw and drawing them
|
400 |
drawHeaderLabel = False
|
401 |
|
402 |
if arguments.makeFancy:
|
403 |
drawHeaderLabel = True
|
404 |
|
405 |
#now that flags are set, draw the appropriate labels
|
406 |
|
407 |
if drawHeaderLabel:
|
408 |
HeaderLabel.Draw()
|
409 |
|
410 |
|
411 |
|
412 |
|
413 |
#drawing the ratio or difference plot if requested
|
414 |
|
415 |
if makeRatioPlots or makeDiffPlots:
|
416 |
Canvas.cd(2)
|
417 |
if makeRatioPlots:
|
418 |
Comparison = ratioHistogram(Histograms[0],Histograms[1])
|
419 |
elif makeDiffPlots:
|
420 |
Comparison = Histograms[0].Clone("diff")
|
421 |
Comparison.Add(Histograms[1],-1)
|
422 |
Comparison.SetTitle("")
|
423 |
Comparison.GetYaxis().SetTitle("hist1-hist2")
|
424 |
Comparison.GetXaxis().SetTitle(xAxisLabel)
|
425 |
Comparison.GetYaxis().CenterTitle()
|
426 |
Comparison.GetYaxis().SetTitleSize(0.1)
|
427 |
Comparison.GetYaxis().SetTitleOffset(0.5)
|
428 |
Comparison.GetXaxis().SetTitleSize(0.15)
|
429 |
Comparison.GetYaxis().SetLabelSize(0.1)
|
430 |
Comparison.GetXaxis().SetLabelSize(0.15)
|
431 |
if makeRatioPlots:
|
432 |
RatioYRange = 1.15
|
433 |
if arguments.ratioYRange:
|
434 |
RatioYRange = float(arguments.ratioYRange)
|
435 |
Comparison.GetYaxis().SetRangeUser(-1*RatioYRange, RatioYRange)
|
436 |
elif makeDiffPlots:
|
437 |
YMax = Comparison.GetMaximum()
|
438 |
YMin = Comparison.GetMinimum()
|
439 |
if YMax <= 0 and YMin <= 0:
|
440 |
Comparison.GetYaxis().SetRangeUser(-1.2*YMin,0)
|
441 |
elif YMax >= 0 and YMin >= 0:
|
442 |
Comparison.GetYaxis().SetRangeUser(0,1.2*YMax)
|
443 |
else: #axis crosses y=0
|
444 |
if abs(YMax) > abs(YMin):
|
445 |
Comparison.GetYaxis().SetRangeUser(-1.2*YMax,1.2*YMax)
|
446 |
else:
|
447 |
Comparison.GetYaxis().SetRangeUser(-1.2*YMin,1.2*YMin)
|
448 |
|
449 |
Comparison.GetYaxis().SetNdivisions(205)
|
450 |
Comparison.Draw()
|
451 |
|
452 |
outputFile.cd()
|
453 |
Canvas.Write()
|
454 |
|
455 |
if arguments.savePDFs:
|
456 |
Canvas.SaveAs("comparison_histograms_pdfs/"+histogramName+".pdf")
|
457 |
|
458 |
|
459 |
##########################################################################################################################################
|
460 |
##########################################################################################################################################
|
461 |
##########################################################################################################################################
|
462 |
|
463 |
|
464 |
#### make output file
|
465 |
outputFileName = "comparison_histograms.root"
|
466 |
if arguments.outputFileName:
|
467 |
outputFileName = arguments.outputFileName
|
468 |
|
469 |
outputFile = TFile(outputFileName, "RECREATE")
|
470 |
|
471 |
first_input = input_sources[0]
|
472 |
|
473 |
#### use the first input file as a template and make comparison versions of all its histograms
|
474 |
testFile = TFile("condor/" + first_input['condor_dir'] + "/" + first_input['dataset'] + ".root")
|
475 |
testFile.cd("OSUAnalysis/" + first_input['channel'])
|
476 |
|
477 |
if arguments.savePDFs:
|
478 |
os.system("rm -rf comparison_histograms_pdfs")
|
479 |
os.system("mkdir comparison_histograms_pdfs")
|
480 |
|
481 |
for key in gDirectory.GetListOfKeys():
|
482 |
if re.match ('TH1', key.GetClassName()): #found a 1D histogram
|
483 |
MakeOneDHist(key.GetName())
|
484 |
|
485 |
testFile.Close()
|
486 |
outputFile.Close()
|