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

# User Rev Content
1 lantonel 1.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()