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