ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/OSUT3Analysis/Configuration/scripts/makeComparisonPlots.py
Revision: 1.2
Committed: Tue Jul 16 14:33:28 2013 UTC (11 years, 9 months ago) by lantonel
Content type: text/x-python
Branch: MAIN
Changes since 1.1: +20 -32 lines
Log Message:
fixed a bug in the ratio-making, and widened the legend (moved the unit area label to the y-axis)

File Contents

# Content
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.SetLineColor(1)
207 ratio.SetMarkerColor(1)
208 ratio.SetLineWidth(2)
209 return ratio
210
211 ##########################################################################################################################################
212 ##########################################################################################################################################
213 ##########################################################################################################################################
214
215
216 def MakeOneDHist(histogramName):
217
218 HeaderLabel = TPaveLabel(header_x_left,header_y_bottom,header_x_right,header_y_top,HeaderText,"NDC")
219 HeaderLabel.SetTextAlign(32)
220 HeaderLabel.SetBorderSize(0)
221 HeaderLabel.SetFillColor(0)
222 HeaderLabel.SetFillStyle(0)
223
224 LumiLabel = TPaveLabel(topLeft_x_left,topLeft_y_bottom,topLeft_x_right,topLeft_y_top,LumiText,"NDC")
225 LumiLabel.SetBorderSize(0)
226 LumiLabel.SetFillColor(0)
227 LumiLabel.SetFillStyle(0)
228
229 NormLabel = TPaveLabel()
230 NormLabel.SetDrawOption("NDC")
231 NormLabel.SetX1NDC(topLeft_x_left)
232 NormLabel.SetX2NDC(topLeft_x_right)
233
234 NormLabel.SetBorderSize(0)
235 NormLabel.SetFillColor(0)
236 NormLabel.SetFillStyle(0)
237
238 NormText = ""
239 if arguments.normalizeToUnitArea:
240 NormText = "Scaled to unit area"
241
242 Legend = TLegend()
243 Legend.SetBorderSize(0)
244 Legend.SetFillColor(0)
245 Legend.SetFillStyle(0)
246
247 Canvas = TCanvas(histogramName)
248 Histograms = []
249 LegendEntries = []
250
251 colorIndex = 0
252
253 for source in input_sources: # loop over different input sources in config file
254 dataset_file = "condor/%s/%s.root" % (source['condor_dir'],source['dataset'])
255 inputFile = TFile(dataset_file)
256 HistogramObj = inputFile.Get("OSUAnalysis/" + source['channel'] + "/" +histogramName)
257 if not HistogramObj:
258 print "WARNING: Could not find histogram " + source['channel'] + "/" + histogramName + " in file " + dataset_file + ". Will skip it and continue."
259 return
260 Histogram = HistogramObj.Clone()
261 Histogram.SetDirectory(0)
262 inputFile.Close()
263 if arguments.rebinFactor:
264 RebinFactor = int(arguments.rebinFactor)
265 #don't rebin histograms which will have less than 5 bins or any gen-matching histograms
266 if Histogram.GetNbinsX() >= RebinFactor*5 and Histogram.GetTitle().find("GenMatch") is -1:
267 Histogram.Rebin(RebinFactor)
268
269 xAxisLabel = Histogram.GetXaxis().GetTitle()
270 unitIndex = xAxisLabel.find("[")
271 yAxisLabel = ""
272 if unitIndex is not -1: #x axis has a unit
273 yAxisLabel = "Entries / " + str(Histogram.GetXaxis().GetBinWidth(1)) + " " + xAxisLabel[unitIndex+1:-1]
274 else:
275 yAxisLabel = "Entries per bin (" + str(Histogram.GetXaxis().GetBinWidth(1)) + " width)"
276 if arguments.normalizeToUnitArea:
277 yAxisLabel = yAxisLabel + " (Unit Area Norm.)"
278
279
280 if not arguments.makeFancy:
281 fullTitle = Histogram.GetTitle()
282 splitTitle = fullTitle.split(":")
283 # print splitTitle
284 histoTitle = splitTitle[1].lstrip(" ")
285 else:
286 histoTitle = ""
287
288 if 'color' in source:
289 Histogram.SetMarkerColor(colors[source['color']])
290 Histogram.SetLineColor(colors[source['color']])
291 else:
292 Histogram.SetMarkerColor(colors[colorList[colorIndex]])
293 Histogram.SetLineColor(colors[colorList[colorIndex]])
294 colorIndex = colorIndex + 1
295 if colorIndex is len(colorList)-1:
296 colorIndex = 0
297
298 markerStyle = 20
299 if 'marker' in source:
300 markerStyle = markers[source['marker']]
301 if 'fill' in source:
302 markerStyle = markerStyle + fills[source['fill']]
303
304 Histogram.SetMarkerStyle(markerStyle)
305
306 Histogram.SetLineWidth(line_width)
307 Histogram.SetFillStyle(0)
308
309 LegendEntries.append(source['legend_entry'])
310 Histograms.append(Histogram)
311
312 ### scaling histograms as per user's specifications
313 for histogram in Histograms:
314 if arguments.normalizeToUnitArea and histogram.Integral() > 0:
315 histogram.Scale(1./histogram.Integral())
316
317
318 ### formatting histograms and adding to legend
319 legendIndex = 0
320 for histogram in Histograms:
321 Legend.AddEntry(histogram,LegendEntries[legendIndex],"LEP")
322 legendIndex = legendIndex+1
323
324 ### finding the maximum value of anything going on the canvas, so we know how to set the y-axis
325 finalMax = 0
326 for histogram in Histograms:
327 currentMax = histogram.GetMaximum() + histogram.GetBinError(histogram.GetMaximumBin())
328 if(currentMax > finalMax):
329 finalMax = currentMax
330 finalMax = 1.5*finalMax
331 if arguments.setYMax:
332 finalMax = float(arguments.setYMax)
333
334 ### Drawing histograms to canvas
335
336 makeRatioPlots = arguments.makeRatioPlots
337 makeDiffPlots = arguments.makeDiffPlots
338
339 yAxisMin = 0.0001
340 if arguments.setYMin:
341 yAxisMin = float(arguments.setYMin)
342
343 if makeRatioPlots or makeDiffPlots:
344 Canvas.SetFillStyle(0)
345 Canvas.Divide(1,2)
346 Canvas.cd(1)
347 gPad.SetPad(0,0.25,1,1)
348 gPad.SetMargin(0.15,0.05,0.01,0.07)
349 gPad.SetFillStyle(0)
350 gPad.Update()
351 gPad.Draw()
352 if arguments.setLogY:
353 gPad.SetLogy()
354 Canvas.cd(2)
355 gPad.SetPad(0,0,1,0.25)
356 #format: gPad.SetMargin(l,r,b,t)
357 gPad.SetMargin(0.15,0.05,0.4,0.01)
358 gPad.SetFillStyle(0)
359 gPad.SetGridy(1)
360 gPad.Update()
361 gPad.Draw()
362
363 Canvas.cd(1)
364
365 histCounter = 0
366 plotting_options = ""
367 if arguments.plot_hist:
368 plotting_options = "HIST"
369
370 for histogram in Histograms:
371 histogram.SetTitle(histoTitle)
372 histogram.Draw(plotting_options)
373 histogram.GetXaxis().SetTitle(xAxisLabel)
374 histogram.GetYaxis().SetTitle(yAxisLabel)
375 histogram.SetMaximum(finalMax)
376 histogram.SetMinimum(yAxisMin)
377 if makeRatioPlots or makeDiffPlots:
378 histogram.GetXaxis().SetLabelSize(0)
379 if histCounter is 0:
380 plotting_options = plotting_options + " SAME"
381 histCounter = histCounter + 1
382
383 #legend coordinates, empirically determined :-)
384
385 x_left = 0.1677852
386 x_right = 0.9647651
387 y_min = 0.6765734
388 y_max = 0.9
389
390 Legend.SetX1NDC(x_left)
391 Legend.SetY1NDC(y_min)
392 Legend.SetX2NDC(x_right)
393 Legend.SetY2NDC(y_max)
394 Legend.Draw()
395
396
397 # Deciding which text labels to draw and drawing them
398 drawHeaderLabel = False
399
400 if arguments.makeFancy:
401 drawHeaderLabel = True
402
403 #now that flags are set, draw the appropriate labels
404
405 if drawHeaderLabel:
406 HeaderLabel.Draw()
407
408
409
410
411 #drawing the ratio or difference plot if requested
412
413 if makeRatioPlots or makeDiffPlots:
414 Canvas.cd(2)
415 if makeRatioPlots:
416 Comparison = ratioHistogram(Histograms[0],Histograms[1])
417 elif makeDiffPlots:
418 Comparison = Histograms[0].Clone("diff")
419 Comparison.Add(Histograms[1],-1)
420 Comparison.SetTitle("")
421 Comparison.GetYaxis().SetTitle("hist1-hist2")
422 Comparison.GetXaxis().SetTitle(xAxisLabel)
423 Comparison.GetYaxis().CenterTitle()
424 Comparison.GetYaxis().SetTitleSize(0.1)
425 Comparison.GetYaxis().SetTitleOffset(0.5)
426 Comparison.GetXaxis().SetTitleSize(0.15)
427 Comparison.GetYaxis().SetLabelSize(0.1)
428 Comparison.GetXaxis().SetLabelSize(0.15)
429 if makeRatioPlots:
430 RatioYRange = 1.15
431 if arguments.ratioYRange:
432 RatioYRange = float(arguments.ratioYRange)
433 Comparison.GetYaxis().SetRangeUser(-1*RatioYRange, RatioYRange)
434 elif makeDiffPlots:
435 YMax = Comparison.GetMaximum()
436 YMin = Comparison.GetMinimum()
437 if YMax <= 0 and YMin <= 0:
438 Comparison.GetYaxis().SetRangeUser(-1.2*YMin,0)
439 elif YMax >= 0 and YMin >= 0:
440 Comparison.GetYaxis().SetRangeUser(0,1.2*YMax)
441 else: #axis crosses y=0
442 if abs(YMax) > abs(YMin):
443 Comparison.GetYaxis().SetRangeUser(-1.2*YMax,1.2*YMax)
444 else:
445 Comparison.GetYaxis().SetRangeUser(-1.2*YMin,1.2*YMin)
446
447 Comparison.GetYaxis().SetNdivisions(205)
448 Comparison.Draw()
449
450 outputFile.cd()
451 Canvas.Write()
452
453 if arguments.savePDFs:
454 Canvas.SaveAs("comparison_histograms_pdfs/"+histogramName+".pdf")
455
456
457 ##########################################################################################################################################
458 ##########################################################################################################################################
459 ##########################################################################################################################################
460
461
462 #### make output file
463 outputFileName = "comparison_histograms.root"
464 if arguments.outputFileName:
465 outputFileName = arguments.outputFileName
466
467 outputFile = TFile(outputFileName, "RECREATE")
468
469 first_input = input_sources[0]
470
471 #### use the first input file as a template and make comparison versions of all its histograms
472 testFile = TFile("condor/" + first_input['condor_dir'] + "/" + first_input['dataset'] + ".root")
473 testFile.cd("OSUAnalysis/" + first_input['channel'])
474
475 if arguments.savePDFs:
476 os.system("rm -rf comparison_histograms_pdfs")
477 os.system("mkdir comparison_histograms_pdfs")
478
479 for key in gDirectory.GetListOfKeys():
480 if re.match ('TH1', key.GetClassName()): #found a 1D histogram
481 MakeOneDHist(key.GetName())
482
483 testFile.Close()
484 outputFile.Close()