NumuCC1PiVars.h
Go to the documentation of this file.
1 #pragma once
2 #include "CAFAna/Core/Var.h"
3 #include "CAFAna/Core/MultiVar.h"
4 
6 
7 
8 //#include "NDAna/nd_shared/NDXSecBPFMuonPID.h"
9 
10 #include <vector>
11 #include <iostream>
12 
13 #include "TVector.h"
14 
17 //#include "NDAna/numucc_1Pi/NumuCC1PiCuts.h"
18 
19 namespace ana
20 {
21 
22  //----- Truth Vars --------
23 
24 
26  {
27  // Kinetic energy of true charged pions
28  //
29  // NB: This returns one value per event, so it is really only
30  // applicable if there is a single pion.
31 
32 
33  int nPrims = nu->prim.size();
34  int nPions = 0;
35 
36  double energy = -5.0;
37 
38  for(int prim_idx = 0; prim_idx < nPrims; prim_idx++){
39  auto &prim = nu->prim[prim_idx];
40  int pdg = prim.pdg;
41 
42  if(abs(pdg) == 211){
43  //if(prim.pdg == 211){
44 
45  energy = prim.p.E;
46  energy -= 0.139;
47  nPions++;
48  }
49  }
50 
51  if(nPions != 1){
52  std::cout << nPions << "\n";
53  std::cout << "FOO \n";
54  abort();
55  }
56 
57 
58  return energy;
59  });
60 
61 
63  {
64  // True kinetic energy of most energetic proton
65 
66  double ret_value = -5;
67 
68  int nPrims = nu->prim.size();
69 
70  double energy = -5.0;
71 
72  for(int prim_idx = 0; prim_idx < nPrims; prim_idx++){
73  auto &prim = nu->prim[prim_idx];
74  int pdg = prim.pdg;
75 
76  //if(prim.pdg == 2212){
77  if(abs(pdg) == 2212){
78  energy = prim.p.E;
79  if(energy > ret_value) ret_value = energy;
80  }// end if(proton)
81 
82  }// end for(prim)
83 
84  ret_value -= 0.938;
85 
86  return ret_value;
87  });
88 
89 
90  const std::map<int, double> pdg_to_mass = {
91  {11, 0.0005},
92  {13, 0.105},
93  {22, 0.0},
94  {111, 0.134},
95  {211, 0.139},
96  {2212, 0.938},
97  };
98 
99  const NuTruthVar kNPions_NT([](const caf::SRNeutrinoProxy* nu)
100  {
101  int nPrims = nu->prim.size();
102 
103  int nPions = 0;
104 
105  for(int prim_idx = 0; prim_idx < nPrims; prim_idx++){
106  auto &prim = nu->prim[prim_idx];
107  int pdg = prim.pdg;
108 
109  if(abs(pdg) == 211){
110  nPions++;
111  }
112  }
113 
114  return nPions;
115 
116  });
117 
118  const NuTruthVar kNProtons_NT([](const caf::SRNeutrinoProxy* nu)
119  {
120  int nPrims = nu->prim.size();
121 
122  int nProtons = 0;
123 
124  for(int prim_idx = 0; prim_idx < nPrims; prim_idx++){
125  auto &prim = nu->prim[prim_idx];
126  int pdg = prim.pdg;
127 
128  if(abs(pdg) == 2212){
129  nProtons++;
130  }
131  }
132 
133  return nProtons;
134 
135  });
136 
137 
138 
139  // ------ Regular Vars -------
140 
141  const Var kTrueE_nu([](const caf::SRProxy* sr){
142  if(sr->mc.nnu < 1) return -5.f;
143  return float(sr->mc.nu[0].E);
144 
145  });
146 
147  const Var kN3dPngs([](const caf::SRProxy* sr){
148 
149  if(sr->vtx.elastic.IsValid == false ) return -5;
150 
151  int nPngs = sr->vtx.elastic.fuzzyk.png.size();
152 
153  return nPngs;
154 
155  });
156 
157 
158  const Var kN2dPngs([](const caf::SRProxy* sr){
159 
160  if(sr->vtx.elastic.IsValid == false ) return -5;
161 
162  int nPngs2d = sr->vtx.elastic.fuzzyk.png2d.size();
163 
164  return nPngs2d;
165 
166  });
167 
168  /*
169  const Var kTruePionTrackEnergy([](const caf::SRProxy* sr)
170  {
171 
172  int nTrks = sr->trk.kalman.ntracks;
173  double energy = -5.0;
174 
175  if(nTrks <= 0) return -5.0;
176 
177  for(int trk_idx = 0; trk_idx < nTrks; trk_idx++){
178  auto &trk = sr->trk.kalman.tracks[trk_idx];
179  if( abs(trk.truth.pdg) == 211 ){
180 
181  energy = trk.truth.p.E;
182  energy -= pion_mass;
183  }
184  }
185  return energy;
186  });
187 
188 
189  const Var kTruePionProngEnergy([](const caf::SRProxy* sr){
190 
191  if(sr->vtx.elastic.IsValid == false) return -5.0;
192 
193  int nPions = 0;
194  int nPngs = sr->vtx.elastic.fuzzyk.png.size();
195  double energy = -5.0;
196 
197  if(nPngs <= 0) return -5.0;
198 
199  for(int png_idx = 0; png_idx < nPngs; png_idx++){
200  auto &png = sr->vtx.elastic.fuzzyk.png[png_idx];
201  if(abs(png.truth.pdg) == 211){
202  // TODO: Do I need to add in motherpdg allowance above? Not sure its valid
203 
204 
205  // Just a sanity check
206  if(nPions != 0){
207  //std::cout << "-----------------WARNING----------------" << std::endl;
208  //std::cout << "Somehow there is more than one pion!! (Prong)" << std::endl;
209  }
210 
211  energy = png.truth.p.E;
212  energy -= pion_mass;
213  nPions++;
214  }
215  }
216 
217  return energy;
218  });
219  */
220 
223 
224 
227 
228 
229  const MultiVar kRun_byPdg(int pdg)
230  // Totally pointless multivar, where a regular var would do
231  // Only necessary because I want it in TH2 vs a multivar
232  {
233  const MultiVar kRetVar([pdg](const caf::SRProxy* sr)
234  {
235  std::vector<double> ret_vec;
236 
237  if(sr->vtx.elastic.IsValid == false) return ret_vec;
238 
239  int nPngs = sr->vtx.elastic.fuzzyk.npng;
240  for(int ipng=0; ipng<nPngs; ipng++){
241  auto &png = sr->vtx.elastic.fuzzyk.png[ipng];
242  if(png.truth.pdg == pdg){
243  ret_vec.push_back((double)sr->hdr.run);
244  }// end if(pdg)
245  }// end for(ipng)
246  return ret_vec;
247  });// end kRetVar
248  return kRetVar;
249  }// end kRun_byPdg
250 
251 
252  const MultiVar kProngEnergyRes_byPdg(int pdg, bool massSub = false,
253  double minGev = 0.0, double maxGev = 100.0)
254  // Energy resolution of prongs with backtraced
255  // PDG == pdg
256  // and minGev < E_true {- m} < maxGev
257  {
258  const MultiVar kProngERes([pdg, massSub, minGev, maxGev](const caf::SRProxy* sr)
259  {
260  std::vector<double> ERes_vec;
261 
262  double E_true;
263  double E_reco;
264  double E_res;
265 
266  std::map<int, double> mass_map{ { 13, 0.1057},
267  {2212, 0.9383},
268  { 211, 0.1396},
269  {-211, 0.1396}
270  };
271 
272  if(sr->vtx.elastic.IsValid == false) return ERes_vec;
273 
274  int nPngs = sr->vtx.elastic.fuzzyk.npng;
275 
276  for(int ipng=0; ipng<nPngs; ipng++){
277 
278  auto &png = sr->vtx.elastic.fuzzyk.png[ipng];
279 
280  if(png.truth.pdg == pdg){
281 
282 
283  double mass = mass_map.at(pdg);
284 
285  E_true = png.truth.p.E;
286 
287  if(massSub == true) E_true -= mass;
288 
289  if( E_true < minGev || E_true > maxGev){
290  ERes_vec.push_back(-5.0);
291  continue;
292  }
293 
294  if(pdg == 13){
295  E_reco = kMuE(sr);
296  }
297  else{
298  E_reco = png.calE;
299  }
300 
301 
302  E_res = (E_reco - E_true)/E_true;
303 
304  //std::cout << "E_res: " << E_res << std::endl;
305  ERes_vec.push_back(E_res);
306 
307  }
308 
309  }// end for ipng
310  return ERes_vec;
311  }); // end kProngERes
312 
313  return kProngERes;
314  }// end kProngEnergyRes_byPdg
315 
316 
318  // same as above, but Var instead of MultiVar
319  {
320  const Var kProngERes([pdg](const caf::SRProxy* sr)
321  {
322 
323  double E_true;
324  double E_reco;
325  double E_res;
326 
327  if(sr->vtx.elastic.IsValid == false) return -5.0;
328 
329  int nPngs = sr->vtx.elastic.fuzzyk.npng;
330 
331  for(int ipng=0; ipng<nPngs; ipng++){
332 
333  auto &png = sr->vtx.elastic.fuzzyk.png[ipng];
334 
335  if(png.truth.pdg == pdg){
336 
337  E_true = png.truth.p.E;
338 
339  if(pdg == 13){
340  E_reco = kMuE(sr);
341  }
342  else{
343  E_reco = png.calE;
344  }
345 
346  E_res = (E_reco - E_true)/E_true;
347  //std::cout << E_res << std::endl;
348  return E_res;
349 
350  }// end if(pdg)
351 
352  }// end for ipng
353  return -5.0;
354  }); // end kProngERes
355 
356  return kProngERes;
357  }// end kProngEnergyRes_byPdg_nonMulti
358 
359 
360 
361  const Var kProngPur_byPdg(int pdg)
362  // purity of most energetic prong matching pdg
363  // NB not checking motherpdg here as it doesn't feel right
364  {
365  const Var kPur([pdg](const caf::SRProxy* sr)
366  {
367  double maxCalE = 0.0;
368  int maxIdx = -5;
369  //double pur;
370 
371  if(sr->vtx.elastic.IsValid == false) return -5.0;
372  int nPngs = sr->vtx.elastic.fuzzyk.npng;
373  for(int ipng=0; ipng<nPngs; ipng++){
374  auto &png = sr->vtx.elastic.fuzzyk.png[ipng];
375 
376  if(abs(png.truth.pdg) != pdg) continue;
377 
378  if(png.calE > maxCalE){
379  maxCalE = png.calE;
380  maxIdx = ipng;
381  }
382  }// end for ipng
383 
384  if(maxIdx < 0 ) return -5.0;
385  double pur = sr->vtx.elastic.fuzzyk.png[maxIdx].truth.pur;
386  return pur;
387  }); // end kPur
388 
389  return kPur;
390  }// end kProngPur_byPdg
391 
392 
393  const Var kProngEff_byPdg(int pdg)
394  // efficiency of most energetic prong matching pdg
395  // NB not checking motherpdg here as it doesn't feel right
396  {
397  const Var kEff([pdg](const caf::SRProxy* sr)
398  {
399  double maxCalE = 0.0;
400  int maxIdx = -5;
401  //double eff;
402 
403 
404 
405  if(sr->vtx.elastic.IsValid == false) return -5.0;
406  int nPngs = sr->vtx.elastic.fuzzyk.npng;
407  for(int ipng=0; ipng<nPngs; ipng++){
408  auto &png = sr->vtx.elastic.fuzzyk.png[ipng];
409 
410  //if(abs(png.truth.pdg) != pdg) continue;
411  if(abs(png.truth.pdg) != pdg ) continue;
412 
413  if(png.calE > maxCalE){
414  maxCalE = png.calE;
415  maxIdx = ipng;
416  }
417  }// end for ipng
418 
419  if(maxIdx < 0 ) return -5.0;
420  double eff = sr->vtx.elastic.fuzzyk.png[maxIdx].truth.eff;
421 
422  return eff;
423  }); // end kEff
424 
425  return kEff;
426  }// end kProngEff_byPdg
427 
428 
429  //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
431  // Energy resolution of tracks with backtraced
432  // PDG == pdg
433  {
434  const MultiVar kTrackERes([pdg](const caf::SRProxy* sr)
435  {
436  std::vector<double> ERes_vec;
437 
438  double E_true;
439  double E_reco;
440  double E_res;
441 
442  int nTrks = sr->trk.kalman.ntracks;
443 
444  if(nTrks <=0 ) return ERes_vec;
445 
446  for(int itrk=0; itrk<nTrks; itrk++){
447 
448  auto &trk = sr->trk.kalman.tracks[itrk];
449 
450  if(trk.truth.pdg == pdg){
451 
452  E_true = trk.truth.p.E;
453 
454  if(pdg == 13){
455  E_reco = kMuE(sr);
456  }
457  else{
458  E_reco = trk.calE;
459  }
460 
461 
462  E_res = (E_reco - E_true)/E_true;
463  ERes_vec.push_back(E_res);
464 
465  }
466 
467  }// end for itrk
468  return ERes_vec;
469  }); // end kTrackERes
470 
471  return kTrackERes;
472  }// end kKalmanEnergyRes_byPdg
473 
474  //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
475 
477  // Shower start in either X, Y, or Z with backtraced
478  // PDG == pdg
479  {
480 
481  assert( (coord == "X") || (coord == "Y") || (coord == "Z") );
482 
483  const MultiVar kShowerStart([pdg, coord](const caf::SRProxy* sr)
484  {
485  std::vector<double> start_vec;
486 
487  if(sr->vtx.elastic.IsValid == false) return start_vec;
488 
489  int nPngs = sr->vtx.elastic.fuzzyk.npng;
490 
491  for(int ipng=0; ipng<nPngs; ipng++){
492 
493  auto &png = sr->vtx.elastic.fuzzyk.png[ipng];
494 
495  if(png.truth.pdg == pdg){
496 
497  const TVector3 start = png.shwlid.start;
498 
499  double ret;
500  if(coord == "X") ret = start.X();
501  else if(coord == "Y") ret = start.Y();
502  else ret = start.Z();
503 
504  start_vec.push_back(ret);
505 
506  }
507 
508  }// end for ipng
509  return start_vec;
510  }); // end kProngERes
511 
512  return kShowerStart;
513  }
514 
515 
516  const MultiVar kShowerStop_byPdg(int pdg, std::string coord)
517  // Shower stop in either X, Y, or Z with backtraced
518  // PDG == pdg
519  {
520 
521  assert( (coord == "X") || (coord == "Y") || (coord == "Z") );
522 
523  const MultiVar kShowerStop([pdg, coord](const caf::SRProxy* sr)
524  {
525  std::vector<double> stop_vec;
526 
527  if(sr->vtx.elastic.IsValid == false) return stop_vec;
528 
529  int nPngs = sr->vtx.elastic.fuzzyk.npng;
530 
531  for(int ipng=0; ipng<nPngs; ipng++){
532 
533  auto &png = sr->vtx.elastic.fuzzyk.png[ipng];
534 
535  if(png.truth.pdg == pdg){
536 
537  const TVector3 stop = png.shwlid.stop;
538 
539  double ret;
540  if(coord == "X") ret = stop.X();
541  else if(coord == "Y") ret = stop.Y();
542  else ret = stop.Z();
543 
544  stop_vec.push_back(ret);
545 
546  }
547 
548  }// end for ipng
549  return stop_vec;
550  }); // end kShowerStop
551 
552  return kShowerStop;
553  }
554 
555 
557  std::string coord,
558  std::string minmax = "mix")
559 
560 
561  // for each prong matching "pdg" we return min / max of
562  // (showerStart, showerStop) in the "coord" dimension
563  //
564  // If we specify "mix" then we return absolute value of distance of shower
565  // from nearest edge in "coord"
566  {
567  assert( (coord == "X") || (coord == "Y") || (coord == "Z") );
568  assert( (minmax == "min") || (minmax == "max") || (minmax == "mix") );
569 
570  const MultiVar kRetVar([pdg, coord, minmax](const caf::SRProxy* sr)
571  {
572 
573  std::vector<double> ret_vec;
574 
575  if(sr->vtx.elastic.IsValid == false) return ret_vec;
576  int nPngs = sr->vtx.elastic.fuzzyk.npng;
577 
578  for(int ipng=0; ipng<nPngs; ipng++){
579 
580  auto &png = sr->vtx.elastic.fuzzyk.png[ipng];
581 
582  if(png.truth.pdg == pdg){
583 
584  const TVector3 start = png.shwlid.start;
585  const TVector3 stop = png.shwlid.stop;
586 
587  //std::cout << "Shower Start Z: " << start.Z() << std::endl;
588  //std::cout << "Shower Stop Z: " << stop.Z() << std::endl;
589 
590  double min_dist ; //minimum distance from "left" edge
591  double max_dist ; //minimum distance from "right" edge
592 
593 
594  double ret;
595  if(coord == "X"){
596  if(minmax == "min"){
597  ret = std::min( start.X(), stop.X() );
598  }
599  else if(minmax == "max"){
600  ret = std::max( start.X(), stop.X() );
601  }
602  else{
603  min_dist = abs( -191 - std::min(start.X(), stop.X()) );
604  max_dist = abs( 192 - std::max(start.X(), stop.X()) );
605 
606  ret = std::min(min_dist, max_dist);
607  }
608  }// end if(X)
609  else if(coord == "Y"){
610  if(minmax == "min"){
611  ret = std::min( start.Y(), stop.Y() );
612  }
613  else if(minmax == "max"){
614  ret = std::max( start.Y(), stop.Y() );
615  }
616  else{
617 
618  min_dist = abs( -187 - std::min(start.Y(), stop.Y()) );
619  max_dist = abs( 194 - std::max(start.Y(), stop.Y()) );
620 
621  ret = std::min(min_dist, max_dist);
622  }
623  }// end if(Y)
624  else if(coord == "Z"){
625  if(minmax == "min"){
626  ret = std::min( start.Z(), stop.Z() );
627  }
628  else if(minmax == "max"){
629  ret = std::max( start.Z(), stop.Z() );
630  }
631  else{
632  min_dist = abs( 0 - std::min(start.Z(), stop.Z()) );
633  max_dist = abs( 1587 - std::max(start.Z(), stop.Z()) );
634 
635  ret = std::min(min_dist, max_dist);
636  }
637  }// end if(Z)
638 
639  ret_vec.push_back(ret);
640 
641  }
642 
643  }// end for ipng
644  return ret_vec;
645 
646  }); // end kRetVar
647 
648  return kRetVar;
649 
650  }// end kShowerMinMax_byPdg
651 
652 
654  std::string coord,
655  std::string minmax = "mix")
656  // same as above, but Var instead of Multivar. This is for purposes of
657  // handscanning investigation to see why there is a jump in energy
658  // resolution at ~440cm in Z
659  {
660  assert( (coord == "X") || (coord == "Y") || (coord == "Z") );
661  assert( (minmax == "min") || (minmax == "max") || (minmax == "mix") );
662 
663  const Var kRetVar([pdg, coord, minmax](const caf::SRProxy* sr)
664  {
665 
666  if(sr->vtx.elastic.IsValid == false) return -5.0;
667  int nPngs = sr->vtx.elastic.fuzzyk.npng;
668 
669  for(int ipng=0; ipng<nPngs; ipng++){
670 
671  auto &png = sr->vtx.elastic.fuzzyk.png[ipng];
672 
673  if(png.truth.pdg == pdg){
674 
675  const TVector3 start = png.shwlid.start;
676  const TVector3 stop = png.shwlid.stop;
677 
678  double min_dist ; //minimum distance from "left" edge
679  double max_dist ; //minimum distance from "right" edge
680 
681 
682  double ret = false; //initialize with value to avoid annoying warning
683  if(coord == "X"){
684  if(minmax == "min"){
685  ret = std::min( start.X(), stop.X() );
686  }
687  else if(minmax == "max"){
688  ret = std::max( start.X(), stop.X() );
689  }
690  else{
691  min_dist = abs( -191 - std::min(start.X(), stop.X()) );
692  max_dist = abs( 192 - std::max(start.X(), stop.X()) );
693 
694  ret = std::min(min_dist, max_dist);
695  }
696  }// end if(X)
697  else if(coord == "Y"){
698  if(minmax == "min"){
699  ret = std::min( start.Y(), stop.Y() );
700  }
701  else if(minmax == "max"){
702  ret = std::max( start.Y(), stop.Y() );
703  }
704  else{
705 
706  min_dist = abs( -187 - std::min(start.Y(), stop.Y()) );
707  max_dist = abs( 194 - std::max(start.Y(), stop.Y()) );
708 
709  ret = std::min(min_dist, max_dist);
710  }
711  }// end if(Y)
712  else if(coord == "Z"){
713  if(minmax == "min"){
714  ret = std::min( start.Z(), stop.Z() );
715  }
716  else if(minmax == "max"){
717  ret = std::max( start.Z(), stop.Z() );
718  }
719  else{
720  min_dist = abs( 0 - std::min(start.Z(), stop.Z()) );
721  max_dist = abs( 1587 - std::max(start.Z(), stop.Z()) );
722 
723  ret = std::min(min_dist, max_dist);
724  }
725  }// end if(Z)
726 
727  //std::cout << ret << std::endl;
728  return ret;
729 
730  }
731 
732  }// end for ipng
733  return -5.0;
734 
735  }); // end kRetVar
736 
737  return kRetVar;
738 
739  }// end kShowerMinMax_byPdg_nonMulti
740 
741 
742 
744  std::string coord,
745  std::string minmax)
746 
747  // same as kShowerMinMax_byPdg but for kalman tracks
748  {
749  assert( (coord == "X") || (coord == "Y") || (coord == "Z") );
750  assert( (minmax == "min") || (minmax == "max") );
751 
752  const MultiVar kRetVar([pdg, coord, minmax](const caf::SRProxy* sr)
753  {
754 
755  std::vector<double> ret_vec;
756 
757  int nTrks = sr->trk.kalman.ntracks;
758 
759  if(nTrks <= 0) return ret_vec;
760 
761 
762  for(int itrk=0; itrk<nTrks; itrk++){
763 
764  auto &trk = sr->trk.kalman.tracks[itrk];
765 
766  if(trk.truth.pdg == pdg){
767 
768  const TVector3 start = trk.start;
769  const TVector3 stop = trk.stop;
770 
771  double ret;
772  if(coord == "X"){
773  if(minmax == "min"){
774  ret = std::min( start.X(), stop.X() );
775  }
776  else{
777  ret = std::max( start.X(), stop.X() );
778  }
779  }// end if(X)
780  else if(coord == "Y"){
781  if(minmax == "min"){
782  ret = std::min( start.Y(), stop.Y() );
783  }
784  else{
785  ret = std::max( start.Y(), stop.Y() );
786  }
787  }// end if(Y)
788  else if(coord == "Z"){
789  if(minmax == "min"){
790  ret = std::min( start.Z(), stop.Z() );
791  }
792  else{
793  ret = std::max( start.Z(), stop.Z() );
794  }
795  }// end if(Z)
796 
797  ret_vec.push_back(ret);
798 
799  }
800 
801  }// end for itrk
802  return ret_vec;
803 
804  }); // end kRetVar
805 
806  return kRetVar;
807 
808  }// end kKalmanMinMax_byPdg
809 
810 
811  /*
812  const Var kScan_BPF([](const caf::SRProxy* sr){
813  // This is just a way to be able to scan TTree variables that I can't
814  // normally due to ROOT limitations on Tree depth
815 
816  if(sr->vtx.elastic.IsValid == false) return -5;
817 
818  int run = sr->hdr.run;
819  int subrun = sr->hdr.subrun;
820  int evt = sr->hdr.evt;
821  int slc = sr->hdr.subevt;
822  int cycle = sr->hdr.cycle;
823 
824  std::cout << std::endl;
825 
826  std::cout << run << " / " << subrun << " / " << evt << " / " << slc
827  << " / " << cycle << std::endl;
828 
829  for(uint png_idx = 0; png_idx < sr->vtx.elastic.fuzzyk.npng; png_idx++){
830 
831  auto &shwlid = sr->vtx.elastic.fuzzyk.png[png_idx].shwlid;
832  const TVector3 shwlid_start = shwlid.start;
833  const TVector3 shwlid_stop = shwlid.stop;
834 
835  auto &bpf = sr->vtx.elastic.fuzzyk.png[png_idx].bpf;
836  int n_bpf = bpf.size();
837  if(n_bpf < 1) return -5;
838 
839  int pdg = bpf[0].truth.pdg;
840  int motherpdg = bpf[0].truth.motherpdg;
841 
842  std::cout << "PDG : " << pdg << std::endl;
843  std::cout << "motherPDG : " << motherpdg << std::endl;
844 
845  std::cout << "shwlid_start : ";
846  printTVector(shwlid_start);
847 
848  std::cout << "shwlid_stop : ";
849  printTVector(shwlid_stop);
850  std::cout << std::endl;
851 
852 
853  for(int bpf_idx = 0; bpf_idx < n_bpf; bpf_idx++){
854  const TVector3 bpf_start = bpf[bpf_idx].start;
855  const TVector3 bpf_stop = bpf[bpf_idx].stop;
856 
857  std::cout << "bpf_pdg : " << bpf[bpf_idx].pdg << std::endl;
858 
859  std::cout << "bpf_start : ";
860  printTVector(bpf_start);
861 
862  std::cout << "bpf_stop : ";
863  printTVector(bpf_stop);
864 
865  std::cout << std::endl;
866 
867 
868  } // end for bpf
869  // std::cout << "--------" << std::endl;
870 
871  } // end for png
872  std::cout << "******************" << std::endl;
873  return 1;
874 
875  });
876  */
877 
878  const Var kTrueMuonProng([](const caf::SRProxy* sr)
879  // Returns idx of true muon prong
880  // TODO what is there is more than one muon prong?
881 
882  {
883  if(sr->vtx.elastic.IsValid == false) return -5;
884 
885  int itrue_muon_png = -5;
886 
887  int nPng = sr->vtx.elastic.fuzzyk.npng;
888 
889  for(int png_idx = 0; png_idx < nPng; png_idx++){
890  auto &bpf = sr->vtx.elastic.fuzzyk.png[png_idx].bpf;
891  //int n_bpf = bpf.size();
892  //if(n_bpf < 1) continue;
893  int pdg;
894  int motherpdg;
895 
896  if(bpf.muon.IsValid == true){
897  pdg = bpf.muon.truth.pdg;
898  motherpdg = bpf.muon.truth.motherpdg;
899  }
900  else if(bpf.pion.IsValid == true){
901  pdg = bpf.pion.truth.pdg;
902  motherpdg = bpf.pion.truth.motherpdg;
903  }
904  else if(bpf.proton.IsValid == true){
905  pdg = bpf.proton.truth.pdg;
906  motherpdg = bpf.proton.truth.motherpdg;
907  }
908  else continue;
909 
910  if(pdg == 13 && motherpdg == 13){
911  if(itrue_muon_png != -5){
912  std::cout << "More than one muon prong! There is a problem here..." << std::endl;
913  return -5;
914  }
915  itrue_muon_png = png_idx;
916  };
917 
918  } // end for png
919  return itrue_muon_png;
920  });
921 
922 
923 
924  const Var kBestMuonPng([](const caf::SRProxy* sr){
925  // Takes longest prong as muon if any prong > 500cm, otherwise takes
926  // prong with highest CVN muon score
927  //
928  // Removing the following as it makes it too confusing{, where muon score
929  // must also be higher than other scores for that prong}
930 
931  if(sr->vtx.elastic.IsValid == false) return -5;
932 
933  int nPng = sr->vtx.elastic.fuzzyk.npng;
934 
935  double maxLen = -5;
936  int maxLen_idx = -5;
937 
938  double maxCVN = -5;
939  int maxCVN_idx = -5;
940 
941  for(int png_idx = 0; png_idx < nPng; ++png_idx){
942  auto &png = sr->vtx.elastic.fuzzyk.png[png_idx];
943 
944  double len = png.len;
945  if(len > maxLen){
946  maxLen = len;
947  maxLen_idx = png_idx;
948  }
949 
950  double cvn = png.spprongcvnpart5label.muonid;
951  if(cvn > maxCVN){
952  maxCVN = cvn;
953  maxCVN_idx = png_idx;
954  }
955  } // end for png
956 
957  if(maxLen > 500) return maxLen_idx;
958 
959  return maxCVN_idx;
960 
961  });
962 
963  //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
964  const Var kBestPionPng([](const caf::SRProxy* sr)
965  {
966  // Returns idx of prong with highest cvn pion score excluding muon prong
967 
968  if(sr->vtx.elastic.IsValid == false) return -5;
969 
970  int nPng = sr->vtx.elastic.fuzzyk.npng;
971 
972  double maxCVN = -5;
973  int maxCVN_idx = -5;
974 
975  for(int png_idx = 0; png_idx < nPng; ++png_idx){
976 
977  //if(png_idx == kBestBPFTrack(sr)) continue;
978  if(png_idx == kBestMuonPng(sr)) continue;
979 
980  auto &png = sr->vtx.elastic.fuzzyk.png[png_idx];
981 
982  double cvn = png.spprongcvnpart5label.pionid;
983  if(cvn > maxCVN ){
984  maxCVN = cvn;
985  maxCVN_idx = png_idx;
986  }
987 
988  } // end for png
989 
990  return maxCVN_idx;
991  });
992 
993 
994  const Var kHighestPionCVN([](const caf::SRProxy* sr)
995  {
996  int piIdx = kBestPionPng(sr);
997  if(piIdx < 0) return -5.;
998 
999  double ret = sr->vtx.elastic.fuzzyk.png[piIdx].spprongcvnpart5label.pionid;
1000  return ret;
1001 
1002  });
1003 
1004 
1005  //~~~~~~~~~~~~~~~~~~~~~~~~~~~
1006 
1007 
1008  const Var kHighestPionCVN_KE([](const caf::SRProxy* sr)
1009  {
1010  int piIdx = kBestPionPng(sr);
1011  if(piIdx < 0) return -5.;
1012 
1013  double ret = sr->vtx.elastic.fuzzyk.png[piIdx].truth.p.E - 0.139;
1014  return ret;
1015  });
1016 
1017 
1018  const Var kHighestPionCVN_Cos([](const caf::SRProxy* sr)
1019  {
1020  int piIdx = kBestPionPng(sr);
1021  if(piIdx < 0) return -5.;
1022 
1023 
1024 
1025  auto &png = sr->vtx.elastic.fuzzyk.png[piIdx];
1026  TVector3 beamDir = NuMIBeamDirection(sr->hdr.det);
1027 
1028  double ret = (beamDir.X() * png.dir.x) +
1029  (beamDir.Y() * png.dir.y) +
1030  (beamDir.Z() * png.dir.z);
1031 
1032  return ret;
1033 
1034  });
1035 
1036 
1037 
1038  const Var kHighestMuonCVN([](const caf::SRProxy* sr)
1039  {
1040  // Returns highest cvn muon score, unless any prong is longer than 500 cm,
1041  // in which case returns 1.0
1042  // (Only prongs )
1043 // TODO should change this to use kBestMuonPng
1044  std::cout << "FOO \n";
1045  std::cout << "BAR \n";
1046 
1047  if(sr->vtx.elastic.IsValid == false) return -5.0;
1048 
1049  int nPng = sr->vtx.elastic.fuzzyk.npng;
1050 
1051  double maxCVN = -5.0;
1052 
1053  for(int png_idx = 0; png_idx < nPng; ++png_idx){
1054 
1055  auto &png = sr->vtx.elastic.fuzzyk.png[png_idx];
1056  double len = png.len;
1057 
1058  std::cout << "len: " << len << "\n";
1059 
1060  if(len > 500.0) return 1.0;
1061 
1062  double cvn = png.spprongcvnpart5label.muonid;
1063  if(cvn > maxCVN ){
1064  maxCVN = cvn;
1065  }
1066 
1067  } // end for png
1068 
1069  std::cout << "maxCVN: " << maxCVN << "\n";
1070 
1071  return maxCVN;
1072  });
1073 
1074 
1075 
1077  {
1078  // returns highest muon CVN if partID matches pdg
1079  const Var khighestmuoncvn([pdg](const caf::SRProxy* sr)
1080  {
1081  if(sr->vtx.elastic.IsValid == false) return -5.0;
1082 
1083  int nPng = sr->vtx.elastic.fuzzyk.npng;
1084  int iBestPng = kBestMuonPng(sr);
1085  if(iBestPng < 0 || iBestPng > nPng) return -5.0;
1086  auto &png = sr->vtx.elastic.fuzzyk.png[iBestPng];
1087 
1088  if(png.truth.pdg != pdg) return -5.0;
1089 
1090  if(png.len > 500.0) return 1.0;
1091 
1092  double score = png.cvnpart.muonid;
1093 
1094  return score;
1095  });
1096 
1097  return khighestmuoncvn;
1098  }
1099 
1100 
1101 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~`
1102 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~`
1103 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~`
1104 
1106  {
1107  // returns muon CVN score of all prongs which match pdg
1108  const MultiVar kmuoncvn([pdg](const caf::SRProxy* sr)
1109  {
1110  std::vector<double> ret;
1111 
1112  if(sr->vtx.elastic.IsValid == false) return ret;
1113 
1114  int nPng = sr->vtx.elastic.fuzzyk.npng;
1115 
1116  for(int iPng=0; iPng < nPng; iPng++){
1117  auto &png = sr->vtx.elastic.fuzzyk.png[iPng];
1118  if(png.truth.pdg != pdg) continue;
1119 
1120  double score;
1121  if(png.len > 500.0){
1122  score = 1.0;
1123  }
1124  else{
1125  score = png.spprongcvnpart5label.muonid;
1126  }
1127  ret.push_back(score);
1128  }//end loop over prongs
1129 
1130  return ret;
1131  });
1132  return kmuoncvn;
1133  }
1134 
1135 
1137  {
1138  // returns electron CVN score of all prongs which match pdg
1139  const MultiVar kelectroncvn([pdg](const caf::SRProxy* sr)
1140  {
1141  std::vector<double> ret;
1142 
1143  if(sr->vtx.elastic.IsValid == false) return ret;
1144 
1145  int nPng = sr->vtx.elastic.fuzzyk.npng;
1146 
1147  for(int iPng=0; iPng < nPng; iPng++){
1148  auto &png = sr->vtx.elastic.fuzzyk.png[iPng];
1149  if(png.truth.pdg != pdg) continue;
1150 
1151  double score = png.spprongcvnpart5label.electronid;
1152 
1153  ret.push_back(score);
1154  }//end loop over prongs
1155 
1156  return ret;
1157  });
1158  return kelectroncvn;
1159  }
1160 
1161 
1163  {
1164  // returns proton CVN score of all prongs which match pdg
1165  const MultiVar kprotoncvn([pdg](const caf::SRProxy* sr)
1166  {
1167  std::vector<double> ret;
1168 
1169  if(sr->vtx.elastic.IsValid == false) return ret;
1170 
1171  int nPng = sr->vtx.elastic.fuzzyk.npng;
1172 
1173  for(int iPng=0; iPng < nPng; iPng++){
1174  auto &png = sr->vtx.elastic.fuzzyk.png[iPng];
1175  if(png.truth.pdg != pdg) continue;
1176 
1177  double score = png.spprongcvnpart5label.protonid;
1178 
1179  ret.push_back(score);
1180  }//end loop over prongs
1181 
1182  return ret;
1183  });
1184  return kprotoncvn;
1185  }
1186 
1187 
1189  {
1190  // returns pion CVN score of all prongs which match pdg
1191  const MultiVar kpioncvn([pdg](const caf::SRProxy* sr)
1192  {
1193  std::vector<double> ret;
1194 
1195  if(sr->vtx.elastic.IsValid == false) return ret;
1196 
1197  int nPng = sr->vtx.elastic.fuzzyk.npng;
1198 
1199  for(int iPng=0; iPng < nPng; iPng++){
1200  auto &png = sr->vtx.elastic.fuzzyk.png[iPng];
1201  if(png.truth.pdg != pdg) continue;
1202 
1203  double score = png.spprongcvnpart5label.pionid;
1204 
1205  ret.push_back(score);
1206  }//end loop over prongs
1207 
1208  return ret;
1209  });
1210  return kpioncvn;
1211  }
1212 
1213 
1215  {
1216  // returns photon CVN score of all prongs which match pdg
1217  const MultiVar kphotoncvn([pdg](const caf::SRProxy* sr)
1218  {
1219  std::vector<double> ret;
1220 
1221  if(sr->vtx.elastic.IsValid == false) return ret;
1222 
1223  int nPng = sr->vtx.elastic.fuzzyk.npng;
1224 
1225  for(int iPng=0; iPng < nPng; iPng++){
1226  auto &png = sr->vtx.elastic.fuzzyk.png[iPng];
1227  if(png.truth.pdg != pdg) continue;
1228 
1229  double score = png.spprongcvnpart5label.photonid;
1230 
1231  ret.push_back(score);
1232  }//end loop over prongs
1233 
1234  return ret;
1235  });
1236  return kphotoncvn;
1237  }
1238 
1239 //****************************************
1240 
1241 
1243  {
1244  // returns muon CVN score of all prongs which match pdg
1245  const MultiVar kmuoncvn([pdg](const caf::SRProxy* sr)
1246  {
1247  std::vector<double> ret;
1248 
1249  if(sr->vtx.elastic.IsValid == false) return ret;
1250 
1251  int nPng = sr->vtx.elastic.fuzzyk.npng;
1252 
1253  for(int iPng=0; iPng < nPng; iPng++){
1254  auto &png = sr->vtx.elastic.fuzzyk.png[iPng];
1255  if(png.truth.pdg != pdg) continue;
1256 
1257  double score;
1258  if(png.len > 500.0){
1259  score = 1.0;
1260  }
1261  else{
1262  score = png.cvnpart.muonid;
1263  }
1264  ret.push_back(score);
1265  }//end loop over prongs
1266 
1267  return ret;
1268  });
1269  return kmuoncvn;
1270  }
1271 
1272 
1274  {
1275  // returns electron CVN score of all prongs which match pdg
1276  const MultiVar kelectroncvn([pdg](const caf::SRProxy* sr)
1277  {
1278  std::vector<double> ret;
1279 
1280  if(sr->vtx.elastic.IsValid == false) return ret;
1281 
1282  int nPng = sr->vtx.elastic.fuzzyk.npng;
1283 
1284  for(int iPng=0; iPng < nPng; iPng++){
1285  auto &png = sr->vtx.elastic.fuzzyk.png[iPng];
1286  if(png.truth.pdg != pdg) continue;
1287 
1288  double score = png.cvnpart.electronid;
1289 
1290  ret.push_back(score);
1291  }//end loop over prongs
1292 
1293  return ret;
1294  });
1295  return kelectroncvn;
1296  }
1297 
1298 
1300  {
1301  // returns proton CVN score of all prongs which match pdg
1302  const MultiVar kprotoncvn([pdg](const caf::SRProxy* sr)
1303  {
1304  std::vector<double> ret;
1305 
1306  if(sr->vtx.elastic.IsValid == false) return ret;
1307 
1308  int nPng = sr->vtx.elastic.fuzzyk.npng;
1309 
1310  for(int iPng=0; iPng < nPng; iPng++){
1311  auto &png = sr->vtx.elastic.fuzzyk.png[iPng];
1312  if(png.truth.pdg != pdg) continue;
1313 
1314  double score = png.cvnpart.protonid;
1315 
1316  ret.push_back(score);
1317  }//end loop over prongs
1318 
1319  return ret;
1320  });
1321  return kprotoncvn;
1322  }
1323 
1324 
1326  {
1327  // returns pion CVN score of all prongs which match pdg
1328  const MultiVar kpioncvn([pdg](const caf::SRProxy* sr)
1329  {
1330  std::vector<double> ret;
1331 
1332  if(sr->vtx.elastic.IsValid == false) return ret;
1333 
1334  int nPng = sr->vtx.elastic.fuzzyk.npng;
1335 
1336  for(int iPng=0; iPng < nPng; iPng++){
1337  auto &png = sr->vtx.elastic.fuzzyk.png[iPng];
1338  if(png.truth.pdg != pdg) continue;
1339 
1340  double score = png.cvnpart.pionid;
1341 
1342  ret.push_back(score);
1343  }//end loop over prongs
1344 
1345  return ret;
1346  });
1347  return kpioncvn;
1348  }
1349 
1350 
1352  {
1353  // returns photon CVN score of all prongs which match pdg
1354  const MultiVar kphotoncvn([pdg](const caf::SRProxy* sr)
1355  {
1356  std::vector<double> ret;
1357 
1358  if(sr->vtx.elastic.IsValid == false) return ret;
1359 
1360  int nPng = sr->vtx.elastic.fuzzyk.npng;
1361 
1362  for(int iPng=0; iPng < nPng; iPng++){
1363  auto &png = sr->vtx.elastic.fuzzyk.png[iPng];
1364  if(png.truth.pdg != pdg) continue;
1365 
1366  double score = png.cvnpart.photonid;
1367 
1368  ret.push_back(score);
1369  }//end loop over prongs
1370 
1371  return ret;
1372  });
1373  return kphotoncvn;
1374  }
1375 
1376 
1377 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~`
1378 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~`
1379 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~`
1380 
1381 
1382 
1383 
1384  const Var kNGoodPionPngs([](const caf::SRProxy* sr)
1385  {
1386  // Returns number of pion prongs other than muon
1387  // with > 6 hits, and pion score > 0.8
1388  //
1389  // (Only prongs )
1390  // NB using 4-view CVN for prod5, as it is all that's available for now
1391 
1392  if(sr->vtx.elastic.IsValid == false) return -5;
1393 
1394  int nPng = sr->vtx.elastic.fuzzyk.npng;
1395 
1396  int nPions = 0;
1397 
1398  for(int png_idx = 0; png_idx < nPng; ++png_idx){
1399 
1400  if(png_idx == kBestMuonPng(sr)) continue;
1401 
1402  auto &png = sr->vtx.elastic.fuzzyk.png[png_idx];
1403 
1404  if(png.nhit <= 6) continue;
1405 
1406 
1407  double cvn = png.cvnpart.pionid;
1408  if(cvn > 0.8 ){
1409  nPions++;
1410  }
1411 
1412  } // end for png
1413 
1414 
1415  //std::cout << "kNGoodPions = " << nPions << std::endl;
1416  return nPions;
1417  });
1418 
1419 
1420 
1421 
1422  const Var kNhit_3d(int pdg)
1423  {
1424  // returns nhits of most energetic 3d prong matching pdg
1425 
1426  const Var knhit_3d([pdg](const caf::SRProxy* sr)
1427  {
1428  double ret_value = -5.0;
1429  double max_energy = -5.0;
1430 
1431 
1432  if(sr->vtx.elastic.IsValid == false) return ret_value;
1433 
1434  int nPngs = sr->vtx.elastic.fuzzyk.npng;
1435  for(int png_idx=0; png_idx<nPngs; png_idx++){
1436 
1437  auto &png = sr->vtx.elastic.fuzzyk.png[png_idx];
1438  if( (abs(png.truth.pdg) == pdg) || (abs(png.truth.motherpdg) == pdg) ){
1439 
1440  if(png.truth.p.E > max_energy){
1441  max_energy = png.truth.p.E;
1442  ret_value = png.nhit;
1443  }
1444 
1445  }
1446 
1447  }// end for png_idx
1448  return ret_value;
1449  }); // end knhit_3d
1450 
1451  return knhit_3d;
1452  }
1453 
1454 
1455  const Var kNhit_2d(int pdg)
1456  {
1457  // returns nhits of most energetic 2d prong matching pdg
1458 
1459  const Var knhit_2d([pdg](const caf::SRProxy* sr)
1460  {
1461  double ret_value = -5.0;
1462  double max_energy = -5.0;
1463 
1464 
1465  if(sr->vtx.elastic.IsValid == false) return ret_value;
1466 
1467  int nPngs = sr->vtx.elastic.fuzzyk.npng2d;
1468  for(int png_idx=0; png_idx<nPngs; png_idx++){
1469 
1470  auto &png = sr->vtx.elastic.fuzzyk.png2d[png_idx];
1471  if(png.truth.pdg == pdg){
1472 
1473  if(png.truth.p.E > max_energy){
1474  max_energy = png.truth.p.E;
1475  ret_value = png.nhit;
1476  }
1477 
1478  }
1479 
1480  }// end for png_idx
1481  return ret_value;
1482  }); // end knhit_2d
1483 
1484  return knhit_2d;
1485  }
1486 
1487 
1488  const Var kTruePngKE(int pdg)
1489  {
1490  // returns True KE of most energetic 3d prong matching pdg
1491 
1492  const Var ktruePngE([pdg](const caf::SRProxy* sr)
1493  {
1494  double max_energy = -5.0;
1495 
1496 
1497  if(sr->vtx.elastic.IsValid == false) return max_energy;
1498 
1499  int nPngs = sr->vtx.elastic.fuzzyk.npng;
1500  for(int png_idx=0; png_idx<nPngs; png_idx++){
1501 
1502  auto &png = sr->vtx.elastic.fuzzyk.png[png_idx];
1503  if(abs(png.truth.pdg) == pdg || abs(png.truth.motherpdg) == pdg){
1504 
1505  if(png.truth.p.E > max_energy){
1506  max_energy = png.truth.p.E;
1507  }
1508 
1509  }
1510 
1511  }// end for png_idx
1512  max_energy -= pdg_to_mass.at(abs(pdg));
1513 
1514  return max_energy;
1515  }); // end ktruepngE
1516 
1517  return ktruePngE;
1518  }
1519 
1520 
1521 
1522  /*
1523  const Var kLongestPionLen([](const caf::SRProxy* sr)
1524  //Returns length of longest pion shower
1525  {
1526 
1527  double maxLen = -5.0;
1528 
1529  int nPngs = sr->vtx.elastic.fuzzyk.npng;
1530  for(int png_idx=0; png_idx<nPngs; png_idx++){
1531 
1532 
1533  auto &png = sr->vtx.elastic.fuzzyk.png[png_idx];
1534  if(png.truth.pdg != 211) continue;
1535 
1536  double len = png.shwlid.len;
1537 
1538  if(len > maxLen) maxLen = len;
1539 
1540  }// end for png_idx
1541  return maxLen;
1542  }); // end kLongestPionLen
1543  */
1544 
1545 
1546 
1547  /*
1548  const Var kBestPionCvnScore([](const caf::SRProxy* sr)
1549  // Returns cvn pion score of kBestPionPng
1550  {
1551  if(!kGoodPionPngExists(sr)) return -5.0;
1552 
1553  int idx = kBestPionPng(sr);
1554  auto &png = sr->vtx.elastic.fuzzyk.png[idx];
1555  return double(png.cvnpart2FlatFluxBal.pionid);
1556 
1557  });
1558 
1559 
1560 
1561  // TODO this doesn't account for motherpdg
1562  const MultiVar knhitByPdg(int pdg)
1563  // Returns values of nhit but only for prongs with true Pdgcode == pdg
1564  // (skips best muon prong)
1565  {
1566  const MultiVar knhit([pdg](const caf::SRProxy* sr)
1567  {
1568  std::vector<double> nhit_vec;
1569 
1570  if(sr->vtx.elastic.IsValid == false) return nhit_vec;
1571 
1572  int nPngs = sr->vtx.elastic.fuzzyk.npng;
1573  for(int png_idx=0; png_idx<nPngs; png_idx++){
1574 
1575  if(png_idx == kBestBPFTrack(sr)) continue;
1576 
1577  auto &png = sr->vtx.elastic.fuzzyk.png[png_idx];
1578  if(png.truth.pdg == pdg){
1579  nhit_vec.push_back(png.nhit);
1580  }
1581 
1582  }// end for png_idx
1583  return nhit_vec;
1584  }); // end knhit
1585 
1586  return knhit;
1587  }// end knhitByPdg
1588 
1589 
1590  const MultiVar kcalEByPdg(int pdg)
1591  // Returns values of calE but only for prongs with true Pdgcode == pdg
1592  // (skips best muon prong)
1593  {
1594  const MultiVar kcalE([pdg](const caf::SRProxy* sr)
1595  {
1596  std::vector<double> calE_vec;
1597 
1598  if(sr->vtx.elastic.IsValid == false) return calE_vec;
1599 
1600  int nPngs = sr->vtx.elastic.fuzzyk.npng;
1601  for(int png_idx=0; png_idx<nPngs; png_idx++){
1602 
1603  if(png_idx == kBestBPFTrack(sr)) continue;
1604 
1605  auto &png = sr->vtx.elastic.fuzzyk.png[png_idx];
1606  if(png.truth.pdg == pdg){
1607  calE_vec.push_back(png.calE);
1608  }
1609 
1610  }// end for png_idx
1611  return calE_vec;
1612  }); // end kcalE
1613 
1614  return kcalE;
1615  }// end kcalEByPdg
1616 
1617 
1618  const MultiVar klenByPdg(int pdg)
1619  // Returns values of len but only for prongs with true Pdgcode == pdg
1620  // (skips best muon prong)
1621  {
1622  const MultiVar klen([pdg](const caf::SRProxy* sr)
1623  {
1624  std::vector<double> len_vec;
1625 
1626  if(sr->vtx.elastic.IsValid == false) return len_vec;
1627 
1628  int nPngs = sr->vtx.elastic.fuzzyk.npng;
1629  for(int png_idx=0; png_idx<nPngs; png_idx++){
1630 
1631  if(png_idx == kBestBPFTrack(sr)) continue;
1632 
1633  auto &png = sr->vtx.elastic.fuzzyk.png[png_idx];
1634  if(png.truth.pdg == pdg){
1635  len_vec.push_back(png.len);
1636  }
1637 
1638  }// end for png_idx
1639  return len_vec;
1640  }); // end klen
1641 
1642  return klen;
1643  }// end klenByPdg
1644 
1645 
1646 
1647  const MultiVar ksimple_dEdxByPdg(int pdg)
1648  // Returns values of simple_dEdx but only for prongs with true Pdgcode == pdg
1649  // (skips best muon prong)
1650  //
1651  // simple_dEdx = calE / len
1652  {
1653  const MultiVar ksimple_dEdx([pdg](const caf::SRProxy* sr)
1654  {
1655  std::vector<double> simple_dEdx_vec;
1656 
1657  if(sr->vtx.elastic.IsValid == false) return simple_dEdx_vec;
1658 
1659  int nPngs = sr->vtx.elastic.fuzzyk.npng;
1660  for(int png_idx=0; png_idx<nPngs; png_idx++){
1661 
1662  if(png_idx == kBestBPFTrack(sr)) continue;
1663 
1664  auto &png = sr->vtx.elastic.fuzzyk.png[png_idx];
1665  if(png.truth.pdg == pdg){
1666 
1667  double simple_dEdx;
1668  if(png.len == 0.0) simple_dEdx = -5.0;
1669  // *100 converts to GeV/m
1670  else simple_dEdx = 100*(png.calE / png.len);
1671  simple_dEdx_vec.push_back(simple_dEdx);
1672  }
1673 
1674  }// end for png_idx
1675  return simple_dEdx_vec;
1676  }); // end ksimple_dEdx
1677 
1678  return ksimple_dEdx;
1679  }// end ksimple_dEdxByPdg
1680 
1681 
1682 
1683  const Var kVtxX([](const caf::SRProxy* sr){
1684 
1685  if(sr->vtx.elastic.IsValid == false ) return -5.;
1686  double vtx_x = sr->vtx.elastic.vtx.x;
1687  return vtx_x;
1688  });
1689 
1690  const Var kVtxY([](const caf::SRProxy* sr){
1691 
1692  if(sr->vtx.elastic.IsValid == false ) return -5.;
1693  double vtx_y = sr->vtx.elastic.vtx.y;
1694  return vtx_y;
1695  });
1696 
1697  const Var kVtxZ([](const caf::SRProxy* sr){
1698 
1699  if(sr->vtx.elastic.IsValid == false ) return -5.;
1700  double vtx_z = sr->vtx.elastic.vtx.z;
1701  return vtx_z;
1702 
1703  });
1704  */
1705 
1706 
1707  const MultiVar kTestMultiVar_A([](const caf::SRProxy* sr)
1708  {
1709  std::vector<double> ret_vec;
1710 
1711  if(sr->vtx.elastic.IsValid == false) return ret_vec;
1712 
1713  int nPngs = sr->vtx.elastic.fuzzyk.npng;
1714 
1715  std::cout << "nPngs: " << nPngs << std::endl;
1716 
1717  for(int ipng=0; ipng<nPngs; ipng++){
1718 
1719  std::cout << (double)ipng << std::endl;
1720 
1721  ret_vec.push_back((double)ipng);
1722  }// end for ipng
1723  return ret_vec;
1724 
1725  }); // end kTestMultiVar_A
1726 
1727 
1728  const MultiVar kProngEnergyRes([](const caf::SRProxy* sr)
1729  {
1730  std::vector<double> ERes_vec;
1731 
1732  double E_true;
1733  double E_reco;
1734  double E_res;
1735 
1736  if(sr->vtx.elastic.IsValid == false) return ERes_vec;
1737 
1738  int nPngs = sr->vtx.elastic.fuzzyk.npng;
1739 
1740  for(int ipng=0; ipng<nPngs; ipng++){
1741 
1742  auto &png = sr->vtx.elastic.fuzzyk.png[ipng];
1743 
1744  E_true = png.truth.p.E;
1745  E_reco = png.calE;
1746 
1747  E_res = (E_reco - E_true)/E_true;
1748  //E_res = (E_true - E_reco)/E_reco;
1749 
1750  ERes_vec.push_back(E_res);
1751 
1752  }// end for ipng
1753  return ERes_vec;
1754  }); // end kProngEnergyRes
1755 
1756 
1757 
1758  //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1759  // Vars used in G4Rwgt example macro
1760  const Var kLongestPionLen([](const caf::SRProxy* sr)
1761  //Returns length of longest pion shower
1762  {
1763 
1764  double maxLen = -5.0;
1765 
1766  int nPngs = sr->vtx.elastic.fuzzyk.npng;
1767  for(int png_idx=0; png_idx<nPngs; png_idx++){
1768 
1769 
1770  auto &png = sr->vtx.elastic.fuzzyk.png[png_idx];
1771  if(png.truth.pdg != 211) continue;
1772 
1773  double len = png.shwlid.len;
1774 
1775  if(len > maxLen) maxLen = len;
1776 
1777  }// end for png_idx
1778  return maxLen;
1779  }); // end kLongestPionLen
1780 
1781 
1782 
1783  Var kCalE([](const caf::SRProxy* sr){
1784 
1785  float calE = sr->slc.calE;
1786  return calE;
1787  });
1788 
1789 
1790 
1791  NuTruthVar kNuE_NT([](const caf::SRNeutrinoProxy* nu){
1792 
1793  float nuE = -5.;
1794  nuE = nu->E;
1795  if(nuE < 0) std::cout << "BAR \n";
1796 
1797  return nuE;
1798  });
1799 
1801 
1802  NuTruthVar kPiE_NT([](const caf::SRNeutrinoProxy* nu)
1803  {
1804 
1805  double maxE = -5.0;
1806 
1807  int nPrim = nu->prim.size();
1808  for(int i=0; i<nPrim; i++){
1809  auto &prim = nu->prim[i];
1810  if(prim.pdg != 211) continue;
1811  double energy = prim.p.E;
1812  if(energy > maxE ) maxE = energy;
1813  }
1814 
1815  if(maxE < 0) std::cout << "FOO \n";
1816 
1817  return maxE;
1818  });
1819 
1820  //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1821 
1822 }
const XML_Char int len
Definition: expat.h:262
caf::Proxy< size_t > npng
Definition: SRProxy.h:2038
const MultiVar kKalmanMinMax_byPdg(int pdg, std::string coord, std::string minmax)
T max(const caf::Proxy< T > &a, T b)
Var kNuE
const MultiVar k4viewElectronCVN_byPDG(int pdg)
caf::Proxy< caf::SRFuzzyK > fuzzyk
Definition: SRProxy.h:2059
_Var< caf::SRNeutrinoProxy > NuTruthVar
Definition: Var.h:9
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const MultiVar kProngEnergyRes_byPdg(int pdg, bool massSub=false, double minGev=0.0, double maxGev=100.0)
const MultiVar kSPProtonCVN_byPDG(int pdg)
const Var kProngPur_byPdg(int pdg)
const Var kTruePionEnergy
const MultiVar kShowerMinMax_byPdg(int pdg, std::string coord, std::string minmax="mix")
caf::Proxy< std::vector< caf::SRProng > > png2d
Definition: SRProxy.h:2044
const MultiVar k4viewProtonCVN_byPDG(int pdg)
const MultiVar kKalmanEnergyRes_byPdg(int pdg)
caf::Proxy< size_t > ntracks
Definition: SRProxy.h:1778
Proxy for caf::SRNeutrino.
Definition: SRProxy.h:510
NuTruthVar kNuE_NT([](const caf::SRNeutrinoProxy *nu){float nuE=-5.;nuE=nu->E;if(nuE< 0) std::cout<< "BAR \n";return nuE;})
caf::Proxy< caf::SRHeader > hdr
Definition: SRProxy.h:2137
const Var kTrueMuonProng([](const caf::SRProxy *sr) {if(sr->vtx.elastic.IsValid==false) return-5;int itrue_muon_png=-5;int nPng=sr->vtx.elastic.fuzzyk.npng;for(int png_idx=0;png_idx< nPng;png_idx++){auto &bpf=sr->vtx.elastic.fuzzyk.png[png_idx].bpf; int pdg;int motherpdg;if(bpf.muon.IsValid==true){pdg=bpf.muon.truth.pdg;motherpdg=bpf.muon.truth.motherpdg;}else if(bpf.pion.IsValid==true){pdg=bpf.pion.truth.pdg;motherpdg=bpf.pion.truth.motherpdg;}else if(bpf.proton.IsValid==true){pdg=bpf.proton.truth.pdg;motherpdg=bpf.proton.truth.motherpdg;}else continue;if(pdg==13 &&motherpdg==13){if(itrue_muon_png!=-5){std::cout<< "More than one muon prong! There is a problem here..."<< std::endl;return-5;}itrue_muon_png=png_idx;};}return itrue_muon_png;})
const Var kN3dPngs([](const caf::SRProxy *sr){if(sr->vtx.elastic.IsValid==false) return-5;int nPngs=sr->vtx.elastic.fuzzyk.png.size();return nPngs;})
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
const MultiVar k4viewPhotonCVN_byPDG(int pdg)
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:618
const Var kHighestPionCVN([](const caf::SRProxy *sr){int piIdx=kBestPionPng(sr);if(piIdx< 0) return-5.;double ret=sr->vtx.elastic.fuzzyk.png[piIdx].spprongcvnpart5label.pionid;return ret;})
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.spprongcvnpart5label.muonid;if(cvn > maxCVN){maxCVN=cvn;maxCVN_idx=png_idx;}}if(maxLen > 500) return maxLen_idx;return maxCVN_idx;})
void abs(TH1 *hist)
caf::Proxy< short int > nnu
Definition: SRProxy.h:617
const Var kTrueProtonEnergy
caf::Proxy< float > E
Definition: SRProxy.h:520
const MultiVar k4viewPionCVN_byPDG(int pdg)
const Var kHighestMuonCVN([](const caf::SRProxy *sr){ std::cout<< "FOO \n";std::cout<< "BAR \n";if(sr->vtx.elastic.IsValid==false) return-5.0;int nPng=sr->vtx.elastic.fuzzyk.npng;double maxCVN=-5.0;for(int png_idx=0;png_idx< nPng;++png_idx){auto &png=sr->vtx.elastic.fuzzyk.png[png_idx];double len=png.len;std::cout<< "len: "<< len<< "\n";if(len > 500.0) return 1.0;double cvn=png.spprongcvnpart5label.muonid;if(cvn > maxCVN){maxCVN=cvn;}}std::cout<< "maxCVN: "<< maxCVN<< "\n";return maxCVN;})
Defines an enumeration for prong classification.
const MultiVar kRun_byPdg(int pdg)
const Var kNPions
caf::Proxy< caf::SRElastic > elastic
Definition: SRProxy.h:2118
const MultiVar kShowerStop_byPdg(int pdg, std::string coord)
const Var kLongestPionLen([](const caf::SRProxy *sr){double maxLen=-5.0;int nPngs=sr->vtx.elastic.fuzzyk.npng;for(int png_idx=0;png_idx< nPngs;png_idx++){auto &png=sr->vtx.elastic.fuzzyk.png[png_idx];if(png.truth.pdg!=211) continue;double len=png.shwlid.len;if(len > maxLen) maxLen=len;}return maxLen;})
Track finder for cosmic rays.
caf::Proxy< std::vector< caf::SRFuzzyKProng > > png
Definition: SRProxy.h:2043
const MultiVar k4viewMuonCVN_byPDG(int pdg)
const Var kHighestPionCVN_KE([](const caf::SRProxy *sr){int piIdx=kBestPionPng(sr);if(piIdx< 0) return-5.;double ret=sr->vtx.elastic.fuzzyk.png[piIdx].truth.p.E-0.139;return ret;})
caf::Proxy< caf::SRTrackBranch > trk
Definition: SRProxy.h:2145
const Var kNhit_2d(int pdg)
const MultiVar kSPPionCVN_byPDG(int pdg)
const MultiVar kSPElectronCVN_byPDG(int pdg)
double energy
Definition: plottest35.C:25
const Var kProngEff_byPdg(int pdg)
caf::StandardRecord * sr
const std::map< int, double > pdg_to_mass
Definition: NumuCC1PiVars.h:90
OStream cout
Definition: OStream.cxx:6
const Var kProngEnergyRes_byPdg_nonMulti(int pdg)
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2138
const NuTruthVar kNProtons_NT([](const caf::SRNeutrinoProxy *nu){int nPrims=nu->prim.size();int nProtons=0;for(int prim_idx=0;prim_idx< nPrims;prim_idx++){auto &prim=nu->prim[prim_idx];int pdg=prim.pdg;if(abs(pdg)==2212){nProtons++;}}return nProtons;})
const NuTruthVar kNPions_NT([](const caf::SRNeutrinoProxy *nu){int nPrims=nu->prim.size();int nPions=0;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++;}}return nPions;})
const Var kShowerMinMax_byPdg_nonMulti(int pdg, std::string coord, std::string minmax="mix")
const MultiVar kTestMultiVar_A([](const caf::SRProxy *sr){std::vector< double > ret_vec;if(sr->vtx.elastic.IsValid==false) return ret_vec;int nPngs=sr->vtx.elastic.fuzzyk.npng;std::cout<< "nPngs: "<< nPngs<< std::endl;for(int ipng=0;ipng< nPngs;ipng++){std::cout<< (double) ipng<< std::endl;ret_vec.push_back((double) ipng);}return ret_vec;})
caf::Proxy< caf::SRKalman > kalman
Definition: SRProxy.h:1797
Var VarFromNuTruthVar(const NuTruthVar &stv, double _default)
Definition: Var.cxx:8
caf::Proxy< bool > IsValid
Definition: SRProxy.h:2058
caf::Proxy< caf::SRSlice > slc
Definition: SRProxy.h:2142
const MultiVar kSPPhotonCVN_byPDG(int pdg)
caf::Proxy< float > calE
Definition: SRProxy.h:1292
const Var kN2dPngs([](const caf::SRProxy *sr){if(sr->vtx.elastic.IsValid==false) return-5;int nPngs2d=sr->vtx.elastic.fuzzyk.png2d.size();return nPngs2d;})
const NuTruthVar kTrueProtonEnergy_NT([](const caf::SRNeutrinoProxy *nu){double ret_value=-5;int nPrims=nu->prim.size();double energy=-5.0;for(int prim_idx=0;prim_idx< nPrims;prim_idx++){auto &prim=nu->prim[prim_idx];int pdg=prim.pdg;if(abs(pdg)==2212){energy=prim.p.E;if(energy > ret_value) ret_value=energy;}}ret_value-=0.938;return ret_value;})
Var kCalE([](const caf::SRProxy *sr){float calE=sr->slc.calE;return calE;})
const Var kNProtons
const MultiVar kProngEnergyRes([](const caf::SRProxy *sr){std::vector< double > ERes_vec;double E_true;double E_reco;double E_res;if(sr->vtx.elastic.IsValid==false) return ERes_vec;int nPngs=sr->vtx.elastic.fuzzyk.npng;for(int ipng=0;ipng< nPngs;ipng++){auto &png=sr->vtx.elastic.fuzzyk.png[ipng];E_true=png.truth.p.E;E_reco=png.calE;E_res=(E_reco-E_true)/E_true;ERes_vec.push_back(E_res);}return ERes_vec;})
assert(nhit_max >=nhit_nbins)
const Var kTrueE_nu([](const caf::SRProxy *sr){if(sr->mc.nnu< 1) return-5.f;return float(sr->mc.nu[0].E);})
TVector3 NuMIBeamDirection(caf::Det_t det)
Average direction of NuMI neutrinos in a given detector This function is a copy of geo::GeometryBase:...
Definition: Utilities.cxx:328
const MultiVar kSPMuonCVN_byPDG(int pdg)
caf::Proxy< std::vector< caf::SRKalmanTrack > > tracks
Definition: SRProxy.h:1780
const Var kBestPionPng([](const caf::SRProxy *sr){if(sr->vtx.elastic.IsValid==false) return-5;int nPng=sr->vtx.elastic.fuzzyk.npng;double maxCVN=-5;int maxCVN_idx=-5;for(int png_idx=0;png_idx< nPng;++png_idx){if(png_idx==kBestMuonPng(sr)) continue;auto &png=sr->vtx.elastic.fuzzyk.png[png_idx];double cvn=png.spprongcvnpart5label.pionid;if(cvn > maxCVN){maxCVN=cvn;maxCVN_idx=png_idx;}}return maxCVN_idx;})
caf::Proxy< size_t > npng2d
Definition: SRProxy.h:2039
T min(const caf::Proxy< T > &a, T b)
caf::Proxy< caf::SRVertexBranch > vtx
Definition: SRProxy.h:2146
const Var kHighestMuonCVN_byPDG(int pdg)
const Var kHighestPionCVN_Cos([](const caf::SRProxy *sr){int piIdx=kBestPionPng(sr);if(piIdx< 0) return-5.;auto &png=sr->vtx.elastic.fuzzyk.png[piIdx];TVector3 beamDir=NuMIBeamDirection(sr->hdr.det);double ret=(beamDir.X()*png.dir.x)+(beamDir.Y()*png.dir.y)+(beamDir.Z()*png.dir.z);return ret;})
double maxE
Definition: plot_hist.C:8
Template for Vars applied to any type of object.
Definition: FwdDeclare.h:12
const Var kNGoodPionPngs([](const caf::SRProxy *sr){ if(sr->vtx.elastic.IsValid==false) return-5;int nPng=sr->vtx.elastic.fuzzyk.npng;int nPions=0;for(int png_idx=0;png_idx< nPng;++png_idx){if(png_idx==kBestMuonPng(sr)) continue;auto &png=sr->vtx.elastic.fuzzyk.png[png_idx];if(png.nhit<=6) continue;double cvn=png.cvnpart.pionid;if(cvn > 0.8){nPions++;}} return nPions;})
NuTruthVar kPiE_NT([](const caf::SRNeutrinoProxy *nu){double maxE=-5.0;int nPrim=nu->prim.size();for(int i=0;i< nPrim;i++){auto &prim=nu->prim[i];if(prim.pdg!=211) continue;double energy=prim.p.E;if(energy > maxE) maxE=energy;}if(maxE< 0) std::cout<< "FOO \n";return maxE;})
const Var kNhit_3d(int pdg)
const MultiVar kShowerStart_byPdg(int pdg, std::string coord)
const Var kTruePngKE(int pdg)
caf::Proxy< std::vector< caf::SRTrueParticle > > prim
Definition: SRProxy.h:555
caf::Proxy< caf::Det_t > det
Definition: SRProxy.h:232
const NuTruthVar kTruePionEnergy_NT([](const caf::SRNeutrinoProxy *nu){ int nPrims=nu->prim.size();int nPions=0;double energy=-5.0;for(int prim_idx=0;prim_idx< nPrims;prim_idx++){auto &prim=nu->prim[prim_idx];int pdg=prim.pdg;if(abs(pdg)==211){energy=prim.p.E;energy-=0.139;nPions++;}}if(nPions!=1){std::cout<< nPions<< "\n";std::cout<< "FOO \n";abort();}return energy;})
enum BeamMode string