prod5_pid_validation.py
Go to the documentation of this file.
1 #!/usr/bin/env python
2 import h5py
3 import os
4 import sys
5 sys.path.append('../..')
6 
7 import numpy as np
8 np.seterr(invalid='ignore', divide='ignore') # ignore warnings about nan and divide by zero
9 import pandas as pd
10 
11 from PandAna.core.core import *
12 from PandAna.cut.analysis_cuts import *
13 
14 import matplotlib
15 matplotlib.use('Agg')
16 matplotlib.rcParams.update({'font.size': 12,
17  'patch.linewidth': 1})
18 
19 import matplotlib.pyplot as plt
20 from matplotlib.lines import Line2D
21 
22 KLP = KL + ['rec.vtx.elastic.fuzzyk.png_idx']
23 
24 def get_files(det, swap_str, horn_current):
25  d = '/lfstev/nnet/R19-02-23-miniprod5/'
26  if 'FD' in det:
27  d = d + 'FD-' + swap_str + "-" + horn_current + '-Eval/'
28  else:
29  d = d + 'ND-Mont-' + horn_current + '/'
30  return [os.path.join(d, f) for f in os.listdir(d) if 'h5caf.h5' in f]
31 
32 # workaround for applying slice cuts such as preselections to prong level spectra
33 def ApplySlcCutsPngDF(png_df, slc_df):
34  png_name = png_df.name
35  slc_name = slc_df.name
36 
37  # reset to slice level index
38  df = png_df.reset_index().set_index(KL)
39 
40  val = df[png_name] # get var column
41  val = val.groupby(level=KL).agg(list) # each row contains list of prong values
42 
43  # slc_df has indices that pass slice cuts. concat with prong column. join='inner' drops rows in prong column that don't pass slice cuts
44  # drop slice var column
45  val = pd.concat([val, slc_df], axis=1, join='inner').drop(slc_name, axis=1)
46 
47  # if no prongs pass prong and slice cuts, return an empty df
48  if val.empty:
49  val = pd.DataFrame([], columns=[png_name])
50 
51  else:
52  # unlist each row
53  val = val[png_name].apply(pd.Series).stack()
54  val = val.reset_index()
55  val.columns = KL + ['idx', png_name]
56  # set new index
57  val = val.set_index(KL + ['idx'])
58 
59  # contains prong dataframe that has slice cuts applied
60  return val[png_name]
61 
62 # selection cuts
63 ## FD check
64 def kIsFarDet(tables):
65  query = tables._files.query
66  if not type(query) is list: query = [query]
67  return 'fardet' in query[0]
68 
69 ############### nue cuts ######################
70 def kNueContain(tables):
71  if kIsFarDet(tables):
72  return kNueProngContainment(tables) & kNueBasicPart(tables)
73  else:
74  return kNueNDContain(tables)
75 
76 kNueContain = Cut(kNueContain)
77 
78 def kNuePresel(tables):
79  if kIsFarDet(tables):
80  return kNueCorePresel(tables)
81  else:
82  return kNueNDPresel(tables)
83 
84 kNuePresel = Cut(kNuePresel)
85 
86 ############### numu cuts #####################
87 # kCCE isn't working yet
88 kNumuNoPIDNoCCEFD = kNumuBasicQuality & kNumuContainFD
89 kNumuNoPIDNoCCEND = kNumuBasicQuality & kNumuContainND
90 def kNumuContain(tables):
91  if kIsFarDet(tables):
92  return kNumuContainFD(tables)
93  else:
94  return kNumuContainND(tables)
95 
96 kNumuContain = Cut(kNumuContain)
97 
98 def kNumuPresel(tables):
99  if kIsFarDet(tables):
100  return kNumuNoPIDNoCCEFD(tables)
101  else:
102  return kNumuNoPIDNoCCEND(tables)
103 
104 kNumuPresel = Cut(kNumuPresel)
105 
106 ############### nus cuts ######################
107 def kNusContain(tables):
108  if kIsFarDet(tables):
109  return kNusFDContain(tables)
110  else:
111  return kNusNDContain(tables)
112 kNusContain = Cut(kNusContain)
113 
114 def kNusPresel(tables):
115  if kIsFarDet(tables):
116  return kNusFDPresel(tables)
117  else:
118  return kNusNDPresel(tables)
119 kNusPresel = Cut(kNusPresel)
120 
121 ############### ORd cuts #####################
122 kCosVeto = kVeto
123 kOrContainment = kNumuContain | kNusContain | kNueContain
124 kOrPreselection = kNumuPresel | kNusPresel | kNuePresel
125 
126 # prong quality cuts
127 kHas2020Score = Cut(lambda tables: tables['rec.vtx.elastic.fuzzyk.png.cvnpart2020fhc']['muonid'] > 0)
128 kHas2018Score = Cut(lambda tables: tables['rec.vtx.elastic.fuzzyk.png.cvnpart']['muonid'] > 0)
129 kProngLength = Cut(lambda tables: tables['rec.vtx.elastic.fuzzyk.png']['len'] < 500)
130 
131 kProngCuts = kProngLength & kHas2020Score & kHas2018Score
132 
133 # particle truth cuts
134 kIsChargedPion = Cut(lambda tables: abs(tables['rec.vtx.elastic.fuzzyk.png.truth']['pdg']) == 211)
135 kIsPhoton = Cut(lambda tables: tables['rec.vtx.elastic.fuzzyk.png.truth']['pdg'] == 22)
136 kIsNeutron = Cut(lambda tables: abs(tables['rec.vtx.elastic.fuzzyk.png.truth']['pdg']) == 2112)
137 kIsProton = Cut(lambda tables: abs(tables['rec.vtx.elastic.fuzzyk.png.truth']['pdg']) == 2212)
138 kIsElectron = Cut(lambda tables: abs(tables['rec.vtx.elastic.fuzzyk.png.truth']['pdg']) == 11)
139 kIsMuon = Cut(lambda tables: abs(tables['rec.vtx.elastic.fuzzyk.png.truth']['pdg']) == 13)
140 
141 # PID vars prod5
142 def kProtonPIDProd5FHC(tables):
143  cvn_png_df = tables[('rec.vtx.elastic.fuzzyk.png.cvnpart2020fhc')]
144  return cvn_png_df['protonid']
145 def kProtonPIDProd5RHC(tables):
146  cvn_png_df = tables[('rec.vtx.elastic.fuzzyk.png.cvnpart2020rhc')]
147  return cvn_png_df['protonid']
148 
149 def kPhotonPIDProd5FHC(tables):
150  cvn_png_df = tables[('rec.vtx.elastic.fuzzyk.png.cvnpart2020fhc')]
151  return cvn_png_df['photonid']
152 def kPhotonPIDProd5RHC(tables):
153  cvn_png_df = tables[('rec.vtx.elastic.fuzzyk.png.cvnpart2020rhc')]
154  return cvn_png_df['photonid']
155 
156 def kPionPIDProd5FHC(tables):
157  cvn_png_df = tables[('rec.vtx.elastic.fuzzyk.png.cvnpart2020fhc')]
158  return cvn_png_df['pionid']
159 def kPionPIDProd5RHC(tables):
160  cvn_png_df = tables[('rec.vtx.elastic.fuzzyk.png.cvnpart2020rhc')]
161  return cvn_png_df['pionid']
162 
164  cvn_png_df = tables[('rec.vtx.elastic.fuzzyk.png.cvnpart2020fhc')]
165  return cvn_png_df['electronid']
167  cvn_png_df = tables[('rec.vtx.elastic.fuzzyk.png.cvnpart2020rhc')]
168  return cvn_png_df['electronid']
169 
170 def kMuonPIDProd5FHC(tables):
171  cvn_png_df = tables[('rec.vtx.elastic.fuzzyk.png.cvnpart2020fhc')]
172  return cvn_png_df['muonid']
173 def kMuonPIDProd5RHC(tables):
174  cvn_png_df = tables[('rec.vtx.elastic.fuzzyk.png.cvnpart2020rhc')]
175  return cvn_png_df['muonid']
176 
177 def kEMPIDProd5FHC(tables):
178  emid_df = kElectronPIDProd5FHC(tables) + kPhotonPIDProd5FHC(tables)
179  emid_df.name = 'emid' # name the series, we'll need it later
180  return emid_df
181 def kEMPIDProd5RHC(tables):
182  emid_df = kElectronPIDProd5RHC(tables) + kPhotonPIDProd5RHC(tables)
183  emid_df.name = 'emid' # name the series, we'll need it later
184  return emid_df
185 
186 
187 # PID vars prod4
188 def kProtonPID(tables):
189  cvn_png_df = tables[('rec.vtx.elastic.fuzzyk.png.cvnpart')]
190  return cvn_png_df['protonid']
191 
192 def kPhotonPID(tables):
193  cvn_png_df = tables[('rec.vtx.elastic.fuzzyk.png.cvnpart')]
194  return cvn_png_df['photonid']
195 
196 def kPionPID(tables):
197  cvn_png_df = tables[('rec.vtx.elastic.fuzzyk.png.cvnpart')]
198  return cvn_png_df['pionid']
199 
200 def kElectronPID(tables):
201  cvn_png_df = tables[('rec.vtx.elastic.fuzzyk.png.cvnpart')]
202  return cvn_png_df['electronid']
203 
204 def kMuonPID(tables):
205  cvn_png_df = tables[('rec.vtx.elastic.fuzzyk.png.cvnpart')]
206  return cvn_png_df['muonid']
207 
208 def kEMPID(tables):
209  emid_df = kElectronPID(tables) + kPhotonPID(tables)
210  emid_df.name = 'emid' # name the series, we'll need it later
211  return emid_df
212 
213 # fhc trained networks
214 kPhotonPIDProd5FHC = Var(kPhotonPIDProd5FHC)
215 kProtonPIDProd5FHC = Var(kProtonPIDProd5FHC)
216 kPionPIDProd5FHC = Var(kPionPIDProd5FHC)
217 kElectronPIDProd5FHC = Var(kElectronPIDProd5FHC)
218 kMuonPIDProd5FHC = Var(kMuonPIDProd5FHC)
219 kEMPIDProd5FHC = Var(kEMPIDProd5FHC)
220 # rhc trained networks
221 kPhotonPIDProd5RHC = Var(kPhotonPIDProd5RHC)
222 kProtonPIDProd5RHC = Var(kProtonPIDProd5RHC)
223 kPionPIDProd5RHC = Var(kPionPIDProd5RHC)
224 kElectronPIDProd5RHC = Var(kElectronPIDProd5RHC)
225 kMuonPIDProd5RHC = Var(kMuonPIDProd5RHC)
226 kEMPIDProd5RHC = Var(kEMPIDProd5RHC)
227 
228 # prod4
229 kPhotonPID = Var(kPhotonPID)
230 kProtonPID = Var(kProtonPID)
231 kPionPID = Var(kPionPID)
232 kElectronPID = Var(kElectronPID)
233 kMuonPID = Var(kMuonPID)
234 kEMPID = Var(kEMPID)
235 
236 
237 cut_levels = {'Veto' : kCosVeto,
238  'Containment' : kCosVeto & kOrContainment,
239  'Preselection' : kCosVeto & kOrContainment & kOrPreselection}
240 
241 datasets = {'FD_Nonswap_FHC' : 'prod_h5_R19-02-23-miniprod5.n_fd_genie_default_nonswap_fhc_nova_v08_full_v1_eval',
242  'FD_Fluxswap_FHC' : 'prod_h5_R19-02-23-miniprod5.n_fd_genie_default_fluxswap_fhc_nova_v08_full_v1_eval',
243  'ND_Nonswap_FHC' : 'junk',
244  'ND_Nonswap_RHC' : 'junk'}
245 
246 pids = ['photonid', 'protonid', 'pionid', 'electronid', 'muonid', 'emid']
247 
248 
249 pid_scores_fhc = {'photonid_prod5fhc' : kPhotonPIDProd5FHC,
250  'protonid_prod5fhc' : kProtonPIDProd5FHC,
251  'pionid_prod5fhc' : kPionPIDProd5FHC,
252  'electronid_prod5fhc' : kElectronPIDProd5FHC,
253  'muonid_prod5fhc' : kMuonPIDProd5FHC}
254 pid_scores_rhc = {'photonid_prod5rhc' : kPhotonPIDProd5RHC,
255  'protonid_prod5rhc' : kProtonPIDProd5RHC,
256  'pionid_prod5rhc' : kPionPIDProd5RHC,
257  'electronid_prod5rhc' : kElectronPIDProd5RHC,
258  'muonid_prod5rhc' : kMuonPIDProd5RHC}
259 
260 pid_scores = {'photonid_prod4' : kPhotonPID,
261  'protonid_prod4' : kProtonPID,
262  'pionid_prod4' : kPionPID,
263  'electronid_prod4' : kElectronPID,
264  'muonid_prod4' : kMuonPID,
265  'emid_prod4' : kEMPID,
266  'photonid_prod5fhc' : kPhotonPIDProd5FHC,
267  'protonid_prod5fhc' : kProtonPIDProd5FHC,
268  'pionid_prod5fhc' : kPionPIDProd5FHC,
269  'electronid_prod5fhc' : kElectronPIDProd5FHC,
270  'muonid_prod5fhc' : kMuonPIDProd5FHC,
271  'emid_prod5fhc' : kEMPIDProd5FHC,
272  'photonid_prod5rhc' : kPhotonPIDProd5RHC,
273  'protonid_prod5rhc' : kProtonPIDProd5RHC,
274  'pionid_prod5rhc' : kPionPIDProd5RHC,
275  'electronid_prod5rhc' : kElectronPIDProd5RHC,
276  'muonid_prod5rhc' : kMuonPIDProd5RHC,
277  'emid_prod5rhc' : kEMPIDProd5RHC}
278 
279 particle_cuts = {'true_pion' : kIsChargedPion,
280  'true_photon' : kIsPhoton,
281  'true_proton' : kIsProton,
282  'true_electron' : kIsElectron,
283  'true_muon' : kIsMuon}
284 
285 particle_colors = {'true_pion' : 'darkgreen',
286  'true_pizero' : 'chartreuse',
287  'true_photon' : 'crimson',
288  'true_neutron' : 'indigo',
289  'true_proton' : 'orangered',
290  'true_electron' : 'blueviolet',
291  'true_muon' : 'cyan',
292  'true_em' : 'green'}
293 
294 def main(wc, limit, stride, offset, spectra_file, make_plots, output, caf):
295 
296  if spectra_file is None:
297  loaders = {}
298  slc_tables = {} # loader for slice level cuts only
299  prong_tables = {} # loader for prong level cuts only
300  for dataset_name, data in datasets.items():
301  if wc:
302  det, swap, horn = dataset_name.split('_')
303  slc_tables[dataset_name] = loader(get_files(det, swap, horn),
304  stride=stride,
305  limit=limit,
306  offset=offset)
307  prong_tables[dataset_name] = loader(get_files(det, swap, horn),
308  stride=stride,
309  limit=limit,
310  offset=offset, index=KLP)
311  loaders[dataset_name] = associate([slc_tables[dataset_name],
312  prong_tables[dataset_name]]) # read same data only once for both loaders
313  else:
314  query = 'defname: {}'.format(data)
315  if limit is not None or stride is not None:
316  query += ' with'
317  if limit is not None:
318  query += ' limit {}'.format(limit)
319  if stride is not None:
320  query += ' stride {}'.format(stride)
321  if offset is not None:
322  query += ' offset {}'.format(offset)
323  print(query)
324  slc_tables[dataset_name] = loader(query)
325  prong_tables[dataset_name] = loader(query, index=KLP)
326  loaders[dataset_name] = associate([slc_tables[dataset_name], prong_tables[dataset_name]])
327 
328 
329  specs = {}
330  slc_specs = {}
331  save_specs = []
332  save_labels = []
333  for cut_name, cut in cut_levels.items():
334  slc_specs[cut_name] = {}
335  specs[cut_name] = {}
336  for dataset_name, data in datasets.items():
337  slc_specs[cut_name][dataset_name] \
338  = spectrum(slc_tables[dataset_name], cut, kCaloE) # use dummy var for slice spectra with slice cuts
339  specs[cut_name][dataset_name] = {}
340  for particle_name, particle_cut in particle_cuts.items():
341  specs[cut_name][dataset_name][particle_name] = {}
342  for var_name, var in pid_scores.items():
343  specs[cut_name][dataset_name][particle_name][var_name] \
344  = spectrum(prong_tables[dataset_name], kProngCuts & particle_cut, var) # apply only particle truth cuts to prong level var
345 
346  for loader_name, load in loaders.items():
347  load.Go() # go, go, go
348 
349  filename = output + '/prod5_pid_validation_spectra'
350  if stride:
351  filename += '_s{}'.format(stride)
352  if limit:
353  filename += '_l{}'.format(limit)
354  if offset:
355  filename += '_o{}'.format(offset)
356  filename += '.hdf5'
357  for cut_name, cut in cut_levels.items():
358  for dataset_name, data in datasets.items():
359  for particle_name, particle_cut in particle_cuts.items():
360  for var_name in pid_scores:
361  # get prong dataframe
362  df_prong = specs[cut_name][dataset_name][particle_name][var_name].df()
363  df_weight = specs[cut_name][dataset_name][particle_name][var_name]._weight
364  # apply slice cuts to prong dataframe
365  df_prong = ApplySlcCutsPngDF(df_prong, slc_specs[cut_name][dataset_name].df())
366  df_weight = ApplySlcCutsPngDF(df_weight, slc_specs[cut_name][dataset_name].df())
367  # reset prong spectrum with new dataframe
368  specs[cut_name][dataset_name][particle_name][var_name]._df = df_prong
369  specs[cut_name][dataset_name][particle_name][var_name]._weight = df_weight
370  # save
371  save_specs.append(specs[cut_name][dataset_name][particle_name][var_name])
372  save_labels.append('{}_{}_{}_{}'.format(cut_name, particle_name, var_name, dataset_name))
373 
374 
375  save_spectra(filename,
376  save_specs,
377  save_labels)
378 
379  else:
380  print('Loading spectra from {}'.format(spectra_file))
381  if caf:
382  print('Loading caf results from {}'.format(caf))
383  caf_file = h5py.File(caf, 'r')
384 
385  specs = {}
386  pid_score_names = list(pid_scores.keys())
387  for cut_name, cut in cut_levels.items():
388  specs[cut_name] = {}
389  for particle_name, particle_cut in particle_cuts.items():
390  specs[cut_name][particle_name] = {}
391  for dataset_name, data in datasets.items():
392  specs[cut_name][particle_name][dataset_name] = {}
393  for var_name in pid_scores:
394  if caf and 'prod4' in var_name:
395  specs[cut_name][particle_name][dataset_name][var_name+'caf'] = {}
396  spec_name = '{}_{}_{}_{}'.format(cut_name, particle_name, var_name, dataset_name)
397  specs[cut_name][particle_name][dataset_name][var_name] = load_spectra(spectra_file,
398  spec_name)
399  if caf and 'prod4' in var_name:
400  specs[cut_name][particle_name][dataset_name][var_name + 'caf'] \
401  = caf_spectra(caf_file,
402  '{}_{}_{}_{}'.format(cut_name,
403  dataset_name,
404  var_name.split('_')[0],
405  particle_name))
406  if caf: caf_file.close()
407 
408 
409  if make_plots or spectra_file:
410  particle_id = ['photonid',
411  'pionid',
412  'protonid',
413  'electronid',
414  'muonid',
415  'emid']
416 
417  network_linestyles = {'prod4': '-',
418  'prod5fhc': '--',
419  'prod5rhc': '-.',
420  'prod4caf': '-'}
421 
422  for cut_name, cut in cut_levels.items():
423  for dataset_name, data in datasets.items():
424  for pid in particle_id:
425  prod4pid = pid + '_prod4'
426  prod5pidfhc = pid + '_prod5fhc'
427  prod5pidrhc = pid + '_prod5rhc'
428  prod4pid_caf = pid + '_prod4caf'
429  networks = [prod4pid, prod5pidfhc, prod5pidrhc]
430  if caf:
431  networks = networks + [prod4pid_caf]
432  for particle_name, particle_cut in particle_cuts.items():
433  plot_spectra = []
434  plot_spectra_styles = []
435  for network in networks:
436  plot_spectra_styles.append({})
437  if 'caf' in network: plot_spectra.append(specs[cut_name][particle_name][dataset_name][network])
438  else: plot_spectra.append(copy_spectrum(specs[cut_name][particle_name][dataset_name][network]))
439 
440  plot_spectra_styles[-1]['color'] = particle_colors[particle_name]
441  plot_spectra_styles[-1]['label'] = network.split('_')[-1]
442  plot_spectra_styles[-1]['linestyle'] = network_linestyles[network.split('_')[-1]]
443  if 'caf' in network:
444  plot_spectra_styles[-1]['color'] = 'springgreen'
445  plot_spectra_styles[-1]['yerr'] = True
446  else:
447  # only plot error bars for the caf hist
448  # try overriding with None
449  plot_spectra_styles[-1]['yerr'] = False
450 
451  plot_ratio(plot_spectra,
452  30,
453  pid,
454  'Prongs',
455  'pid_plots_technote/{}_{}_{}_{}'.format(cut_name, dataset_name, pid, particle_name),
456  logy=True,
457  spec_styles=plot_spectra_styles)
458 
460  def __init__(self, input_file, group):
461  self.contents = input_file[group + '/contents'][()]
462  self.edges = input_file[group + '/edges'][()]
463  self.pot = input_file[group + '/pot'][()]
464 
465  def histogram(self, bins=None, range=None, POT=None):
466  if not POT: POT = self.pot
467  return (self.contents * POT / self.pot, self.edges)
468 
469  def integral(self, POT=None):
470  if not POT: POT = self.pot
471  return self.contents.sum() * POT / self.pot
472 
473 def plot(spec, nbins, xlabel, ylabel, name):
474  fig, ax = plt.subplots()
475  n, bins = spec.histogram(bins=nbins, range=(0, 1))
476  ax.hist(bins[:-1], bins, weights=n, histtype='step')
477  ax.set_xlabel(xlabel)
478  ax.set_ylabel(ylabel)
479  plt.legend(loc='upper right')
480  fig.savefig(name + '.png')
481  plt.close()
482 
483 def event_count_report(label, specs, networks):
484  integrals = [spec.integral() for spec in specs]
485  nominal = specs[0].integral()
486  print('{} | '.format(label), end='')
487  for spec, label, integral in zip(specs, networks, integrals):
488  perdiff = (nominal - integral) / nominal
489  print(' {}: {} ({:.2f}) |'.format(label, integral, perdiff), end='')
490  print('')
491 
492 def plot_ratio(specs, nbins, xlabel, ylabel, plot_name,
493  logy = True, spec_styles = None, additional_handles=None):
494  # assume the first entry is the one to take ratios wrt to
495  if not type(specs) is list: specs = [specs]
496  if spec_styles is not None:
497  if not type(spec_styles) is list: spec_styles = [spec_styles]
498  else:
499  spec_styles = [{} for _ in range(len(specs))]
500  assert len(specs) == len(spec_styles), 'Each spectrum must have a style if any are provided'
501 
502  pot = specs[0].POT()
503  fig, ax = plt.subplots(2, sharex=True, gridspec_kw = {'hspace': 0, 'height_ratios': [3, 1]})
504 
505  one_x = np.linspace(0, 1)
506  one_y = [1 for _ in range(len(one_x))]
507  ax[1].plot(one_x, one_y, color='k', linestyle='dashed')
508  for spec, spec_style in zip(specs, spec_styles):
509  plot_ratio_err = spec_style['yerr']
510  spec_style.pop('yerr')
511  n, bins = spec.histogram(bins=nbins, range=(0, 1), POT=pot)
512  nratio, binsratio = ratio(spec, specs[0], nbins, pot)
513 
514  try:
515  # matplotlib doesn't like plotting histograms with all weights = 0
516  ax[0].hist(bins[:-1], bins, weights=n, histtype='step',
517  **spec_style)
518  ax[1].hist(binsratio[:-1], binsratio, weights=nratio, histtype='step',
519  **spec_style)
520  except:
521  print('{} is empty. skipping...'.format(plot_name))
522 
523  if plot_ratio_err:
524  bin_centers = (bins[:-1] + bins[1:])/2
525  nominal, nominal_bins = specs[0].histogram(bins=nbins, range=(0, 1))
526  err = np.sqrt(1/n + 1/nominal) * nominal / n
527  ax[1].errorbar(bin_centers, nratio, yerr=err, ecolor='k', fmt='k.')
528 
529 
530 
531  ax[1].set_xlabel(xlabel)
532  ax[1].set_ylabel('Ratio to prod4')
533  ax[1].set_xlim([0, 1])
534  ax[0].set_xlim([0, 1])
535  ax[0].set_ylabel(ylabel)
536  ax[0].set_yscale('log')
537 
538  # configure the legend on the main figure
539  handles, labels = ax[0].get_legend_handles_labels()
540  line_handles = [Line2D([], [], c=h.get_edgecolor(), ls=h.get_linestyle()) for h in handles]
541  if additional_handles is not None:
542  line_handles = line_handles + additional_handles['handles']
543  labels = labels + additional_handles['labels']
544  ax[0].legend(handles=line_handles, labels=labels, loc='upper center')
545 
546  plt.tight_layout()
547  plt.subplots_adjust(top=0.93)
548  fig.suptitle(plot_name.split('/')[-1].replace('_', ' '))
549  fig.savefig(plot_name + '.png')
550  print('Created {}.png'.format(plot_name))
551  plt.close()
552 
553 
554 def ratio(spec1, spec2, nbins, pot, binrange=(0, 1)):
555  h1, bins1 = spec1.histogram(bins=nbins, range=binrange, POT=pot)
556  h2, bins2 = spec2.histogram(bins=nbins, range=binrange, POT=pot)
557  ratio = np.nan_to_num(h1) / np.nan_to_num(h2)
558  return np.nan_to_num(ratio), bins1
559 
560 if __name__ == '__main__':
561  import argparse
562 
563  parser = argparse.ArgumentParser('Plot prod5 and prod4 PID distributions')
564  parser.add_argument('--limit', '-l', default=None, type=int,
565  help='Limit the number of files to process')
566  parser.add_argument('--stride', '-s', default=1, type=int,
567  help='Set stride')
568  parser.add_argument('--offset', '-o', default=0, type=int,
569  help='Set file offset')
570  parser.add_argument('--spectra_file', '-f', default=None,
571  help='Only plot spectra previously filled by this application')
572  parser.add_argument('--wc', action='store_true',
573  help='Flag for running on WC')
574  parser.add_argument('--make_plots', action='store_true',
575  help='Make plots. Always true if providing a spectra file')
576  parser.add_argument('--output', default='.',
577  help='Output directory')
578  parser.add_argument('--caf', default=None,
579  help='Results of the Prod4 network from analysis cafs')
580  args = parser.parse_args()
581 
582  main(args.wc, args.limit, args.stride, args.offset,
583  args.spectra_file, args.make_plots, args.output,
584  args.caf)
585 
def apply(command)
Definition: g4zmq.py:38
def main(wc, limit, stride, offset, spectra_file, make_plots, output, caf)
const Cut kNusNDPresel
Definition: NusCuts.h:69
void abs(TH1 *hist)
def get_files(det, swap_str, horn_current)
def plot(spec, nbins, xlabel, ylabel, name)
const Cut kNumuContainFD([](const caf::SRProxy *sr){ std::pair< int, int > planes=calcFirstLastLivePlane(sr->slc.firstplane, std::bitset< 14 >(sr->hdr.dibmask));int planestofront=sr->slc.firstplane-planes.first;int planestoback=planes.second-sr->slc.lastplane;return( sr->slc.ncellsfromedge > 1 &&planestofront > 1 &&planestoback > 1 &&sr->sel.contain.kalfwdcell > 10 &&sr->sel.contain.kalbakcell > 10 &&sr->sel.contain.cosfwdcell > 0 &&sr->sel.contain.cosbakcell > 0);})
Definition: NumuCuts.h:20
const Cut kNusFDContain([](const caf::SRProxy *sr){const caf::SRNueCosRejProxy &cr=sr->sel.nuecosrej;if(std::min(cr.starteast, cr.stopeast)< 10) return false;if(std::min(cr.startwest, cr.stopwest)< 10) return false;if(std::min(cr.starttop, cr.stoptop) < 10) return false;if(std::min(cr.startbottom, cr.stopbottom)< 10) return false;if(std::min(cr.startfront, cr.stopfront)< 10) return false;if(std::min(cr.startback, cr.stopback)< 10) return false;return true;})
Containment variable for NC events from docdb 14241.
Definition: NusCuts.h:24
def ratio(spec1, spec2, nbins, pot, binrange=(0, 1))
def load_spectra(fname, groups)
Definition: core.py:164
const Cut kNumuContainND([](const caf::SRProxy *sr){return( sr->trk.kalman.ntracks > sr->trk.kalman.idxremid &&sr->slc.ncellsfromedge > 1 &&sr->slc.firstplane > 1 &&sr->slc.lastplane< 212 &&sr->trk.kalman.tracks[0].start.Z()< 1150 &&( sr->trk.kalman.tracks[0].stop.Z()< 1275 ||sr->sel.contain.kalyposattrans< 55) &&( sr->energy.numu.ndhadcalcatE +sr->energy.numu.ndhadcaltranE)< 0.03 &&sr->sel.contain.kalfwdcellnd > 4 &&sr->sel.contain.kalbakcellnd > 8);})
Definition: NumuCuts.h:22
kNumuNoPIDNoCCEFD
numu cuts ##################### kCCE isn&#39;t working yet
const ana::Cut kNueBasicPart
Definition: prod4_pid.C:60
def histogram(self, bins=None, range=None, POT=None)
def plot_ratio(specs, nbins, xlabel, ylabel, plot_name, logy=True, spec_styles=None, additional_handles=None)
bool print
std::string format(const int32_t &value, const int &ndigits=8)
Definition: HexUtils.cpp:14
std::vector< double > POT
const Cut kNusNDContain([](const caf::SRProxy *sr){const caf::SRNueCosRejProxy &cr=sr->sel.nuecosrej;if(std::min(cr.starteast, cr.stopeast)< 25) return false;if(std::min(cr.startwest, cr.stopwest)< 25) return false;if(std::min(cr.starttop, cr.stoptop) < 25) return false;if(std::min(cr.startbottom, cr.stopbottom)< 25) return false;if(std::min(cr.startfront, cr.stopfront)< 25) return false;if(std::min(cr.startback, cr.stopback)< 25) return false;return true;})
Containment variable for NC events from docdb 15242.
Definition: NusCuts.h:58
const Cut kNusFDPresel
Definition: NusCuts.h:38
def save_spectra(fname, spectra, groups)
Definition: core.py:116