ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/OSUT3Analysis/Configuration/scripts/makePlots.py
Revision: 1.59
Committed: Tue Sep 3 09:01:53 2013 UTC (11 years, 7 months ago) by lantonel
Content type: text/x-python
Branch: MAIN
CVS Tags: HEAD
Changes since 1.58: +84 -13 lines
Log Message:
changed the formatting of 2D histograms, also, don't skip an entire canvas if one of the histograms is missing, just draw it without that one

File Contents

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