ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/OSUT3Analysis/Configuration/scripts/makePlots.py
Revision: 1.15
Committed: Fri Mar 8 16:34:27 2013 UTC (12 years, 1 month ago) by lantonel
Content type: text/x-python
Branch: MAIN
Changes since 1.14: +59 -3 lines
Log Message:
added ability to make data/MC ratio plots with the -r option

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