demo_prong.py
Go to the documentation of this file.
1 import os
2 import sys
3 
4 import numpy as np
5 import matplotlib.pyplot as plt
6 import pandas as pd
7 
8 import PandAna.core as pa
9 from PandAna.core.core import KL, KLN, KLS
10 from PandAna.cut.analysis_cuts import *
11 from PandAna.var.analysis_vars import *
12 
13 # Demo script for analysing neutron prongs on ND RHC data, similar to Miranda's caf level analysis
14 
15 KLP = KL + ['rec.vtx.elastic.fuzzyk.png_idx']
16 
17 # demo slice cuts
18 def kNDRockFilter(tables):
19  df = tables['rec.vtx.elastic']
20  ret = (df['vtx.x'] > -180) & \
21  (df['vtx.x'] < 180) & \
22  (df['vtx.y'] > -180) & \
23  (df['vtx.y'] < 180) & \
24  (df['vtx.z'] > 20)
25  ret.name = 'rock'
26  return ret
27 kNDRockFilter = Cut(kNDRockFilter)
28 
29 # Prong Cuts that need slice indices
30 def kProngDispl(tables):
31  dfvtx = tables['rec.vtx.elastic'][['vtx.x', 'vtx.y', 'vtx.z']]
32  dfpng = tables['rec.vtx.elastic.fuzzyk.png'][['start.x', 'start.y', 'start.z', 'rec.vtx.elastic.fuzzyk.png_idx']]
33 
34  xdiff = dfvtx['vtx.x'] - dfpng['start.x']
35  ydiff = dfvtx['vtx.y'] - dfpng['start.y']
36  zdiff = dfvtx['vtx.z'] - dfpng['start.z']
37  length = np.sqrt(xdiff*xdiff + ydiff*ydiff + zdiff*zdiff)
38  length.name = 'disp'
39  idx = dfpng['rec.vtx.elastic.fuzzyk.png_idx']
40  ret = pd.concat([length, idx], axis=1)
41  ret = ret.reset_index().set_index(KLP)
42 
43  return (ret['disp'] > 20).where(ret['disp'] > 20)
44 
45 
46 # prong cuts that don't need slice indices
47 kProngNhitsCut = Cut(lambda tables: tables['rec.vtx.elastic.fuzzyk.png']['nhit'] < 6)
48 kPhotonID2018Cut = Cut(lambda tables: tables['rec.vtx.elastic.fuzzyk.png.cvnpart']['photonid'] < 0.6)
49 kIsNeutron = Cut(lambda tables: abs(tables['rec.vtx.elastic.fuzzyk.png.truth']['motherpdg']) == 2112)
50 
51 def kIgnoreMuon2018(tables):
52  dfpng = tables['rec.vtx.elastic.fuzzyk.png.cvnpart']['muonid']
53  dflen = tables['rec.vtx.elastic.fuzzyk.png']['len']
54  if not dfpng.empty:
55  dflen = dflen.reset_index().set_index(KL)
56  dfpng = dfpng.reset_index().set_index(KL)
57  maxid = dfpng['muonid'].groupby(level=KL).agg(np.max)
58  ret = pd.concat([dfpng['muonid'] - maxid, dfpng['rec.vtx.elastic.fuzzyk.png_idx'], dflen['len']], axis=1)
59  ret = ret.reset_index().set_index(KLP)
60  return ((ret['muonid'] < 0) & (ret['len'] <= 500))
61  else:
62  return dfpng
63 kIgnoreMuon2018 = Cut(kIgnoreMuon2018)
64 
65 # kIgnoreMuon2018 has to come first here before other cuts.
66 # Otherwise the maximum muon id won't be taken over all prongs in the slice but just on prongs that pass some cuts which would have removed muons
67 kProngCuts = kIgnoreMuon2018 & kIsNeutron & kProngNhitsCut & kPhotonID2018Cut
68 
69 # prong level var
70 kProtonID = Var(lambda tables: abs(tables['rec.vtx.elastic.fuzzyk.png.cvnpart']['protonid']))
71 
72 limit = 100
73 filelist = '/pnfs/nova/persistent/users/karlwarb/HDF5-Training-19-02-26/ND-ProngCVN-RHC/*.h5'
74 tablesProng = pa.loader(filelist,limit=limit, index=KLP)
75 tablesSlc = pa.loader(filelist, limit=limit)
76 # associate loads the same data for tablesProng and tablesSlc just once
77 a = pa.associate([tablesProng, tablesSlc])
78 
79 # get dataframe for kProngDispl first
80 specDummySlc = pa.spectrum(tablesSlc, kNDRockFilter & kNumuContainND, kProngDispl)
81 # get dataframe for prong level var
82 specProng = pa.spectrum(tablesProng, kProngCuts, kProtonID)
83 
84 # go, go, go
85 a.Go()
86 
87 prongdispldf = specDummySlc.df()
88 vardf = specProng.df()
89 print(prongdispldf.head())
90 print(vardf.head())
91 # this accomplishes two things at once
92 # prongdispldf has all prong displacements passing slice level cuts and displacement > 20
93 # doing pd.concat([prongdispldf, vardf], axis=1, join='inner') therefore only selects rows of vardf that exist in prongdispldf when we use join='inner'
94 finaldf = pd.concat([prongdispldf, vardf], axis=1, join='inner')[vardf.name]
95 
96 # create a spectrum object out of finaldf
97 finalspec = pa.filledSpectrum(finaldf, specProng.POT())
98 
99 # make the prong level var plot
100 n, bins = finalspec.histogram(bins=25, range=(0,1))
101 
102 print('Selected ' + str(n.sum()) + ' events from ' + str(finalspec.POT()) + ' POT.')
103 
104 plt.hist(bins[:-1], bins=bins, weights=n, histtype='step', color='blue', label='ND RHC')
105 plt.xlabel('Proton ID 2018')
106 plt.ylabel('Prongs')
107 
108 plt.legend(loc='best')
109 
110 plt.show()
def kProngDispl(tables)
Definition: demo_prong.py:30
void abs(TH1 *hist)
bool print