plot_validation_datamc_2018.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 import ROOT
9 
10 gLeaked = []
11 def Caption(fName, axisTitle, decomp):
12  beamMode = 'neutrino'
13  if 'rhc' in fName:
14  beamMode = 'antineutrino'
15 
16  lcvn = 'CVNeLow' in fName
17  hcvn = 'CVNeHigh' in fName
18  cut = 'preselection cuts'
19  if lcvn:
20  cut = 'low CVN cut'
21  if hcvn:
22  cut = 'high CVN cut'
23 
24  decompStr = ''
25  if decomp:
26  if 'rhc' in fName:
27  decompType = 'proportional'
28  else:
29  decompType = 'combo'
30  decompStr = 'Solid line indicates simulation has been corrected by {} decomposition.'.format(decompType)
31 
32  caption = '{} in {} beam mode with {}. Top panel shows data (markers) and simulation (colors) spectra. {} Bottom panel shows the ratio of data to simulation.'.format(axisTitle, beamMode, cut, decompStr)
33 
34  txtName = fName.replace('pdf', 'txt')
35  textFile = open(txtName, 'w')
36  textFile.write(caption)
37  textFile.close()
38 
39 def New(cons, *args):
40  ret = cons(*args)
41  gLeaked.append(ret)
42  return ret
43 
44 def Clone(obj):
45  ret = obj.Clone()
46  gLeaked.append(ret)
47  return ret
48 
49 def POTStr(pot):
50  return '%.3g#times10^{20} POT' % (pot/1e20)
51 
52 def DeriveScalePow(spect, pot):
53  """Next smallest power of 1000 than maximum bin"""
54  if spect.hist.GetMaximum() == 0: return 0
55  return 3*(int(math.log10(spect.hist.GetMaximum()*pot/spect.pot))/3)
56 
57 def One():
58  g = New(TGraph)
59  g.SetPoint(0, -1e10, 1)
60  g.SetPoint(1, 1e10, 1)
61  g.SetLineWidth(2)
62  g.SetLineStyle(7)
63  g.Draw('same')
64 
66  tex = New(TLatex, .9, .95, "NOvA Preliminary")
67  tex.SetTextColor(kBlue)
68  tex.SetNDC()
69  tex.SetTextSize(2/30.)
70  tex.SetTextAlign(32)
71  tex.Draw()
72 
73 def BeamModeLabel(beamMode = 'fhc'):
74  txt = 'Neutrino Mode'
75  if 'rhc' in beamMode:
76  txt = 'Antineutrino Mode'
77  lb = New(TLatex, .1, .95, txt)
78  lb.SetNDC()
79  lb.SetTextAlign(13)
80  lb.SetTextColor(kGray+3)
81  lb.SetTextSize(2/45.)
82  lb.Draw()
83 
84 class Spectrum:
85  def __init__(self, dirfile):
86  self.hist = dirfile.Get('hist')
87  self.pot = dirfile.Get('pot').Integral()
88 
89  label = dirfile.Get('label0').GetString().Data()
90  self.hist.GetXaxis().SetTitle(label)
91  self.hist.GetXaxis().CenterTitle()
92  self.hist.GetYaxis().CenterTitle()
93 
94  # Originally made with 200 bins. Way too many
95  if 'shw_cos' in dirfile.GetName():
96  self.hist.Rebin(4)
97 
98 
99 
100  def Clone(self):
101  return copy.copy(self)
102 
103  def Draw(self, pot, opt = '', scalePow = 0, split = False):
104  hc = Clone(self.hist)
105 
106  if split:
107  hc.GetXaxis().SetLabelSize(0)
108  hc.GetXaxis().SetTitle('')
109 
110  hc.Scale(pot/self.pot)
111  hc.Scale(1./pow(10, scalePow))
112  title = 'Events / '+POTStr(pot)
113  if scalePow != 0: title = '10^{'+str(scalePow)+'} '+title
114  hc.GetYaxis().SetTitle(title)
115 
116  # make some room for the legend
117  minY = hc.GetMinimum()
118  maxY = hc.GetMaximum()
119  # root bug breaks log plotting if minY <= 0
120  hc.GetYaxis().SetRangeUser(minY+0.00001, 1.35*maxY)
121  hc.Draw(opt)
122 
123  def __add__(self, other):
124  ret = copy.copy(self)
125  ret.hist = self.hist.Clone()
126  ret.hist.Add(other.hist, other.pot/ret.pot)
127  return ret
128 
129  def SetLineStyle(self, lineStyle):
130  self.hist.SetLineStyle(lineStyle)
131 
132 
134  def __init__(self, dirfile, basename):
135  (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']]
136 
137  self.tot = self.nue + self.numu + self.antinue + self.antinumu + self.nc
138 
139  # Analysis/Style.h
140  self.tot.hist.SetLineColor(kRed)
141  self.nue.hist.SetLineColor(kMagenta-2)
142  self.numu.hist.SetLineColor(kGreen+2)
143  self.antinue.hist.SetLineColor(kMagenta)
144  self.antinumu.hist.SetLineColor(kGreen-3)
145  self.nc.hist.SetLineColor(kAzure)
146 
147  def Draw(self, pot, opt = '', scalePow = 0, split = False, lStyle = 1):
148  for h in self.tot, self.nue, self.numu, self.antinue, self.antinumu, self.nc:
149  h.SetLineStyle(lStyle)
150  h.Draw(pot, opt, scalePow, split)
151 
152 
153 class Ratio:
154  def __init__(self, numSpect, denomSpect):
155  self.hist = numSpect.hist.Clone()
156  self.hist.Divide(denomSpect.hist)
157  self.hist.Scale(denomSpect.pot / numSpect.pot)
158 
159  def Draw(self, opt='ep', color=1):
160  self.hist.SetMarkerColor(color)
161  self.hist.SetLineColor(color)
162  self.hist.GetYaxis().SetTitle('Data / MC')
163  self.hist.GetYaxis().SetRangeUser(.5, 1.5)
164  self.hist.Draw(opt)
165 
166 
167 # Ported from CAFAna/Analysis/Plots.cxx
168 def SplitCanvas(ysplit):
169  if gPad is None: New(TCanvas)
170  p1 = New(TPad, '', '', 0, 0, 1, 1)
171  p2 = New(TPad, '', '', 0, 0, 1, 1)
172 
173  p1.SetBottomMargin(ysplit)
174  p2.SetTopMargin(1-ysplit)
175 
176  # Draw p1 second since it's often the more important one, that the user
177  # would prefer to be able to interact with.
178  for p in [p2, p1]:
179  p.SetFillStyle(0)
180  p.Draw()
181 
182  return (p1, p2)
183 
184 
185 def LegendSplit(spec, decomp = False, ratio = False):
186  y1 = .56
187  y2 = .85
188  if ratio:
189  y1 = .46
190 
191  # Make some room for the legend
192  minY = spec.hist.GetMinimum()
193  maxY = spec.hist.GetMaximinum()
194  spec.hist.GetYaxis().SetRangeUser(minY, 1.3*maxY)
195 
196  leg = New(TLegend, .6, y1, .95, y2)
197  leg.SetNColumns(2)
198  leg.SetBorderSize(1)
199  h = New(TH1F, '', '', 1, 0, 1)
200  leg.SetTextSize(.7*h.GetXaxis().GetTitleSize())
201  h.SetMarkerStyle(kFullCircle)
202  leg.AddEntry(Clone(h), 'ND data', 'lep') #1,1
203 
204  h.SetLineColor(kAzure)
205  leg.AddEntry(Clone(h), 'NC', 'l') #2,1
206 
207  h.SetLineColor(kRed)
208  leg.AddEntry(Clone(h), 'Total MC', 'l') #1,2
209 
210  if decomp:
211  h.SetLineStyle(7)
212  h.SetLineColor(kBlack)
213  leg.AddEntry(Clone(h), 'Uncorr. MC', 'l') #2,2
214 
215  h.SetLineStyle(1)
216  h.SetLineColor(kGreen+2)
217  leg.AddEntry(Clone(h), '#nu_{#mu} CC', 'l') #1,3
218 
219  h.SetLineColor(kGreen-3)
220  leg.AddEntry(Clone(h), '#bar{#nu_{#mu}} CC', 'l') #2,3
221 
222  h.SetLineColor(kMagenta-2)
223  leg.AddEntry(Clone(h), '#nu_{e} CC', 'l') #1,4
224 
225  h.SetLineColor(kMagenta)
226  leg.AddEntry(Clone(h), '#bar{#nu_{e}} CC', 'l') #2,4
227 
228  h.SetMarkerColor(1)
229  h.SetLineColor(1)
230  if decomp:
231  h.SetMarkerColor(17)
232  h.SetLineColor(17)
233  if ratio:
234  leg.AddEntry(Clone(h), '#splitline{Ratio to}{Uncorr. MC}') #1,5
235 
236  h.SetMarkerColor(1)
237  h.SetLineColor(1)
238  if decomp and ratio:
239  leg.AddEntry(Clone(h), '#splitline{Ratio to}{Decomp MC}') #2,5
240  leg.Draw()
241 
242 
243 
244 def Legend(isDecomp = False, isSplit = False):
245  leg = New(TLegend, .78, .4, .95, .89)
246  h = New(TH1F, '', '', 1, 0, 1)
247 
248  leg.SetMargin(.22)
249  leg.SetBorderSize(1)
250 
251  h.SetMarkerStyle(kFullCircle)
252  leg.AddEntry(Clone(h), 'ND data', 'lep')
253  h.SetLineColor(kRed)
254  leg.AddEntry(Clone(h), 'Total MC', 'l')
255  h.SetLineColor(kAzure)
256  leg.AddEntry(Clone(h), 'NC', 'l')
257  h.SetLineColor(kGreen+2)
258  leg.AddEntry(Clone(h), '#nu_{#mu} CC', 'l')
259  h.SetLineColor(kMagenta-2)
260  leg.AddEntry(Clone(h), '#nu_{e} CC', 'l')
261  h.SetLineColor(kGreen-3)
262  leg.AddEntry(Clone(h), '#bar{#nu}_{#mu} CC', 'l')
263  h.SetLineColor(kMagenta)
264  leg.AddEntry(Clone(h), '#bar{#nu}_{e} CC', 'l')
265  if decomp:
266  h.SetLineStyle(7)
267  h.SetLineColor(kBlack)
268  leg.AddEntry(Clone(h), 'Uncorr. MC', 'l')
269  leg.SetTextSize(1.2*leg.GetTextSize())
270  else:
271  leg.SetTextSize(0.8*h.GetXaxis().GetTitleSize())
272  leg.Draw()
273 
274  if isSplit:
275  RatioLegend(isDecomp)
276 
277 # Add plots where the legend would fit better on the
278 # left side
279 def IsLeftAlign(fname):
280  needsLeft = ['id_lid', 'shw_cos_z',
281  'shw_cos_numi', 'id_cvne',
282  'cosangle', 'cosdang',
283  'nplanestofront', 'shwhitfrac',
284  'shwMaxy', 'cvn2d'
285  'cvn_mode', 'lid2d']
286  for fig in needsLeft:
287  if fig.lower() in fname.lower():
288  return True
289 
290  return False
291 
292 def ColorLegend(isDecomp = False, isSplit = False, leftAlign = True):
293  # top legend for line/marker descriptors
294  if not leftAlign: topleg = New(TLegend, .4, .8, .9, .88)
295  else: topleg = New(TLegend, .11, .8, .61, .88)
296  topleg.SetNColumns(3)
297  topleg.SetFillStyle(0)
298 
299  # bottom legend for interaction colors
300  if not leftAlign: leg = New(TLegend, .7, .6, .875, .8)
301  else: leg = New(TLegend, .1, .6, .285, .8)
302  leg.SetNColumns(2)
303  leg.SetMargin(.22)
304  leg.SetFillStyle(0)
305  leg.SetTextAlign(33)
306  if leftAlign:
307  leg.SetTextAlign(13)
308 
309  # dummy histogram
310  h = New(TH1F, '', '', 1, 0, 1)
311 
312  h.SetLineStyle(7)
313  # if not decomp, we don't have uncorrected vs corrected
314  # spectra to distinguish between. Set text color white
315  # so the other entries stay right justified
316  if not leftAlign:
317  if isDecomp:
318  topleg.AddEntry(Clone(h), 'Uncorr. MC', 'l')
319  else:
320  h.SetLineColor(kWhite)
321  blank = topleg.AddEntry(Clone(h), ' ', 'l')
322  blank.SetTextColor(kWhite)
323 
324  h.SetLineStyle(1)
325  h.SetLineColor(kRed)
326  mc = topleg.AddEntry(Clone(h), 'Total MC', 'l')
327 
328  h.SetMarkerStyle(kFullCircle)
329  h.SetLineColor(kBlack)
330  d = topleg.AddEntry(Clone(h), 'ND data', 'lep')
331 
332  else:
333  h.SetMarkerStyle(kFullCircle)
334  h.SetLineColor(kBlack)
335  d = topleg.AddEntry(Clone(h), 'ND data', 'lep')
336 
337  h.SetLineStyle(1)
338  h.SetLineColor(kRed)
339  mc = topleg.AddEntry(Clone(h), 'Total MC', 'l')
340 
341  if isDecomp:
342  h.SetLineStyle(7)
343  h.SetLineColor(kBlack)
344  topleg.AddEntry(Clone(h), 'Uncorr. MC', 'l')
345  else:
346  h.SetLineColor(kWhite)
347  blank = topleg.AddEntry(Clone(h), ' ', 'l')
348  blank.SetTextColor(kWhite)
349 
350 
351  # build the bottom legend
352  h.SetLineColor(kWhite)
353  l2 = leg.AddEntry(Clone(h), '#nu_{#mu} CC', 'l')
354  l2.SetTextColor(kGreen+2)
355 
356  l4 = leg.AddEntry(Clone(h), '#bar{#nu}_{#mu} CC', 'l')
357  l4.SetTextColor(kGreen-3)
358 
359  l3 = leg.AddEntry(Clone(h), '#nu_{e} CC', 'l')
360  l3.SetTextColor(kMagenta-2)
361 
362  l5 = leg.AddEntry(Clone(h), '#bar{#nu}_{e} CC', 'l')
363  l5.SetTextColor(kMagenta)
364 
365  # pad to make NC right justified
366  if not leftAlign:
367  blank = leg.AddEntry(Clone(h), '', 'l')
368  blank.SetTextColor(kWhite)
369  l1 = leg.AddEntry(Clone(h), 'NC', 'l')
370  l1.SetTextColor(kAzure)
371  else:
372  l1 = leg.AddEntry(Clone(h), 'NC', 'l')
373  l1.SetTextColor(kAzure)
374  blank = leg.AddEntry(Clone(h), '', 'l')
375  blank.SetTextColor(kWhite)
376 
377 
378  topleg.SetTextSize(1.3*topleg.GetTextSize())
379  leg.Draw()
380  topleg.Draw()
381 
382  if isSplit:
383  RatioLegend(isDecomp)
384 
385 
386 def TopLegend(isDecomp = False, isSplit = False):
387  leg = New(TLegend, .15, .75, .85, .85)
388  h = New(TH1F, '', '', 1, 0, 1)
389 
390  leg.SetMargin(.22)
391 # leg.SetBorderSize(1)
392  leg.SetNColumns(8)
393 
394  h.SetMarkerStyle(kFullCircle)
395  leg.AddEntry(Clone(h), 'ND data', 'lep')
396  h.SetLineColor(kRed)
397  leg.AddEntry(Clone(h), 'Total MC', 'l')
398  h.SetLineColor(kAzure)
399  leg.AddEntry(Clone(h), 'NC', 'l')
400  h.SetLineColor(kGreen+2)
401  leg.AddEntry(Clone(h), '#nu_{#mu} CC', 'l')
402  h.SetLineColor(kMagenta-2)
403  leg.AddEntry(Clone(h), '#nu_{e} CC', 'l')
404  h.SetLineColor(kGreen-3)
405  leg.AddEntry(Clone(h), '#bar{#nu}_{#mu} CC', 'l')
406  h.SetLineColor(kMagenta)
407  leg.AddEntry(Clone(h), '#bar{#nu}_{e} CC', 'l')
408  if decomp:
409  h.SetLineStyle(7)
410  h.SetLineColor(kBlack)
411  leg.AddEntry(Clone(h), 'Uncorr. MC', 'l')
412  leg.SetTextSize(1.2*leg.GetTextSize())
413  else:
414  leg.SetTextSize(0.8*h.GetXaxis().GetTitleSize())
415  leg.Draw()
416 
417  if isSplit:
418  RatioLegend(isDecomp)
419 
420 
421 def RatioLegend(isDecomp = False):
422  x1 = .12
423  x2 = .32
424  if isDecomp:
425  x2 = .47
426  l = New(TLegend, x1, .13, x2, .18)
427  l.SetBorderSize(1)
428  l.SetNColumns(2)
429  h = New(TH1F, '', '', 1, 0, 1)
430  h.SetMarkerStyle(kFullCircle)
431  h.SetMarkerColor(1)
432  h.SetLineColor(1)
433  if isDecomp:
434  h.SetMarkerColor(17)
435  h.SetLineColor(17)
436  l.AddEntry(Clone(h), 'Uncorr. MC', 'lep')
437 
438  if isDecomp:
439  h.SetMarkerColor(1)
440  h.SetLineColor(1)
441  l.AddEntry(Clone(h), 'Decomp MC', 'lep')
442  l.Draw()
443 
444 
445 parser = argparse.ArgumentParser()
446 parser.add_argument('-b', '--batch', action = 'store_true',
447  help = 'batch mode, no graphics output')
448 parser.add_argument('beamMode', help = 'Beam mode: rhc or fhc')
449 parser.add_argument('dataFile', help = 'Validation histogram file for data')
450 parser.add_argument('mcFile', help = 'Validation histogram file for MC')
451 parser.add_argument('decompFile', help = 'Validation histogram file for MC Decomp')
452 opts = parser.parse_args()
453 
454 # Screws up the argument parsing if it's any earlier
455 from ROOT import *
456 
457 # Style the plots nicely
458 gApplication.ExecuteFile('$SRT_PUBLIC_CONTEXT/Utilities/rootlogon.C')
459 gROOT.ForceStyle()
460 
461 if opts.batch: gROOT.SetBatch(True)
462 beamMode = opts.beamMode
463 fdata = TFile(opts.dataFile)
464 fmc = TFile(opts.mcFile)
465 fdecomp = TFile(opts.decompFile)
466 
467 for k in fdata.GetListOfKeys():
468  print k.GetName()
469  isLeft = IsLeftAlign(k.GetName())
470 
471  New(TCanvas)
472  for decomp in [True, False]:
473  decompStr = '';
474  lineStyle = 1 # line style for uncorrected MC
475  if decomp:
476  decompStr = '_decomp'
477  lineStyle = 7
478 
479  # Draw plots with no ratios
480  data = Spectrum(k.ReadObj())
481  if data.hist.Integral() == 0: continue
482 
483  data.hist.SetMarkerStyle(kFullCircle)
484  scalePow = DeriveScalePow(data, data.pot)
485  data.Draw(data.pot, 'ep', scalePow)
486 
487  if decomp:
488  decomp = Prediction(fdecomp,k.GetName())
489  decomp.Draw(data.pot, 'hist same', scalePow)
490 
491  pred = Prediction(fmc, k.GetName())
492  pred.Draw(data.pot, 'hist same', scalePow, lStyle = lineStyle)
493 
494  data.Draw(data.pot, 'ep same', scalePow) # data should be on top
495  fname = re.sub(r'\{group=Cut,cat=._(.*)\}', r'_\1', k.GetName())
496  Preliminary()
497  BeamModeLabel(beamMode)
498  ColorLegend(decomp, leftAlign = isLeft)
499 
500  gPad.SetLogy(True)
501  gPad.Print('{}_plots/{}_{}{}_log.pdf'.format(beamMode, beamMode, fname, decompStr))
502  gPad.SetLogy(False)
503  gPad.Print('{}_plots/{}_{}{}.pdf'.format(beamMode, beamMode, fname, decompStr))
504 
505  # Draw only datamc ratios
506  New(TCanvas)
507  r = Ratio(Clone(data), Clone(pred.tot))
508  r.Draw('ep')
509  One()
510  BeamModeLabel(beamMode)
511  Preliminary()
512  RatioLegend(decomp)
513  if decomp:
514  rdecomp = Ratio(Clone(data), Clone(decomp.tot))
515  rdecomp.Draw('ep same')
516  r.Draw('ep same', color=17) # uncorrected should be on top
517  gPad.Print('{}_plots/{}_{}{}_ratio.pdf'.format(beamMode, beamMode, fname, decompStr))
518 
519  # Draw split canvas with, spectra and ratio
520  c = New(TCanvas)
521 
522  # De-conflict the y-axis labels in the split plot case
523  data.hist.GetYaxis().SetTitleSize(.9*data.hist.GetYaxis().GetTitleSize())
524  r.hist.GetYaxis().SetTitleSize(.9*r.hist.GetYaxis().GetTitleSize())
525 
526  (p1, p2) = SplitCanvas(.35)
527  p1.cd()
528  data.Draw(data.pot, 'ep', scalePow, True)
529  pred.Draw(data.pot, 'hist same', scalePow, lStyle = lineStyle)
530  if decomp:
531  decomp.Draw(data.pot, 'hist same', scalePow)
532  data.Draw(data.pot, 'ep same', scalePow)
533  ColorLegend(isDecomp = decomp, isSplit = True, leftAlign = isLeft)
534  p2.cd()
535  if decomp:
536  rdecomp.Draw('ep')
537  r.Draw('ep same', color = 17)
538  else:
539  r.Draw('ep')
540 
541  One()
542  Preliminary()
543  BeamModeLabel(beamMode)
544 
545  c.cd()
546  p1.SetLogy(True)
547  gPad.Print('{}_plots/{}_{}{}_split_log.pdf'.format(beamMode, beamMode, fname, decompStr))
548  p1.SetLogy(False)
549 
550  splitName ='{}_plots/{}_{}{}_split.pdf'.format(beamMode, beamMode, fname, decompStr)
551  gPad.Print(splitName)
552 
553  Caption(splitName, data.hist.GetXaxis().GetTitle(), decomp)
554 
555 
556  # Don't flood the screen with plots
557  if not gROOT.IsBatch(): break
558 
559 
560 if not opts.batch:
561  # Just want the console to hang around so the user can look at the plots
562  # without them disappearing. Ctrl-D or "exit()" to quit.
563  print 'Ctrl-D to quit'
564  code.interact(local = locals(), banner='')
double Integral(const Spectrum &s, const double pot, int cvnregion)
def Legend(isDecomp=False, isSplit=False)
constexpr T pow(T x)
Definition: pow.h:75
gargamelle SetTitle("Gargamelle #nu_{e} CC data")
def __init__(self, numSpect, denomSpect)
def TopLegend(isDecomp=False, isSplit=False)
std::string format(const int32_t &value, const int &ndigits=8)
Definition: HexUtils.cpp:14
procfile open("FD_BRL_v0.txt")
string GetString(xmlDocPtr xml_doc, string node_path)
def Caption(fName, axisTitle, decomp)
def Draw(self, pot, opt='', scalePow=0, split=False)
def ColorLegend(isDecomp=False, isSplit=False, leftAlign=True)
def LegendSplit(spec, decomp=False, ratio=False)
def Draw(self, pot, opt='', scalePow=0, split=False, lStyle=1)