TrueProngVars.cxx
Go to the documentation of this file.
2 
4 
8 
10 
11 namespace ana
12 {
13 
14  // ==================== //
15  //
16  // Helper functions and variables for choosing prongs based on truth, for training the numuCC BPF energy estimator:
17  //
18  // ==================== //
19 
20  bool ProngTruthCut(unsigned int pIdx, const caf::SRProxy* sr,
21  std::vector<int> pdgCodes,
22  double eff, double pur) {
23 
24  bool pass = false;
25 
26  // Sanity checks first to prevent segfaults...
27  if(!sr->vtx.elastic.IsValid) return pass;
28  if(sr->vtx.elastic.fuzzyk.npng == 0) return pass;
29  if(pIdx >= sr->vtx.elastic.fuzzyk.npng) return pass;
30 
31  for(unsigned int i = 0; i < pdgCodes.size(); ++i) {
32  if(abs(sr->vtx.elastic.fuzzyk.png[pIdx].truth.pdg) == pdgCodes[i] &&
33  sr->vtx.elastic.fuzzyk.png[pIdx].truth.eff >= eff &&
34  sr->vtx.elastic.fuzzyk.png[pIdx].truth.pur >= pur) pass = true;
35  }
36 
37  return pass;
38 
39  }
40 
41  // to compare true vs. BPF track energy
42  bool hasBPFE(unsigned int pIdx, const caf::SRProxy* sr){
43  // Sanity checks first to prevent segfaults...
44  if(!sr->vtx.elastic.IsValid) return false;
45  if(sr->vtx.elastic.fuzzyk.npng == 0) return false;
46  if(pIdx >= sr->vtx.elastic.fuzzyk.npng) return false;
47 
48  auto &png = sr->vtx.elastic.fuzzyk.png[pIdx];
49 
50  if( abs(png.truth.pdg) == 13 && png.bpf.muon.IsValid && png.bpf.muon.energy >= 0.0) return true;
51  if( abs(png.truth.pdg) == 211 && png.bpf.pion.IsValid && png.bpf.pion.energy >= 0.0) return true;
52  if( abs(png.truth.pdg) == 2212 && png.bpf.proton.IsValid && png.bpf.proton.energy >= 0.0) return true;
53 
54  return false;
55  }
56 
57  //
58  // Define helper functions for computing energies:
59  //
60 
61 
62 // flagged for removal
63  /// getMuonEnergy
64  /// \brief: function to compute the muon energy
65  ///
66  double getMuonEnergyByTruth(TSpline3 spline_calE, TSpline3 spline_BPFE, const caf::SRProxy* sr) {
67 
68  double muE = 0.0;
69  int muIdx = -5;
70 
71  for(unsigned int p = 0; p < sr->vtx.elastic.fuzzyk.png.size(); ++p) {
72  if(ProngTruthCut(p, sr, {13}, 0.6, 0.6)) muIdx = p;
73  }
74 
75  if(muIdx==-5) return muE;
76 
77  // get the muon calE and BPFE
78  // As usual, here is a hardcoded assumption about there being only one vertex...
79  double muBPFE = -1.0;
80  double muCalE = sr->vtx.elastic.fuzzyk.png[muIdx].calE; // need to subtrack overlapE in the loop below
81 
82  // Assuming we've passed the quality cuts, take the value from the muon track.
83  if(sr->vtx.elastic.fuzzyk.png[muIdx].bpf.muon.IsValid) {
84  muBPFE = sr->vtx.elastic.fuzzyk.png[muIdx].bpf.muon.energy;
85  muCalE = muCalE - sr->vtx.elastic.fuzzyk.png[muIdx].bpf.muon.overlapE; // subtract overlapping E from the muon track (to be put in the garbage bin)
86  }
87 
88  // adjust both E values with spline fits
89  double newCalE = spline_calE.Eval(muCalE);
90  double newBPFE = spline_BPFE.Eval(muBPFE);
91 
92  // Finally, combine newBPFE and newCalE.
93  // If no BPF muon track existed, fall back on calE only.
94  if(muBPFE < 0.0) muE = newCalE;
95  else muE = kMuonFitWeight * newBPFE + (1.0 - kMuonFitWeight) * newCalE;
96 
97  return muE;
98 
99  } // end getMuonEnergy
100 
101 
102 
103  /// getEMShowerEnergy
104  /// \brief: function to compute the EM shower energy
105  ///
106  double getEMShowerEnergyByTruth(TSpline3 spline_calE, const caf::SRProxy* sr) {
107 
108  double E = 0.0;
109  double calE = 0.0;
110 
111  for(unsigned int p = 0; p < sr->vtx.elastic.fuzzyk.png.size(); ++p) {
112  if(ProngTruthCut(p, sr, {11, 22, 111}, 0.5, 0.5)) calE += sr->vtx.elastic.fuzzyk.png[p].weightedCalE;
113  } // end loop over 3D prongs
114 
115  // if there are no selected EM prongs, return 0.0
116  if(calE==0.0) return calE;
117 
118  // adjust calE with spline fits
119  E = spline_calE.Eval(calE);
120 
121  return E;
122 
123  } // end getEMShowerEnergy
124 
125 
126 
127  /// getHadronEnergy
128  /// \brief: function to compute the hadron energy
129  ///
130  double getHadronEnergyByTruth(TSpline3 spline_calE, const caf::SRProxy* sr) {
131 
132  double E = 0.0;
133  double calE = 0.0;
134 
135  // Again with the assumtion about only having one vertex (and below too...)
136  for(unsigned int p = 0; p < sr->vtx.elastic.fuzzyk.png.size(); ++p) {
137  if(ProngTruthCut(p, sr, {211}, 0.5, 0.5) || ProngTruthCut(p, sr, {2212}, 0.5, 0.5)) calE += sr->vtx.elastic.fuzzyk.png[p].weightedCalE;
138  } // end loop over 3D prongs
139 
140  // if there are no selected hadron prongs, return 0.0
141  if(calE==0.0) return calE;
142 
143  // adjust calE with spline fits
144  E = spline_calE.Eval(calE);
145 
146  return E;
147 
148  } // end getHadronEnergy
149 
150  /// takeOutTrash
151  /// \brief: function to compute v3 of the garbage bin.
152  /// this takes out muons, EM showers, and general hadrons, everything else in the garbage bin
153  ///
154  double takeOutTrashByTruth(TSpline3 spline_calE, const caf::SRProxy* sr) {
155 
156  double energy = 0.0;
157  double calE = 0.0;
158 
159  // loop over all 3D prongs, skip the muon or if EM shower, sum weighted calE
160  // Again with the assumption about only having one vertex
161  for(unsigned int p = 0; p < sr->vtx.elastic.fuzzyk.png.size(); ++p) {
162  if(ProngTruthCut(p, sr, {13}, 0.6, 0.6)) continue; // skip the muon
163  if(ProngTruthCut(p, sr, {11, 22, 111}, 0.5, 0.5)) continue; // skip if a good shower
164  if(ProngTruthCut(p, sr, {211}, 0.5, 0.5)) continue; // skip if a good pion
165  if(ProngTruthCut(p, sr, {2212}, 0.5, 0.5)) continue; // skip if a good proton
166 
167  calE += sr->vtx.elastic.fuzzyk.png[p].weightedCalE;
168  } // end loop over 3D prongs
169 
170  // loop over all 2D prongs, sum weighted calE
171  for(unsigned int p = 0; p < sr->vtx.elastic.fuzzyk.png2d.size(); ++p) {
172  if(sr->vtx.elastic.fuzzyk.png2d[p].weightedCalE>0.0) calE += sr->vtx.elastic.fuzzyk.png2d[p].weightedCalE;
173  } // end loop over 2D prongs
174 
175  // add in the orphaned hit calE
176  if(sr->vtx.elastic.fuzzyk.orphCalE>0.0) calE += sr->vtx.elastic.fuzzyk.orphCalE;
177 
178  // add in the overlapping energy from the muon track
179  int muIdx = -5;
180  for(unsigned int p = 0; p < sr->vtx.elastic.fuzzyk.png.size(); ++p) {
181  if(ProngTruthCut(p, sr, {13}, 0.6, 0.6)) muIdx = p;
182  }
183 
184  double overlapE = 0.0;
185  if(muIdx!=-5){
186  // Assuming we've passed the quality cuts, take the value from the muon track.
187  if(sr->vtx.elastic.fuzzyk.png[muIdx].bpf.muon.IsValid) {
188  overlapE = sr->vtx.elastic.fuzzyk.png[muIdx].bpf.muon.overlapE;
189  }
190  }
191  calE += overlapE;
192 
193  // convert calE with the spline
194  if(calE==0.0) return calE;
195  energy = spline_calE.Eval(calE);
196 
197  return energy;
198 
199  } // end takeOutTrash
200 
201 
202 
203 
204  const Var kTrueMuonTrueE([](const caf::SRProxy *sr)
205  {
206  double E = 0.0;
207  for(unsigned int pIdx = 0; pIdx < sr->vtx.elastic.fuzzyk.npng; pIdx++) {
208  if(ProngTruthCut(pIdx, sr, {13}, 0.6, 0.6)) {
209  E = sr->vtx.elastic.fuzzyk.png[pIdx].truth.p.E;
210  }
211  }
212 
213  // if there are no prongs that pass selection, return -5.0
214  // prevents inclusion in true vs. reco plots
215  // prevents segfaults in test Eres plots
216  if(E==0.0) return -5.0;
217  return E;
218 
219  });
220 
221  const MultiVar kTrueProtonTrueKE([](const caf::SRProxy *sr)
222  {
223  std::vector<double> E;
224  for(unsigned int pIdx = 0; pIdx < sr->vtx.elastic.fuzzyk.npng; pIdx++) {
225  if(ProngTruthCut(pIdx, sr, {2212}, 0.5, 0.5)) {
226  E.emplace_back(sr->vtx.elastic.fuzzyk.png[pIdx].truth.p.E-0.938272);
227  }
228  }
229  return E;
230  });
231 
232  const MultiVar kTruePionTrueE([](const caf::SRProxy *sr)
233  {
234  std::vector<double> E;
235  for(unsigned int pIdx = 0; pIdx < sr->vtx.elastic.fuzzyk.npng; pIdx++) {
236  if(ProngTruthCut(pIdx, sr, {211}, 0.5, 0.5)) {
237  E.emplace_back(sr->vtx.elastic.fuzzyk.png[pIdx].truth.p.E);
238  }
239  }
240  return E;
241  });
242 
243  const Var kTrueMuonwBPFTrueE([](const caf::SRProxy *sr)
244  {
245  double E = 0.0;
246  for(unsigned int pIdx = 0; pIdx < sr->vtx.elastic.fuzzyk.npng; pIdx++) {
247  if(ProngTruthCut(pIdx, sr, {13}, 0.6, 0.6) && hasBPFE(pIdx, sr)){
248  E = sr->vtx.elastic.fuzzyk.png[pIdx].truth.p.E;
249  }
250  }
251 
252  // if there are no prongs that pass selection, return -5.0
253  // prevents inclusion in true vs. reco plots
254  // prevents segfaults in test Eres plots
255  if(E==0.0) return -5.0;
256  return E;
257  });
258 
259  const MultiVar kTrueProtonwBPFTrueKE([](const caf::SRProxy *sr)
260  {
261  std::vector<double> E;
262  for(unsigned int pIdx = 0; pIdx < sr->vtx.elastic.fuzzyk.npng; pIdx++) {
263  if(ProngTruthCut(pIdx, sr, {2212}, 0.5, 0.5) && hasBPFE(pIdx, sr))
264  E.emplace_back(sr->vtx.elastic.fuzzyk.png[pIdx].truth.p.E-0.938272);
265  }
266  return E;
267  });
268 
269  const MultiVar kTruePionwBPFTrueE([](const caf::SRProxy *sr)
270  {
271  std::vector<double> E;
272  for(unsigned int pIdx = 0; pIdx < sr->vtx.elastic.fuzzyk.npng; pIdx++) {
273  if(ProngTruthCut(pIdx, sr, {211}, 0.5, 0.5) && hasBPFE(pIdx, sr)) {
274  E.emplace_back(sr->vtx.elastic.fuzzyk.png[pIdx].truth.p.E);
275  }
276  }
277  return E;
278  });
279 
280  const MultiVar kTrueEMTrueE([](const caf::SRProxy *sr)
281  {
282  std::vector<double> E;
283  for(unsigned int pIdx = 0; pIdx < sr->vtx.elastic.fuzzyk.npng; pIdx++) {
284  if(ProngTruthCut(pIdx, sr, {11, 22, 111}, 0.5, 0.5)) {
285  E.emplace_back(sr->vtx.elastic.fuzzyk.png[pIdx].truth.p.E);
286  }
287  }
288  return E;
289  });
290 
291  const MultiVar kTrueHadTrueE([](const caf::SRProxy *sr)
292  {
293  std::vector<double> E;
294  for(unsigned int pIdx = 0; pIdx < sr->vtx.elastic.fuzzyk.npng; pIdx++) {
295  if(ProngTruthCut(pIdx, sr, {2212}, 0.5, 0.5)) {
296  E.emplace_back(sr->vtx.elastic.fuzzyk.png[pIdx].truth.p.E-0.938272);
297  }
298  else if(ProngTruthCut(pIdx, sr, {211}, 0.5, 0.5)){
299  E.emplace_back(sr->vtx.elastic.fuzzyk.png[pIdx].truth.p.E);
300  }
301  }
302  return E;
303  });
304 
305  const Var kTrueMuonCalE([](const caf::SRProxy *sr){
306  double E = 0.0;
307  for(unsigned int pIdx = 0; pIdx < sr->vtx.elastic.fuzzyk.npng; pIdx++) {
308  if(ProngTruthCut(pIdx, sr, {13}, 0.6, 0.6)) {
309  E = sr->vtx.elastic.fuzzyk.png[pIdx].weightedCalE;
310  }
311  }
312  if(E==0.0) return -5.0;
313  return E;
314  });
315 
316  const MultiVar kTrueProtonCalE([](const caf::SRProxy *sr){
317  std::vector<double> E;
318  for(unsigned int pIdx = 0; pIdx < sr->vtx.elastic.fuzzyk.npng; pIdx++) {
319  if(ProngTruthCut(pIdx, sr, {2212}, 0.5, 0.5)) {
320  E.emplace_back(sr->vtx.elastic.fuzzyk.png[pIdx].weightedCalE);
321  }
322  }
323  return E;
324  });
325 
326  const MultiVar kTruePionCalE([](const caf::SRProxy *sr){
327  std::vector<double> E;
328  for(unsigned int pIdx = 0; pIdx < sr->vtx.elastic.fuzzyk.npng; pIdx++) {
329  if(ProngTruthCut(pIdx, sr, {211}, 0.5, 0.5)) {
330  E.emplace_back(sr->vtx.elastic.fuzzyk.png[pIdx].weightedCalE);
331  }
332  }
333  return E;
334  });
335 
336  const MultiVar kTrueEMCalE([](const caf::SRProxy *sr){
337  std::vector<double> E;
338  for(unsigned int pIdx = 0; pIdx < sr->vtx.elastic.fuzzyk.npng; pIdx++) {
339  if(ProngTruthCut(pIdx, sr, {11, 22, 111}, 0.5, 0.5)){
340  E.emplace_back(sr->vtx.elastic.fuzzyk.png[pIdx].weightedCalE);
341  }
342  }
343  return E;
344  });
345 
346  const MultiVar kTrueHadCalE([](const caf::SRProxy *sr){
347  std::vector<double> E;
348  for(unsigned int pIdx = 0; pIdx < sr->vtx.elastic.fuzzyk.npng; pIdx++) {
349  if(ProngTruthCut(pIdx, sr, {2212, 211}, 0.5, 0.5)) {
350  E.emplace_back(sr->vtx.elastic.fuzzyk.png[pIdx].weightedCalE);
351  }
352  }
353  return E;
354  });
355 
356  const Var kTrueMuonBPFE([](const caf::SRProxy *sr){
357  double E = 0.0;
358  for(unsigned int pIdx = 0; pIdx < sr->vtx.elastic.fuzzyk.npng; pIdx++) {
359  if(!ProngTruthCut(pIdx, sr, {13}, 0.6, 0.6) || !hasBPFE(pIdx,sr)) continue;
360  auto &png = sr->vtx.elastic.fuzzyk.png[pIdx];
361  if(png.bpf.muon.IsValid) {
362  E = png.bpf.muon.energy;
363  }
364  }
365 
366  // if there are no prongs that pass selection, return -5.0
367  // prevents inclusion in true vs. reco plots
368  if(E==0.0) return -5.0;
369  return E;
370  });
371 
372  const MultiVar kTrueProtonBPFE([](const caf::SRProxy *sr){
373  std::vector<double> E;
374  for(unsigned int pIdx = 0; pIdx < sr->vtx.elastic.fuzzyk.npng; pIdx++) {
375  if(!ProngTruthCut(pIdx, sr, {2212}, 0.5, 0.5) || !hasBPFE(pIdx,sr)) continue;
376  auto &png = sr->vtx.elastic.fuzzyk.png[pIdx];
377  if(png.bpf.proton.IsValid) {
378  E.emplace_back(png.bpf.proton.energy-0.938272);
379  }
380  }
381  return E;
382  });
383 
384  const MultiVar kTruePionBPFE([](const caf::SRProxy *sr){
385  std::vector<double> E;
386  for(unsigned int pIdx = 0; pIdx < sr->vtx.elastic.fuzzyk.npng; pIdx++) {
387  if(!ProngTruthCut(pIdx, sr, {211}, 0.5, 0.5) || !hasBPFE(pIdx,sr)) continue;
388  auto &png = sr->vtx.elastic.fuzzyk.png[pIdx];
389  if(png.bpf.pion.IsValid) {
390  E.emplace_back(png.bpf.pion.energy);
391  }
392  }
393  return E;
394  });
395 
396  const Var kTrueMuonwBPFCalE([](const caf::SRProxy *sr){
397  double E = 0.0;
398  for(unsigned int pIdx = 0; pIdx < sr->vtx.elastic.fuzzyk.npng; pIdx++) {
399  if(!ProngTruthCut(pIdx, sr, {13}, 0.6, 0.6) || !hasBPFE(pIdx,sr)) continue;
400  auto &png = sr->vtx.elastic.fuzzyk.png[pIdx];
401  if(png.bpf.muon.IsValid) {
402  E = png.calE - png.bpf.muon.overlapE;
403  }
404  }
405 
406  // if there are no prongs that pass selection, return -5.0
407  // prevents inclusion in true vs. reco plots
408  if(E==0.0) return -5.0;
409  return E;
410  });
411 
412  const MultiVar kTrueProtonwBPFCalE([](const caf::SRProxy *sr){
413  std::vector<double> E;
414  for(unsigned int pIdx = 0; pIdx < sr->vtx.elastic.fuzzyk.npng; pIdx++) {
415  if(!ProngTruthCut(pIdx, sr, {2212}, 0.5, 0.5) || !hasBPFE(pIdx,sr)) continue;
416  auto &png = sr->vtx.elastic.fuzzyk.png[pIdx];
417  if(png.bpf.proton.IsValid) {
418  E.emplace_back(png.weightedCalE);
419  }
420  }
421  return E;
422  });
423 
424 
425  const MultiVar kTruePionwBPFCalE([](const caf::SRProxy *sr){
426  std::vector<double> E;
427  for(unsigned int pIdx = 0; pIdx < sr->vtx.elastic.fuzzyk.npng; pIdx++) {
428  if(!ProngTruthCut(pIdx, sr, {211}, 0.5, 0.5) || !hasBPFE(pIdx,sr)) continue;
429  auto &png = sr->vtx.elastic.fuzzyk.png[pIdx];
430  if(png.bpf.pion.IsValid) {
431  E.emplace_back(png.weightedCalE);
432  }
433  }
434  return E;
435  });
436 
437 
438  // garbage bin -- everything but selected muon, selected EM, selected had
439  const Var kTrueOtherCalE([](const caf::SRProxy *sr)
440  {
441  if(!sr->vtx.elastic.IsValid) return -5.0;
442  double calE = 0;
443 
444  unsigned int muIdx = (unsigned int)kCVNMuonIdx(sr);
445 
446 // We are using the selection where all 3D prongs are categorized as either EM or hadronic. If we revisit this selection we must update this as well.
447 
448 // for(unsigned int pIdx = 0; pIdx < sr->vtx.elastic.fuzzyk.npng; pIdx++){
449 // if(ProngTruthCut(pIdx, sr, {13}, 0.6, 0.6) || ProngTruthCut(pIdx, sr, {2212}, 0.5, 0.5)
450 // || ProngTruthCut(pIdx, sr, {211}, 0.5, 0.5) || IsEMShowerByTruth(pIdx,sr)) continue;
451 // if(pIdx == muIdx) continue; // skip the muon
452 // if(IsNumuCCEMShowerByCVN(pIdx,sr)) continue; // skip if a good shower
453 // if(IsNumuCCPionByCVN(pIdx,sr)) continue; // skip if a good pion
454 // if(IsNumuCCProtonByCVN(pIdx,sr)) continue; // skip if a good proton
455 //
456 // calE += sr->vtx.elastic.fuzzyk.png[pIdx].weightedCalE;
457 // }
458 
459  for(unsigned int pIdx = 0; pIdx < sr->vtx.elastic.fuzzyk.png2d.size(); pIdx++){
460  auto &png = sr->vtx.elastic.fuzzyk.png2d[pIdx];
461  calE += png.weightedCalE;
462  }
463 
464  calE += sr->vtx.elastic.fuzzyk.orphCalE;
465 
466  double overlapE = 0.0;
467 
468  if(sr->vtx.elastic.fuzzyk.png[muIdx].bpf.muon.IsValid) {
469  overlapE = sr->vtx.elastic.fuzzyk.png[muIdx].bpf.muon.overlapE;
470  }
471  calE += overlapE;
472  return calE;
473  });
474 
475 
476  const Var kTrueGarbageE([](const caf::SRProxy *sr) {
477  if(!sr->vtx.elastic.IsValid) return -5.0;
478  double energy = 0.0;
479 
480  int run = sr->hdr.run;
481  caf::Det_t det = sr->hdr.det;
482 
483  // split by detector
484  if(det == caf::kNEARDET) {
485 
486  // get the period
487  unsigned int period = PeriodFromRunND(run);
488 
489  // get the muon energy
490  energy += getMuonEnergy(allNDMuonFits_calE[period], allNDMuonFits_BPFE[period], sr);
491 
492  // get the EM shower energy
493  energy += getEMShowerEnergy(allNDEMFits[period], sr);
494 
495  // get the hadron energy
496  energy += getHadronEnergy(allNDHadFits[period], sr);
497 
498  }
499 
500  if(det == caf::kFARDET) {
501 
502  // get the period
503  unsigned int period = PeriodFromRunFD(run);
504  // get the muon energy
505  energy += getMuonEnergy(allFDMuonFits_calE[period], allFDMuonFits_BPFE[period], sr);
506 
507  // get the EM shower energy
508  energy += getEMShowerEnergy(allFDEMFits[period], sr);
509 
510  // get the hadron energy
511  energy += getHadronEnergy(allFDHadFits[period], sr);
512 
513  }
514  return (sr->mc.nu[0].E - energy);
515 
516  });
517 
518 
519  const Var kBPFMuonEByTruth([](const caf::SRProxy *sr) {
520 
521  double muE = 0.0;
522 
523  int run = sr->hdr.run;
524  caf::Det_t det = sr->hdr.det;
525 
526  // split by detector
527  if(det == caf::kNEARDET) {
528  unsigned int period = PeriodFromRunND(run);
529  muE = getMuonEnergyByTruth(allNDMuonFits_calE[period], allNDMuonFits_BPFE[period], sr);
530  }
531 
532  if(det == caf::kFARDET) {
533  unsigned int period = PeriodFromRunFD(run);
534  muE = getMuonEnergyByTruth(allFDMuonFits_calE[period], allFDMuonFits_BPFE[period], sr);
535  }
536 
537  return muE;
538 
539  }); // end kBPFMuonE
540 
541 
542  const Var kEMEByTruth([](const caf::SRProxy *sr){
543 
544  double EME = 0.0;
545 
546  int run = sr->hdr.run;
547  caf::Det_t det = sr->hdr.det;
548 
549  // split by detector
550  if(det == caf::kNEARDET) {
551  unsigned int period = PeriodFromRunND(run);
552  EME = getEMShowerEnergyByTruth(allNDEMFits[period], sr);
553  }
554 
555  if(det == caf::kFARDET) {
556  unsigned int period = PeriodFromRunFD(run);
557  EME = getEMShowerEnergyByTruth(allFDEMFits[period], sr);
558  }
559 
560  return EME;
561 
562  });
563 
564 
565  const Var kHadEByTruth([](const caf::SRProxy *sr){
566 
567  double hadE = 0.0;
568 
569  int run = sr->hdr.run;
570  caf::Det_t det = sr->hdr.det;
571 
572  // split by detector
573  if(det == caf::kNEARDET) {
574  unsigned int period = PeriodFromRunND(run);
575  hadE = getHadronEnergyByTruth(allNDHadFits[period], sr);
576  }
577 
578  if(det == caf::kFARDET) {
579  unsigned int period = PeriodFromRunFD(run);
580  hadE = getHadronEnergyByTruth(allFDHadFits[period], sr);
581  }
582 
583  return hadE;
584 
585  });
586 
587 
588  const Var kGarbageEByTruth([](const caf::SRProxy *sr){
589  double garE = 0.0;
590 
591  int run = sr->hdr.run;
592  caf::Det_t det = sr->hdr.det;
593 
594  // split by detector
595  if(det == caf::kNEARDET) {
596  unsigned int period = PeriodFromRunND(run);
597  garE = takeOutTrashByTruth(allNDGarbageFits[period], sr);
598  }
599 
600  if(det == caf::kFARDET) {
601  unsigned int period = PeriodFromRunFD(run);
602  garE = takeOutTrashByTruth(allFDGarbageFits[period], sr);
603  }
604 
605  return garE;
606 
607  });
608 
609 
610  const Var kBPFSliceEByTruth([](const caf::SRProxy *sr) {
611 
612  double energy = 0.0;
613 
614  int run = sr->hdr.run;
615  caf::Det_t det = sr->hdr.det;
616 
617  // split by detector
618  if(det == caf::kNEARDET) {
619 
620  // get the period
621  unsigned int period = PeriodFromRunND(run);
622 
623  // get the muon energy
624  energy += getMuonEnergyByTruth(allNDMuonFits_calE[period], allNDMuonFits_BPFE[period], sr);
625 
626  // get the EM shower energy
627  energy += getEMShowerEnergyByTruth(allNDEMFits[period], sr);
628 
629  // get the hadron energy
630  energy += getHadronEnergyByTruth(allNDHadFits[period], sr);
631 
632  // get the garbage bin energy
633  energy += takeOutTrashByTruth(allNDGarbageFits[period], sr);
634 
635  }
636 
637  if(det == caf::kFARDET) {
638 
639  // get the period
640  unsigned int period = PeriodFromRunFD(run);
641 
642  // get the muon energy
643  energy += getMuonEnergyByTruth(allFDMuonFits_calE[period], allFDMuonFits_BPFE[period], sr);
644 
645  // get the EM shower energy
646  energy += getEMShowerEnergyByTruth(allFDEMFits[period], sr);
647 
648  // get the hadron energy
649  energy += getHadronEnergyByTruth(allFDHadFits[period], sr);
650 
651  // get the garbage bin energy
652  energy += takeOutTrashByTruth(allFDGarbageFits[period], sr);
653 
654  }
655 
656  return energy;
657 
658  }); // end kBPFSliceE
659 
660 
661 }
caf::Proxy< size_t > npng
Definition: SRProxy.h:2038
Near Detector underground.
Definition: SREnums.h:10
Det_t
Which NOvA detector?
Definition: SREnums.h:7
caf::Proxy< caf::SRFuzzyK > fuzzyk
Definition: SRProxy.h:2059
Far Detector at Ash River.
Definition: SREnums.h:11
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const Var kTrueMuonBPFE([](const caf::SRProxy *sr){double E=0.0;for(unsigned int pIdx=0;pIdx< sr->vtx.elastic.fuzzyk.npng;pIdx++){if(!ProngTruthCut(pIdx, sr,{13}, 0.6, 0.6)||!hasBPFE(pIdx, sr)) continue;auto &png=sr->vtx.elastic.fuzzyk.png[pIdx];if(png.bpf.muon.IsValid){E=png.bpf.muon.energy;}} if(E==0.0) return-5.0;return E;})
Definition: TrueProngVars.h:48
double getHadronEnergyByTruth(TSpline3 spline_calE, const caf::SRProxy *sr)
: function to compute the hadron energy
double takeOutTrashByTruth(TSpline3 spline_calE, const caf::SRProxy *sr)
: function to compute v3 of the garbage bin. this takes out muons, EM showers, and general hadrons...
caf::Proxy< std::vector< caf::SRProng > > png2d
Definition: SRProxy.h:2044
const std::vector< TSpline3 > allFDMuonFits_BPFE
caf::Proxy< caf::SRHeader > hdr
Definition: SRProxy.h:2137
const char * p
Definition: xmltok.h:285
const std::vector< TSpline3 > allNDHadFits
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
const MultiVar kTrueProtonCalE([](const caf::SRProxy *sr){std::vector< double > E;for(unsigned int pIdx=0;pIdx< sr->vtx.elastic.fuzzyk.npng;pIdx++){if(ProngTruthCut(pIdx, sr,{2212}, 0.5, 0.5)){E.emplace_back(sr->vtx.elastic.fuzzyk.png[pIdx].weightedCalE);}}return E;})
Definition: TrueProngVars.h:43
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:618
const MultiVar kTruePionCalE([](const caf::SRProxy *sr){std::vector< double > E;for(unsigned int pIdx=0;pIdx< sr->vtx.elastic.fuzzyk.npng;pIdx++){if(ProngTruthCut(pIdx, sr,{211}, 0.5, 0.5)){E.emplace_back(sr->vtx.elastic.fuzzyk.png[pIdx].weightedCalE);}}return E;})
Definition: TrueProngVars.h:44
const MultiVar kTrueProtonwBPFTrueKE([](const caf::SRProxy *sr){std::vector< double > E;for(unsigned int pIdx=0;pIdx< sr->vtx.elastic.fuzzyk.npng;pIdx++){if(ProngTruthCut(pIdx, sr,{2212}, 0.5, 0.5)&&hasBPFE(pIdx, sr)) E.emplace_back(sr->vtx.elastic.fuzzyk.png[pIdx].truth.p.E-0.938272);}return E;})
Definition: TrueProngVars.h:35
void abs(TH1 *hist)
const MultiVar kTrueProtonBPFE([](const caf::SRProxy *sr){std::vector< double > E;for(unsigned int pIdx=0;pIdx< sr->vtx.elastic.fuzzyk.npng;pIdx++){if(!ProngTruthCut(pIdx, sr,{2212}, 0.5, 0.5)||!hasBPFE(pIdx, sr)) continue;auto &png=sr->vtx.elastic.fuzzyk.png[pIdx];if(png.bpf.proton.IsValid){E.emplace_back(png.bpf.proton.energy-0.938272);}}return E;})
Definition: TrueProngVars.h:49
const Var kCVNMuonIdx([](const caf::SRProxy *sr){float longest_idx=-5.0;float longest_len=-5.0;if(sr->vtx.elastic.IsValid){for(unsigned int png_idx=0;png_idx< sr->vtx.elastic.fuzzyk.png.size();png_idx++){auto &png=sr->vtx.elastic.fuzzyk.png[png_idx];if(png.len > longest_len){longest_len=png.len;longest_idx=png_idx;}}} if(longest_len > 500.0) return longest_idx;float best_idx=-5.0;float best_score=-5.0;if(sr->vtx.elastic.IsValid){for(unsigned int png_idx=0;png_idx< sr->vtx.elastic.fuzzyk.png.size();png_idx++){auto &png=sr->vtx.elastic.fuzzyk.png[png_idx];if(png.cvnpart.muonid > best_score){best_score=png.cvnpart.muonid;best_idx=png_idx;}}}return best_idx;})
: Prong index of best muon prong by CVN score & length
Definition: CVNProngVars.h:19
const MultiVar kTruePionwBPFCalE([](const caf::SRProxy *sr){std::vector< double > E;for(unsigned int pIdx=0;pIdx< sr->vtx.elastic.fuzzyk.npng;pIdx++){if(!ProngTruthCut(pIdx, sr,{211}, 0.5, 0.5)||!hasBPFE(pIdx, sr)) continue;auto &png=sr->vtx.elastic.fuzzyk.png[pIdx];if(png.bpf.pion.IsValid){E.emplace_back(png.weightedCalE);}}return E;})
Definition: TrueProngVars.h:56
const MultiVar kTruePionBPFE([](const caf::SRProxy *sr){std::vector< double > E;for(unsigned int pIdx=0;pIdx< sr->vtx.elastic.fuzzyk.npng;pIdx++){if(!ProngTruthCut(pIdx, sr,{211}, 0.5, 0.5)||!hasBPFE(pIdx, sr)) continue;auto &png=sr->vtx.elastic.fuzzyk.png[pIdx];if(png.bpf.pion.IsValid){E.emplace_back(png.bpf.pion.energy);}}return E;})
Definition: TrueProngVars.h:50
const Var kTrueOtherCalE([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-5.0;double calE=0;unsigned int muIdx=(unsigned int) kCVNMuonIdx(sr); for(unsigned int pIdx=0;pIdx< sr->vtx.elastic.fuzzyk.png2d.size();pIdx++){auto &png=sr->vtx.elastic.fuzzyk.png2d[pIdx];calE+=png.weightedCalE;}calE+=sr->vtx.elastic.fuzzyk.orphCalE;double overlapE=0.0;if(sr->vtx.elastic.fuzzyk.png[muIdx].bpf.muon.IsValid){overlapE=sr->vtx.elastic.fuzzyk.png[muIdx].bpf.muon.overlapE;}calE+=overlapE;return calE;})
Definition: TrueProngVars.h:59
caf::Proxy< unsigned int > run
Definition: SRProxy.h:248
unsigned int PeriodFromRunFD(int run)
: Helper function to compute the period given a run number for the FD.
Definition: RunPeriods.cxx:15
const std::vector< TSpline3 > allFDMuonFits_calE
caf::Proxy< caf::SRElastic > elastic
Definition: SRProxy.h:2118
caf::Proxy< std::vector< caf::SRFuzzyKProng > > png
Definition: SRProxy.h:2043
const Var kTrueMuonTrueE([](const caf::SRProxy *sr){double E=0.0;for(unsigned int pIdx=0;pIdx< sr->vtx.elastic.fuzzyk.npng;pIdx++){if(ProngTruthCut(pIdx, sr,{13}, 0.6, 0.6)){E=sr->vtx.elastic.fuzzyk.png[pIdx].truth.p.E;}} if(E==0.0) return-5.0;return E;})
Definition: TrueProngVars.h:27
const MultiVar kTrueHadCalE([](const caf::SRProxy *sr){std::vector< double > E;for(unsigned int pIdx=0;pIdx< sr->vtx.elastic.fuzzyk.npng;pIdx++){if(ProngTruthCut(pIdx, sr,{2212, 211}, 0.5, 0.5)){E.emplace_back(sr->vtx.elastic.fuzzyk.png[pIdx].weightedCalE);}}return E;})
Definition: TrueProngVars.h:46
double getMuonEnergyByTruth(TSpline3 spline_calE, TSpline3 spline_BPFE, const caf::SRProxy *sr)
: function to compute the muon energy
const std::vector< TSpline3 > allNDMuonFits_calE
const std::vector< TSpline3 > allNDEMFits
Float_t E
Definition: plot.C:20
const MultiVar kTrueProtonTrueKE([](const caf::SRProxy *sr){std::vector< double > E;for(unsigned int pIdx=0;pIdx< sr->vtx.elastic.fuzzyk.npng;pIdx++){if(ProngTruthCut(pIdx, sr,{2212}, 0.5, 0.5)){E.emplace_back(sr->vtx.elastic.fuzzyk.png[pIdx].truth.p.E-0.938272);}}return E;})
Definition: TrueProngVars.h:29
static const double kMuonFitWeight
double energy
Definition: plottest35.C:25
const std::vector< TSpline3 > allFDHadFits
const Var kEMEByTruth([](const caf::SRProxy *sr){double EME=0.0;int run=sr->hdr.run;caf::Det_t det=sr->hdr.det;if(det==caf::kNEARDET){unsigned int period=PeriodFromRunND(run);EME=getEMShowerEnergyByTruth(allNDEMFits[period], sr);}if(det==caf::kFARDET){unsigned int period=PeriodFromRunFD(run);EME=getEMShowerEnergyByTruth(allFDEMFits[period], sr);}return EME;})
caf::StandardRecord * sr
const Var kHadEByTruth([](const caf::SRProxy *sr){double hadE=0.0;int run=sr->hdr.run;caf::Det_t det=sr->hdr.det;if(det==caf::kNEARDET){unsigned int period=PeriodFromRunND(run);hadE=getHadronEnergyByTruth(allNDHadFits[period], sr);}if(det==caf::kFARDET){unsigned int period=PeriodFromRunFD(run);hadE=getHadronEnergyByTruth(allFDHadFits[period], sr);}return hadE;})
const std::vector< TSpline3 > allFDGarbageFits
const Var kTrueMuonwBPFTrueE([](const caf::SRProxy *sr){double E=0.0;for(unsigned int pIdx=0;pIdx< sr->vtx.elastic.fuzzyk.npng;pIdx++){if(ProngTruthCut(pIdx, sr,{13}, 0.6, 0.6)&&hasBPFE(pIdx, sr)){E=sr->vtx.elastic.fuzzyk.png[pIdx].truth.p.E;}} if(E==0.0) return-5.0;return E;})
Definition: TrueProngVars.h:33
Definition: run.py:1
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2138
double getEMShowerEnergyByTruth(TSpline3 spline_calE, const caf::SRProxy *sr)
: function to compute the EM shower energy
const Var kBPFSliceEByTruth([](const caf::SRProxy *sr){double energy=0.0;int run=sr->hdr.run;caf::Det_t det=sr->hdr.det;if(det==caf::kNEARDET){unsigned int period=PeriodFromRunND(run);energy+=getMuonEnergyByTruth(allNDMuonFits_calE[period], allNDMuonFits_BPFE[period], sr);energy+=getEMShowerEnergyByTruth(allNDEMFits[period], sr);energy+=getHadronEnergyByTruth(allNDHadFits[period], sr);energy+=takeOutTrashByTruth(allNDGarbageFits[period], sr);}if(det==caf::kFARDET){unsigned int period=PeriodFromRunFD(run);energy+=getMuonEnergyByTruth(allFDMuonFits_calE[period], allFDMuonFits_BPFE[period], sr);energy+=getEMShowerEnergyByTruth(allFDEMFits[period], sr);energy+=getHadronEnergyByTruth(allFDHadFits[period], sr);energy+=takeOutTrashByTruth(allFDGarbageFits[period], sr);}return energy;})
: Main Var for reconstructed neutrino energy.
caf::Proxy< float > orphCalE
Definition: SRProxy.h:2042
caf::Proxy< bool > IsValid
Definition: SRProxy.h:2058
const std::vector< TSpline3 > allNDGarbageFits
const Var kTrueMuonwBPFCalE([](const caf::SRProxy *sr){double E=0.0;for(unsigned int pIdx=0;pIdx< sr->vtx.elastic.fuzzyk.npng;pIdx++){if(!ProngTruthCut(pIdx, sr,{13}, 0.6, 0.6)||!hasBPFE(pIdx, sr)) continue;auto &png=sr->vtx.elastic.fuzzyk.png[pIdx];if(png.bpf.muon.IsValid){E=png.calE-png.bpf.muon.overlapE;}} if(E==0.0) return-5.0;return E;})
Definition: TrueProngVars.h:52
bool ProngTruthCut(unsigned int pIdx, const caf::SRProxy *sr, std::vector< int > pdgCodes, double eff, double pur)
: Returns true if the prong matches the truth criteria.
const Var kGarbageEByTruth([](const caf::SRProxy *sr){double garE=0.0;int run=sr->hdr.run;caf::Det_t det=sr->hdr.det;if(det==caf::kNEARDET){unsigned int period=PeriodFromRunND(run);garE=takeOutTrashByTruth(allNDGarbageFits[period], sr);}if(det==caf::kFARDET){unsigned int period=PeriodFromRunFD(run);garE=takeOutTrashByTruth(allFDGarbageFits[period], sr);}return garE;})
const MultiVar kTrueEMTrueE([](const caf::SRProxy *sr){std::vector< double > E;for(unsigned int pIdx=0;pIdx< sr->vtx.elastic.fuzzyk.npng;pIdx++){if(ProngTruthCut(pIdx, sr,{11, 22, 111}, 0.5, 0.5)){E.emplace_back(sr->vtx.elastic.fuzzyk.png[pIdx].truth.p.E);}}return E;})
Definition: TrueProngVars.h:38
const Var kTrueGarbageE([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-5.0;double energy=0.0;int run=sr->hdr.run;caf::Det_t det=sr->hdr.det;if(det==caf::kNEARDET){unsigned int period=PeriodFromRunND(run);energy+=getMuonEnergy(allNDMuonFits_calE[period], allNDMuonFits_BPFE[period], sr);energy+=getEMShowerEnergy(allNDEMFits[period], sr);energy+=getHadronEnergy(allNDHadFits[period], sr);}if(det==caf::kFARDET){unsigned int period=PeriodFromRunFD(run);energy+=getMuonEnergy(allFDMuonFits_calE[period], allFDMuonFits_BPFE[period], sr);energy+=getEMShowerEnergy(allFDEMFits[period], sr);energy+=getHadronEnergy(allFDHadFits[period], sr);}return(sr->mc.nu[0].E-energy);})
Definition: TrueProngVars.h:94
const MultiVar kTrueEMCalE([](const caf::SRProxy *sr){std::vector< double > E;for(unsigned int pIdx=0;pIdx< sr->vtx.elastic.fuzzyk.npng;pIdx++){if(ProngTruthCut(pIdx, sr,{11, 22, 111}, 0.5, 0.5)){E.emplace_back(sr->vtx.elastic.fuzzyk.png[pIdx].weightedCalE);}}return E;})
Definition: TrueProngVars.h:45
const std::vector< TSpline3 > allNDMuonFits_BPFE
const MultiVar kTruePionwBPFTrueE([](const caf::SRProxy *sr){std::vector< double > E;for(unsigned int pIdx=0;pIdx< sr->vtx.elastic.fuzzyk.npng;pIdx++){if(ProngTruthCut(pIdx, sr,{211}, 0.5, 0.5)&&hasBPFE(pIdx, sr)){E.emplace_back(sr->vtx.elastic.fuzzyk.png[pIdx].truth.p.E);}}return E;})
Definition: TrueProngVars.h:37
caf::Proxy< caf::SRVertexBranch > vtx
Definition: SRProxy.h:2146
const Var kBPFMuonEByTruth([](const caf::SRProxy *sr){double muE=0.0;int run=sr->hdr.run;caf::Det_t det=sr->hdr.det;if(det==caf::kNEARDET){unsigned int period=PeriodFromRunND(run);muE=getMuonEnergyByTruth(allNDMuonFits_calE[period], allNDMuonFits_BPFE[period], sr);}if(det==caf::kFARDET){unsigned int period=PeriodFromRunFD(run);muE=getMuonEnergyByTruth(allFDMuonFits_calE[period], allFDMuonFits_BPFE[period], sr);}return muE;})
: Returns the recoE of the true selected muon.
Definition: TrueProngVars.h:99
const MultiVar kTrueProtonwBPFCalE([](const caf::SRProxy *sr){std::vector< double > E;for(unsigned int pIdx=0;pIdx< sr->vtx.elastic.fuzzyk.npng;pIdx++){if(!ProngTruthCut(pIdx, sr,{2212}, 0.5, 0.5)||!hasBPFE(pIdx, sr)) continue;auto &png=sr->vtx.elastic.fuzzyk.png[pIdx];if(png.bpf.proton.IsValid){E.emplace_back(png.weightedCalE);}}return E;})
Definition: TrueProngVars.h:54
const std::vector< TSpline3 > allFDEMFits
bool hasBPFE(unsigned int pIdx, const caf::SRProxy *sr)
const MultiVar kTrueHadTrueE([](const caf::SRProxy *sr){std::vector< double > E;for(unsigned int pIdx=0;pIdx< sr->vtx.elastic.fuzzyk.npng;pIdx++){if(ProngTruthCut(pIdx, sr,{2212}, 0.5, 0.5)){E.emplace_back(sr->vtx.elastic.fuzzyk.png[pIdx].truth.p.E-0.938272);}else if(ProngTruthCut(pIdx, sr,{211}, 0.5, 0.5)){E.emplace_back(sr->vtx.elastic.fuzzyk.png[pIdx].truth.p.E);}}return E;})
Definition: TrueProngVars.h:39
const MultiVar kTruePionTrueE([](const caf::SRProxy *sr){std::vector< double > E;for(unsigned int pIdx=0;pIdx< sr->vtx.elastic.fuzzyk.npng;pIdx++){if(ProngTruthCut(pIdx, sr,{211}, 0.5, 0.5)){E.emplace_back(sr->vtx.elastic.fuzzyk.png[pIdx].truth.p.E);}}return E;})
Definition: TrueProngVars.h:31
unsigned int PeriodFromRunND(int run)
: Helper function to compute the period given a run number for the ND.
Definition: RunPeriods.cxx:34
caf::Proxy< caf::Det_t > det
Definition: SRProxy.h:232
const Var kTrueMuonCalE([](const caf::SRProxy *sr){double E=0.0;for(unsigned int pIdx=0;pIdx< sr->vtx.elastic.fuzzyk.npng;pIdx++){if(ProngTruthCut(pIdx, sr,{13}, 0.6, 0.6)){E=sr->vtx.elastic.fuzzyk.png[pIdx].weightedCalE;}}if(E==0.0) return-5.0;return E;})
Definition: TrueProngVars.h:41