ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/OSUT3Analysis/Configuration/scripts/makeEfficiencyPlots.py
Revision: 1.2
Committed: Fri Aug 2 11:27:46 2013 UTC (11 years, 9 months ago) by lantonel
Content type: text/x-python
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +3 -3 lines
Log Message:
fixed y 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 lantonel 1.2 from ROOT import TFile, gROOT, gStyle, gDirectory, TStyle, TH1F, TCanvas, TString, TLegend, TLegendEntry, TIter, TKey, TPaveLabel, gPad
66 lantonel 1.1
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     NumHistogramObj = inputFile.Get("OSUAnalysis/" + source['num_channel'] + "/" +histogramName)
258     DenHistogramObj = inputFile.Get("OSUAnalysis/" + source['den_channel'] + "/" +histogramName)
259     if not NumHistogramObj:
260     print "WARNING: Could not find histogram " + source['num_channel'] + "/" + histogramName + " in file " + dataset_file + ". Will skip it and continue."
261     return
262     if not DenHistogramObj:
263     print "WARNING: Could not find histogram " + source['den_channel'] + "/" + histogramName + " in file " + dataset_file + ". Will skip it and continue."
264     return
265    
266     Histogram = NumHistogramObj.Clone()
267     Histogram.SetDirectory(0)
268     DenHistogram = DenHistogramObj.Clone()
269     DenHistogram.SetDirectory(0)
270     inputFile.Close()
271    
272     if arguments.rebinFactor:
273     RebinFactor = int(arguments.rebinFactor)
274     #don't rebin histograms which will have less than 5 bins or any gen-matching histograms
275     if Histogram.GetNbinsX() >= RebinFactor*5 and Histogram.GetTitle().find("GenMatch") is -1:
276     Histogram.Rebin(RebinFactor)
277     DenHistogram.Rebin(RebinFactor)
278    
279     Histogram.Divide(DenHistogram)
280    
281     xAxisLabel = Histogram.GetXaxis().GetTitle()
282     unitBeginIndex = xAxisLabel.find("[")
283     unitEndIndex = xAxisLabel.find("]")
284    
285     if unitBeginIndex is not -1 and unitEndIndex is not -1: #x axis has a unit
286 lantonel 1.2 yAxisLabel = "#epsilon_{ " + cutName + "} (" + str(Histogram.GetXaxis().GetBinWidth(1)) + " " + xAxisLabel[unitBeginIndex+1:unitEndIndex] + " width)"
287 lantonel 1.1 else:
288 lantonel 1.2 yAxisLabel = "#epsilon_{ " + cutName + "} (" + str(Histogram.GetXaxis().GetBinWidth(1)) + " width)"
289 lantonel 1.1 if arguments.normalizeToUnitArea:
290     yAxisLabel = yAxisLabel + " (Unit Area Norm.)"
291    
292    
293     if not arguments.makeFancy:
294     fullTitle = Histogram.GetTitle()
295     splitTitle = fullTitle.split(":")
296     # print splitTitle
297     histoTitle = splitTitle[1].lstrip(" ")
298     else:
299     histoTitle = ""
300    
301     if 'color' in source:
302     Histogram.SetMarkerColor(colors[source['color']])
303     Histogram.SetLineColor(colors[source['color']])
304     else:
305     Histogram.SetMarkerColor(colors[colorList[colorIndex]])
306     Histogram.SetLineColor(colors[colorList[colorIndex]])
307     colorIndex = colorIndex + 1
308     if colorIndex is len(colorList):
309     colorIndex = 0
310    
311     markerStyle = 20
312     if 'marker' in source:
313     markerStyle = markers[source['marker']]
314     if 'fill' in source:
315     markerStyle = markerStyle + fills[source['fill']]
316    
317     Histogram.SetMarkerStyle(markerStyle)
318    
319     Histogram.SetLineWidth(line_width)
320     Histogram.SetFillStyle(0)
321    
322     LegendEntries.append(source['legend_entry'])
323     Histograms.append(Histogram)
324    
325     ### scaling histograms as per user's specifications
326     for histogram in Histograms:
327     if arguments.normalizeToUnitArea and histogram.Integral() > 0:
328     histogram.Scale(1./histogram.Integral())
329    
330    
331     ### formatting histograms and adding to legend
332     legendIndex = 0
333     for histogram in Histograms:
334     Legend.AddEntry(histogram,LegendEntries[legendIndex],"LEP")
335     legendIndex = legendIndex+1
336    
337     ### finding the maximum value of anything going on the canvas, so we know how to set the y-axis
338     finalMax = 0
339     for histogram in Histograms:
340     currentMax = histogram.GetMaximum() + histogram.GetBinError(histogram.GetMaximumBin())
341     if(currentMax > finalMax):
342     finalMax = currentMax
343     finalMax = 1.5*finalMax
344     if finalMax is 0:
345     finalMax = 1
346     if arguments.setYMax:
347     finalMax = float(arguments.setYMax)
348    
349     ### Drawing histograms to canvas
350    
351     makeRatioPlots = arguments.makeRatioPlots
352     makeDiffPlots = arguments.makeDiffPlots
353    
354     yAxisMin = 0.0001
355     if arguments.setYMin:
356     yAxisMin = float(arguments.setYMin)
357    
358     if makeRatioPlots or makeDiffPlots:
359     Canvas.SetFillStyle(0)
360     Canvas.Divide(1,2)
361     Canvas.cd(1)
362     gPad.SetPad(0,0.25,1,1)
363     gPad.SetMargin(0.15,0.05,0.01,0.07)
364     gPad.SetFillStyle(0)
365     gPad.Update()
366     gPad.Draw()
367     if arguments.setLogY:
368     gPad.SetLogy()
369     Canvas.cd(2)
370     gPad.SetPad(0,0,1,0.25)
371     #format: gPad.SetMargin(l,r,b,t)
372     gPad.SetMargin(0.15,0.05,0.4,0.01)
373     gPad.SetFillStyle(0)
374     gPad.SetGridy(1)
375     gPad.Update()
376     gPad.Draw()
377    
378     Canvas.cd(1)
379    
380     histCounter = 0
381     plotting_options = ""
382     if arguments.plot_hist:
383     plotting_options = "HIST"
384    
385     for histogram in Histograms:
386     histogram.SetTitle(histoTitle)
387     histogram.Draw(plotting_options)
388     histogram.GetXaxis().SetTitle(xAxisLabel)
389     histogram.GetYaxis().SetTitle(yAxisLabel)
390     histogram.SetMaximum(finalMax)
391     histogram.SetMinimum(yAxisMin)
392     if makeRatioPlots or makeDiffPlots:
393     histogram.GetXaxis().SetLabelSize(0)
394     if histCounter is 0:
395     plotting_options = plotting_options + " SAME"
396     histCounter = histCounter + 1
397    
398     #legend coordinates, empirically determined :-)
399    
400     x_left = 0.1677852
401     x_right = 0.9647651
402     y_min = 0.6765734
403     y_max = 0.9
404    
405     Legend.SetX1NDC(x_left)
406     Legend.SetY1NDC(y_min)
407     Legend.SetX2NDC(x_right)
408     Legend.SetY2NDC(y_max)
409     Legend.Draw()
410    
411    
412     # Deciding which text labels to draw and drawing them
413     drawHeaderLabel = False
414    
415     if arguments.makeFancy:
416     drawHeaderLabel = True
417    
418     #now that flags are set, draw the appropriate labels
419    
420     if drawHeaderLabel:
421     HeaderLabel.Draw()
422    
423    
424    
425    
426     #drawing the ratio or difference plot if requested
427    
428     if makeRatioPlots or makeDiffPlots:
429     Canvas.cd(2)
430     if makeRatioPlots:
431     Comparison = ratioHistogram(Histograms[0],Histograms[1])
432     elif makeDiffPlots:
433     Comparison = Histograms[0].Clone("diff")
434     Comparison.Add(Histograms[1],-1)
435     Comparison.SetTitle("")
436     Comparison.GetYaxis().SetTitle("hist1-hist2")
437     Comparison.GetXaxis().SetTitle(xAxisLabel)
438     Comparison.GetYaxis().CenterTitle()
439     Comparison.GetYaxis().SetTitleSize(0.1)
440     Comparison.GetYaxis().SetTitleOffset(0.5)
441     Comparison.GetXaxis().SetTitleSize(0.15)
442     Comparison.GetYaxis().SetLabelSize(0.1)
443     Comparison.GetXaxis().SetLabelSize(0.15)
444     if makeRatioPlots:
445     RatioYRange = 1.15
446     if arguments.ratioYRange:
447     RatioYRange = float(arguments.ratioYRange)
448     Comparison.GetYaxis().SetRangeUser(-1*RatioYRange, RatioYRange)
449     elif makeDiffPlots:
450     YMax = Comparison.GetMaximum()
451     YMin = Comparison.GetMinimum()
452     if YMax <= 0 and YMin <= 0:
453     Comparison.GetYaxis().SetRangeUser(-1.2*YMin,0)
454     elif YMax >= 0 and YMin >= 0:
455     Comparison.GetYaxis().SetRangeUser(0,1.2*YMax)
456     else: #axis crosses y=0
457     if abs(YMax) > abs(YMin):
458     Comparison.GetYaxis().SetRangeUser(-1.2*YMax,1.2*YMax)
459     else:
460     Comparison.GetYaxis().SetRangeUser(-1.2*YMin,1.2*YMin)
461    
462     Comparison.GetYaxis().SetNdivisions(205)
463     Comparison.Draw()
464    
465     outputFile.cd()
466     Canvas.Write()
467    
468     if arguments.savePDFs:
469     Canvas.SaveAs("efficiency_histograms_pdfs/"+histogramName+".pdf")
470    
471    
472     ##########################################################################################################################################
473     ##########################################################################################################################################
474     ##########################################################################################################################################
475    
476    
477     #### make output file
478     outputFileName = "efficiency_histograms.root"
479     if arguments.outputFileName:
480     outputFileName = arguments.outputFileName
481    
482     outputFile = TFile(outputFileName, "RECREATE")
483    
484     first_input = input_sources[0]
485    
486     #### use the first input file as a template and make efficiency versions of all its histograms
487     testFile = TFile("condor/" + first_input['condor_dir'] + "/" + first_input['dataset'] + ".root")
488     testFile.cd("OSUAnalysis/" + first_input['num_channel'])
489    
490     if arguments.savePDFs:
491     os.system("rm -rf efficiency_histograms_pdfs")
492     os.system("mkdir efficiency_histograms_pdfs")
493    
494     for key in gDirectory.GetListOfKeys():
495     if re.match ('TH1', key.GetClassName()): #found a 1D histogram
496     MakeOneDHist(key.GetName())
497    
498     testFile.Close()
499     outputFile.Close()