ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/OSUT3Analysis/Configuration/scripts/makePlots.py
Revision: 1.47
Committed: Tue Jun 25 15:51:32 2013 UTC (11 years, 10 months ago) by lantonel
Content type: text/x-python
Branch: MAIN
Changes since 1.46: +26 -15 lines
Log Message:
added y axis labels and cleaned up the formatting a bit

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
12
13
14 ### parse the command-line options
15
16 parser = OptionParser()
17 parser = set_commandline_arguments(parser)
18
19 parser.add_option("-f", "--fancy", action="store_true", dest="makeFancy", default=False,
20 help="removes the title and replaces it with the official CMS plot heading")
21 parser.add_option("--ylog", action="store_true", dest="setLogY", default=False,
22 help="Set logarithmic scale on vertical axis on all plots")
23 parser.add_option("--ymin", dest="setYMin",
24 help="Minimum of y axis")
25 parser.add_option("--ymax", dest="setYMax",
26 help="Maximum of y axis")
27
28 (arguments, args) = parser.parse_args()
29
30
31 if arguments.localConfig:
32 sys.path.append(os.getcwd())
33 exec("from " + arguments.localConfig.rstrip('.py') + " import *")
34
35 #### deal with conflicting arguments
36 if arguments.normalizeToData and arguments.normalizeToUnitArea:
37 print "Conflicting normalizations requsted, will normalize to unit area"
38 arguments.normalizeToData = False
39 if arguments.normalizeToData and arguments.noStack:
40 print "You have asked to scale non-stacked backgrounds to data. This is a very strange request. Will normalize to unit area instead"
41 arguments.normalizeToData = False
42 arguments.normalizeToUnitArea = True
43 if arguments.makeRatioPlots and arguments.makeDiffPlots:
44 print "You have requested both ratio and difference plots. Will make just ratio plots instead"
45 arguments.makeRatioPlots = False
46
47
48 from ROOT import TFile, gROOT, gStyle, gDirectory, TStyle, THStack, TH1F, TCanvas, TString, TLegend, TLegendEntry, THStack, TIter, TKey, TPaveLabel, gPad
49
50
51 ### setting ROOT options so our plots will look awesome and everyone will love us
52
53 gROOT.SetBatch()
54 gStyle.SetOptStat(0)
55 gStyle.SetCanvasBorderMode(0)
56 gStyle.SetPadBorderMode(0)
57 gStyle.SetPadColor(0)
58 gStyle.SetCanvasColor(0)
59 gStyle.SetTextFont(42)
60 gStyle.SetCanvasDefH(600)
61 gStyle.SetCanvasDefW(600)
62 gStyle.SetCanvasDefX(0)
63 gStyle.SetCanvasDefY(0)
64 gStyle.SetPadTopMargin(0.07)
65 gStyle.SetPadBottomMargin(0.13)
66 gStyle.SetPadLeftMargin(0.15)
67 gStyle.SetPadRightMargin(0.05)
68 gStyle.SetTitleColor(1, "XYZ")
69 gStyle.SetTitleFont(42, "XYZ")
70 gStyle.SetTitleSize(0.04, "XYZ")
71 gStyle.SetTitleXOffset(1.1)
72 gStyle.SetTitleYOffset(2)
73 gStyle.SetTextAlign(12)
74 gStyle.SetLabelColor(1, "XYZ")
75 gStyle.SetLabelFont(42, "XYZ")
76 gStyle.SetLabelOffset(0.007, "XYZ")
77 gStyle.SetLabelSize(0.04, "XYZ")
78 gStyle.SetAxisColor(1, "XYZ")
79 gStyle.SetStripDecimals(True)
80 gStyle.SetTickLength(0.03, "XYZ")
81 gStyle.SetNdivisions(510, "XYZ")
82 gStyle.SetPadTickX(1)
83 gStyle.SetPadTickY(1)
84 gROOT.ForceStyle()
85
86
87 #set the text for the luminosity label
88 if(intLumi < 1000.):
89 LumiText = "L_{int} = " + str(intLumi) + " pb^{-1}"
90 LumiText = "L_{int} = " + str.format('{0:.1f}', LumiInPb) + " pb^{-1}"
91 else:
92 LumiInFb = intLumi/1000.
93 LumiText = "L_{int} = " + str.format('{0:.1f}', LumiInFb) + " fb^{-1}"
94
95 #bestest place for lumi. label, in top left corner
96 topLeft_x_left = 0.1375839
97 topLeft_x_right = 0.4580537
98 topLeft_y_bottom = 0.8479021
99 topLeft_y_top = 0.9475524
100 topLeft_y_offset = 0.035
101
102 #set the text for the fancy heading
103 HeaderText = "CMS Preliminary: " + LumiText + " at #sqrt{s} = 8 TeV"
104
105 #position for header
106 header_x_left = 0.2181208
107 header_x_right = 0.9562937
108 header_y_bottom = 0.9479866
109 header_y_top = 0.9947552
110
111
112
113 ##########################################################################################################################################
114 ##########################################################################################################################################
115 ##########################################################################################################################################
116
117 def plainTextString(inputString):
118 inputString.replace(' ','_')
119 inputString.replace('<','lt')
120 inputString.replace('>','gt')
121 inputString.replace('(','')
122 inputString.replace(')','')
123 inputString.replace('=','eq')
124
125 return inputString
126
127 ##########################################################################################################################################
128 ##########################################################################################################################################
129 ##########################################################################################################################################
130
131 # some fancy-ass code from Andrzej Zuranski to merge bins in the ratio plot until the error goes below some threshold
132 def ratioHistogram( dataHist, mcHist, relErrMax=0.10):
133
134 def groupR(group):
135 Data,MC = [float(sum(hist.GetBinContent(i) for i in group)) for hist in [dataHist,mcHist]]
136 return (Data-MC)/MC if MC else 0
137
138 def groupErr(group):
139 Data,MC = [float(sum(hist.GetBinContent(i) for i in group)) for hist in [dataHist,mcHist]]
140 dataErr2,mcErr2 = [sum(hist.GetBinError(i)**2 for i in group) for hist in [dataHist,mcHist]]
141 return abs(math.sqrt( (dataErr2+mcErr2)/(Data-MC)**2 + mcErr2/MC**2 ) * (Data-MC)/MC) if Data and MC else 0
142
143 def regroup(groups):
144 err,iG = max( (groupErr(g),groups.index(g)) for g in groups )
145 if err < relErrMax or len(groups)<3 : return groups
146 iH = max( [iG-1,iG+1], key = lambda i: groupErr(groups[i]) if 0<=i<len(groups) else -1 )
147 iLo,iHi = sorted([iG,iH])
148 return regroup(groups[:iLo] + [groups[iLo]+groups[iHi]] + groups[iHi+1:])
149
150 #don't rebin the histograms of the number of a given object (except for the pileup ones)
151 if ((dataHist.GetName().find("num") is not -1 and dataHist.GetName().find("Primaryvertexs") is -1) or
152 dataHist.GetName().find("CutFlow") is not -1 or
153 dataHist.GetName().find("GenMatch") is not -1):
154 ratio = dataHist.Clone()
155 ratio.Add(mcHist,-1)
156 ratio.Divide(mcHist)
157 ratio.SetTitle("")
158 else:
159 groups = regroup( [(i,) for i in range(1,1+dataHist.GetNbinsX())] )
160 ratio = TH1F("ratio","",len(groups), array('d', [dataHist.GetBinLowEdge(min(g)) for g in groups ] + [dataHist.GetXaxis().GetBinUpEdge(dataHist.GetNbinsX())]) )
161 for i,g in enumerate(groups) :
162 ratio.SetBinContent(i+1,groupR(g))
163 ratio.SetBinError(i+1,groupErr(g))
164
165 ratio.GetYaxis().SetTitle("#frac{Data-MC}{MC}")
166 ratio.SetLineColor(1)
167 ratio.SetLineWidth(2)
168 return ratio
169
170 ##########################################################################################################################################
171 ##########################################################################################################################################
172 ##########################################################################################################################################
173
174
175 def MakeOneDHist(pathToDir,histogramName):
176
177
178 numBgMCSamples = 0
179 numDataSamples = 0
180 numSignalSamples = 0
181
182 Stack = THStack("stack",histogramName)
183
184 HeaderLabel = TPaveLabel(header_x_left,header_y_bottom,header_x_right,header_y_top,HeaderText,"NDC")
185 HeaderLabel.SetTextAlign(32)
186 HeaderLabel.SetBorderSize(0)
187 HeaderLabel.SetFillColor(0)
188 HeaderLabel.SetFillStyle(0)
189
190 LumiLabel = TPaveLabel(topLeft_x_left,topLeft_y_bottom,topLeft_x_right,topLeft_y_top,LumiText,"NDC")
191 LumiLabel.SetBorderSize(0)
192 LumiLabel.SetFillColor(0)
193 LumiLabel.SetFillStyle(0)
194
195 NormLabel = TPaveLabel()
196 NormLabel.SetDrawOption("NDC")
197 NormLabel.SetX1NDC(topLeft_x_left)
198 NormLabel.SetX2NDC(topLeft_x_right)
199
200 NormLabel.SetBorderSize(0)
201 NormLabel.SetFillColor(0)
202 NormLabel.SetFillStyle(0)
203
204 NormText = ""
205 if arguments.normalizeToUnitArea:
206 NormText = "Scaled to unit area"
207 elif arguments.normalizeToData:
208 NormText = "MC scaled to data"
209 NormLabel.SetLabel(NormText)
210
211
212 BgMCLegend = TLegend()
213 BgTitle = BgMCLegend.AddEntry(0, "Data & Bkgd. MC", "H")
214 BgTitle.SetTextAlign(22)
215 BgTitle.SetTextFont(62)
216 BgMCLegend.SetBorderSize(0)
217 BgMCLegend.SetFillColor(0)
218 BgMCLegend.SetFillStyle(0)
219 SignalMCLegend = TLegend()
220 SignalTitle = SignalMCLegend.AddEntry(0, "Signal MC", "H")
221 SignalTitle.SetTextAlign(22)
222 SignalTitle.SetTextFont(62)
223 SignalMCLegend.SetBorderSize(0)
224 SignalMCLegend.SetFillColor(0)
225 SignalMCLegend.SetFillStyle(0)
226
227 outputFile.cd(pathToDir)
228 Canvas = TCanvas(histogramName)
229 BgMCHistograms = []
230 BgMCLegendEntries = []
231 SignalMCHistograms = []
232 SignalMCLegendEntries = []
233 DataHistograms = []
234 DataLegendEntries = []
235
236
237 backgroundIntegral = 0
238 dataIntegral = 0
239 scaleFactor = 1
240
241 for sample in processed_datasets: # loop over different samples as listed in configurationOptions.py
242 dataset_file = "%s/%s.root" % (condor_dir,sample)
243 inputFile = TFile(dataset_file)
244 Histogram = inputFile.Get(pathToDir+"/"+histogramName).Clone()
245 Histogram.SetDirectory(0)
246 inputFile.Close()
247 if arguments.rebinFactor:
248 RebinFactor = int(arguments.rebinFactor)
249 #don't rebin histograms which will have less than 5 bins or any gen-matching histograms
250 if Histogram.GetNbinsX() >= RebinFactor*5 and Histogram.GetTitle().find("GenMatch") is -1:
251 Histogram.Rebin(RebinFactor)
252
253 xAxisLabel = Histogram.GetXaxis().GetTitle()
254 unitIndex = xAxisLabel.find("[")
255 if unitIndex is not -1: #x axis has a unit
256 yAxisLabel = "Entries / " + str(Histogram.GetXaxis().GetBinWidth(1)) + " " + xAxisLabel[unitIndex+1:-1]
257 else:
258 yAxisLabel = "A.U."
259
260 if not arguments.makeFancy:
261 histoTitle = Histogram.GetTitle()
262 else:
263 histoTitle = ""
264
265
266 legLabel = labels[sample]
267 if (arguments.printYields):
268 yieldHist = Histogram.Integral()
269 legLabel = legLabel + " (%.1f)" % yieldHist
270
271 if( types[sample] == "bgMC"):
272
273 numBgMCSamples += 1
274 backgroundIntegral += Histogram.Integral()
275
276 Histogram.SetLineStyle(1)
277 if(arguments.noStack):
278 Histogram.SetFillStyle(0)
279 Histogram.SetLineColor(colors[sample])
280 Histogram.SetLineWidth(2)
281 else:
282 Histogram.SetFillStyle(1001)
283 Histogram.SetFillColor(colors[sample])
284 Histogram.SetLineColor(1)
285 Histogram.SetLineWidth(1)
286
287 BgMCLegendEntries.append(legLabel)
288 BgMCHistograms.append(Histogram)
289
290
291 elif( types[sample] == "signalMC"):
292
293 numSignalSamples += 1
294
295 Histogram.SetFillStyle(0)
296 Histogram.SetLineColor(colors[sample])
297 Histogram.SetLineStyle(1)
298 Histogram.SetLineWidth(2)
299 if(arguments.normalizeToUnitArea and Histogram.Integral() > 0):
300 Histogram.Scale(1./Histogram.Integral())
301
302 SignalMCLegendEntries.append(legLabel)
303 SignalMCHistograms.append(Histogram)
304
305 elif( types[sample] == "data"):
306
307 numDataSamples += 1
308 dataIntegral += Histogram.Integral()
309
310 Histogram.SetMarkerStyle(20)
311 Histogram.SetMarkerSize(0.8)
312 Histogram.SetFillStyle(0)
313 Histogram.SetLineColor(colors[sample])
314 Histogram.SetLineStyle(1)
315 Histogram.SetLineWidth(2)
316 if(arguments.normalizeToUnitArea and Histogram.Integral() > 0):
317 Histogram.Scale(1./Histogram.Integral())
318
319 DataLegendEntries.append(legLabel)
320 DataHistograms.append(Histogram)
321
322 #scaling histograms as per user's specifications
323 if dataIntegral > 0 and backgroundIntegral > 0:
324 scaleFactor = dataIntegral/backgroundIntegral
325 for bgMCHist in BgMCHistograms:
326 if arguments.normalizeToData:
327 bgMCHist.Scale(scaleFactor)
328
329 if arguments.normalizeToUnitArea and not arguments.noStack and backgroundIntegral > 0:
330 bgMCHist.Scale(1./backgroundIntegral)
331 elif arguments.normalizeToUnitArea and arguments.noStack and bgMCHist.Integral() > 0:
332 bgMCHist.Scale(1./bgMCHist.Integral())
333
334 if not arguments.noStack:
335 Stack.Add(bgMCHist)
336
337
338
339 ### formatting data histograms and adding to legend
340 legendIndex = 0
341 for Histogram in DataHistograms:
342 BgMCLegend.AddEntry(Histogram,DataLegendEntries[legendIndex],"LEP")
343 legendIndex = legendIndex+1
344
345
346 ### creating the histogram to represent the statistical errors on the stack
347 if numBgMCSamples is not 0 and not arguments.noStack:
348 ErrorHisto = BgMCHistograms[0].Clone("errors")
349 ErrorHisto.SetFillStyle(3001)
350 ErrorHisto.SetFillColor(13)
351 ErrorHisto.SetLineWidth(0)
352 BgMCLegend.AddEntry(ErrorHisto,"Stat. Errors","F")
353 for Histogram in BgMCHistograms:
354 if Histogram is not BgMCHistograms[0]:
355 ErrorHisto.Add(Histogram)
356
357
358 ### formatting bgMC histograms and adding to legend
359 legendIndex = numBgMCSamples-1
360 for Histogram in reversed(BgMCHistograms):
361 if(arguments.noStack):
362 BgMCLegend.AddEntry(Histogram,BgMCLegendEntries[legendIndex],"L")
363 else:
364 BgMCLegend.AddEntry(Histogram,BgMCLegendEntries[legendIndex],"F")
365 legendIndex = legendIndex-1
366
367
368 ### formatting signalMC histograms and adding to legend
369 legendIndex = 0
370 for Histogram in SignalMCHistograms:
371 SignalMCLegend.AddEntry(Histogram,SignalMCLegendEntries[legendIndex],"L")
372 legendIndex = legendIndex+1
373
374
375 ### finding the maximum value of anything going on the canvas, so we know how to set the y-axis
376 finalMax = 0
377 if numBgMCSamples is not 0 and not arguments.noStack:
378 finalMax = ErrorHisto.GetMaximum() + ErrorHisto.GetBinError(ErrorHisto.GetMaximumBin())
379 else:
380 for bgMCHist in BgMCHistograms:
381 if(bgMCHist.GetMaximum() > finalMax):
382 finalMax = bgMCHist.GetMaximum()
383 for signalMCHist in SignalMCHistograms:
384 if(signalMCHist.GetMaximum() > finalMax):
385 finalMax = signalMCHist.GetMaximum()
386 for dataHist in DataHistograms:
387 if(dataHist.GetMaximum() > finalMax):
388 finalMax = dataHist.GetMaximum() + dataHist.GetBinError(dataHist.GetMaximumBin())
389 finalMax = 1.15*finalMax
390 if arguments.setYMax:
391 finalMax = float(arguments.setYMax)
392
393
394 ### Drawing histograms to canvas
395
396 outputFile.cd(pathToDir)
397
398 makeRatioPlots = arguments.makeRatioPlots
399 makeDiffPlots = arguments.makeDiffPlots
400
401 yAxisMin = 0.0001
402 if arguments.setYMin:
403 yAxisMin = float(arguments.setYMin)
404
405 if numBgMCSamples is 0 or numDataSamples is not 1:
406 makeRatioPlots = False
407 makeDiffPlots = False
408 if makeRatioPlots or makeDiffPlots:
409 Canvas.SetFillStyle(0)
410 Canvas.Divide(1,2)
411 Canvas.cd(1)
412 gPad.SetPad(0,0.25,1,1)
413 gPad.SetMargin(0.15,0.05,0.01,0.07)
414 gPad.SetFillStyle(0)
415 gPad.Update()
416 gPad.Draw()
417 if arguments.setLogY:
418 gPad.SetLogy()
419 Canvas.cd(2)
420 gPad.SetPad(0,0,1,0.25)
421 #format: gPad.SetMargin(l,r,b,t)
422 gPad.SetMargin(0.15,0.05,0.4,0.01)
423 gPad.SetFillStyle(0)
424 gPad.SetGridy(1)
425 gPad.Update()
426 gPad.Draw()
427
428 Canvas.cd(1)
429
430 if numBgMCSamples is not 0: # the first thing to draw to the canvas is a bgMC sample
431
432 if not arguments.noStack: # draw unstacked background samples
433 Stack.SetTitle(histoTitle)
434 Stack.Draw("HIST")
435 Stack.GetXaxis().SetTitle(xAxisLabel)
436 Stack.GetYaxis().SetTitle(yAxisLabel)
437 Stack.SetMaximum(finalMax)
438 Stack.SetMinimum(yAxisMin)
439 if makeRatioPlots or makeDiffPlots:
440 Stack.GetHistogram().GetXaxis().SetLabelSize(0)
441 #draw shaded error bands
442 ErrorHisto.Draw("A E2 SAME")
443
444 else: #draw the stacked backgrounds
445 BgMCHistograms[0].SetTitle(histoTitle)
446 BgMCHistograms[0].Draw("HIST")
447 BgMCHistograms[0].GetXaxis().SetTitle(xAxisLabel)
448 BgMCHistograms[0].GetYaxis().SetTitle(yAxisLabel)
449 BgMCHistograms[0].SetMaximum(finalMax)
450 BgMCHistograms[0].SetMinimum(yAxisMin)
451 for bgMCHist in BgMCHistograms:
452 bgMCHist.Draw("A HIST SAME")
453
454 for signalMCHist in SignalMCHistograms:
455 signalMCHist.Draw("A HIST SAME")
456 for dataHist in DataHistograms:
457 dataHist.Draw("A E X0 SAME")
458
459
460 elif numSignalSamples is not 0: # the first thing to draw to the canvas is a signalMC sample
461 SignalMCHistograms[0].SetTitle(histoTitle)
462 SignalMCHistograms[0].Draw("HIST")
463 SignalMCHistograms[0].GetXaxis().SetTitle(xAxisLabel)
464 SignalMCHistograms[0].GetYaxis().SetTitle(yAxisLabel)
465 SignalMCHistograms[0].SetMaximum(finalMax)
466 SignalMCHistograms[0].SetMinimum(yAxisMin)
467
468 for signalMCHist in SignalMCHistograms:
469 if(signalMCHist is not SignalMCHistograms[0]):
470 signalMCHist.Draw("A HIST SAME")
471 for dataHist in DataHistograms:
472 dataHist.Draw("A E X0 SAME")
473
474
475 elif(numDataSamples is not 0): # the first thing to draw to the canvas is a data sample
476 DataHistograms[0].SetTitle(histoTitle)
477 DataHistograms[0].Draw("E")
478 DataHistograms[0].GetXaxis().SetTitle(xAxisLabel)
479 DataHistograms[0].GetYaxis().SetTitle(yAxisLabel)
480 DataHistograms[0].SetMaximum(finalMax)
481 DataHistograms[0].SetMinimum(yAxisMin)
482 for dataHist in DataHistograms:
483 if(dataHist is not DataHistograms[0]):
484 dataHist.Draw("A E X0 SAME")
485
486
487
488 #legend coordinates, empirically determined :-)
489 x_left = 0.6761745
490 x_right = 0.9328859
491 x_width = x_right - x_left
492 y_max = 0.9
493 entry_height = 0.05
494
495 if(numBgMCSamples is not 0 or numDataSamples is not 0): #then draw the data & bgMC legend
496
497 numExtraEntries = 1 # count the legend title
498 BgMCLegend.SetX1NDC(x_left)
499 if numBgMCSamples > 0:
500 numExtraEntries = numExtraEntries + 1 # count the stat. errors entry
501
502 BgMCLegend.SetY1NDC(y_max-entry_height*(numExtraEntries+numBgMCSamples+numDataSamples))
503 BgMCLegend.SetX2NDC(x_right)
504 BgMCLegend.SetY2NDC(y_max)
505 BgMCLegend.Draw()
506
507 if(numSignalSamples is not 0): #then draw the signalMC legend to the left of the other one
508 SignalMCLegend.SetX1NDC(x_left-x_width)
509 SignalMCLegend.SetY1NDC(y_max-entry_height*(1+numSignalSamples)) # add one for the title
510 SignalMCLegend.SetX2NDC(x_left)
511 SignalMCLegend.SetY2NDC(y_max)
512 SignalMCLegend.Draw()
513
514 elif numSignalSamples is not 0: #draw the signalMC legend in the upper right corner
515 SignalMCLegend.SetX1NDC(x_left)
516 SignalMCLegend.SetY1NDC(y_max-entry_height*(1+numSignalSamples)) # add one for the title
517 SignalMCLegend.SetX2NDC(x_right)
518 SignalMCLegend.SetY2NDC(y_max)
519 SignalMCLegend.Draw()
520
521
522 # Deciding which text labels to draw and drawing them
523 drawLumiLabel = False
524 drawNormLabel = False
525 offsetNormLabel = False
526 drawHeaderLabel = False
527
528 if not arguments.normalizeToUnitArea or numDataSamples > 0: #don't draw the lumi label if there's no data and it's scaled to unit area
529 drawLumiLabel = True
530 #move the normalization label down before drawing if we drew the lumi. label
531 offsetNormLabel = True
532 if arguments.normalizeToUnitArea or arguments.normalizeToData:
533 drawNormLabel = True
534 if arguments.makeFancy:
535 drawHeaderLabel = True
536 drawLumiLabel = False
537
538 #now that flags are set, draw the appropriate labels
539
540 if drawLumiLabel:
541 LumiLabel.Draw()
542
543 if drawNormLabel:
544 if offsetNormLabel:
545 NormLabel.SetY1NDC(topLeft_y_bottom-topLeft_y_offset)
546 NormLabel.SetY2NDC(topLeft_y_top-topLeft_y_offset)
547 else:
548 NormLabel.SetY1NDC(topLeft_y_bottom)
549 NormLabel.SetY2NDC(topLeft_y_top)
550 NormLabel.Draw()
551
552 if drawHeaderLabel:
553 HeaderLabel.Draw()
554
555
556
557
558 #drawing the ratio or difference plot if requested
559
560 if makeRatioPlots or makeDiffPlots:
561 Canvas.cd(2)
562 BgSum = Stack.GetStack().Last()
563 if makeRatioPlots:
564 Comparison = ratioHistogram(DataHistograms[0],BgSum)
565 elif makeDiffPlots:
566 Comparison = DataHistograms[0].Clone("diff")
567 Comparison.Add(BgSum,-1)
568 Comparison.SetTitle("")
569 Comparison.GetYaxis().SetTitle("Data-MC")
570 Comparison.GetXaxis().SetTitle(xAxisLabel)
571 Comparison.GetYaxis().CenterTitle()
572 Comparison.GetYaxis().SetTitleSize(0.1)
573 Comparison.GetYaxis().SetTitleOffset(0.5)
574 Comparison.GetXaxis().SetTitleSize(0.15)
575 Comparison.GetYaxis().SetLabelSize(0.1)
576 Comparison.GetXaxis().SetLabelSize(0.15)
577 if makeRatioPlots:
578 RatioYRange = 1.15
579 if arguments.ratioYRange:
580 RatioYRange = float(arguments.ratioYRange)
581 Comparison.GetYaxis().SetRangeUser(-1*RatioYRange, RatioYRange)
582 elif makeDiffPlots:
583 YMax = Comparison.GetMaximum()
584 YMin = Comparison.GetMinimum()
585 if YMax <= 0 and YMin <= 0:
586 Comparison.GetYaxis().SetRangeUser(-1.2*YMin,0)
587 elif YMax >= 0 and YMin >= 0:
588 Comparison.GetYaxis().SetRangeUser(0,1.2*YMax)
589 else: #axis crosses y=0
590 if abs(YMax) > abs(YMin):
591 Comparison.GetYaxis().SetRangeUser(-1.2*YMax,1.2*YMax)
592 else:
593 Comparison.GetYaxis().SetRangeUser(-1.2*YMin,1.2*YMin)
594
595 Comparison.GetYaxis().SetNdivisions(205)
596 Comparison.Draw()
597
598 Canvas.Write()
599 if arguments.savePDFs:
600 pathToDirString = plainTextString(pathToDir)
601 Canvas.SaveAs(condor_dir+"/stacked_histograms_pdfs/"+pathToDirString+"/"+histogramName+".pdf")
602
603
604 ##########################################################################################################################################
605 ##########################################################################################################################################
606 ##########################################################################################################################################
607
608 def MakeTwoDHist(pathToDir,histogramName):
609 numBgMCSamples = 0
610 numDataSamples = 0
611 numSignalSamples = 0
612
613 LumiLabel = TPaveLabel(0.1,0.8,0.34,0.9,LumiText,"NDC")
614 LumiLabel.SetBorderSize(0)
615 LumiLabel.SetFillColor(0)
616 LumiLabel.SetFillStyle(0)
617
618 BgMCLegend = TLegend(0.76,0.65,0.99,0.9)
619 BgMCLegend.AddEntry (0, "Data & Bkgd. MC", "H").SetTextFont (62)
620 BgMCLegend.SetBorderSize(0)
621 BgMCLegend.SetFillColor(0)
622 BgMCLegend.SetFillStyle(0)
623 SignalMCLegend = TLegend(0.76,0.135,0.99,0.377)
624 SignalMCLegend.AddEntry (0, "Signal MC", "H").SetTextFont (62)
625 SignalMCLegend.SetBorderSize(0)
626 SignalMCLegend.SetFillColor(0)
627 SignalMCLegend.SetFillStyle(0)
628
629 outputFile.cd(pathToDir)
630 Canvas = TCanvas(histogramName)
631 Canvas.SetRightMargin(0.2413793);
632 BgMCHistograms = []
633 SignalMCHistograms = []
634 DataHistograms = []
635
636 for sample in processed_datasets: # loop over different samples as listed in configurationOptions.py
637 dataset_file = "%s/%s.root" % (condor_dir,sample)
638 inputFile = TFile(dataset_file)
639 Histogram = inputFile.Get(pathToDir+"/"+histogramName).Clone()
640 Histogram.SetDirectory(0)
641 RebinFactor = int(arguments.rebinFactor)
642 if arguments.rebinFactor and Histogram.GetNbinsX() >= RebinFactor*10 and Histogram.GetNbinsY() >= RebinFactor*10:
643 Histogram.Rebin2D(RebinFactor)
644 inputFile.Close()
645 xAxisLabel = Histogram.GetXaxis().GetTitle()
646 yAxisLabel = Histogram.GetYaxis().GetTitle()
647 if not arguments.makeFancy:
648 histoTitle = Histogram.GetTitle()
649 else:
650 histoTitle = ""
651
652 if( types[sample] == "bgMC"):
653
654 numBgMCSamples += 1
655 Histogram.SetMarkerColor(colors[sample])
656 Histogram.SetFillColor(colors[sample])
657 BgMCLegend.AddEntry(Histogram,labels[sample],"F").SetTextFont (42)
658 BgMCHistograms.append(Histogram)
659
660 elif( types[sample] == "signalMC"):
661
662 numSignalSamples += 1
663 Histogram.SetMarkerColor(colors[sample])
664 Histogram.SetFillColor(colors[sample])
665 SignalMCLegend.AddEntry(Histogram,labels[sample],"F").SetTextFont (42)
666 SignalMCHistograms.append(Histogram)
667
668 elif( types[sample] == "data"):
669
670 numDataSamples += 1
671 Histogram.SetMarkerColor(colors[sample])
672 Histogram.SetFillColor(colors[sample])
673 BgMCLegend.AddEntry(Histogram,labels[sample],"F").SetTextFont (42)
674 DataHistograms.append(Histogram)
675
676
677 outputFile.cd(pathToDir)
678
679 if(numBgMCSamples is not 0):
680 BgMCHistograms[0].SetTitle(histoTitle)
681 BgMCHistograms[0].GetXaxis().SetTitle(xAxisLabel)
682 BgMCHistograms[0].GetYaxis().SetTitle(yAxisLabel)
683 BgMCHistograms[0].Draw()
684 for signalMCHist in SignalMCHistograms:
685 signalMCHist.Draw("SAME")
686 for dataHist in DataHistograms:
687 dataHist.Draw("SAME")
688
689 elif(numSignalSamples is not 0):
690 SignalMCHistograms[0].SetTitle(histoTitle)
691 SignalMCHistograms[0].Draw()
692 SignalMCHistograms[0].GetXaxis().SetTitle(xAxisLabel)
693 SignalMCHistograms[0].GetYaxis().SetTitle(yAxisLabel)
694 for signalMCHist in SignalMCHistograms:
695 if(signalMCHist is not SignalMCHistograms[0]):
696 signalMCHist.Draw("SAME")
697 for dataHist in DataHistograms:
698 dataHist.Draw("SAME")
699
700 elif(numDataSamples is not 0):
701 DataHistograms[0].SetTitle(histoTitle)
702 DataHistograms[0].GetXaxis().SetTitle(xAxisLabel)
703 DataHistograms[0].GetYaxis().SetTitle(yAxisLabel)
704 DataHistograms[0].Draw()
705 for dataHist in DataHistograms:
706 if(dataHist is not DataHistograms[0]):
707 dataHist.Draw("SAME")
708
709
710 if(numBgMCSamples is not 0 or numDataSamples is not 0):
711 BgMCLegend.Draw()
712 if(numSignalSamples is not 0):
713 SignalMCLegend.Draw()
714 if not arguments.normalizeToUnitArea or numDataSamples > 0:
715 LumiLabel.Draw()
716
717 Canvas.Write()
718
719
720
721
722 ##########################################################################################################################################
723 ##########################################################################################################################################
724 ##########################################################################################################################################
725
726 processed_datasets = []
727
728 condor_dir = set_condor_output_dir(arguments)
729
730 #### check which input datasets have valid output files
731 for sample in datasets:
732 fileName = condor_dir + "/" + sample + ".root"
733 if not os.path.exists(fileName):
734 continue
735 testFile = TFile(fileName)
736 if testFile.IsZombie() or not testFile.GetNkeys():
737 continue
738 processed_datasets.append(sample)
739
740 if len(processed_datasets) is 0:
741 sys.exit("No datasets have been processed")
742
743
744 #### make output file
745 outputFileName = "stacked_histograms.root"
746 if arguments.outputFileName:
747 outputFileName = arguments.outputFileName
748
749 outputFile = TFile(condor_dir + "/" + outputFileName, "RECREATE")
750
751
752
753 #### use the first input file as a template and make stacked versions of all its histograms
754 inputFile = TFile(condor_dir + "/" + processed_datasets[0] + ".root")
755 inputFile.cd()
756 outputFile.cd()
757
758 if arguments.savePDFs:
759 os.system("rm -rf %s/stacked_histograms_pdfs" % (condor_dir))
760 os.system("mkdir %s/stacked_histograms_pdfs" % (condor_dir))
761
762
763 #get root directory in the first layer, generally "OSUAnalysis"
764 for key in inputFile.GetListOfKeys():
765 if (key.GetClassName() != "TDirectoryFile"):
766 continue
767 rootDirectory = key.GetName()
768 outputFile.mkdir(rootDirectory)
769 if arguments.savePDFs:
770 os.system("mkdir %s/stacked_histograms_pdfs/%s" % (condor_dir,plainTextString(rootDirectory)))
771
772 #cd to root directory and look for histograms
773 inputFile.cd(rootDirectory)
774 for key2 in gDirectory.GetListOfKeys():
775
776 if re.match ('TH1', key2.GetClassName()): # found a 1-D histogram
777 MakeOneDHist(rootDirectory,key2.GetName())
778 elif re.match ('TH2', key2.GetClassName()) and arguments.draw2DPlots: # found a 2-D histogram
779 MakeTwoDHist(rootDirectory,key2.GetName())
780
781 elif (key2.GetClassName() == "TDirectoryFile"): # found a directory, cd there and look for histograms
782 level2Directory = rootDirectory+"/"+key2.GetName()
783
784 #make a corresponding directory in the output file
785 outputFile.cd(rootDirectory)
786 gDirectory.mkdir(key2.GetName())
787 if arguments.savePDFs:
788 os.system("mkdir %s/stacked_histograms_pdfs/%s" % (condor_dir,plainTextString(level2Directory)))
789
790 #####################################################
791 ### This layer is typically the "channels" layer ###
792 #####################################################
793
794 inputFile.cd(level2Directory)
795 for key3 in gDirectory.GetListOfKeys():
796 if re.match ('TH1', key3.GetClassName()): # found a 1-D histogram
797 MakeOneDHist(level2Directory,key3.GetName())
798 elif re.match ('TH2', key3.GetClassName()) and arguments.draw2DPlots: # found a 2-D histogram
799 MakeTwoDHist(level2Directory,key3.GetName())
800
801 elif (key3.GetClassName() == "TDirectoryFile"): # found a directory, cd there and look for histograms
802 level3Directory = level2Directory+"/"+key3.GetName()
803
804 #make a corresponding directory in the output file
805 outputFile.cd(level2Directory)
806 gDirectory.mkdir(key3.GetName())
807 if arguments.savePDFs:
808 os.system("mkdir %s/stacked_histograms_pdfs/%s" % (condor_dir,plainTextString(level3Directory)))
809
810 #################################################
811 ### This layer is typically the "cuts" layer ###
812 #################################################
813
814 inputFile.cd(level3Directory)
815 for key3 in gDirectory.GetListOfKeys():
816 if re.match ('TH1', key3.GetClassName()): # found a 1-D histogram
817 MakeOneDHist(level3Directory,key3.GetName())
818 elif re.match ('TH2', key3.GetClassName()) and arguments.draw2DPlots: # found a 2-D histogram
819 MakeTwoDHist(level3Directory,key3.GetName())
820
821
822 outputFile.Close()