ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/OSUT3Analysis/Configuration/scripts/makeComparisonPlots.py
Revision: 1.3
Committed: Fri Aug 2 11:25:53 2013 UTC (11 years, 9 months ago) by lantonel
Content type: text/x-python
Branch: MAIN
CVS Tags: HEAD
Changes since 1.2: +7 -5 lines
Log Message:
fixed axis label

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.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()