ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/OSUT3Analysis/Configuration/scripts/makePlots.py
Revision: 1.28
Committed: Tue Apr 9 21:42:18 2013 UTC (12 years ago) by ahart
Content type: text/x-python
Branch: MAIN
Changes since 1.27: +18 -17 lines
Log Message:
Perform the full weighting at merging instead of splitting it between the merging and plotting steps.

File Contents

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