plot_validation_datamc.py
Go to the documentation of this file.
1 #!/usr/bin/env python
2 
3 import argparse
4 import code
5 import copy
6 import re
7 import math
8 
9 gLeaked = []
10 def New(cons, *args):
11  ret = cons(*args)
12  gLeaked.append(ret)
13  return ret
14 
15 def Clone(obj):
16  ret = obj.Clone()
17  gLeaked.append(ret)
18  return ret
19 
20 def POTStr(pot):
21  return '%.3g#times10^{20} POT' % (pot/1e20)
22 
23 def DeriveScalePow(spect, pot):
24  """Next smallest power of 1000 than maximum bin"""
25  if spect.hist.GetMaximum() == 0: return 0
26  return 3*(int(math.log10(spect.hist.GetMaximum()*pot/spect.pot))/3)
27 
28 def One():
29  g = New(TGraph)
30  g.SetPoint(0, -1e10, 1)
31  g.SetPoint(1, 1e10, 1)
32  g.SetLineWidth(2)
33  g.SetLineStyle(7)
34  g.Draw('same')
35 
36 class Spectrum:
37  def __init__(self, dirfile):
38  self.hist = dirfile.Get('hist')
39  self.pot = dirfile.Get('pot').Integral()
40 
41  label = dirfile.Get('label0').GetString().Data()
42  self.hist.GetXaxis().SetTitle(label)
43  self.hist.GetXaxis().CenterTitle()
44  self.hist.GetYaxis().CenterTitle()
45 
46  def Draw(self, pot, opt = '', scalePow = 0, split = False):
47  hc = Clone(self.hist)
48 
49  if split:
50  hc.GetXaxis().SetLabelSize(0)
51  hc.GetXaxis().SetTitle('')
52 
53  hc.Scale(pot/self.pot)
54  hc.Scale(1./pow(10, scalePow))
55  title = 'Events / '+POTStr(pot)
56  if scalePow != 0: title = '10^{'+str(scalePow)+'} '+title
57  hc.GetYaxis().SetTitle(title)
58  hc.Draw(opt)
59 
60  def __add__(self, other):
61  ret = copy.copy(self)
62  ret.hist = self.hist.Clone()
63  ret.hist.Add(other.hist, other.pot/ret.pot)
64  return ret
65 
66 
67 class Prediction:
68  def __init__(self, dirfile, basename):
69  (self.nue, self.numu, self.antinue, self.antinumu, self.nc) = [Spectrum(dirfile.GetDirectory(basename+'{group=Component,cat='+c+'}')) for c in ['0_nue', '1_numu', '3_antinue', '4_antinumu', '2_nc']]
70 
71  self.tot = self.nue + self.numu + self.antinue + self.antinumu + self.nc
72 
73  # Analysis/Style.h
74  self.tot.hist.SetLineColor(kRed)
75  self.nue.hist.SetLineColor(kMagenta)
76  self.numu.hist.SetLineColor(kGreen+2)
77  self.antinue.hist.SetLineColor(kMagenta-2)
78  self.antinumu.hist.SetLineColor(kGreen+1)
79  self.nc.hist.SetLineColor(kAzure)
80 
81  self.antinue.hist.SetLineStyle(7)
82  self.antinumu.hist.SetLineStyle(7)
83 
84  def Draw(self, pot, opt = '', scalePow = 0, split = False):
85  for h in self.tot, self.nue, self.numu, self.antinue, self.antinumu, self.nc:
86  h.Draw(pot, opt, scalePow, split)
87 
88 
89 class Ratio:
90  def __init__(self, numSpect, denomSpect):
91  self.hist = numSpect.hist.Clone()
92  self.hist.Divide(denomSpect.hist)
93  self.hist.Scale(denomSpect.pot / numSpect.pot)
94 
95  def Draw(self):
96  self.hist.GetYaxis().SetTitle('Data / MC')
97  self.hist.GetYaxis().SetRangeUser(.5, 1.5)
98  self.hist.Draw('ep')
99 
100 
101 # Ported from CAFAna/Analysis/Plots.cxx
102 def SplitCanvas(ysplit):
103  if gPad is None: New(TCanvas)
104  p1 = New(TPad, '', '', 0, 0, 1, 1)
105  p2 = New(TPad, '', '', 0, 0, 1, 1)
106 
107  p1.SetBottomMargin(ysplit)
108  p2.SetTopMargin(1-ysplit)
109 
110  # Draw p1 second since it's often the more important one, that the user
111  # would prefer to be able to interact with.
112  for p in [p2, p1]:
113  p.SetFillStyle(0)
114  p.Draw()
115 
116  return (p1, p2)
117 
118 
119 
120 def Legend():
121  leg = New(TLegend, .4, .25, .6, .75)
122  h = New(TH1F, '', '', 1, 0, 1)
123  leg.SetTextSize(h.GetXaxis().GetTitleSize())
124  h.SetMarkerStyle(kFullCircle)
125  leg.AddEntry(Clone(h), 'ND data', 'lep')
126  h.SetLineColor(kRed)
127  leg.AddEntry(Clone(h), 'Total MC', 'l')
128  h.SetLineColor(kAzure)
129  leg.AddEntry(Clone(h), 'NC', 'l')
130  h.SetLineColor(kGreen+2)
131  leg.AddEntry(Clone(h), '#nu_{#mu} CC', 'l')
132  h.SetLineColor(kMagenta)
133  leg.AddEntry(Clone(h), '#nu_{e} CC', 'l')
134  h.SetLineColor(kGreen+1)
135  h.SetLineStyle(7)
136  leg.AddEntry(Clone(h), '#bar{#nu_{#mu}} CC', 'l')
137  h.SetLineColor(kMagenta-2)
138  leg.AddEntry(Clone(h), '#bar{#nu_{e}} CC', 'l')
139  leg.Draw()
140 
141 
142 parser = argparse.ArgumentParser()
143 parser.add_argument('-b', '--batch', action = 'store_true',
144  help = 'batch mode, no graphics output')
145 parser.add_argument('dataFile', help = 'Validation histogram file for data')
146 parser.add_argument('mcFile', help = 'Validation histogram file for MC')
147 opts = parser.parse_args()
148 
149 # Screws up the argument parsing if it's any earlier
150 from ROOT import *
151 
152 # Style the plots nicely
153 gApplication.ExecuteFile('$SRT_PUBLIC_CONTEXT/Utilities/rootlogon.C')
154 gROOT.ForceStyle()
155 
156 if opts.batch: gROOT.SetBatch(True)
157 
158 Legend()
159 gPad.Print('plots/legend.pdf')
160 
161 fdata = TFile(opts.dataFile)
162 fmc = TFile(opts.mcFile)
163 
164 for k in fdata.GetListOfKeys():
165  print k.GetName()
166 
167  New(TCanvas)
168 
169  data = Spectrum(k.ReadObj())
170  if data.hist.Integral() == 0: continue
171  data.hist.SetMarkerStyle(kFullCircle)
172  scalePow = DeriveScalePow(data, data.pot)
173  data.Draw(data.pot, 'ep', scalePow)
174 
175  pred = Prediction(fmc, k.GetName())
176  pred.Draw(data.pot, 'hist same', scalePow)
177 
178  data.Draw(data.pot, 'ep same', scalePow) # data should be on top
179 
180  fname = re.sub(r'\{group=Cut,cat=._(.*)\}', r'_\1', k.GetName())
181 
182  gPad.SetLogy(True)
183  gPad.Print('plots/'+fname+'_log.pdf')
184  gPad.SetLogy(False)
185  gPad.Print('plots/'+fname+'.pdf')
186 
187  New(TCanvas)
188  r = Ratio(data, pred.tot)
189  r.Draw()
190  One()
191 
192  gPad.Print('plots/'+fname+'_ratio.pdf')
193 
194  # De-conflict the y-axis labels in the split plot case
195  data.hist.GetYaxis().SetTitleSize(.9*data.hist.GetYaxis().GetTitleSize())
196  r.hist.GetYaxis().SetTitleSize(.9*r.hist.GetYaxis().GetTitleSize())
197 
198  c = New(TCanvas)
199  (p1, p2) = SplitCanvas(.35)
200  p1.cd()
201  data.Draw(data.pot, 'ep', scalePow, True)
202  pred.Draw(data.pot, 'hist same', scalePow)
203  data.Draw(data.pot, 'ep same', scalePow)
204  p2.cd()
205  r.Draw()
206  One()
207  c.cd()
208 
209  p1.SetLogy(True)
210  gPad.Print('plots/'+fname+'_split_log.pdf')
211  p1.SetLogy(False)
212  gPad.Print('plots/'+fname+'_split.pdf')
213 
214  # Don't flood the screen with plots
215  if not gROOT.IsBatch(): break
216 
217 
218 if not opts.batch:
219  # Just want the console to hang around so the user can look at the plots
220  # without them disappearing. Ctrl-D or "exit()" to quit.
221  print 'Ctrl-D to quit'
222  code.interact(local = locals(), banner='')
double Integral(const Spectrum &s, const double pot, int cvnregion)
constexpr T pow(T x)
Definition: pow.h:75
def Draw(self, pot, opt='', scalePow=0, split=False)
def Draw(self, pot, opt='', scalePow=0, split=False)
gargamelle SetTitle("Gargamelle #nu_{e} CC data")
def __init__(self, numSpect, denomSpect)
string GetString(xmlDocPtr xml_doc, string node_path)
def __init__(self, dirfile, basename)