analysis_vars.py
Go to the documentation of this file.
1 import numpy as np
2 import pandas as pd
3 from PandAna.core.core import Var, KL
4 from PandAna.utils import *
5 
6 from PandAna.var.numuE_utils import *
7 
8 ###################################################################################
9 #
10 # General Vars
11 #
12 ###################################################################################
13 
14 # cvnProd3Train is messed up in the current iteration of files.
15 # NOTE (from mbaird): Would really like to rename these 'kCVNe17'...
16 kCVNe = Var(lambda tables: tables['rec.sel.cvn2017']['nueid'])
17 kCVNm = Var(lambda tables: tables['rec.sel.cvn2017']['numuid'])
18 kCVNnc = Var(lambda tables: tables['rec.sel.cvn2017']['ncid'])
19 
20 kNHit = Var(lambda tables: tables['rec.slc']['nhit'])
21 
22 kRHC = Var(lambda tables: tables['rec.spill']['isRHC'])
23 kDetID = Var(lambda tables: tables['rec.hdr']['det'])
24 
25 kMeantime = Var(lambda tables: tables['rec.slc']['meantime'])
26 
27 # used in the 2020 analysis
28 kCVNe20 = Var(lambda tables: tables['rec.sel.cvnloosepreselptp']['nueid'])
29 kCVNm20 = Var(lambda tables: tables['rec.sel.cvnloosepreselptp']['numuid'])
30 kCVNnc20 = Var(lambda tables: tables['rec.sel.cvnloosepreselptp']['ncid'])
31 kCVNcos20 = Var(lambda tables: tables['rec.sel.cvnloosepreselptp']['cosmicid'])
32 
33 
34 
35 ###################################################################################
36 #
37 # Nue Vars
38 #
39 ###################################################################################
40 
41 def kLongestProng(tables):
42  df = tables['rec.vtx.elastic.fuzzyk.png']['len']
43  return df.groupby(level=KL).agg(np.max)
44 kLongestProng = Var(kLongestProng)
45 
46 kDistAllTop = Var(lambda tables: tables['rec.sel.nuecosrej']['distallpngtop'])
47 kDistAllBottom = Var(lambda tables: tables['rec.sel.nuecosrej']['distallpngbottom'])
48 kDistAllWest = Var(lambda tables: tables['rec.sel.nuecosrej']['distallpngwest'])
49 kDistAllEast = Var(lambda tables: tables['rec.sel.nuecosrej']['distallpngeast'])
50 kDistAllBack = Var(lambda tables: tables['rec.sel.nuecosrej']['distallpngback'])
51 kDistAllFront = Var(lambda tables: tables['rec.sel.nuecosrej']['distallpngfront'])
52 
53 kHitsPerPlane = Var(lambda tables: tables['rec.sel.nuecosrej']['hitsperplane'])
54 
55 kPtP = Var(lambda tables: tables['rec.sel.nuecosrej']['partptp'])
56 
57 def kMaxY(tables):
58  df = tables['rec.vtx.elastic.fuzzyk.png.shwlid']
59  df = df[['start.y','stop.y']].max(axis=1)
60  return df.groupby(level=KL).agg(np.max)
61 kMaxY = Var(kMaxY)
62 
63 kSparsenessAsymm = Var(lambda tables: tables['rec.sel.nuecosrej']['sparsenessasymm'])
64 
65 kCaloE = Var(lambda tables: tables['rec.slc']['calE'])
66 
67 def kNueCalibrationCorr(det, ismc, run):
68  if (det != detector.kFD): return 1.
69  if not ismc: return 0.9949
70  if run < 20753: return 0.9949/0.9844
71  return 1.
72 
73 def kEMEnergy(tables):
74  lng_png = kLongestProng(tables)
75  isRHC = kRHC(tables)
76 
77  shwlid_df = tables['rec.vtx.elastic.fuzzyk.png.shwlid']
78  prim_png = shwlid_df['calE'][(shwlid_df['rec.vtx.elastic.fuzzyk.png_idx']==0)]
79 
80 
81  png_df = tables['rec.vtx.elastic.fuzzyk.png']
82  cvn_png_df = tables['rec.vtx.elastic.fuzzyk.png.cvnpart']
83 
84  cvn_em_pid_df = cvn_png_df[['photonid', \
85  'pizeroid', \
86  'electronid']].sum(axis=1)
87  cvn_had_pid_df = cvn_png_df[['pionid', \
88  'protonid', \
89  'neutronid',\
90  'otherid', \
91  'muonid']].sum(axis=1)
92 
93  cvn_em_calE = shwlid_df['calE'].where( \
94  (cvn_em_pid_df > 0) & \
95  (cvn_em_pid_df >= cvn_had_pid_df), \
96  0).groupby(level=KL).agg(np.sum)
97 
98  if isRHC.agg(np.all):
99  cvn_em_calE[cvn_em_calE == 0] = prim_png[cvn_em_calE == 0]
100  else:
101  cvn_em_calE[lng_png >= 500] = prim_png[lng_png >= 500]
102 
103  hdr_df = tables['rec.hdr']
104  det = hdr_df['det']
105  ismc = tables['rec.hdr']['ismc']
106  runs = hdr_df.assign(run=hdr_df.index.get_level_values('run'))['run']
107  df = pd.concat([runs, det, ismc], axis=1).dropna()
108  scale = df.apply(lambda x: kNueCalibrationCorr(x['det'], x['ismc'], x['run']), axis=1, result_type='reduce')
109 
110  cvn_em_calE *= scale
111  cvn_em_calE.name = 'calE'
112 
113  return cvn_em_calE
114 kEMEnergy = Var(kEMEnergy)
115 
116 def kHadEnergy(tables):
117  EMEnergy = kEMEnergy(tables)
118 
119  hdr_df = tables['rec.hdr']
120  det = hdr_df['det']
121  ismc = tables['rec.hdr']['ismc']
122  runs = hdr_df.assign(run=hdr_df.index.get_level_values('run'))['run']
123  df = pd.concat([runs, det, ismc], axis=1).dropna()
124  scale = df.apply(lambda x: kNueCalibrationCorr(x['det'], x['ismc'], x['run']), axis=1, result_type='reduce')
125 
126  calE = tables['rec.slc']['calE']*scale
127  calE.name = 'calE'
128 
129  HadEnergy = calE - EMEnergy
130  return HadEnergy.where(HadEnergy > 0, 0)
131 kHadEnergy = Var(kHadEnergy)
132 
133 def kNueEnergy(tables):
134  EMEnergy = kEMEnergy(tables)
135  HadEnergy = kHadEnergy(tables)
136  isRHC = kRHC(tables)
137 
138  p0 = 0.0
139  p1 = 1.00756
140  p2 = 1.07093
141  p3 = 0.0
142  p4 = 1.28608e-02
143  p5 = 2.27129e-01
144  norm = 0.0501206
145  if isRHC.agg(np.all):
146  p0 = 0.0
147  p1 = 0.980479
148  p2 = 1.45170
149  p3 = 0.0
150  p4 = -5.82609e-03
151  p5 = -2.27599e-01
152  norm = 0.001766
153 
154  NueEnergy = 1./(1+norm)*(HadEnergy*HadEnergy*p5 + \
155  EMEnergy*EMEnergy*p4 + \
156  EMEnergy*HadEnergy*p3 + \
157  HadEnergy*p2 + \
158  EMEnergy*p1 + p0)
159  return NueEnergy.where((HadEnergy >= 0) & (EMEnergy >= 0), -5.)
160 kNueEnergy = Var(kNueEnergy)
161 
162 ###################################################################################
163 #
164 # Numu Vars
165 #
166 ###################################################################################
167 
168 
169 def kCosNumi(tables):
170  df = tables['rec.trk.kalman.tracks']
171  # Primary kalman track only
172  df = df[df['rec.trk.kalman.tracks_idx']==0]
173  KalDir = df[['dir.x','dir.y', 'dir.z']]
174 
175  # Use separate beam dir for each detector
176  KalDirFD = KalDir[kDetID(tables)==detector.kFD]
177  KalDirND = KalDir[kDetID(tables)==detector.kND]
178 
179  CosFD = KalDirFD.mul(BeamDirFD, axis=1).sum(axis=1)
180  CosND = KalDirND.mul(BeamDirND, axis=1).sum(axis=1)
181 
182  return pd.concat([CosFD, CosND])
183 kCosNumi = Var(kCosNumi)
184 
185 def kNumuMuE(tables):
186  det = kDetID(tables)
187  isRHC = kRHC(tables)
188 
189  hdr_df = tables['rec.hdr']
190  runs = hdr_df.assign(run=hdr_df.index.get_level_values('run'))['run']
191  ntracks = (tables['rec.trk.kalman']['ntracks'] != 0)
192 
193  if (det == detector.kFD).agg(np.all) and not det.empty:
194  ismc = tables['rec.hdr']['ismc']
195 
196  trks = tables['rec.trk.kalman.tracks']
197  trklen = trks['len'][(trks['rec.trk.kalman.tracks_idx'] == 0)]/100
198  if not ismc.agg(np.all) and not ismc.empty:
199  trklen = trklen*0.9957
200 
201  df = pd.concat([runs, det, isRHC, trklen], axis=1).dropna()
202  muE = df.apply(lambda x: \
203  GetSpline(x[0], x[1], x[2], "muon")(x[3]), axis = 1)
204  return muE.where(ntracks, -5.)
205 
206  else:
207  trklenact = tables['rec.energy.numu']['ndtrklenact']/100.
208  trklencat = tables['rec.energy.numu']['ndtrklencat']/100.
209  trkcalactE = tables['rec.energy.numu']['ndtrkcalactE']
210  trkcaltranE = tables['rec.energy.numu']['ndtrkcaltranE']
211 
212  df = pd.concat([runs, det, isRHC, trklenact, trklencat], axis=1)
213  muE = df.apply(lambda x: \
214  GetSpline(x['run'], x['det'], x['isRHC'], "act")(x['ndtrklenact']) + \
215  GetSpline(x['run'], x['det'], x['isRHC'], "cat")(x['ndtrklencat']), axis = 1)
216  muE[(trkcalactE == 0.) & (trkcaltranE == 0.)] = -5.
217  return muE.where(ntracks, -5.)
218 kNumuMuE = Var(kNumuMuE)
219 
220 def kNumuHadE(tables):
221  det = kDetID(tables)
222  isRHC = kRHC(tables)
223 
224  hdr_df = tables['rec.hdr']
225  runs = hdr_df.assign(run=hdr_df.index.get_level_values('run'))['run']
226 
227  hadvisE = tables['rec.energy.numu']['hadtrkE'] + tables['rec.energy.numu']['hadcalE']
228  hadvisE.name = 'hadvisE'
229  ntracks = (tables['rec.trk.kalman']['ntracks'] != 0)
230 
231  if (det == detector.kFD).agg(np.all) and not det.empty:
232  ismc = tables['rec.hdr']['ismc']
233  if not ismc.agg(np.all) and not ismc.empty:
234  periods = runs.apply(lambda x: GetPeriod(x, det))
235  hadvisE[periods > 2] = 0.9949*hadvisE[periods > 2]
236  hadvisE[periods <= 2] = 0.9844*hadvisE[periods <= 2]
237 
238  df = pd.concat([runs, det, isRHC, hadvisE], axis=1)
239  hadE = df.apply(lambda x: \
240  GetSpline(x['run'], x['det'], x['isRHC'], "had")(x['hadvisE']), axis = 1)
241  return hadE.where(ntracks, -5.)
242 kNumuHadE = Var(kNumuHadE)
243 
244 kNumuE = kNumuMuE + kNumuHadE
245 kCCE = kNumuE
246 
247 ###################################################################################
248 #
249 # Nus Vars
250 #
251 ###################################################################################
252 
253 def NusScale(isRHC, det):
254  if (isRHC and det == detector.kFD): return 1.18
255  if (isRHC and det == detector.kND): return 1.15
256  if (not isRHC and det == detector.kFD): return 1.2
257  if (not isRHC and det == detector.kND): return 1.11
258  else: return -5.
259 
260 def kNusEnergy(tables):
261  det = kDetID(tables)
262  isRHC = kRHC(tables)
263  cale = kCaloE(tables)
264  df = pd.concat([det, isRHC, cale], axis=1)
265  return df.apply(lambda x: \
266  NusScale(x['isRHC'], x['det'])*x['calE'], axis = 1)
267 kNusEnergy = Var(kNusEnergy)
268 
269 kClosestSlcTime = Var(lambda tables: tables['rec.slc']['closestslicetime'])
270 kClosestSlcMinDist = Var(lambda tables: tables['rec.slc']['closestslicemindist'])
271 kClosestSlcMinTop = Var(lambda tables: tables['rec.slc']['closestsliceminfromtop'])
def NusScale(isRHC, det)
Nus Vars.
def GetSpline(run, det, isRHC, comp)
Definition: numuE_utils.py:128
def kNueCalibrationCorr(det, ismc, run)
Double_t sum
Definition: plot.C:31
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
def GetPeriod(run, det)
Definition: misc.py:17