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

# 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 lantonel 1.2 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 lantonel 1.1
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 lantonel 1.3 ratio.GetXaxis().SetLabelOffset(0.03)
207 lantonel 1.1 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 lantonel 1.3 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 lantonel 1.1 else:
277     yAxisLabel = "Entries per bin (" + str(Histogram.GetXaxis().GetBinWidth(1)) + " width)"
278 lantonel 1.2 if arguments.normalizeToUnitArea:
279     yAxisLabel = yAxisLabel + " (Unit Area Norm.)"
280    
281 lantonel 1.1
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 lantonel 1.3 if colorIndex is len(colorList):
298 lantonel 1.1 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 lantonel 1.2 if arguments.normalizeToUnitArea and histogram.Integral() > 0:
317 lantonel 1.1 histogram.Scale(1./histogram.Integral())
318    
319 lantonel 1.2
320 lantonel 1.1 ### 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 lantonel 1.2 currentMax = histogram.GetMaximum() + histogram.GetBinError(histogram.GetMaximumBin())
330     if(currentMax > finalMax):
331     finalMax = currentMax
332     finalMax = 1.5*finalMax
333 lantonel 1.1 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 lantonel 1.2
387     x_left = 0.1677852
388     x_right = 0.9647651
389     y_min = 0.6765734
390 lantonel 1.1 y_max = 0.9
391    
392     Legend.SetX1NDC(x_left)
393 lantonel 1.2 Legend.SetY1NDC(y_min)
394 lantonel 1.1 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()