All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MakePlots.py
Go to the documentation of this file.
1 ## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2 #
3 # Plotting Outputs
4 #
5 # (from covariance and sensitivity
6 # calculations)
7 #
8 ## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
9 
10 import ROOT
11 from ROOT import TFile, TCanvas, TH2D, TH1D, TGraph, TGraph2D, TStyle, TLegend, THStack, TPad
12 import argparse
13 import os
14 import math
15 
16 
17 ## Main function
18 ## ~~~~~~~~~~~~~
19 
20 def make_count_plot(args):
21 
22  # Count plots
23 
24  countfile = TFile(args.cntfile)
25 
26  dets = []
27  for key in countfile.GetListOfKeys():
28  key = key.GetName()
29  if key[:key.find("_")] not in dets:
30  dets.append(key[:key.find("_")])
31 
32  samples = []
33  for key in countfile.GetListOfKeys():
34  key = key.GetName()
35  if key[key.find("_")+1:key.find("_")+1+key[key.find("_")+1:].find("_")] not in samples:
36  samples.append(key[key.find("_")+1:key.find("_")+1+key[key.find("_")+1:].find("_")])
37 
38  samplename = {"numu": "#nu_{#mu}", "nue": "#nu_{e}"}
39  canvases = []; stacks = []; ct_hists = []; bkg_hists = []; legends = []
40  for s, sample in enumerate(samples):
41 
42  canvases.append(TCanvas(sample+"_canvas"))
43  canvases[s].Divide(len(dets), 1)
44 
45  ct_hists.append([])
46  bkg_hists.append([])
47  stacks.append([])
48 
49  legends.append([])
50 
51  for d, det in enumerate(dets):
52 
53  ct_hists[s].append(countfile.Get(det+"_"+sample+"_cts"))
54  ct_hists[s][d].SetLineColor(38)
55  ct_hists[s][d].SetFillColor(38)
56 
57  bkg_hists[s].append(countfile.Get(det+"_"+sample+"_bkg"))
58  bkg_hists[s][d].SetLineColor(30)
59  bkg_hists[s][d].SetFillColor(30)
60 
61  for b in range(ct_hists[s][d].GetNbinsX()):
62  ct_hists[s][d].SetBinContent(1+b, ct_hists[s][d].GetBinContent(1+b) / ct_hists[s][d].GetBinWidth(1+b))
63  bkg_hists[s][d].SetBinContent(1+b, bkg_hists[s][d].GetBinContent(1+b) / bkg_hists[s][d].GetBinWidth(1+b))
64 
65  stacks[s].append(THStack(det+"_"+sample+"_stack", det+"; Counts; Energy (GeV)"))
66  stacks[s][d].Add(bkg_hists[s][d])
67  stacks[s][d].Add(ct_hists[s][d])
68 
69  legends[s].append(TLegend())
70  legends[s][d].AddEntry(ct_hists[s][d], 'CC -> #mu + X', 'f')
71  legends[s][d].AddEntry(bkg_hists[s][d], 'NC -> #pi^{#pm} + X', 'f')
72 
73  canvases[s].cd(d+1)
74  stacks[s][d].Draw("hist")
75  legends[s][d].Draw()
76 
77  # NOTE: This will have to be changed when nue samples are added; this assumed a simple
78  # count + background histogram stack, whereas nue plots in the proposal have the
79  # background broken up into many different 'types'
80 
81  canvases[s].SaveAs(args.outdir + sample + "_selection.pdf")
82 
83 
84 def plot_cov_output(args):
85 
86  # Covariance, fractional covariance and correlation
87 
88  covfile = TFile(args.covfile)
89 
90  covcanvas = TCanvas()
91 
92  gStyle = TStyle()
93  gStyle.SetPadLeftMargin(0.25); gStyle.SetPadRightMargin(0.15)
94  gStyle.SetPalette(56)
95 
96  for matname in ('cov', 'fcov', 'corr'):
97 
98  mat = covfile.Get(matname)
99 
100  mat.GetYaxis().LabelsOption('v') # doesn't work...
101  mat.GetYaxis().SetLabelSize(0.07)
102  mat.GetXaxis().SetLabelSize(0.07)
103 
104  mat.SetTitleSize(0.3, 't') # doesn't work...
105 
106  if matname == 'corr':
107  mat.GetZaxis().SetRangeUser(-0.4, 1)
108  mat.SetTitle("Flux Correlation Matrix")
109 
110  mat.Draw("colz")
111  mat.SetStats(False)
112  covcanvas.SetLeftMargin(0.25)
113  covcanvas.Update()
114  covcanvas.SaveAs(args.outdir + matname + "_plot.pdf")
115 
116 
117 def plot_chi2_output(args):
118 
119  # Chi squareds
120 
121  chi2file = TFile(args.chifile)
122 
123  gStyle = TStyle()
124 
125  chi2 = chi2file.Get('chisq')
126 
127  chi2canvas = TCanvas()
128 
129  chi2.SetTitle('#chi^{2}; log_{10}(sin^{2}(2#theta)); log_{10}(#Delta m^{2}); #chi^{2}');
130  gStyle.SetPalette(1)
131  chi2.Draw('surf1')
132  chi2canvas.SaveAs(args.outdir + "chisq.pdf")
133 
134 
135  # Contours
136  contcanvas = TCanvas('cont_canvas', '', 1020, 990)
137  gStyle.SetPadLeftMargin(0.15); gStyle.SetPadRightMargin(0.15)
138 
139  colours = [30, 38, 46]
140  contours = [chi2file.Get('90pct'),
141  chi2file.Get('3sigma'),
142  chi2file.Get('5sigma')]
143 
144  for g in range(len(contours)):
145 
146  contours[g].SetMarkerStyle(20)
147  contours[g].SetMarkerSize(0.25)
148  contours[g].SetMarkerColor(colours[g])
149  contours[g].SetLineColor(colours[g])
150 
151  gr_range = TGraph()
152  gr_range.SetPoint(0, 0.001, 0.01)
153  gr_range.SetPoint(1, 1, 100)
154  gr_range.SetMarkerColor(0)
155 
156  bestfit = TGraph()
157  bestfit.SetPoint(0, 0.062, 1.7)
158  bestfit.SetMarkerStyle(29)
159  bestfit.SetMarkerSize(1.6)
160  bestfit.SetMarkerColor(40)
161 
162  gr_range.SetTitle('SBN Sensitivity; sin^{2}(2#theta); #Delta m^{2} (eV^{2})')
163 
164  legend = TLegend()
165  legend.AddEntry(contours[0], '90% CL', 'l')
166  legend.AddEntry(contours[1], '3#sigma CL', 'l')
167  legend.AddEntry(contours[2], '5#sigma CL', 'l')
168  legend.AddEntry(bestfit, 'Best Fit Point', 'p')
169 
170  contcanvas.SetLogy()
171  contcanvas.SetLogx()
172 
173  gr_range.Draw('AP')
174  gr_range.GetXaxis().SetRangeUser(0.001, 1)
175  gr_range.GetYaxis().SetRangeUser(0.01, 100)
176 
177  contours[0].Draw('P same')
178  contours[1].Draw('P same')
179  contours[2].Draw('P same')
180 
181  legend.Draw()
182  bestfit.Draw('P same')
183 
184  contcanvas.SetLeftMargin(0.15)
185  contcanvas.Update()
186  contcanvas.SaveAs(args.outdir+'Sensitivity.pdf')
187 
188 
189 def compare_w_proposal(args):
190 
191  chi2file = TFile(args.chifile)
192 
193  gStyle = TStyle()
194  gStyle.SetPadLeftMargin(0.15); gStyle.SetPadRightMargin(0.15)
195 
196  colours = [30, 38, 46]
197  contours = [chi2file.Get('90pct'),
198  chi2file.Get('3sigma'),
199  chi2file.Get('5sigma')]
200 
201  propcontours = []
202  contournames = ['90pct', '3s', '5s']
203  contourtitles = ['90% Confidence Level', '3#sigma Confidence Level', '5#sigma Confidence Level']
204 
205  gr_range = TGraph()
206  gr_range.SetPoint(0, 0.001, 0.01)
207  gr_range.SetPoint(1, 1, 100)
208  gr_range.SetMarkerColor(0)
209 
210  bestfit = TGraph()
211  bestfit.SetPoint(0, 0.062, 1.7)
212  bestfit.SetMarkerStyle(29)
213  bestfit.SetMarkerSize(1.6)
214  bestfit.SetMarkerColor(40)
215 
216  print("contours has length " + str(len(contours)))
217 
218  for i in range(len(contours)):
219 
220  x = []
221  y = []
222 
223  with open(args.compdir+'numu'+contournames[i]+'.txt') as f:
224  for l, line in enumerate(f):
225  if l == 0: continue
226  x.append(float(line.split(', ')[0]))
227  y.append(float(line.split(', ')[1].replace("\n", "")))
228 
229  propcontours.append(TGraph())
230  for j in range(len(x)):
231  propcontours[i].SetPoint(j, x[j], y[j])
232 
233  tempcanvas = TCanvas('temp_canvas', '', 1020, 990)
234 
235  templegend = TLegend()
236  templegend.AddEntry(contours[i], 'Our contour', 'l')
237  templegend.AddEntry(propcontours[i], 'From proposal', 'l')
238  templegend.AddEntry(bestfit, 'Best Fit Point', 'p')
239 
240  tempcanvas.SetLogy()
241  tempcanvas.SetLogx()
242 
243  gr_range.SetTitle(contourtitles[i]+' Comparison; sin^{2}(2#theta); #Delta m^{2} (eV^{2})')
244 
245  gr_range.Draw('AP')
246  gr_range.GetXaxis().SetRangeUser(0.001, 1)
247  gr_range.GetYaxis().SetRangeUser(0.01, 100)
248 
249  for lst in (contours, propcontours):
250  lst[i].SetMarkerStyle(20)
251  lst[i].SetMarkerSize(0.25)
252  lst[i].SetMarkerColor(colours[i] if lst == contours else 1)
253  lst[i].SetLineColor(colours[i] if lst == contours else 1)
254 
255  contours[i].Draw('P same')
256  propcontours[i].Draw('P same')
257 
258  templegend.Draw()
259  bestfit.Draw('P same')
260 
261  tempcanvas.SaveAs(args.outdir+contournames[i]+'_comparison.pdf')
262 
263 
264 def dm2_chi2_slice(args):
265 
266  dm2s = [float(dm2) for dm2 in args.dm2list.split(",")]
267 
268  chi2file = TFile(args.chifile)
269  chi2 = chi2file.Get('chisq')
270 
271  min_sin = chi2.GetXmin(); max_sin = chi2.GetXmax()
272  min_dm2 = chi2.GetYmin(); max_dm2 = chi2.GetYmax()
273 
274  chi2_vals = [float(chi) for chi in chi2.GetZ()]
275  NP = int(math.sqrt(len(chi2_vals)))
276 
277  reusable_canvas = TCanvas()
278  print(dm2s)
279  print("min_sin = " + str(min_sin) + " and max_sin = " + str(max_sin))
280  print("min_dm2 = " + str(min_dm2) + " and max_dm2 + " + str(max_dm2))
281  print("NP + " + str(NP))
282  slices = []
283  for m, dm2 in enumerate(dm2s):
284 
285  jval = int((math.log10(dm2) - min_dm2)*(NP-1)/(max_dm2-min_dm2))
286  tempslice = [chi2_vals[i*NP + jval] for i in range(NP)]
287 
288  slices.append(TGraph())
289  for c, chi in enumerate(tempslice):
290  slices[m].SetPoint(c, (min_sin + c*(max_sin - min_sin)/(NP-1)), chi)
291 
292  if args.chilims:
293  chilims = [float(lim) for lim in args.chilims.split(",")]
294  if len(chilims) == 1:
295  slices[m].GetYaxis().SetRangeUser(0, chilims[0])
296  else:
297  slices[m].GetYaxis().SetRangeUser(chilims[0], chilims[1])
298 
299  slices[m].SetTitle("#chi^{2} @ #Delta m^{2} = " + str(dm2) + "; log_{10}(sin^{2}(2#theta)); #chi^{2}")
300  slices[m].Draw("AP")
301 
302  reusable_canvas.SaveAs(args.outdir + "dm2_"+str(dm2).replace(".", "-")+"_slice.pdf")
303 
304 def sin_chi2_slice(args):
305 
306  sins = [float(sin) for sin in args.sinlist.split(",")]
307 
308  chi2file = TFile(args.chifile)
309  chi2 = chi2file.Get('chisq')
310 
311  min_sin = chi2.GetXmin(); max_sin = chi2.GetXmax()
312  min_dm2 = chi2.GetYmin(); max_dm2 = chi2.GetYmax()
313 
314  chi2_vals = [float(chi) for chi in chi2.GetZ()]
315  NP = int(math.sqrt(len(chi2_vals)))
316 
317  reusable_canvas = TCanvas()
318 
319  slices = []
320  for s, sin in enumerate(sins):
321 
322  ival = int((math.log10(sin) - min_sin)*(NP-1)/(max_sin-min_sin))
323  tempslice = [chi2_vals[ival*NP + j] for j in range(NP)]
324 
325  slices.append(TGraph())
326  for c, chi in enumerate(tempslice):
327  slices[s].SetPoint(c, 10**(min_dm2 + c*(max_dm2 - min_dm2)/(NP-1)), chi)
328 
329  if args.chilims:
330  chilims = [float(lim) for lim in args.chilims.split(",")]
331  if len(chilims) == 1:
332  slices[s].GetYaxis().SetRangeUser(0, chilims[0])
333  else:
334  slices[s].GetYaxis().SetRangeUser(chilims[0], chilims[1])
335 
336  slices[s].SetTitle("#chi^{2} @ #sin^{2}(2#theta) = " + str(sin) + "; log_{10}(#Delta m^{2}); #chi^{2}")
337  slices[s].Draw("AP")
338  reusable_canvas.SaveAs(args.outdir + "sin_"+str(sin).replace(".", "-")+"_slice.pdf")
339 
340 
341 if __name__ == "__main__":
342 
343  buildpath = os.environ['SBN_LIB_DIR']
344  if not buildpath:
345  print('ERROR: SBNDDAQ_ANALYSIS_BUILD_PATH not set')
346  sys.exit()
347 
348  ROOT.gROOT.ProcessLine('.L ' + buildpath + '/libsbnanalysis_Event.so')
349  ROOT.gROOT.ProcessLine('.L ' + buildpath + '/libsbnanalysis_SBNOsc_classes.so')
350 
351  parser = argparse.ArgumentParser()
352  parser.add_argument("-chi", "--chifile", default = False)
353  parser.add_argument("-cov", "--covfile", default = False)
354  parser.add_argument("-cts", "--cntfile", default = False)
355  parser.add_argument("-o", "--outdir", required = True)
356  parser.add_argument("-comp", "--compdir", default = False)
357  parser.add_argument("-dm2slice", "--dm2list", default = False)
358  parser.add_argument("-sinslice", "--sinlist", default = False)
359  parser.add_argument("-chilims", "--chilims", default = False)
360 
361  args = parser.parse_args()
362 
363  if args.cntfile: make_count_plot(args)
364  if args.covfile: plot_cov_output(args)
365  if args.chifile: plot_chi2_output(args)
366  if args.compdir: compare_w_proposal(args)
367  if args.dm2list: dm2_chi2_slice(args)
368  if args.sinlist: sin_chi2_slice(args)
369 
370  #if parser.parse_args().compare:
371  # compare_w_proposal(parser_parse_args())
372  # with open('filename') as f:
373  # for line in f:
374  # data = [float(x) for x in line.split(",")]
375 
376 
377 
378 
print
do one_file $F done echo for F in find $TOP name CMakeLists txt print
Definition: update_sources.sh:94
MakePlots.dm2_chi2_slice
def dm2_chi2_slice
Definition: MakePlots.py:264
MakePlots.sin_chi2_slice
def sin_chi2_slice
Definition: MakePlots.py:304
util::enumerate
auto enumerate(Iterables &&...iterables)
Range-for loop helper tracking the number of iteration.
Definition: enumerate.h:69
MakePlots.compare_w_proposal
def compare_w_proposal
Definition: MakePlots.py:189
MakePlots.make_count_plot
def make_count_plot
Main function ~~~~~~~~~~~~~.
Definition: MakePlots.py:20
MakePlots.plot_cov_output
def plot_cov_output
Definition: MakePlots.py:84
MakePlots.plot_chi2_output
def plot_chi2_output
Definition: MakePlots.py:117
cd
do cd
Definition: getStatistics.sh:62
open
open(RACETRACK) or die("Could not open file $RACETRACK for writing")
Generated by   doxygen 1.8.5