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