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

# 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, TH1F, TCanvas, TString, TLegend, TLegendEntry, 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 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 yAxisLabel = "#epsilon_{ " + cutName + "} (" + str(Histogram.GetXaxis().GetBinWidth(1)) + " " + xAxisLabel[unitBeginIndex+1:unitEndIndex] + " width)"
287 else:
288 yAxisLabel = "#epsilon_{ " + cutName + "} (" + str(Histogram.GetXaxis().GetBinWidth(1)) + " width)"
289 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()