analysis_cuts.py
Go to the documentation of this file.
1 import pandas as pd
2 import numpy as np
3 
4 from PandAna.core.core import Cut, KL
5 from PandAna.utils.enums import *
6 from PandAna.var.analysis_vars import *
7 
8 kIsFD = kDetID == detector.kFD
9 
10 ###################################################################################
11 #
12 # Basic Cuts
13 #
14 ###################################################################################
15 
16 # Basic cosrej
17 kVeto = Cut(lambda tables: tables['rec.sel.veto']['keep'] == 1)
18 
19 # Basic Reco Cuts
20 kHasVtx = Cut(lambda tables: tables['rec.vtx']['nelastic'] > 0)
21 kHasPng = Cut(lambda tables: tables['rec.vtx.elastic.fuzzyk']['npng'] > 0)
22 
23 ###################################################################################
24 #
25 # Nue Cuts
26 #
27 ###################################################################################
28 
29 #Data Quality
30 
32  mask = l[0]
33 
34  fp = l[1]
35  fpmin = fp
36  fpmax = fp
37 
38  lp = l[2]
39  lpmin = lp
40  lpmax = lp
41 
42  for i in range(fp, 14, 1):
43  if mask[13-i] == '0':
44  break
45  else:
46  fpmax = i
47 
48  for i in range(fp, -1, -1):
49  if mask[13-i] == '0':
50  break
51  else:
52  fpmin = i
53 
54  for i in range(lp, 14, 1):
55  if mask[13-i] == '0':
56  break
57  else:
58  lpmax = i
59 
60  for i in range(lp, -1, -1):
61  if mask[13-i] == '0':
62  break
63  else:
64  lpmin = i
65  return (fpmin==lpmin) & (fpmax==lpmax) & (lpmax-fpmin+1>=4)
66 
67 def kNueApplyMask(tables):
68  mask = tables['rec.hdr']['dibmask']
69  fp = tables['rec.slc']['firstplane']//64
70  lp = tables['rec.slc']['lastplane']//64
71  df = mask.apply(lambda x: bin(x)[2:].zfill(14))
72  df = pd.concat([df,fp,lp],axis=1)
73  return df.apply(kDibMaskHelper, axis=1)
74 kNueApplyMask = Cut(kNueApplyMask)
75 
76 kNueDQ = kHasVtx & kHasPng & (kHitsPerPlane < 8)
77 
78 kNueBasicPart = kVeto & kIsFD & kNueDQ & kNueApplyMask
79 
80 # Presel
81 kNuePresel = (kNueEnergy > 1) & (kNueEnergy < 4) & \
82  (kNHit > 30) & (kNHit < 150) & \
83  (kLongestProng > 100) & (kLongestProng < 500)
84 
85 kNueProngContainment = (kDistAllTop > 63) & (kDistAllBottom > 12) & \
86  (kDistAllEast > 12) & (kDistAllWest > 12) & \
87  (kDistAllFront > 18) & (kDistAllBack > 18)
88 
89 kNueBackwardCut = ((kDistAllBack < 200) & (kSparsenessAsymm < -0.1)) | (kDistAllBack >= 200)
90 
91 kNuePtPCut = (kPtP < 0.58) | ((kPtP >= 0.58) & (kPtP < 0.8) & (kMaxY < 590)) | ((kPtP >= 0.8) & (kMaxY < 350))
92 
93 kNueCorePart = kNueProngContainment & kNueBackwardCut & kNuePtPCut & kNuePresel
94 
95 kNueCorePresel = kNueBasicPart & kNueCorePart
96 
97 # PID Selections
98 kNueCVNFHC = 0.84
99 kNueCVNRHC = 0.89
100 
101 def kNueCVNCut(tables):
102  df = kCVNe(tables)
103  dfRHC = df[kRHC(tables)==1] >= kNueCVNRHC
104  dfFHC = df[kRHC(tables)!=1] >= kNueCVNFHC
105 
106  return pd.concat([dfRHC, dfFHC])
107 kNueCVNCut = Cut(kNueCVNCut)
108 
109 # Full FD Selection
110 kNueFD = kNueCVNCut & kNueCorePresel
111 
112 def kNueNDFiducial(tables):
113  check = tables['rec.vtx.elastic']['rec.vtx.elastic_idx'] == 0
114  df = tables['rec.vtx.elastic'][check]
115  return (df['vtx.x'] > -100) & \
116  (df['vtx.x'] < 160) & \
117  (df['vtx.y'] > -160) & \
118  (df['vtx.y'] < 100) & \
119  (df['vtx.z'] > 150) & \
120  (df['vtx.z'] < 900)
121 kNueNDFiducial = Cut(kNueNDFiducial)
122 
123 def kNueNDContain(tables):
124  df = tables['rec.vtx.elastic.fuzzyk.png.shwlid']
125  df_trans = df[['start.y','stop.y', 'start.x', 'stop.x']]
126  df_long = df[['start.z', 'stop.z']]
127 
128  return ((df_trans.min(axis=1) > -170) & (df_trans.max(axis=1) < 170) & \
129  (df_long.min(axis=1) > 100) & (df_long.max(axis=1) < 1225)).groupby(level=KL).agg(np.all)
130 kNueNDContain = Cut(kNueNDContain)
131 
132 kNueNDFrontPlanes = Cut(lambda tables: tables['rec.sel.contain']['nplanestofront'] > 6)
133 
134 kNueNDNHits = (kNHit >= 20) & (kNHit <= 200)
135 
136 kNueNDEnergy = (kNueEnergy < 4.5)
137 
138 kNueNDProngLength = (kLongestProng > 100) & (kLongestProng < 500)
139 
140 kNueNDPresel = kNueDQ & kNueNDFiducial & kNueNDContain & kNueNDFrontPlanes & \
141  kNueNDNHits & kNueNDEnergy & kNueNDProngLength
142 
143 kNueNDCVNSsb = kNueNDPresel & kNueCVNCut
144 
145 
146 
147 ###################################################################################
148 #
149 # Numu Cuts
150 #
151 ###################################################################################
152 
153 def kNumuBasicQuality(tables):
154  df_numutrkcce=tables['rec.energy.numu']['trkccE']
155  df_remid=tables['rec.sel.remid']['pid']
156  df_nhit=tables['rec.slc']['nhit']
157  df_ncontplanes=tables['rec.slc']['ncontplanes']
158  df_cosmicntracks=tables['rec.trk.cosmic']['ntracks']
159  return(df_numutrkcce > 0) &\
160  (df_remid > 0) &\
161  (df_nhit > 20) &\
162  (df_ncontplanes > 4) &\
163  (df_cosmicntracks > 0)
164 kNumuBasicQuality = Cut(kNumuBasicQuality)
165 
166 kNumuQuality = kNumuBasicQuality & (kCCE < 5.)
167 
168 # FD
169 
170 kNumuProngsContainFD = (kDistAllTop > 60) & (kDistAllBottom > 12) & (kDistAllEast > 16) & \
171  (kDistAllWest > 12) & (kDistAllFront > 18) & (kDistAllBack > 18)
172 
174  mask = l[0]
175 
176  fd = l[1]//64
177  ld = l[2]//64
178 
179  dmin = 0
180  dmax = 13
181 
182  for i in range(fd, 14, 1):
183  if mask[13-i] == '0':
184  break
185  else:
186  dmax = i
187 
188  for i in range(fd, -1, -1):
189  if mask[13-i] == '0':
190  break
191  else:
192  dmin = i
193 
194  return ((l[1]-64*dmin) > 1) & ((64*(dmax+1)-l[2]-1) > 1)
195 
197  mask = tables['rec.hdr']['dibmask']
198  fp = tables['rec.slc']['firstplane']
199  lp = tables['rec.slc']['lastplane']
200  df = mask.apply(lambda x: bin(x)[2:].zfill(14))
201  df = pd.concat([df,fp,lp],axis=1)
202  df = df.apply(kNumuDibMaskHelper, axis=1, result_type='reduce')
203 
204  df_containkalfwdcell = tables['rec.sel.contain']['kalfwdcell'] > 6
205  df_containkalbakcell = tables['rec.sel.contain']['kalbakcell'] > 6
206  df_containcosfwdcell = tables['rec.sel.contain']['cosfwdcell'] > 0
207  df_containcosbakcell = tables['rec.sel.contain']['cosbakcell'] > 7
208 
209  return df & df_containkalfwdcell & df_containkalbakcell & \
210  df_containcosfwdcell & df_containkalbakcell
211 kNumuOptimizedContainFD = Cut(kNumuOptimizedContainFD)
212 
213 kNumuContainFD = kNumuProngsContainFD & kNumuOptimizedContainFD
214 
215 kNumuNoPIDFD = kNumuQuality & kNumuContainFD
216 # ND
217 def kNumuContainND(tables):
218  check = tables['rec.vtx.elastic']['rec.vtx.elastic_idx'] == 0
219  shw_df = tables['rec.vtx.elastic.fuzzyk.png.shwlid'][check]
220  shw_df_trans = shw_df[['start.y','stop.y', 'start.x', 'stop.x']]
221  shw_df_long = shw_df[['start.z', 'stop.z']]
222  no_shw = (tables['rec.vtx.elastic.fuzzyk']['nshwlid'] == 0)
223 
224  shw_contain = ((shw_df_trans.min(axis=1) >= -180.) & (shw_df_trans.max(axis=1) <= 180.) & \
225  (shw_df_long.min(axis=1) >= 20.) & (shw_df_long.max(axis=1) <= 1525.)).groupby(level=KL).agg(np.all)
226  shw_contain = (shw_contain | no_shw)
227 
228  trk_df = tables['rec.trk.kalman.tracks'][['start.z', 'stop.z', 'rec.trk.kalman.tracks_idx']]
229  kalman_contain = False
230  if not trk_df.empty:
231  trk_df = trk_df.apply(lambda x: (x['rec.trk.kalman.tracks_idx'] == 0) | ((x['start.z'] <= 1275) & (x['stop.z'] <= 1275)), axis = 1)
232  kalman_contain = trk_df.groupby(level=KL).agg(np.all)
233 
234  df_ntracks = tables['rec.trk.kalman']['ntracks']
235  df_remid = tables['rec.trk.kalman']['idxremid']
236  df_firstplane = tables['rec.slc']['firstplane']
237  df_lastplane = tables['rec.slc']['lastplane']
238 
239  first_trk = tables['rec.trk.kalman.tracks']['rec.trk.kalman.tracks_idx'] == 0
240  df_startz = tables['rec.trk.kalman.tracks'][first_trk]['start.z']
241  df_stopz = tables['rec.trk.kalman.tracks'][first_trk]['stop.z']
242 
243  df_containkalposttrans = tables['rec.sel.contain']['kalyposattrans']
244  df_containkalfwdcellnd = tables['rec.sel.contain']['kalfwdcellnd']
245  df_containkalbakcellnd = tables['rec.sel.contain']['kalbakcellnd']
246 
247  return (df_ntracks > df_remid) &\
248  (df_firstplane > 1) &\
249  (df_lastplane < 212) &\
250  (df_containkalfwdcellnd > 5) &\
251  (df_containkalbakcellnd > 10) &\
252  (df_startz < 1100 ) & (( df_containkalposttrans < 55) | (df_stopz < 1275) ) &\
253  shw_contain &\
254  kalman_contain
255 
256 kNumuContainND = Cut(kNumuContainND)
257 
258 kNumuNCRej = Cut(lambda tables: tables['rec.sel.remid']['pid'] > 0.75)
259 
260 kNumuNoPIDND = kNumuQuality & kNumuContainND
261 
262 ###################################################################################
263 #
264 # Nus Cuts
265 #
266 ###################################################################################
267 
268 
269 # FD
270 
271 kNusFDContain = (kDistAllTop > 100) & (kDistAllBottom > 10) & \
272  (kDistAllEast > 50) & (kDistAllWest > 50) & \
273  (kDistAllFront > 50) & (kDistAllBack > 50)
274 
275 kNusContPlanes = Cut(lambda tables: tables['rec.slc']['ncontplanes'] > 2)
276 
277 kNusEventQuality = kHasVtx & kHasPng & \
278  (kHitsPerPlane < 8) & kNusContPlanes
279 
280 kNusFDPresel = kNueApplyMask & kVeto & kNusEventQuality & kNusFDContain
281 
282 kNusBackwardCut = ((kDistAllBack < 200) & (kSparsenessAsymm < -0.1)) | (kDistAllBack >= 200)
283 
284 kNusEnergyCut = (kNusEnergy >= 0.5) & (kNusEnergy <= 20.)
285 
286 kNusSlcTimeGap = (kClosestSlcTime > -150.) & (kClosestSlcTime < 50.)
287 kNusSlcDist = (kClosestSlcMinTop < 100) & (kClosestSlcMinDist < 500)
288 kNusShwPtp = ((kMaxY > 580) & (kPtP > 0.2)) | ((kMaxY > 540) & (kPtP > 0.4))
289 
290 # Nus Cosrej Cuts use TMVA trained BDT natively
291 kNusNoPIDFD = (kNusFDPresel & kNusBackwardCut) & (~(kNusSlcTimeGap & kNusSlcDist)) & \
292  (~kNusShwPtp) & kNusEnergyCut
293 
294 
295 # ND
296 
297 def kNusNDFiducial(tables):
298  check = tables['rec.vtx.elastic']['rec.vtx.elastic_idx'] == 0
299  df = tables['rec.vtx.elastic'][check]
300  return (df['vtx.x'] > -100) & \
301  (df['vtx.x'] < 100) & \
302  (df['vtx.y'] > -100) & \
303  (df['vtx.y'] < 100) & \
304  (df['vtx.z'] > 150) & \
305  (df['vtx.z'] < 1000)
306 kNusNDFiducial = Cut(kNusNDFiducial)
307 
308 kNusNDContain = (kDistAllTop > 25) & (kDistAllBottom > 25) & \
309  (kDistAllEast > 25) & (kDistAllWest > 25) & \
310  (kDistAllFront > 25) & (kDistAllBack > 25)
311 
312 kNusNDPresel = kNusEventQuality & kNusNDFiducial & kNusNDContain
313 kNusNoPIDND = kNusNDPresel & (kPtP <= 0.8) & kNusEnergyCut
314 
315 ###################################################################################
316 #
317 # OR'd cuts for Near and Far
318 #
319 ###################################################################################
320 # FD check. Hacky but it works with standard file name conventions
321 def kIsFarDet(tables):
322  query = tables._files.query
323  if not type(query) is list: query = [query]
324  return 'fardet' in query[0]
325 
326 ############### nue cuts ######################
327 def kNueContain(tables):
328  if kIsFarDet(tables):
329  return kNueProngContainment(tables) & kNueBasicPart(tables)
330  else:
331  return kNueNDContain(tables)
332 
333 kNueContain = Cut(kNueContain)
334 
335 def kNuePresel(tables):
336  if kIsFarDet(tables):
337  return kNueCorePresel(tables)
338  else:
339  return kNueNDPresel(tables)
340 
341 kNuePresel = Cut(kNuePresel)
342 
343 ############### numu cuts #####################
344 # kCCE isn't working yet
345 kNumuNoPIDNoCCEFD = kNumuBasicQuality & kNumuContainFD
346 kNumuNoPIDNoCCEND = kNumuBasicQuality & kNumuContainND
347 def kNumuContain(tables):
348  if kIsFarDet(tables):
349  return kNumuContainFD(tables)
350  else:
351  return kNumuContainND(tables)
352 
353 kNumuContain = Cut(kNumuContain)
354 
355 def kNumuPresel(tables):
356  if kIsFarDet(tables):
357  return kNumuNoPIDNoCCEFD(tables)
358  else:
359  return kNumuNoPIDNoCCEND(tables)
360 
361 kNumuPresel = Cut(kNumuPresel)
362 
363 ############### nus cuts ######################
364 def kNusContain(tables):
365  if kIsFarDet(tables):
366  return kNusFDContain(tables)
367  else:
368  return kNusNDContain(tables)
369 kNusContain = Cut(kNusContain)
370 
371 def kNusPresel(tables):
372  if kIsFarDet(tables):
373  return kNusFDPresel(tables)
374  else:
375  return kNusNDPresel(tables)
376 kNusPresel = Cut(kNusPresel)
377 
378 ############### ORd cuts #####################
379 kCosVeto = kVeto
380 kOrContainment = kNumuContain | kNusContain | kNueContain
381 kOrPreselection = kNumuPresel | kNusPresel | kNuePresel
382 
const Var kCVNe
PID
Definition: Vars.cxx:35
tuple kNusFDContain
Nus Cuts.
def kDibMaskHelper(l)
Nue Cuts.
kNumuNoPIDNoCCEFD
numu cuts ##################### kCCE isn&#39;t working yet
float bin[41]
Definition: plottest35.C:14
def kIsFarDet(tables)
OR&#39;d cuts for Near and Far.