ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/OSUT3Analysis/Configuration/scripts/makePlots.py
Revision: 1.31
Committed: Mon Apr 15 13:48:09 2013 UTC (12 years ago) by lantonel
Content type: text/x-python
Branch: MAIN
Changes since 1.30: +32 -16 lines
Log Message:
added shaded regions to represent stat. error bars on the background stack

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 parser = OptionParser()
13 parser = set_commandline_arguments(parser)
14 (arguments, args) = parser.parse_args()
15
16 if arguments.localConfig:
17 sys.path.append(os.getcwd())
18 exec("from " + arguments.localConfig.rstrip('.py') + " import *")
19
20
21 outputFileName = "stacked_histograms.root"
22 if arguments.outputFileName:
23 outputFileName = arguments.outputFileName
24
25 condor_dir = set_condor_output_dir(arguments)
26
27
28
29 #### deal with conflicting arguments
30 if arguments.normalizeToData and arguments.normalizeToUnitArea:
31 print "Conflicting normalizations requsted, will normalize to unit area"
32 arguments.normalizeToData = False
33 if arguments.normalizeToData and arguments.noStack:
34 print "You have asked to scale non-stacked backgrounds to data. This is a very strange request. Will normalize to unit area instead"
35 arguments.normalizeToData = False
36 arguments.normalizeToUnitArea = True
37 if arguments.makeRatioPlots and arguments.makeDiffPlots:
38 print "You have requested both ratio and difference plots. Will make just ratio plots instead"
39 arguments.makeRatioPlots = False
40
41 from ROOT import TFile, gROOT, gStyle, gDirectory, TStyle, THStack, TH1F, TCanvas, TString, TLegend, TArrow, THStack, TIter, TKey, TPaveLabel, gPad
42
43 gROOT.SetBatch()
44 gStyle.SetOptStat(0)
45 gStyle.SetCanvasBorderMode(0)
46 gStyle.SetPadBorderMode(0)
47 gStyle.SetPadColor(0)
48 gStyle.SetCanvasColor(0)
49 gStyle.SetTextFont(42)
50 gROOT.ForceStyle()
51 outputFile = TFile(condor_dir + "/" + outputFileName, "RECREATE")
52
53 channels = []
54 processed_datasets = []
55
56 #### check which input datasets have valid output files
57 for sample in datasets:
58 fileName = condor_dir + "/" + sample + ".root"
59 if not os.path.exists(fileName):
60 continue
61 testFile = TFile(fileName)
62 if testFile.IsZombie() or not testFile.GetNkeys():
63 continue
64 processed_datasets.append(sample)
65
66 if len(processed_datasets) is 0:
67 sys.exit("No datasets have been processed")
68
69 #### open first input file and re-make its directory structure in the output file
70 testFile = TFile(condor_dir + "/" + processed_datasets[0] + ".root")
71 testFile.cd()
72 for key in testFile.GetListOfKeys():
73 if (key.GetClassName() != "TDirectoryFile"):
74 continue
75 outputFile.cd()
76 outputFile.mkdir(key.GetName())
77 rootDirectory = key.GetName()
78
79 testFile.cd(key.GetName())
80 for key2 in gDirectory.GetListOfKeys():
81 if (key2.GetClassName() != "TDirectoryFile"):
82 continue
83 outputFile.cd(key.GetName())
84 gDirectory.mkdir(key2.GetName())
85 channels.append(key2.GetName())
86
87
88 ## weight = intLumi / 10000.0
89 ## for dataset in processed_datasets:
90 ## dataset_file = "%s/%s.root" % (condor_dir,dataset)
91 ## fin = TFile (dataset_file)
92 ## flags = fin.Get ("flags")
93 ## noWeights = flags and flags.GetBinContent (1)
94 ## fin.Close ()
95
96 ## if types[dataset] != "data" and not noWeights:
97 ## os.system("mergeTFileServiceHistograms -i %s -o %s -w %g" % (dataset_file, dataset_file + "_tmp", weight))
98 ## else:
99 ## os.system("mergeTFileServiceHistograms -i %s -o %s -w %g" % (dataset_file, dataset_file + "_tmp", 1.0))
100
101 for channel in channels: # loop over final states, which each have their own directory
102
103 testFile.cd(rootDirectory+"/"+channel)
104
105 for key in gDirectory.GetListOfKeys(): # loop over histograms in the current directory
106 histogramName = key.GetName()
107
108 if re.match ('TH1', key.GetClassName()): # plot a 1-D histogram
109
110 numBgMCSamples = 0
111 numDataSamples = 0
112 numSignalSamples = 0
113
114 Stack = THStack("stack",histogramName)
115
116 if(intLumi < 1000.):
117 LumiText = "L_{int} = " + str(intLumi) + " pb^{-1}"
118 else:
119 getcontext().prec = 2
120 LumiInFb = intLumi/1000.
121 LumiText = "L_{int} = " + str(LumiInFb) + " fb^{-1}"
122
123 LumiLabel = TPaveLabel(0.1,0.8,0.34,0.9,LumiText,"NDC")
124 LumiLabel.SetBorderSize(0)
125 LumiLabel.SetFillColor(0)
126 LumiLabel.SetFillStyle(0)
127
128
129 BgMCLegend = TLegend(0.70,0.65,0.94,0.89)
130 BgMCLegend.AddEntry (0, "Data & Bkgd. MC", "H").SetTextFont (62)
131 BgMCLegend.SetBorderSize(0)
132 BgMCLegend.SetFillColor(0)
133 BgMCLegend.SetFillStyle(0)
134 SignalMCLegend = TLegend(0.45,0.65,0.70,0.89)
135 SignalMCLegend.AddEntry (0, "Signal MC", "H").SetTextFont (62)
136 SignalMCLegend.SetBorderSize(0)
137 SignalMCLegend.SetFillColor(0)
138 SignalMCLegend.SetFillStyle(0)
139
140 outputFile.cd(rootDirectory+"/"+channel)
141 Canvas = TCanvas(histogramName)
142 BgMCHistograms = []
143 SignalMCHistograms = []
144 DataHistograms = []
145
146 backgroundIntegral = 0
147 dataIntegral = 0
148 scaleFactor = 1
149
150 for sample in processed_datasets: # loop over different samples as listed in configurationOptions.py
151 dataset_file = "%s/%s.root" % (condor_dir,sample)
152 inputFile = TFile(dataset_file)
153 Histogram = inputFile.Get(rootDirectory+"/"+channel+"/"+histogramName).Clone()
154 Histogram.SetDirectory(0)
155 if arguments.rebinFactor:
156 RebinFactor = int(arguments.rebinFactor)
157 if Histogram.GetNbinsX() >= RebinFactor*10:
158 Histogram.Rebin(RebinFactor)
159 inputFile.Close()
160 xAxisLabel = Histogram.GetXaxis().GetTitle()
161 histoTitle = Histogram.GetTitle()
162
163 legLabel = labels[sample]
164 if (arguments.printYields):
165 yieldHist = Histogram.Integral()
166 legLabel = legLabel + " (%.1f)" % yieldHist
167
168 if( types[sample] == "bgMC"):
169
170 numBgMCSamples += 1
171
172 if(arguments.noStack):
173 Histogram.SetFillStyle(0)
174 Histogram.SetLineColor(colors[sample])
175 Histogram.SetLineWidth(2)
176 BgMCLegend.AddEntry(Histogram,legLabel,"L").SetTextFont (42)
177 else:
178 Histogram.SetFillStyle(1001)
179 Histogram.SetFillColor(colors[sample])
180 Histogram.SetLineColor(1)
181 Histogram.SetLineWidth(1)
182 BgMCLegend.AddEntry(Histogram,legLabel,"F").SetTextFont (42)
183
184 Histogram.SetLineStyle(1)
185
186 backgroundIntegral += Histogram.Integral()
187
188 BgMCHistograms.append(Histogram)
189
190
191 elif( types[sample] == "signalMC"):
192
193 numSignalSamples += 1
194
195 Histogram.SetFillStyle(0)
196 Histogram.SetLineColor(colors[sample])
197 Histogram.SetLineStyle(1)
198 Histogram.SetLineWidth(2)
199 if(arguments.normalizeToUnitArea and Histogram.Integral() > 0):
200 Histogram.Scale(1./Histogram.Integral())
201 SignalMCLegend.AddEntry(Histogram,legLabel,"L").SetTextFont (42)
202 SignalMCHistograms.append(Histogram)
203
204 elif( types[sample] == "data"):
205
206 numDataSamples += 1
207
208 Histogram.SetFillStyle(0)
209 Histogram.SetLineColor(colors[sample])
210 Histogram.SetLineStyle(1)
211 Histogram.SetLineWidth(2)
212 if(arguments.normalizeToUnitArea and Histogram.Integral() > 0):
213 Histogram.Scale(1./Histogram.Integral())
214
215 dataIntegral += Histogram.Integral()
216
217 BgMCLegend.AddEntry(Histogram,legLabel,"LEP").SetTextFont (42)
218 DataHistograms.append(Histogram)
219
220 if dataIntegral > 0 and backgroundIntegral > 0:
221 scaleFactor = dataIntegral/backgroundIntegral
222 for bgMCHist in BgMCHistograms:
223 if arguments.normalizeToData:
224 bgMCHist.Scale(scaleFactor)
225 if arguments.normalizeToUnitArea and not arguments.noStack and backgroundIntegral > 0:
226 bgMCHist.Scale(1./backgroundIntegral)
227 elif arguments.normalizeToUnitArea and arguments.noStack and bgMCHist.Integral() > 0:
228 bgMCHist.Scale(1./bgMCHist.Integral())
229 if not arguments.noStack:
230 Stack.Add(bgMCHist)
231
232
233 finalMax = 0
234 if not arguments.noStack:
235 finalMax = Stack.GetMaximum()
236 else:
237 for bgMCHist in BgMCHistograms:
238 if(bgMCHist.GetMaximum() > finalMax):
239 finalMax = bgMCHist.GetMaximum()
240 for signalMCHist in SignalMCHistograms:
241 if(signalMCHist.GetMaximum() > finalMax):
242 finalMax = signalMCHist.GetMaximum()
243 for dataHist in DataHistograms:
244 if(dataHist.GetMaximum() > finalMax):
245 finalMax = dataHist.GetMaximum()
246
247
248
249 if len(DataHistograms) is 1:
250 dataIntegral += DataHistograms[0].Integral()
251
252
253 ### Drawing histograms to canvas
254
255 if numBgMCSamples is not 0:
256 ErrorHisto = BgMCHistograms[0].Clone("errors")
257 ErrorHisto.SetFillStyle(3001)
258 ErrorHisto.SetFillColor(13)
259 ErrorHisto.SetLineWidth(0)
260 for Histogram in BgMCHistograms:
261 if Histogram is not BgMCHistograms[0]:
262 ErrorHisto.Add(Histogram)
263
264 outputFile.cd(rootDirectory+"/"+channel)
265
266 makeRatioPlots = arguments.makeRatioPlots
267 makeDiffPlots = arguments.makeDiffPlots
268
269 if numBgMCSamples is 0 or numDataSamples is not 1:
270 makeRatioPlots = False
271 makeDiffPlots = False
272 if makeRatioPlots or makeDiffPlots:
273 Canvas.SetFillStyle(0)
274 Canvas.Divide(1,2)
275 Canvas.cd(1)
276 gPad.SetPad(0.01,0.25,0.99,0.99)
277 gPad.SetMargin(0.1,0.05,0.02,0.07)
278 gPad.SetFillStyle(0)
279 gPad.Update()
280 gPad.Draw()
281 Canvas.cd(2)
282 gPad.SetPad(0.01,0.01,0.99,0.25)
283 #format: gPad.SetMargin(l,r,b,t)
284 gPad.SetMargin(0.1,0.05,0.4,0.02)
285 gPad.SetFillStyle(0)
286 gPad.SetGridy(1)
287 gPad.Update()
288 gPad.Draw()
289
290 Canvas.cd(1)
291
292 if numBgMCSamples is not 0:
293
294 if not arguments.noStack:
295 Stack.SetTitle(histoTitle)
296 Stack.Draw("HIST")
297 Stack.GetXaxis().SetTitle(xAxisLabel)
298 Stack.SetMaximum(1.1*finalMax)
299 Stack.SetMinimum(0.0001)
300 if makeRatioPlots or makeDiffPlots:
301 Stack.GetHistogram().GetXaxis().SetLabelSize(0)
302 #draw shaded error bands
303 ErrorHisto.Draw("E2 SAME")
304
305
306 else:
307 BgMCHistograms[0].SetTitle(histoTitle)
308 BgMCHistograms[0].Draw("HIST")
309 BgMCHistograms[0].GetXaxis().SetTitle(xAxisLabel)
310 BgMCHistograms[0].SetMaximum(1.1*finalMax)
311 BgMCHistograms[0].SetMinimum(0.0001)
312 for bgMCHist in BgMCHistograms:
313 bgMCHist.Draw("HIST SAME")
314 for signalMCHist in SignalMCHistograms:
315 signalMCHist.Draw("HIST SAME")
316 for dataHist in DataHistograms:
317 dataHist.Draw("E SAME")
318
319 elif numSignalSamples is not 0:
320 SignalMCHistograms[0].SetTitle(histoTitle)
321 SignalMCHistograms[0].Draw("HIST")
322 SignalMCHistograms[0].GetXaxis().SetTitle(xAxisLabel)
323 SignalMCHistograms[0].SetMaximum(1.1*finalMax)
324 SignalMCHistograms[0].SetMinimum(0.0001)
325 for signalMCHist in SignalMCHistograms:
326 if(signalMCHist is not SignalMCHistograms[0]):
327 signalMCHist.Draw("HIST SAME")
328 for dataHist in DataHistograms:
329 dataHist.Draw("E SAME")
330
331 elif(numDataSamples is not 0):
332 DataHistograms[0].SetTitle(histoTitle)
333 DataHistograms[0].Draw("E")
334 DataHistograms[0].GetXaxis().SetTitle(xAxisLabel)
335 DataHistograms[0].SetMaximum(1.1*finalMax)
336 DataHistograms[0].SetMinimum(0.0001)
337 for dataHist in DataHistograms:
338 if(dataHist is not DataHistograms[0]):
339 dataHist.Draw("E SAME")
340
341
342 if(numBgMCSamples is not 0 or numDataSamples is not 0):
343 BgMCLegend.Draw()
344 if(numSignalSamples is not 0):
345 SignalMCLegend.Draw()
346
347 if not arguments.normalizeToUnitArea or numDataSamples > 0:
348 LumiLabel.Draw()
349 if arguments.normalizeToData and numBgMCSamples > 0 and numDataSamples > 0:
350 NormLabel = TPaveLabel(0.1,0.75,0.35,0.85,"MC scaled to data","NDC")
351 NormLabel.SetBorderSize(0)
352 NormLabel.SetFillColor(0)
353 NormLabel.SetFillStyle(0)
354 NormLabel.Draw()
355 elif arguments.normalizeToUnitArea:
356 NormLabel = TPaveLabel(0.1,0.75,0.35,0.85,"Scaled to unit area","NDC")
357 NormLabel.SetBorderSize(0)
358 NormLabel.SetFillColor(0)
359 NormLabel.SetFillStyle(0)
360 NormLabel.Draw()
361
362
363 if makeRatioPlots or makeDiffPlots:
364 Canvas.cd(2)
365 BgSum = Stack.GetStack().Last()
366 Comparison = DataHistograms[0].Clone()
367 Comparison.Add(BgSum,-1)
368 if not makeDiffPlots:
369 Comparison.Divide(BgSum)
370 Comparison.SetTitle("")
371 Comparison.GetXaxis().SetTitle(xAxisLabel)
372 if makeRatioPlots:
373 Comparison.GetYaxis().SetTitle("#frac{Data-MC}{MC}")
374 elif makeDiffPlots:
375 Comparison.GetYaxis().SetTitle("Data-MC")
376 Comparison.GetYaxis().CenterTitle()
377 Comparison.GetYaxis().SetTitleSize(0.1)
378 Comparison.GetYaxis().SetTitleOffset(0.35)
379 Comparison.GetXaxis().SetTitleSize(0.15)
380 Comparison.GetYaxis().SetLabelSize(0.1)
381 Comparison.GetXaxis().SetLabelSize(0.15)
382 if makeRatioPlots:
383 Comparison.GetYaxis().SetRangeUser(-1,1)
384 elif makeDiffPlots:
385 YMax = Comparison.GetMaximum()
386 YMin = Comparison.GetMinimum()
387 if YMax <= 0 and YMin <= 0:
388 Comparison.GetYaxis().SetRangeUser(-1.2*YMin,0)
389 elif YMax >= 0 and YMin >= 0:
390 Comparison.GetYaxis().SetRangeUser(0,1.2*YMax)
391 else: #axis crosses y=0
392 if abs(YMax) > abs(YMin):
393 Comparison.GetYaxis().SetRangeUser(-1.2*YMax,1.2*YMax)
394 else:
395 Comparison.GetYaxis().SetRangeUser(-1.2*YMin,1.2*YMin)
396
397 Comparison.GetYaxis().SetNdivisions(205)
398 Comparison.Draw()
399
400 Canvas.Write()
401
402
403 if re.match ('TH2', key.GetClassName()) and arguments.draw2DPlots: # plot a 2-D histogram
404
405 numBgMCSamples = 0
406 numDataSamples = 0
407 numSignalSamples = 0
408
409 if(intLumi < 1000.):
410 LumiText = "L_{int} = " + str(intLumi) + " pb^{-1}"
411 else:
412 getcontext().prec = 2
413 LumiInFb = intLumi/1000.
414 LumiText = "L_{int} = " + str(LumiInFb) + " fb^{-1}"
415
416 LumiLabel = TPaveLabel(0.1,0.8,0.34,0.9,LumiText,"NDC")
417 LumiLabel.SetBorderSize(0)
418 LumiLabel.SetFillColor(0)
419 LumiLabel.SetFillStyle(0)
420
421 BgMCLegend = TLegend(0.76,0.65,0.99,0.9)
422 BgMCLegend.AddEntry (0, "Data & Bkgd. MC", "H").SetTextFont (62)
423 BgMCLegend.SetBorderSize(0)
424 BgMCLegend.SetFillColor(0)
425 BgMCLegend.SetFillStyle(0)
426 SignalMCLegend = TLegend(0.76,0.135,0.99,0.377)
427 SignalMCLegend.AddEntry (0, "Signal MC", "H").SetTextFont (62)
428 SignalMCLegend.SetBorderSize(0)
429 SignalMCLegend.SetFillColor(0)
430 SignalMCLegend.SetFillStyle(0)
431
432 outputFile.cd(rootDirectory+"/"+channel)
433 Canvas = TCanvas(histogramName)
434 Canvas.SetRightMargin(0.2413793);
435 BgMCHistograms = []
436 SignalMCHistograms = []
437 DataHistograms = []
438
439 for sample in processed_datasets: # loop over different samples as listed in configurationOptions.py
440 dataset_file = "%s/%s.root" % (condor_dir,sample)
441 inputFile = TFile(dataset_file)
442 Histogram = inputFile.Get(rootDirectory+"/"+channel+"/"+histogramName).Clone()
443 Histogram.SetDirectory(0)
444 RebinFactor = int(arguments.rebinFactor)
445 if arguments.rebinFactor and Histogram.GetNbinsX() >= RebinFactor*10 and Histogram.GetNbinsY() >= RebinFactor*10:
446 Histogram.Rebin2D(RebinFactor)
447 inputFile.Close()
448 xAxisLabel = Histogram.GetXaxis().GetTitle()
449 yAxisLabel = Histogram.GetYaxis().GetTitle()
450 histoTitle = Histogram.GetTitle()
451
452 if( types[sample] == "bgMC"):
453
454 numBgMCSamples += 1
455 Histogram.SetMarkerColor(colors[sample])
456 Histogram.SetFillColor(colors[sample])
457 BgMCLegend.AddEntry(Histogram,labels[sample],"F").SetTextFont (42)
458 BgMCHistograms.append(Histogram)
459
460 elif( types[sample] == "signalMC"):
461
462 numSignalSamples += 1
463 Histogram.SetMarkerColor(colors[sample])
464 Histogram.SetFillColor(colors[sample])
465 SignalMCLegend.AddEntry(Histogram,labels[sample],"F").SetTextFont (42)
466 SignalMCHistograms.append(Histogram)
467
468 elif( types[sample] == "data"):
469
470 numDataSamples += 1
471 Histogram.SetMarkerColor(colors[sample])
472 Histogram.SetFillColor(colors[sample])
473 BgMCLegend.AddEntry(Histogram,labels[sample],"F").SetTextFont (42)
474 DataHistograms.append(Histogram)
475
476
477 outputFile.cd(rootDirectory+"/"+channel)
478
479 if(numBgMCSamples is not 0):
480 BgMCHistograms[0].SetTitle(histoTitle)
481 BgMCHistograms[0].GetXaxis().SetTitle(xAxisLabel)
482 BgMCHistograms[0].GetYaxis().SetTitle(yAxisLabel)
483 BgMCHistograms[0].Draw()
484 for signalMCHist in SignalMCHistograms:
485 signalMCHist.Draw("SAME")
486 for dataHist in DataHistograms:
487 dataHist.Draw("SAME")
488
489 elif(numSignalSamples is not 0):
490 SignalMCHistograms[0].SetTitle(histoTitle)
491 SignalMCHistograms[0].Draw()
492 SignalMCHistograms[0].GetXaxis().SetTitle(xAxisLabel)
493 SignalMCHistograms[0].GetYaxis().SetTitle(yAxisLabel)
494 for signalMCHist in SignalMCHistograms:
495 if(signalMCHist is not SignalMCHistograms[0]):
496 signalMCHist.Draw("SAME")
497 for dataHist in DataHistograms:
498 dataHist.Draw("SAME")
499
500 elif(numDataSamples is not 0):
501 DataHistograms[0].SetTitle(histoTitle)
502 DataHistograms[0].GetXaxis().SetTitle(xAxisLabel)
503 DataHistograms[0].GetYaxis().SetTitle(yAxisLabel)
504 DataHistograms[0].Draw()
505 for dataHist in DataHistograms:
506 if(dataHist is not DataHistograms[0]):
507 dataHist.Draw("SAME")
508
509
510 if(numBgMCSamples is not 0 or numDataSamples is not 0):
511 BgMCLegend.Draw()
512 if(numSignalSamples is not 0):
513 SignalMCLegend.Draw()
514 if not arguments.normalizeToUnitArea or numDataSamples > 0:
515 LumiLabel.Draw()
516
517 Canvas.Write()
518
519
520
521 #for dataset in processed_datasets:
522 # dataset_file = "%s/%s.root_tmp" % (condor_dir,dataset)
523 # os.remove(dataset_file)
524
525 outputFile.Close()