ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/OSUT3Analysis/Configuration/scripts/makePlots.py
Revision: 1.27
Committed: Thu Mar 28 10:45:05 2013 UTC (12 years, 1 month ago) by lantonel
Content type: text/x-python
Branch: MAIN
Changes since 1.26: +4 -3 lines
Log Message:
fixed bug in rebinning 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 (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 if types[dataset] != "data" and not noWeights:
95 os.system("mergeTFileServiceHistograms -i %s -o %s -w %g" % (dataset_file, dataset_file + "_tmp", weight))
96 else:
97 os.system("mergeTFileServiceHistograms -i %s -o %s -w %g" % (dataset_file, dataset_file + "_tmp", 1.0))
98
99 for channel in channels: # loop over final states, which each have their own directory
100
101 testFile.cd(rootDirectory+"/"+channel)
102
103 for key in gDirectory.GetListOfKeys(): # loop over histograms in the current directory
104 histogramName = key.GetName()
105
106 if re.match ('TH1', key.GetClassName()): # plot a 1-D histogram
107
108 numBgMCSamples = 0
109 numDataSamples = 0
110 numSignalSamples = 0
111
112 Stack = THStack("stack",histogramName)
113
114 if(intLumi < 1000.):
115 LumiText = "L_{int} = " + str(intLumi) + " pb^{-1}"
116 else:
117 getcontext().prec = 2
118 LumiInFb = intLumi/1000.
119 LumiText = "L_{int} = " + str(LumiInFb) + " fb^{-1}"
120
121 LumiLabel = TPaveLabel(0.1,0.8,0.34,0.9,LumiText,"NDC")
122 LumiLabel.SetBorderSize(0)
123 LumiLabel.SetFillColor(0)
124 LumiLabel.SetFillStyle(0)
125
126 BgMCLegend = TLegend(0.70,0.65,0.94,0.89, "Data & Bkgd. MC")
127 BgMCLegend.SetBorderSize(0)
128 BgMCLegend.SetFillColor(0)
129 BgMCLegend.SetFillStyle(0)
130 SignalMCLegend = TLegend(0.45,0.65,0.70,0.89,"Signal MC")
131 SignalMCLegend.SetBorderSize(0)
132 SignalMCLegend.SetFillColor(0)
133 SignalMCLegend.SetFillStyle(0)
134
135 outputFile.cd(rootDirectory+"/"+channel)
136 Canvas = TCanvas(histogramName)
137 BgMCHistograms = []
138 SignalMCHistograms = []
139 DataHistograms = []
140
141 backgroundIntegral = 0
142 dataIntegral = 0
143 scaleFactor = 1
144
145 for sample in processed_datasets: # loop over different samples as listed in configurationOptions.py
146 dataset_file = "%s/%s.root_tmp" % (condor_dir,sample)
147 inputFile = TFile(dataset_file)
148 Histogram = inputFile.Get(rootDirectory+"/"+channel+"/"+histogramName).Clone()
149 Histogram.SetDirectory(0)
150 if arguments.rebinFactor:
151 RebinFactor = int(arguments.rebinFactor)
152 if Histogram.GetNbinsX() >= RebinFactor*10:
153 Histogram.Rebin(RebinFactor)
154 inputFile.Close()
155 xAxisLabel = Histogram.GetXaxis().GetTitle()
156 histoTitle = Histogram.GetTitle()
157
158 if( types[sample] == "bgMC"):
159
160 numBgMCSamples += 1
161
162 if(arguments.noStack):
163 Histogram.SetFillStyle(0)
164 Histogram.SetLineColor(colors[sample])
165 Histogram.SetLineWidth(2)
166 BgMCLegend.AddEntry(Histogram,labels[sample],"L")
167 else:
168 Histogram.SetFillStyle(1001)
169 Histogram.SetFillColor(colors[sample])
170 Histogram.SetLineColor(1)
171 Histogram.SetLineWidth(1)
172 BgMCLegend.AddEntry(Histogram,labels[sample],"F")
173 Histogram.SetLineStyle(1)
174
175 backgroundIntegral += Histogram.Integral()
176
177 BgMCHistograms.append(Histogram)
178
179 elif( types[sample] == "signalMC"):
180
181 numSignalSamples += 1
182
183 Histogram.SetFillStyle(0)
184 Histogram.SetLineColor(colors[sample])
185 Histogram.SetLineStyle(1)
186 Histogram.SetLineWidth(2)
187 if(arguments.normalizeToUnitArea and Histogram.Integral() > 0):
188 Histogram.Scale(1./Histogram.Integral())
189 SignalMCLegend.AddEntry(Histogram,labels[sample],"L")
190 SignalMCHistograms.append(Histogram)
191
192 elif( types[sample] == "data"):
193
194 numDataSamples += 1
195
196 Histogram.SetFillStyle(0)
197 Histogram.SetLineColor(colors[sample])
198 Histogram.SetLineStyle(1)
199 Histogram.SetLineWidth(2)
200 if(arguments.normalizeToUnitArea and Histogram.Integral() > 0):
201 Histogram.Scale(1./Histogram.Integral())
202
203 dataIntegral += Histogram.Integral()
204
205 BgMCLegend.AddEntry(Histogram,labels[sample],"LEP")
206 DataHistograms.append(Histogram)
207
208 if dataIntegral > 0 and backgroundIntegral > 0:
209 scaleFactor = dataIntegral/backgroundIntegral
210 for bgMCHist in BgMCHistograms:
211 if arguments.normalizeToData:
212 bgMCHist.Scale(scaleFactor)
213 if arguments.normalizeToUnitArea and not arguments.noStack and backgroundIntegral > 0:
214 bgMCHist.Scale(1./backgroundIntegral)
215 elif arguments.normalizeToUnitArea and arguments.noStack and bgMCHist.Integral() > 0:
216 bgMCHist.Scale(1./bgMCHist.Integral())
217 if not arguments.noStack:
218 Stack.Add(bgMCHist)
219
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 if not arguments.noStack:
275 Stack.SetTitle(histoTitle)
276 Stack.Draw("HIST")
277 Stack.GetXaxis().SetTitle(xAxisLabel)
278 Stack.SetMaximum(1.1*finalMax)
279 Stack.SetMinimum(0.0001)
280 if makeRatioPlots or makeDiffPlots:
281 Stack.GetHistogram().GetXaxis().SetLabelSize(0)
282 else:
283 BgMCHistograms[0].SetTitle(histoTitle)
284 BgMCHistograms[0].Draw("HIST")
285 BgMCHistograms[0].GetXaxis().SetTitle(xAxisLabel)
286 BgMCHistograms[0].SetMaximum(1.1*finalMax)
287 BgMCHistograms[0].SetMinimum(0.0001)
288 for bgMCHist in BgMCHistograms:
289 bgMCHist.Draw("HIST SAME")
290 for signalMCHist in SignalMCHistograms:
291 signalMCHist.Draw("HIST SAME")
292 for dataHist in DataHistograms:
293 dataHist.Draw("E SAME")
294
295 elif(numSignalSamples is not 0):
296 SignalMCHistograms[0].SetTitle(histoTitle)
297 SignalMCHistograms[0].Draw("HIST")
298 SignalMCHistograms[0].GetXaxis().SetTitle(xAxisLabel)
299 SignalMCHistograms[0].SetMaximum(1.1*finalMax)
300 SignalMCHistograms[0].SetMinimum(0.0001)
301 for signalMCHist in SignalMCHistograms:
302 if(signalMCHist is not SignalMCHistograms[0]):
303 signalMCHist.Draw("HIST SAME")
304 for dataHist in DataHistograms:
305 dataHist.Draw("E SAME")
306
307 elif(numDataSamples is not 0):
308 DataHistograms[0].SetTitle(histoTitle)
309 DataHistograms[0].Draw("E")
310 DataHistograms[0].GetXaxis().SetTitle(xAxisLabel)
311 DataHistograms[0].SetMaximum(1.1*finalMax)
312 DataHistograms[0].SetMinimum(0.0001)
313 for dataHist in DataHistograms:
314 if(dataHist is not DataHistograms[0]):
315 dataHist.Draw("E SAME")
316
317
318 if(numBgMCSamples is not 0 or numDataSamples is not 0):
319 BgMCLegend.Draw()
320 if(numSignalSamples is not 0):
321 SignalMCLegend.Draw()
322
323 if not arguments.normalizeToUnitArea or numDataSamples > 0:
324 LumiLabel.Draw()
325 if arguments.normalizeToData and numBgMCSamples > 0 and numDataSamples > 0:
326 NormLabel = TPaveLabel(0.1,0.75,0.35,0.85,"MC scaled to data","NDC")
327 NormLabel.SetBorderSize(0)
328 NormLabel.SetFillColor(0)
329 NormLabel.SetFillStyle(0)
330 NormLabel.Draw()
331 elif arguments.normalizeToUnitArea:
332 NormLabel = TPaveLabel(0.1,0.75,0.35,0.85,"Scaled to unit area","NDC")
333 NormLabel.SetBorderSize(0)
334 NormLabel.SetFillColor(0)
335 NormLabel.SetFillStyle(0)
336 NormLabel.Draw()
337
338
339 if makeRatioPlots or makeDiffPlots:
340 Canvas.cd(2)
341 BgSum = Stack.GetStack().Last()
342 Comparison = DataHistograms[0].Clone()
343 Comparison.Add(BgSum,-1)
344 if not makeDiffPlots:
345 Comparison.Divide(BgSum)
346 Comparison.SetTitle("")
347 Comparison.GetXaxis().SetTitle(xAxisLabel)
348 if makeRatioPlots:
349 Comparison.GetYaxis().SetTitle("#frac{Data-MC}{MC}")
350 elif makeDiffPlots:
351 Comparison.GetYaxis().SetTitle("Data-MC")
352 Comparison.GetYaxis().CenterTitle()
353 Comparison.GetYaxis().SetTitleSize(0.1)
354 Comparison.GetYaxis().SetTitleOffset(0.35)
355 Comparison.GetXaxis().SetTitleSize(0.15)
356 Comparison.GetYaxis().SetLabelSize(0.1)
357 Comparison.GetXaxis().SetLabelSize(0.15)
358 if makeRatioPlots:
359 Comparison.GetYaxis().SetRangeUser(-1,1)
360 elif makeDiffPlots:
361 YMax = Comparison.GetMaximum()
362 YMin = Comparison.GetMinimum()
363 if YMax <= 0 and YMin <= 0:
364 Comparison.GetYaxis().SetRangeUser(-1.2*YMin,0)
365 elif YMax >= 0 and YMin >= 0:
366 Comparison.GetYaxis().SetRangeUser(0,1.2*YMax)
367 else: #axis crosses y=0
368 if abs(YMax) > abs(YMin):
369 Comparison.GetYaxis().SetRangeUser(-1.2*YMax,1.2*YMax)
370 else:
371 Comparison.GetYaxis().SetRangeUser(-1.2*YMin,1.2*YMin)
372
373 Comparison.GetYaxis().SetNdivisions(205)
374 Comparison.Draw()
375
376 Canvas.Write()
377
378
379 if re.match ('TH2', key.GetClassName()) and arguments.draw2DPlots: # plot a 2-D histogram
380
381 numBgMCSamples = 0
382 numDataSamples = 0
383 numSignalSamples = 0
384
385 if(intLumi < 1000.):
386 LumiText = "L_{int} = " + str(intLumi) + " pb^{-1}"
387 else:
388 getcontext().prec = 2
389 LumiInFb = intLumi/1000.
390 LumiText = "L_{int} = " + str(LumiInFb) + " fb^{-1}"
391
392 LumiLabel = TPaveLabel(0.1,0.8,0.34,0.9,LumiText,"NDC")
393 LumiLabel.SetBorderSize(0)
394 LumiLabel.SetFillColor(0)
395 LumiLabel.SetFillStyle(0)
396
397 BgMCLegend = TLegend(0.76,0.65,0.99,0.9, "Data & Bkgd. MC")
398 BgMCLegend.SetBorderSize(0)
399 BgMCLegend.SetFillColor(0)
400 BgMCLegend.SetFillStyle(0)
401 SignalMCLegend = TLegend(0.76,0.135,0.99,0.377,"Signal MC")
402 SignalMCLegend.SetBorderSize(0)
403 SignalMCLegend.SetFillColor(0)
404 SignalMCLegend.SetFillStyle(0)
405
406 outputFile.cd(rootDirectory+"/"+channel)
407 Canvas = TCanvas(histogramName)
408 Canvas.SetRightMargin(0.2413793);
409 BgMCHistograms = []
410 SignalMCHistograms = []
411 DataHistograms = []
412
413 for sample in processed_datasets: # loop over different samples as listed in configurationOptions.py
414 dataset_file = "%s/%s.root_tmp" % (condor_dir,sample)
415 inputFile = TFile(dataset_file)
416 Histogram = inputFile.Get(rootDirectory+"/"+channel+"/"+histogramName).Clone()
417 Histogram.SetDirectory(0)
418 RebinFactor = int(arguments.rebinFactor)
419 if arguments.rebinFactor and Histogram.GetNbinsX() >= RebinFactor*10 and Histogram.GetNbinsY() >= RebinFactor*10:
420 Histogram.Rebin2D(RebinFactor)
421 inputFile.Close()
422 xAxisLabel = Histogram.GetXaxis().GetTitle()
423 yAxisLabel = Histogram.GetYaxis().GetTitle()
424 histoTitle = Histogram.GetTitle()
425
426 if( types[sample] == "bgMC"):
427
428 numBgMCSamples += 1
429 Histogram.SetMarkerColor(colors[sample])
430 Histogram.SetFillColor(colors[sample])
431 BgMCLegend.AddEntry(Histogram,labels[sample],"F")
432 BgMCHistograms.append(Histogram)
433
434 elif( types[sample] == "signalMC"):
435
436 numSignalSamples += 1
437 Histogram.SetMarkerColor(colors[sample])
438 Histogram.SetFillColor(colors[sample])
439 SignalMCLegend.AddEntry(Histogram,labels[sample],"F")
440 SignalMCHistograms.append(Histogram)
441
442 elif( types[sample] == "data"):
443
444 numDataSamples += 1
445 Histogram.SetMarkerColor(colors[sample])
446 Histogram.SetFillColor(colors[sample])
447 BgMCLegend.AddEntry(Histogram,labels[sample],"F")
448 DataHistograms.append(Histogram)
449
450
451 outputFile.cd(rootDirectory+"/"+channel)
452
453 if(numBgMCSamples is not 0):
454 BgMCHistograms[0].SetTitle(histoTitle)
455 BgMCHistograms[0].GetXaxis().SetTitle(xAxisLabel)
456 BgMCHistograms[0].GetYaxis().SetTitle(yAxisLabel)
457 BgMCHistograms[0].Draw()
458 for signalMCHist in SignalMCHistograms:
459 signalMCHist.Draw("SAME")
460 for dataHist in DataHistograms:
461 dataHist.Draw("SAME")
462
463 elif(numSignalSamples is not 0):
464 SignalMCHistograms[0].SetTitle(histoTitle)
465 SignalMCHistograms[0].Draw()
466 SignalMCHistograms[0].GetXaxis().SetTitle(xAxisLabel)
467 SignalMCHistograms[0].GetYaxis().SetTitle(yAxisLabel)
468 for signalMCHist in SignalMCHistograms:
469 if(signalMCHist is not SignalMCHistograms[0]):
470 signalMCHist.Draw("SAME")
471 for dataHist in DataHistograms:
472 dataHist.Draw("SAME")
473
474 elif(numDataSamples is not 0):
475 DataHistograms[0].SetTitle(histoTitle)
476 DataHistograms[0].GetXaxis().SetTitle(xAxisLabel)
477 DataHistograms[0].GetYaxis().SetTitle(yAxisLabel)
478 DataHistograms[0].Draw()
479 for dataHist in DataHistograms:
480 if(dataHist is not DataHistograms[0]):
481 dataHist.Draw("SAME")
482
483
484 if(numBgMCSamples is not 0 or numDataSamples is not 0):
485 BgMCLegend.Draw()
486 if(numSignalSamples is not 0):
487 SignalMCLegend.Draw()
488 if not arguments.normalizeToUnitArea or numDataSamples > 0:
489 LumiLabel.Draw()
490
491 Canvas.Write()
492
493
494
495 for dataset in processed_datasets:
496 dataset_file = "%s/%s.root_tmp" % (condor_dir,dataset)
497 os.remove(dataset_file)
498
499 outputFile.Close()