ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/OSUT3Analysis/Configuration/scripts/makeComparisonPlots.py
Revision: 1.1
Committed: Fri Jul 12 14:23:47 2013 UTC (11 years, 9 months ago) by lantonel
Content type: text/x-python
Branch: MAIN
CVS Tags: V02-03-02
Log Message:
first commit of a script very similar to makePlots.py, but for plotting several channels in the same canvas instead of several datasets.  sample config file in AnaTools/test/sampleComparisonPlotConfig.py

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