ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/OSUT3Analysis/Configuration/scripts/makePlots.py
Revision: 1.57
Committed: Fri Aug 2 10:59:01 2013 UTC (11 years, 9 months ago) by wulsin
Content type: text/x-python
Branch: MAIN
Changes since 1.56: +139 -109 lines
Log Message:
Do not include data in blinded histograms specified in histsToBlind list in localOptions file

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