NumuCC1PiCuts.h
Go to the documentation of this file.
1 #pragma once
2 #include "CAFAna/Core/Cut.h"
4 #include "CAFAna/Vars/Vars.h"
5 //#include "3FlavorAna/Cuts/NumuCuts.h"
6 
7 #include "TVector3.h"
8 
9 //#include "StandardRecord/Proxy/SRProxy.h"
11 //#include "NDAna/numucc_1Pi/NumuCC1PiCuts.h"
12 //#include "NDAna/numucc_1Pi/NumuCC1PiVars.h"
13 
14 namespace ana
15 {
16 
17  // Fiducial volume limits
18  const TVector3 my_vtxmin(-150, -140, 150);
19  const TVector3 my_vtxmax( 150, 150, 1050);
20 
21 
22  //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
23 
24  //------NuTruthCuts-----------------
25 
27  {
28  int nPions = 0;
29  int nMuons = 0;
30 
31  int nPrims = nu->prim.size();
32  for(int prim_idx = 0; prim_idx < nPrims; prim_idx++){
33  auto &prim = nu->prim[prim_idx];
34  int pdg = prim.pdg;
35 
36  if(abs(pdg) == 211) ++nPions;
37  else if (pdg == 13) ++nMuons;
38 
39  }// end for prim
40  bool is_numucc = nu->iscc && nu->pdg == 14;
41 
42  bool ret = is_numucc && (nPions == 1) && (nMuons == 1) ;
43 
44  /*
45  if(ret) std::cout << "good" << std::endl;
46  else std::cout << "not good" << std::endl;
47  */
48  return ret;
49 
50  });
51 
52 
54  {
55  bool is_numucc = nu->iscc && nu->pdg == 14;
56 
57  return is_numucc;
58  });
59 
60 
62  {
63 
64 
65  return (nu->vtx.X() < my_vtxmax.X() &&
66  nu->vtx.X() > my_vtxmin.X() &&
67  nu->vtx.Y() < my_vtxmax.Y() &&
68  nu->vtx.Y() > my_vtxmin.Y() &&
69  nu->vtx.Z() < my_vtxmax.Z() &&
70  nu->vtx.Z() > my_vtxmin.Z() );
71 
72  });
73 
74 
76  {
77 
78  bool fid =
79  (nu->vtx.X() < 192 &&
80  nu->vtx.X() > -191 &&
81  nu->vtx.Y() < 194 &&
82  nu->vtx.Y() > -187 &&
83  nu->vtx.Z() < 1270 &&
84  nu->vtx.Z() > 0 );
85 
86  /*
87  if(fid) std::cout << "fid" << std::endl;
88  else std::cout << "not fid" << std::endl;
89  */
90 
91  return fid;
92 
93  // TODO here I've extended limits to *active* detector limits
94  // to use to determine fiducucial and containment cuts.
95  // I hope excluding muon-catcher here is ok.
96  });
97 
98 
99 
100  //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
101 
104 
105  //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
106 
107  //------Regular Cuts---------------------
108 
109  // Same as kTrueFidSig_NT but additionally require event formed slice
110  const Cut kTrueFidSig = CutFromNuTruthCut(kTrueFidSig_NT);
111 
112 
113  const Cut kFiducialRecoVtx([](const caf::SRProxy* sr)
114  //Reco vtx in fiducial volume
115  {
116 
117  if(sr->vtx.elastic.IsValid == false) return false;
118 
119  bool isFid = (sr->vtx.elastic.vtx.x < my_vtxmax.X() &&
120  sr->vtx.elastic.vtx.x > my_vtxmin.X() &&
121  sr->vtx.elastic.vtx.y < my_vtxmax.Y() &&
122  sr->vtx.elastic.vtx.y > my_vtxmin.Y() &&
123  sr->vtx.elastic.vtx.z < my_vtxmax.Z() &&
124  sr->vtx.elastic.vtx.z > my_vtxmin.Z() );
125 
126  return isFid;
127  });
128 
129  const Cut kRecoDetector([](const caf::SRProxy* sr)
130  //Reco vtx in fiducial volume
131  {
132 
133  if(sr->vtx.elastic.IsValid == false) return false;
134 
135  bool isFid = (sr->vtx.elastic.vtx.x < 192 &&
136  sr->vtx.elastic.vtx.x > -191 &&
137  sr->vtx.elastic.vtx.y < 194 &&
138  sr->vtx.elastic.vtx.y > -187 &&
139  sr->vtx.elastic.vtx.z < 1270 &&
140  sr->vtx.elastic.vtx.z > 0 );
141 
142  return isFid;
143  });
144 
145 
146 
147 
148  const Cut kMyProngQuality([](const caf::SRProxy* sr)
149  {
150 
151  if(sr->vtx.elastic.IsValid == false) return false;
152 
153  int nPngs = sr->vtx.elastic.fuzzyk.npng;
154 
155  bool isQuality = ( nPngs > 0 &&
156  sr->slc.nhit > 20 &&
157  sr->slc.ncontplanes > 4);
158 
159  return isQuality;
160 
161  });
162 
163 
164  const Cut kMaxNPngs(const int max)
165  {
166 
167  const Cut kNProngs([max](const caf::SRProxy* sr)
168  {
169  if(sr->vtx.elastic.IsValid == false) return false;
170  int nPngs = sr->vtx.elastic.fuzzyk.npng;
171  if(nPngs > max) return false;
172  else return true;
173  });
174 
175  return kNProngs;
176  }
177 
178  /*
179  const Cut kGoodMuonPngExists([](const caf::SRProxy* sr)
180  {
181  int muon_idx = kBestBPFTrack(sr);
182  float muon_id = kBPFMuonID(sr); // highest muonID
183 
184  int nPngs = sr->vtx.elastic.fuzzyk.npng;
185 
186  if(muon_idx < 0 || muon_idx >= nPngs) return false;
187  if(muon_id < 0.24) return false; // TODO this cut needs to be optimized
188  else return true;
189  });
190 
191 
192  const Cut kGoodPionPngExists([](const caf::SRProxy* sr)
193  {
194  int pion_idx = kBestPionPng(sr);
195  int nPngs = sr->vtx.elastic.fuzzyk.npng;
196 
197  if(pion_idx < 0 || pion_idx >= nPngs) return false;
198  else return true;
199  });
200  */
201 
202  const Cut kMinPngCalE(const double min)
203  // requires that ALL prongs have at least <min> GeV calE
204  {
205 
206  const Cut kPngsCalE([min](const caf::SRProxy* sr)
207  {
208  if(sr->vtx.elastic.IsValid == false) return false;
209  int nPngs = sr->vtx.elastic.fuzzyk.npng;
210 
211  for(int i=0; i<nPngs; ++i){
212  double calE = sr->vtx.elastic.fuzzyk.png[i].calE;
213  if(calE < min) return false;
214  }
215  return true;
216  });
217 
218  return kPngsCalE;
219  }
220 
221 
222 
223  const Cut kMinPngNhit(const int min)
224  // requires that ALL prongs have at least <min> nHit
225  {
226 
227  const Cut kPngsNhit([min](const caf::SRProxy* sr)
228  {
229  if(sr->vtx.elastic.IsValid == false) return false;
230  int nPngs = sr->vtx.elastic.fuzzyk.npng;
231 
232  for(int i=0; i<nPngs; ++i){
233  int nHit = sr->vtx.elastic.fuzzyk.png[i].nhit;
234  if(nHit < min) return false;
235  }
236  return true;
237  });
238 
239  return kPngsNhit;
240  }
241 
242 
243  const Cut kAnyPionProngs([](const caf::SRProxy* sr){
244  // Are there any 3d prongs who trace back to a pion or pion daughter
245  if(sr->vtx.elastic.IsValid == false) return false;
246 
247  bool anyPions = false;
248  int nPngs = sr->vtx.elastic.fuzzyk.png.size();
249 
250  for(int png_idx = 0; png_idx < nPngs; png_idx++){
251  auto &png = sr->vtx.elastic.fuzzyk.png[png_idx];
252  if( abs(png.truth.pdg) == 211 || abs(png.truth.motherpdg) == 211 ){
253  anyPions = true;
254  }
255  }
256  return anyPions;
257  });
258 
259  const Cut kAnyPionProngs2d([](const caf::SRProxy* sr){
260  // Are there any 2d prongs who trace back to a pion or pion daughter
261  if(sr->vtx.elastic.IsValid == false) return false;
262 
263  bool anyPions2d = false;
264  int nPngs2d = sr->vtx.elastic.fuzzyk.png2d.size();
265 
266  for(int png_idx = 0; png_idx < nPngs2d; png_idx++){
267  auto &png2d = sr->vtx.elastic.fuzzyk.png2d[png_idx];
268  if( abs(png2d.truth.pdg) == 211 || abs(png2d.truth.motherpdg) == 211 ){
269  anyPions2d = true;
270  }
271  }
272  return anyPions2d;
273  });
274 
275 
276 
277  const Cut kAnyProtonProngs([](const caf::SRProxy* sr){
278  // Are there any 3d prongs who trace back to a proton or proton daughter
279  if(sr->vtx.elastic.IsValid == false) return false;
280 
281  bool anyProtons = false;
282  int nPngs = sr->vtx.elastic.fuzzyk.png.size();
283 
284  for(int png_idx = 0; png_idx < nPngs; png_idx++){
285  auto &png = sr->vtx.elastic.fuzzyk.png[png_idx];
286  if( abs(png.truth.pdg) == 2212 || abs(png.truth.motherpdg) == 2212 ){
287  anyProtons = true;
288  }
289  }
290  return anyProtons;
291  });
292 
293  const Cut kAnyProtonProngs2d([](const caf::SRProxy* sr){
294  // Are there any 2d prongs who trace back to a proton or proton daughter
295  if(sr->vtx.elastic.IsValid == false) return false;
296 
297  bool anyProtons2d = false;
298  int nPngs2d = sr->vtx.elastic.fuzzyk.png2d.size();
299 
300  for(int png_idx = 0; png_idx < nPngs2d; png_idx++){
301  auto &png2d = sr->vtx.elastic.fuzzyk.png2d[png_idx];
302  if( abs(png2d.truth.pdg) == 2212 || abs(png2d.truth.motherpdg) == 2212 ){
303  anyProtons2d = true;
304  }
305  }
306  return anyProtons2d;
307  });
308 
309  //---------------
310 
311  const Cut kAnyPionProngs_noMother([](const caf::SRProxy* sr){
312  // removed check of motherpdg
313  if(sr->vtx.elastic.IsValid == false) return false;
314 
315  bool anyPions = false;
316  int nPngs = sr->vtx.elastic.fuzzyk.png.size();
317 
318  for(int png_idx = 0; png_idx < nPngs; png_idx++){
319  auto &png = sr->vtx.elastic.fuzzyk.png[png_idx];
320  if( abs(png.truth.pdg) == 211 ){
321  anyPions = true;
322  }
323  }
324  return anyPions;
325  });
326 
327  const Cut kAnyPionProngs2d_noMother([](const caf::SRProxy* sr){
328 
329  if(sr->vtx.elastic.IsValid == false) return false;
330 
331  bool anyPions2d = false;
332  int nPngs2d = sr->vtx.elastic.fuzzyk.png2d.size();
333 
334  for(int png_idx = 0; png_idx < nPngs2d; png_idx++){
335  auto &png2d = sr->vtx.elastic.fuzzyk.png2d[png_idx];
336  if( abs(png2d.truth.pdg) == 211 ){
337  anyPions2d = true;
338  }
339  }
340  return anyPions2d;
341  });
342 
343  //--------------
344 
345  const Cut kAnyPionTracks([](const caf::SRProxy* sr){
346 
347  int nTrks = sr->trk.kalman.tracks.size();
348 
349  if(nTrks < 1) return false;
350 
351  bool anyPions = false;
352 
353  for(int trk_idx = 0; trk_idx < nTrks; trk_idx++){
354  auto &trk = sr->trk.kalman.tracks[trk_idx];
355  if( abs(trk.truth.pdg) == 211 || abs(trk.truth.motherpdg) == 211 ){
356  anyPions = true;
357  }
358  }
359  return anyPions;
360  });
361 
362  //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
363  /*
364  Cut kTrueNumuCCPiExc = CutFromNuTruthCut(kIsNumuCC1Pi_NT);
365  */
366  //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
367 
368  // TODO - below should switch to using kTrueFidSig to bring
369  // in line with denom of efficiency.
370  // Currently no requirement on true vtx.
371 
372 
373  // All prongs contained
374  //
375  // Gonna keep this here for posterity, but I'm not sure this is a
376  // preselection. Surely presel shoudln't have TrueSignal in it. Maybe I just
377  // need to rename this cut.
378  /*
379  const Cut kMyPngPresel = kTrueFidSig &&
380  kFiducialRecoVtx &&
381  kMyProngQuality &&
382  kSimpleProngContainment;
383 
384 
385  const Cut kMyPngPresel_tight = kTrueFidSig &&
386  kFiducialRecoVtx &&
387  kMyProngQuality &&
388  kSimpleProngContainment_tight;
389 
390  //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
391  const Cut kRecoFidSig = kTrueNumuCCPiExc &&
392  kFiducialRecoVtx;
393 
394 
395  const Cut kAnyPionProngs_anyView = kAnyPionProngs || kAnyPionProngs2d;
396 
397  const Cut kAnyPionProngs_anyView_noMother = kAnyPionProngs_noMother ||
398  kAnyPionProngs2d_noMother;
399  */
400  //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
401 
402 
403 
404  //--------------------
405 
406  const Cut kProngContainment([](const caf::SRProxy* sr)
407 
408  // Reimplementing kNumuTightContainND but without any dependance on
409  // kalman tracks
410  {
411 
412  if( sr->vtx.elastic.IsValid == false ) return false;
413 
414  //std::cout << "FOO 1" << std::endl;
415 
416  int iMuonPng = kBestMuonPng(sr);
417 
418  // Reconstructed showers all contained
419  // (I would like to understand showerLID better or else use
420  // something else here)
421  for( unsigned int i = 0; i < sr->vtx.elastic.fuzzyk.nshwlid; ++i ) {
422 
423  const TVector3 start = sr->vtx.elastic.fuzzyk.png[i].shwlid.start;
424  const TVector3 stop = sr->vtx.elastic.fuzzyk.png[i].shwlid.stop;
425 
426  if( std::min( start.X(), stop.X() ) < -180.0 ) return false;
427  if( std::max( start.X(), stop.X() ) > 180.0 ) return false;
428  if( std::min( start.Y(), stop.Y() ) < -180.0 ) return false;
429  if( std::max( start.Y(), stop.Y() ) > 180.0 ) return false;
430  if( std::min( start.Z(), stop.Z() ) < 20.0 ) return false;
431  if( std::max( start.Z(), stop.Z() ) > 1525.0 ) return false;
432  }
433 
434  //std::cout << "FOO 2" << std::endl;
435 
436  // only primary muon BPF track present in muon catcher
437 
438  int nPngs = sr->vtx.elastic.fuzzyk.npng;
439  if(nPngs < 1) return false;
440 
441  for( int i = 0; i < nPngs; ++i ) {
442 
443  if( i == iMuonPng ) continue;
444 
445  // TODO should I be using shwlid here or BPF?
446  // Shwlid may be too tight a constraint, but if we use BPF
447  // we require that all prongs form BPF track. Also we need to
448  // choose BPF hypothesis. Should look at shwlid in evd to assess
449  // this
450  else if( sr->vtx.elastic.fuzzyk.png[i].shwlid.start.Z() > 1275 ||
451  sr->vtx.elastic.fuzzyk.png[i].shwlid.stop.Z() > 1275 )
452  return false;
453  }
454 
455  //std::cout << "FOO 3" << std::endl;
456 
457  // Removing this, under the assumption that the user
458  // will apply a cut ensuring there is a good muon prong
459  // e.g. kHighestMuonCVN > 0.3
460  /*
461  // If there is no muon prong, we are done (cont.)
462  if(!kGoodMuonPngExists(sr)) return true; // |
463  // |
464  // ______________________________________|
465  // |
466  // V
467  // (cont.) Otherwise we need to check muon prong
468  */
469 
470  if(iMuonPng < 0 || iMuonPng > nPngs) return false;
471 
472  auto &muon_prong = sr->vtx.elastic.fuzzyk.png[iMuonPng];
473 
474  // TODO this doesn't seem right - if muon doesn't form
475  // BPF track we say event is uncontained??
476  if(muon_prong.bpf.muon.IsValid == false) return false;
477 
478  auto &muon_bpf = muon_prong.bpf.muon;
479 
480  //std::cout << "FOO 4" << std::endl;
481 
482  bool ret = ( ( muon_bpf.stop.z < 1275 || muon_bpf.trkyposattrans < 55 )
483  // ^^^ air gap
484  && muon_bpf.trkfwdcellnd > 5
485  && muon_bpf.trkbakcellnd > 10);
486 
487 
488  //if(ret == true ) std::cout << "FOO 5" << std::endl;
489 
490  return ret;
491  });
492 
493  //-----------------------------
494 
495 
496  const Cut kSimpleProngContainment([](const caf::SRProxy* sr)
497 
498  // Very simple contaiment cut, with no assumption about where the muon is.
499  //
500  // N.B. - as a result all particles are allow in muon catcher, and there is
501  // no air gap checking
502  {
503 
504  if( sr->vtx.elastic.IsValid == false ) return false;
505 
506  // Reconstructed showers all contained
507  // (I would like to understand showerLID better or else use
508  // something else here)
509  for( unsigned int i = 0; i < sr->vtx.elastic.fuzzyk.nshwlid; ++i ) {
510 
511  const TVector3 start = sr->vtx.elastic.fuzzyk.png[i].shwlid.start;
512  const TVector3 stop = sr->vtx.elastic.fuzzyk.png[i].shwlid.stop;
513 
514  if( std::min( start.X(), stop.X() ) < -180.0 ) return false;
515  if( std::max( start.X(), stop.X() ) > 180.0 ) return false;
516  if( std::min( start.Y(), stop.Y() ) < -180.0 ) return false;
517  if( std::max( start.Y(), stop.Y() ) > 180.0 ) return false;
518  if( std::min( start.Z(), stop.Z() ) < 20.0 ) return false;
519  if( std::max( start.Z(), stop.Z() ) > 1525.0 ) return false;
520  }
521 
522  return true;
523  });
524 
525 
526 
528 
529  // Same as above but we require no particles in muon catcher
530  {
531 
532  if( sr->vtx.elastic.IsValid == false ) return false;
533 
534  // Reconstructed showers all contained
535  // (I would like to understand showerLID better or else use
536  // something else here)
537  for( unsigned int i = 0; i < sr->vtx.elastic.fuzzyk.nshwlid; ++i ) {
538 
539  const TVector3 start = sr->vtx.elastic.fuzzyk.png[i].shwlid.start;
540  const TVector3 stop = sr->vtx.elastic.fuzzyk.png[i].shwlid.stop;
541 
542  if( std::min( start.X(), stop.X() ) < -180.0 ) return false;
543  if( std::max( start.X(), stop.X() ) > 180.0 ) return false;
544  if( std::min( start.Y(), stop.Y() ) < -180.0 ) return false;
545  if( std::max( start.Y(), stop.Y() ) > 180.0 ) return false;
546  if( std::min( start.Z(), stop.Z() ) < 20.0 ) return false;
547  if( std::max( start.Z(), stop.Z() ) > 1275.0 ) return false;
548  }
549 
550  return true;
551  });
552 
553  //--------------------------------
554 
555 
556 
558 
559  // Very simple contaiment cut, with no assumption about where the muon is.
560  //
561  // N.B. - as a result all particles are allow in muon catcher, and there is
562  // no air gap checking
563  {
564 
565  if( sr->vtx.elastic.IsValid == false ) return false;
566 
567  // Reconstructed showers all contained
568  // (I would like to understand showerLID better or else use
569  // something else here)
570  for( unsigned int i = 0; i < sr->vtx.elastic.fuzzyk.nshwlid; ++i ) {
571 
572  const TVector3 start = sr->vtx.elastic.fuzzyk.png[i].shwlid.start;
573  const TVector3 stop = sr->vtx.elastic.fuzzyk.png[i].shwlid.stop;
574 
575  if( std::min( start.X(), stop.X() ) < -191.0 ) return false;
576  if( std::max( start.X(), stop.X() ) > 194.0 ) return false;
577  if( std::min( start.Y(), stop.Y() ) < -187.0 ) return false;
578  if( std::max( start.Y(), stop.Y() ) > 194.0 ) return false;
579  if( std::min( start.Z(), stop.Z() ) < 0.0 ) return false;
580  if( std::max( start.Z(), stop.Z() ) > 1587.0 ) return false;
581  }
582 
583  return true;
584  });
585 
586 
587 
588  //--------------------------------
589 
590  const Cut kTrue([](const caf::SRProxy* sr)
591  {
592  return true;
593  });
594 
595  /*
596  const Cut kBestPionPng_truePdg(int pdg)
597  // True if true PdgCode == pdg for BestPionPng
598  //
599  // Allows me to look at distribution of pion CVN score for best pion,
600  // broken down by true pdg.
601  {
602  const Cut kTruePdg([pdg](const caf::SRProxy* sr)
603  {
604  if(!kGoodPionPngExists(sr)) return false;
605 
606  int idx = kBestPionPng(sr);
607  auto &png = sr->vtx.elastic.fuzzyk.png[idx];
608 
609  if(png.truth.pdg == pdg) return true;
610  else return false;
611  });
612 
613  return kTruePdg;
614  }
615 
616 
617  const Cut kBestPionPng_calE_slice(double low, double high)
618  // True if for best pion prong: low <= calE < high
619  {
620 
621  const Cut kCalE([low,high](const caf::SRProxy* sr)
622  {
623  if(!kGoodPionPngExists(sr)) return false;
624 
625  int pion_idx = kBestPionPng(sr);
626  auto &png = sr->vtx.elastic.fuzzyk.png[pion_idx];
627 
628  double calE = png.calE;
629 
630  if(calE >= low && calE < high) return true;
631  else return false;
632 
633  });
634 
635  return kCalE;
636  }
637  */
638 
639 
640 
641  const Cut kMuonCVN_correct([](const caf::SRProxy* sr)
642  {
643  // True if prong with highest muon cvn score is true muon
644  // or if longest prong is true muon (if any prongs > 500cm)
645 
646  if(sr->vtx.elastic.IsValid == false) return false;
647 
648  int nPng = sr->vtx.elastic.fuzzyk.npng;
649 
650  double maxCVN = -5.0;
651  double maxLen = -5.0;
652  int ibest_muon = -5;
653 
654  bool evalCVN = true;
655 
656  for(int png_idx = 0; png_idx < nPng; ++png_idx){
657 
658  auto &png = sr->vtx.elastic.fuzzyk.png[png_idx];
659  double len = png.len;
660 
661  if(len > 500.0 && len > maxLen){
662  maxLen = len;
663  ibest_muon = png_idx;
664  evalCVN = false;
665  }
666 
667  if(evalCVN){
668  double cvn = png.cvnpart.muonid;
669  if(cvn > maxCVN ){
670  maxCVN = cvn;
671  ibest_muon = png_idx;
672  }
673  }
674 
675  } // end for png
676 
677  if(ibest_muon < 0 || ibest_muon > nPng) return false;
678 
679  auto &bestmuonpng = sr->vtx.elastic.fuzzyk.png[ibest_muon];
680  if(bestmuonpng.truth.pdg == 13 || bestmuonpng.truth.motherpdg == 13){
681  return true;
682  }
683 
684  return false;
685  });
686 
687 
688  /*
689  const Cut kPre_Presel = kFiducialRecoVtx &&
690  kProngContainment &&
691  kMyProngQuality &&
692  kMinPngCalE(0.150);
693 
694  const Cut kPresel_3prongs = kPre_Presel &&
695  (kN3dPngs == 3);
696 
697 
698  const Cut kPresel_2prongs = kPre_Presel &&
699  (kN3dPngs == 2);
700 
701  const Cut kNumuCC_sel = kPresel &&
702  kGoodMuonPngExists &&
703  kGoodPionPngExists;
704  // As it stands, kGoodPionPngExists just enforces
705  // that there is more than one prong
706  */
707 
708  // Selection {kFiducialRecoVtx, kSimpleProngContainement, kMyProngQuality,
709  // kN3dPngs == 2, kMinPngNhit(5)}
710  //
711  // Signal {CutFromNuTruthCut(kIsNumuCC1Pi_NT && kTrueFiducial_NT)}
712  }
713 
caf::Proxy< size_t > npng
Definition: SRProxy.h:2037
caf::Proxy< unsigned int > nshwlid
Definition: SRProxy.h:2039
T max(const caf::Proxy< T > &a, T b)
const Cut kAnyPionTracks([](const caf::SRProxy *sr){int nTrks=sr->trk.kalman.tracks.size();if(nTrks< 1) return false;bool anyPions=false;for(int trk_idx=0;trk_idx< nTrks;trk_idx++){auto &trk=sr->trk.kalman.tracks[trk_idx];if(abs(trk.truth.pdg)==211||abs(trk.truth.motherpdg)==211){anyPions=true;}}return anyPions;})
const NuTruthCut kTrueFidSig_NT
caf::Proxy< caf::SRFuzzyK > fuzzyk
Definition: SRProxy.h:2058
const Cut kMyProngQuality([](const caf::SRProxy *sr){if(sr->vtx.elastic.IsValid==false) return false;int nPngs=sr->vtx.elastic.fuzzyk.npng;bool isQuality=(nPngs > 0 && sr->slc.nhit > 20 && sr->slc.ncontplanes > 4);return isQuality;})
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
caf::Proxy< std::vector< caf::SRProng > > png2d
Definition: SRProxy.h:2043
const Cut kRecoDetector([](const caf::SRProxy *sr){if(sr->vtx.elastic.IsValid==false) return false;bool isFid=(sr->vtx.elastic.vtx.x< 192 && sr->vtx.elastic.vtx.x >-191 && sr->vtx.elastic.vtx.y< 194 && sr->vtx.elastic.vtx.y >-187 && sr->vtx.elastic.vtx.z< 1270 && sr->vtx.elastic.vtx.z > 0);return isFid;})
const Cut kAnyPionProngs2d_noMother([](const caf::SRProxy *sr){if(sr->vtx.elastic.IsValid==false) return false;bool anyPions2d=false;int nPngs2d=sr->vtx.elastic.fuzzyk.png2d.size();for(int png_idx=0;png_idx< nPngs2d;png_idx++){auto &png2d=sr->vtx.elastic.fuzzyk.png2d[png_idx];if(abs(png2d.truth.pdg)==211){anyPions2d=true;}}return anyPions2d;})
const NuTruthCut kTrueFiducial_NT([](const caf::SRNeutrinoProxy *nu){return(nu->vtx.X()< my_vtxmax.X()&& nu->vtx.X() > my_vtxmin.X()&& nu->vtx.Y()< my_vtxmax.Y()&& nu->vtx.Y() > my_vtxmin.Y()&& nu->vtx.Z()< my_vtxmax.Z()&& nu->vtx.Z() > my_vtxmin.Z());})
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2125
caf::Proxy< unsigned int > ncontplanes
Definition: SRProxy.h:1313
void abs(TH1 *hist)
const Cut kAnyProtonProngs2d([](const caf::SRProxy *sr){if(sr->vtx.elastic.IsValid==false) return false;bool anyProtons2d=false;int nPngs2d=sr->vtx.elastic.fuzzyk.png2d.size();for(int png_idx=0;png_idx< nPngs2d;png_idx++){auto &png2d=sr->vtx.elastic.fuzzyk.png2d[png_idx];if(abs(png2d.truth.pdg)==2212||abs(png2d.truth.motherpdg)==2212){anyProtons2d=true;}}return anyProtons2d;})
const Cut kSimpleProngContainment_loose([](const caf::SRProxy *sr) {if(sr->vtx.elastic.IsValid==false) return false; for(unsigned int i=0;i< sr->vtx.elastic.fuzzyk.nshwlid;++i){const TVector3 start=sr->vtx.elastic.fuzzyk.png[i].shwlid.start;const TVector3 stop=sr->vtx.elastic.fuzzyk.png[i].shwlid.stop;if(std::min(start.X(), stop.X())< -191.0) return false;if(std::max(start.X(), stop.X()) > 194.0) return false;if(std::min(start.Y(), stop.Y())< -187.0) return false;if(std::max(start.Y(), stop.Y()) > 194.0) return false;if(std::min(start.Z(), stop.Z())< 0.0) return false;if(std::max(start.Z(), stop.Z()) > 1587.0) return false;}return true;})
const Cut kMinPngCalE(const double min)
_Cut< caf::SRNeutrinoProxy > NuTruthCut
Cut designed to be used over the nuTree, ie all neutrinos, not just those that got slices...
Definition: Cut.h:104
const NuTruthCut kTrueDetector_NT([](const caf::SRNeutrinoProxy *nu){bool fid=(nu->vtx.X()< 192 && nu->vtx.X() >-191 && nu->vtx.Y()< 194 && nu->vtx.Y() >-187 && nu->vtx.Z()< 1270 && nu->vtx.Z() > 0);return fid; })
caf::Proxy< caf::SRElastic > elastic
Definition: SRProxy.h:2117
Track finder for cosmic rays.
caf::Proxy< std::vector< caf::SRFuzzyKProng > > png
Definition: SRProxy.h:2042
const Cut kAnyPionProngs([](const caf::SRProxy *sr){if(sr->vtx.elastic.IsValid==false) return false;bool anyPions=false;int nPngs=sr->vtx.elastic.fuzzyk.png.size();for(int png_idx=0;png_idx< nPngs;png_idx++){auto &png=sr->vtx.elastic.fuzzyk.png[png_idx];if(abs(png.truth.pdg)==211||abs(png.truth.motherpdg)==211){anyPions=true;}}return anyPions;})
caf::Proxy< unsigned int > nhit
Definition: SRProxy.h:1314
caf::Proxy< float > z
Definition: SRProxy.h:107
caf::Proxy< caf::SRTrackBranch > trk
Definition: SRProxy.h:2144
caf::Proxy< float > x
Definition: SRProxy.h:105
const Var kBestMuonPng([](const caf::SRProxy *sr){ if(sr->vtx.elastic.IsValid==false) return-5;int nPng=sr->vtx.elastic.fuzzyk.npng;double maxLen=-5;int maxLen_idx=-5;double maxCVN=-5;int maxCVN_idx=-5;for(int png_idx=0;png_idx< nPng;++png_idx){auto &png=sr->vtx.elastic.fuzzyk.png[png_idx];double len=png.len;if(len > maxLen){maxLen=len;maxLen_idx=png_idx;}double cvn=png.cvnpart.muonid;if(cvn > maxCVN){maxCVN=cvn;maxCVN_idx=png_idx;}}if(maxLen > 500) return maxLen_idx;return maxCVN_idx;})
const NuTruthCut kMy_IsNumuCC_NT([](const caf::SRNeutrinoProxy *nu){bool is_numucc=nu->iscc &&nu->pdg==14;return is_numucc;})
const Cut kAnyPionProngs_noMother([](const caf::SRProxy *sr){if(sr->vtx.elastic.IsValid==false) return false;bool anyPions=false;int nPngs=sr->vtx.elastic.fuzzyk.png.size();for(int png_idx=0;png_idx< nPngs;png_idx++){auto &png=sr->vtx.elastic.fuzzyk.png[png_idx];if(abs(png.truth.pdg)==211){anyPions=true;}}return anyPions;})
const Cut kMaxNPngs(const int max)
_Cut< caf::SRProxy > Cut
Representation of a cut (selection) to be applied to a caf::StandardRecord object.
Definition: Cut.h:96
const Cut kSimpleProngContainment_tight([](const caf::SRProxy *sr){if(sr->vtx.elastic.IsValid==false) return false; for(unsigned int i=0;i< sr->vtx.elastic.fuzzyk.nshwlid;++i){const TVector3 start=sr->vtx.elastic.fuzzyk.png[i].shwlid.start;const TVector3 stop=sr->vtx.elastic.fuzzyk.png[i].shwlid.stop;if(std::min(start.X(), stop.X())< -180.0) return false;if(std::max(start.X(), stop.X()) > 180.0) return false;if(std::min(start.Y(), stop.Y())< -180.0) return false;if(std::max(start.Y(), stop.Y()) > 180.0) return false;if(std::min(start.Z(), stop.Z())< 20.0) return false;if(std::max(start.Z(), stop.Z()) > 1275.0) return false;}return true;})
caf::StandardRecord * sr
const TVector3 my_vtxmin(-150,-140, 150)
const Cut kTrue([](const caf::SRProxy *sr){return true;})
const NuTruthCut kIsNumuCC1Pi_NT([](const caf::SRNeutrinoProxy *nu){int nPions=0;int nMuons=0;int nPrims=nu->prim.size();for(int prim_idx=0;prim_idx< nPrims;prim_idx++){auto &prim=nu->prim[prim_idx];int pdg=prim.pdg;if(abs(pdg)==211)++nPions;else if(pdg==13)++nMuons;}bool is_numucc=nu->iscc &&nu->pdg==14;bool ret=is_numucc &&(nPions==1)&&(nMuons==1);return ret;})
const Var kNProngs([](const caf::SRProxy *sr){return(sr->vtx.elastic[0].fuzzyk.npng);})
const Cut kMuonCVN_correct([](const caf::SRProxy *sr){ if(sr->vtx.elastic.IsValid==false) return false;int nPng=sr->vtx.elastic.fuzzyk.npng;double maxCVN=-5.0;double maxLen=-5.0;int ibest_muon=-5;bool evalCVN=true;for(int png_idx=0;png_idx< nPng;++png_idx){auto &png=sr->vtx.elastic.fuzzyk.png[png_idx];double len=png.len;if(len > 500.0 &&len > maxLen){maxLen=len;ibest_muon=png_idx;evalCVN=false;}if(evalCVN){double cvn=png.cvnpart.muonid;if(cvn > maxCVN){maxCVN=cvn;ibest_muon=png_idx;}}}if(ibest_muon< 0||ibest_muon > nPng) return false;auto &bestmuonpng=sr->vtx.elastic.fuzzyk.png[ibest_muon];if(bestmuonpng.truth.pdg==13||bestmuonpng.truth.motherpdg==13){return true;}return false;})
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
caf::Proxy< caf::SRKalman > kalman
Definition: SRProxy.h:1796
const Cut kTrueFidSig
const Cut kSimpleProngContainment([](const caf::SRProxy *sr) {if(sr->vtx.elastic.IsValid==false) return false; for(unsigned int i=0;i< sr->vtx.elastic.fuzzyk.nshwlid;++i){const TVector3 start=sr->vtx.elastic.fuzzyk.png[i].shwlid.start;const TVector3 stop=sr->vtx.elastic.fuzzyk.png[i].shwlid.stop;if(std::min(start.X(), stop.X())< -180.0) return false;if(std::max(start.X(), stop.X()) > 180.0) return false;if(std::min(start.Y(), stop.Y())< -180.0) return false;if(std::max(start.Y(), stop.Y()) > 180.0) return false;if(std::min(start.Z(), stop.Z())< 20.0) return false;if(std::max(start.Z(), stop.Z()) > 1525.0) return false;}return true;})
caf::Proxy< bool > IsValid
Definition: SRProxy.h:2057
caf::Proxy< caf::SRSlice > slc
Definition: SRProxy.h:2141
const Cut kFiducialRecoVtx([](const caf::SRProxy *sr){if(sr->vtx.elastic.IsValid==false) return false;bool isFid=(sr->vtx.elastic.vtx.x< my_vtxmax.X()&& sr->vtx.elastic.vtx.x > my_vtxmin.X()&& sr->vtx.elastic.vtx.y< my_vtxmax.Y()&& sr->vtx.elastic.vtx.y > my_vtxmin.Y()&& sr->vtx.elastic.vtx.z< my_vtxmax.Z()&& sr->vtx.elastic.vtx.z > my_vtxmin.Z());return isFid;})
caf::Proxy< float > y
Definition: SRProxy.h:106
const Cut kAnyProtonProngs([](const caf::SRProxy *sr){if(sr->vtx.elastic.IsValid==false) return false;bool anyProtons=false;int nPngs=sr->vtx.elastic.fuzzyk.png.size();for(int png_idx=0;png_idx< nPngs;png_idx++){auto &png=sr->vtx.elastic.fuzzyk.png[png_idx];if(abs(png.truth.pdg)==2212||abs(png.truth.motherpdg)==2212){anyProtons=true;}}return anyProtons;})
const TVector3 my_vtxmax(150, 150, 1050)
caf::Proxy< caf::SRVector3D > vtx
Definition: SRProxy.h:2072
const Cut kMinPngNhit(const int min)
caf::Proxy< std::vector< caf::SRKalmanTrack > > tracks
Definition: SRProxy.h:1779
T min(const caf::Proxy< T > &a, T b)
Template for Cut and SpillCut.
Definition: Cut.h:15
caf::Proxy< caf::SRVertexBranch > vtx
Definition: SRProxy.h:2145
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
const Cut kProngContainment([](const caf::SRProxy *sr) {if(sr->vtx.elastic.IsValid==false) return false;int iMuonPng=kBestMuonPng(sr); for(unsigned int i=0;i< sr->vtx.elastic.fuzzyk.nshwlid;++i){const TVector3 start=sr->vtx.elastic.fuzzyk.png[i].shwlid.start;const TVector3 stop=sr->vtx.elastic.fuzzyk.png[i].shwlid.stop;if(std::min(start.X(), stop.X())< -180.0) return false;if(std::max(start.X(), stop.X()) > 180.0) return false;if(std::min(start.Y(), stop.Y())< -180.0) return false;if(std::max(start.Y(), stop.Y()) > 180.0) return false;if(std::min(start.Z(), stop.Z())< 20.0) return false;if(std::max(start.Z(), stop.Z()) > 1525.0) return false;} int nPngs=sr->vtx.elastic.fuzzyk.npng;if(nPngs< 1) return false;for(int i=0;i< nPngs;++i){if(i==iMuonPng) continue; else if(sr->vtx.elastic.fuzzyk.png[i].shwlid.start.Z() > 1275|| sr->vtx.elastic.fuzzyk.png[i].shwlid.stop.Z() > 1275) return false;} if(iMuonPng< 0||iMuonPng > nPngs) return false;auto &muon_prong=sr->vtx.elastic.fuzzyk.png[iMuonPng]; if(muon_prong.bpf.muon.IsValid==false) return false;auto &muon_bpf=muon_prong.bpf.muon;bool ret=((muon_bpf.stop.z< 1275||muon_bpf.trkyposattrans< 55) &&muon_bpf.trkfwdcellnd > 5 &&muon_bpf.trkbakcellnd > 10);return ret;})
Cut CutFromNuTruthCut(const NuTruthCut &stc)
Definition: Cut.cxx:7
const Cut kAnyPionProngs2d([](const caf::SRProxy *sr){if(sr->vtx.elastic.IsValid==false) return false;bool anyPions2d=false;int nPngs2d=sr->vtx.elastic.fuzzyk.png2d.size();for(int png_idx=0;png_idx< nPngs2d;png_idx++){auto &png2d=sr->vtx.elastic.fuzzyk.png2d[png_idx];if(abs(png2d.truth.pdg)==211||abs(png2d.truth.motherpdg)==211){anyPions2d=true;}}return anyPions2d;})