tute_pid_validation.py
Go to the documentation of this file.
1 from PandAna import *
2 from PandAna.var.analysis_vars import *
3 from PandAna.cut.analysis_cuts import *
4 from PandAna.utils.enums import *
5 
6 import os
7 import sys
8 import numpy as np
9 import pandas as pd
10 import argparse
11 import re
12 
13 # your definition to run over
14 definition = 'def_snapshot karlwarb-MLTute2019-HDF5_R19-11-18-prod5reco.n_fluxswap'
15 
16 # define PID vars as two column dataframes for two CVNs
17 var = ['rec.sel.cvnloosepreselptp', 'rec.sel.cvnoldpresel']
18 def kNewCVNe(tables):
19  df = []
20  for v in var:
21  df.append(tables[v]['nueid'])
22  ret = pd.concat(df, axis=1)
23  ret.columns = [re.sub(r'rec.sel.', '', v) for v in var]
24  return ret
25 kNewCVNe = Var(kNewCVNe)
26 
27 def kNewCVNm(tables):
28  df = []
29  for v in var:
30  df.append(tables[v]['numuid'])
31  ret = pd.concat(df, axis=1)
32  ret.columns = [re.sub(r'rec.sel.', '', v) for v in var]
33  return ret
34 kNewCVNm = Var(kNewCVNm)
35 
36 def kNewCVNnc(tables):
37  df = []
38  for v in var:
39  df.append(tables[v]['ncid'])
40  ret = pd.concat(df, axis=1)
41  ret.columns = [re.sub(r'rec.sel.', '', v) for v in var]
42  return ret
43 kNewCVNnc = Var(kNewCVNnc)
44 
45 cvn_vars = {}
46 cvn_vars['nueid'] = kNewCVNe
47 cvn_vars['numuid'] = kNewCVNm
48 cvn_vars['ncid'] = kNewCVNnc
49 
50 # define a few more truth vars
51 kInelasticity = Var(lambda tables: tables['rec.mc.nu']['y'])
52 kTrueE = Var(lambda tables: tables['rec.mc.nu']['E'])
53 
54 # dumb oscillation weights to apply to spectra
55 def kWoscDumb(tables, weight):
56  sel = (tables['rec.mc']['nnu'] == 1)
57  weight[sel] *= tables['rec.mc.nu']['woscdumb']
58  return weight
59 
60 # now define some truth cuts
61 kPDG = Var(lambda tables: tables['rec.mc.nu']['pdg'])
62 kPDGAbs = Var(lambda tables: abs(tables['rec.mc.nu']['pdg']))
63 kMode = Var(lambda tables: tables['rec.mc.nu']['mode'])
64 kIsCC = Cut(lambda tables: tables['rec.mc.nu']['iscc'] == 1)
65 kIsNC = Cut(lambda tables: tables['rec.mc.nu']['iscc'] == 0)
66 
67 # function to return signal and background cuts depending on which analysis PID we're looking at
68 def SigCuts(idx='nueid'):
69  cut = Cut(lambda tables: tables['rec.mc']['nnu'] == 1)
70  if 'nue' in idx:
71  cut = cut & (kPDGAbs == 12) & kIsCC
72  if 'numu' in idx:
73  cut = cut & (kPDGAbs == 14) & kIsCC
74  if 'nc' in idx:
75  cut = cut & (~kIsCC)
76  return cut
77 
78 def BkgCuts(idx='nueid'):
79  cut = Cut(lambda tables: tables['rec.mc']['nnu'] == 1)
80  if 'nue' in idx:
81  cut = cut & ((kPDGAbs != 12) | ~kIsCC)
82  if 'numu' in idx:
83  cut = cut & ((kPDGAbs != 14) | ~kIsCC)
84  if 'nc' in idx:
85  cut = cut & (~kIsNC)
86  return cut
87 
88 # some basic preselections for each analysis
89 kNueFDSel = kNueProngContainment & kVeto
90 
91 kNumuFDSel = kNumuBasicQuality & kNumuContainFD
92 
93 kNusFDSel = kNusFDContain & kVeto
94 
95 ids = ['nueid', 'numuid', 'ncid']
96 
97 sel_cuts = {}
98 sel_cuts['nueid'] = kNueFDSel
99 sel_cuts['numuid'] = kNumuFDSel
100 sel_cuts['ncid'] = kNusFDSel
101 
102 truth_cuts = {}
103 for idx in ids:
104  truth_cuts[idx] = {}
105  truth_cuts[idx]['sig'] = SigCuts(idx)
106  truth_cuts[idx]['bkg'] = BkgCuts(idx)
107 
108 # also split sample by mode
109 def ModeCuts(mode=-1):
110  cut = Cut(lambda tables: tables['rec.mc']['nnu'] == 1)
111  if mode >= 0:
112  cut = cut & (kMode == mode)
113  return cut
114 
115 mode_cuts = {}
116 mode_cuts['kAll'] = ModeCuts(-1)
117 mode_cuts['kQE'] = ModeCuts(mode.kQE)
118 mode_cuts['kRes'] = ModeCuts(mode.kRes)
119 mode_cuts['kDIS'] = ModeCuts(mode.kDIS)
120 
121 # some helpful arguments for definition
122 parser = argparse.ArgumentParser(
123  description='Save CVN distributions for various truth cuts')
124 parser.add_argument(
125  '--limit', help=('Limit of files to run: '),
126  type=int, default=None)
127 parser.add_argument(
128  '--stride', help=('stride of files to run: '),
129  type=int, default=1)
130 parser.add_argument(
131  '--offset', help=('offset of files to run: '),
132  type=int, default=0)
133 parser.add_argument(
134  '--plot', help=('Run plotting part of the code '),
135  action='store_true', default=False)
136 
137 args = parser.parse_args()
138 
139 filename = 'cvn_truth_dist.h5'
140 # initialize loader
141 tables = loader(definition, limit=args.limit, stride=args.stride, offset=args.offset)
142 
143 # create PID, trueE and trueY spectra for each cut above and apply dumb oscillations
144 specs = {}
145 for selid, sel_cut in sel_cuts.items():
146  for ptype, pcut in truth_cuts[selid].items():
147  for mtype, mcut in mode_cuts.items():
148  cut = sel_cut & pcut & mcut
149  specid = 'cvn_%s_%s_%s' % (selid, ptype, mtype)
150  specs[specid] = spectrum(tables, cut, cvn_vars[selid], weight=kWoscDumb)
151  trueEid = 'trueE_%s_%s_%s' % (selid, ptype, mtype)
152  specs[trueEid] = spectrum(tables, cut, kTrueE, weight=kWoscDumb)
153  trueYid = 'trueY_%s_%s_%s' % (selid, ptype, mtype)
154  specs[trueYid] = spectrum(tables, cut, kInelasticity, weight=kWoscDumb)
155 
156 if not args.plot:
157  # Run!
158  tables.Go()
159  # save spectra to HDF5 output
160  save_tree(filename, specs.values(), specs.keys())
161 
162 # plot
163 else:
164  from tute_plots import *
165  # load from file
166  specs = load_tree(filename, specs.keys())
167  combined = {}
168  for selid, sel_cut in sel_cuts.items():
169  for ptype, pcut in truth_cuts[selid].items():
170  for mtype, mcut in mode_cuts.items():
171  cid = '%s_%s_%s' % (selid, ptype, mtype)
172  cvnid = 'cvn_%s_%s_%s' % (selid, ptype, mtype)
173  trueEid = 'trueE_%s_%s_%s' % (selid, ptype, mtype)
174  trueYid = 'trueY_%s_%s_%s' % (selid, ptype, mtype)
175 
176  # combine dataframes
177  combined[cid] = specs[cvnid] & specs[trueEid] & specs[trueYid]
178 
179  # Make some plots
180  MakeCVNDistPlot(combined['nueid_sig_kAll'], combined['nueid_bkg_kAll'], title='nueid', folder='plots')
181  MakeDeltaCVNPlots(combined['nueid_sig_kAll'], title='nueid', folder='plots')
182  MakeDeltaCVNPlots(combined['nueid_sig_kAll'], title='nueid', slc={'E':[0, 1, 2, 3, 4, 5]},folder='plots')
183  MakeDeltaCVNPlots(combined['nueid_sig_kAll'], title='nueid', slc={'y':[0., 0.2, 0.4, 0.6, 0.8]},folder='plots')
def MakeCVNDistPlot(sig, bkg, title='test', folder='.')
Definition: tute_plots.py:86
void abs(TH1 *hist)
def save_tree(fname, spectra, groups, attrs=True)
Definition: core.py:133
def MakeDeltaCVNPlots(sig, title='test', color='blue', slc={}, folder='.')
Definition: tute_plots.py:15
def load_tree(fname, groups, attrs=True)
Definition: core.py:183