NumubarCCpi0_Vars.cxx
Go to the documentation of this file.
3 
4 #include <iostream>
5 namespace ana
6 {
7  namespace numubarccpi0
8  {
9 
10  //====================//
11  // Truth vars //
12  //====================//
13 
14  const NuTruthVar kTruthTag_NT([](const caf::SRNeutrinoProxy* truth)
15  {
16  bool AntiNuSig = kAntiNuSig_NT(truth);
17  bool NuSig = kNuSig_NT(truth);
18  bool SecoPi0 = kNumuCC_NT(truth) && kSecoPi0_NT(truth);
19  bool NoPi0 = kNumuCC_NT(truth) && kNoPi0_NT(truth);
20  bool NC = kNC_NT(truth);
21  bool Other = kIsOther_NT(truth);
22  if(AntiNuSig) return 1;
23  else if(NuSig) return 2;
24  else if(SecoPi0) return 3;
25  else if(NoPi0) return 4;
26  else if(NC) return 5;
27  else if(Other) return 6;
28  else return 0;
29  });
30 
31  // Number of true pi0s
32  //=========================================================================
33  const NuTruthVar kNPrimaryPi0_NT([](const caf::SRNeutrinoProxy* truth)
34  {
35  int nprims = truth->prim.size();
36  int countpi0 = 0;
37  for(int i = 0; i < nprims; ++i){
38  if(truth->prim[i].pdg == 111) countpi0++;
39  }
40  return countpi0;
41  });
43  {
44  int nprims = truth->prim.size();
45  int countpi0 = 0;
46  for(int i = 0; i < nprims; ++i){
47  auto& daughters = truth->prim[i].daughterlist;
48  if(daughters.empty()) continue;
49  for(const auto& pdg : daughters){
50  if(pdg == 111) countpi0++;
51  }
52  }
53  return countpi0;
54  });
55  //=========================================================================
56  // True Neutrino Interaction
57  const NuTruthVar kTrueNuE_NT([](const caf::SRNeutrinoProxy* nu)
58  {
59  return nu->E;
60  });
61 
62  const NuTruthVar kTrueNuQ2_NT([](const caf::SRNeutrinoProxy* nu)
63  {
64  return nu->q2;
65  });
66 
67  const NuTruthVar kTrueNuW_NT([](const caf::SRNeutrinoProxy* nu)
68  {
69  float val = nu->W2;
70  if(val < 0) return (double)-5.0;
71  return (double)sqrt(nu->W2);
72  });
73 
74 
75 
76 
77  //True Muon Kinematic Variables
78  //=========================================================================
80  ([](const caf::SRNeutrinoProxy* nu)
81  {
82  float ke = -5;
83  if(abs(nu->pdg) != 14 || !nu->iscc)
84  return ke;
85  int nprims = nu->prim.size();
86  for(int iprim = 0; iprim < nprims; iprim++){
87  if(abs(nu->prim[iprim].pdg) == 13){
88  double E= nu->prim[iprim].p.T();
89  ke = E-MuonMass();
90  }
91  }//end loop over primaries
92  return ke;
93  });
94 
96  ([](const caf::SRNeutrinoProxy* truth)
97  {
98  float pmag = -5;
99  if(abs(truth->pdg) != 14 || !truth->iscc)
100  return pmag;
101  int nprims = truth->prim.size();
102  for(int iprim = 0; iprim < nprims; iprim++){
103  if(abs(truth->prim[iprim].pdg) == 13){
104  pmag = truth->prim[iprim].p.Mag();
105  }
106  }
107 
108  return pmag;
109  });
110 
111 
113  ([](const caf::SRNeutrinoProxy* truth)
114  {
115  if(abs(truth->pdg) != 14 || !truth->iscc)
116  return -5.0;
117  int nprims = truth->prim.size();
118  for(int iprim = 0; iprim < nprims; iprim++){
119  if(abs(truth->prim[iprim].pdg) == 13){
120  TVector3 mudir = truth->prim[iprim].p.Vect();
121  TVector3 beamdir = NuMIBeamDirection(caf::kNEARDET);
122  return mudir.Unit().Dot(beamdir.Unit());
123  }
124  }//end loop over primaries
125  return -5.0;
126  });
127 
129  ([](const caf::SRNeutrinoProxy* truth)
130  {
131  if(abs(truth->pdg) != 14 || !truth->iscc)
132  return -5.0;
133  int nprims = truth->prim.size();
134  for(int iprim = 0; iprim < nprims; iprim++){
135  if(abs(truth->prim[iprim].pdg) == 13){
136  TVector3 mudir = truth->prim[iprim].p.Vect();
137  TVector3 beamdir = NuMIBeamDirection(caf::kNEARDET);
138  return mudir.Unit().Angle(beamdir.Unit());
139  }
140  }//end loop over primaries
141  return -5.0;
142  });
143 
145  ([](const caf::SRNeutrinoProxy* truth)
146  {
147  float TrueMuDirZ_Rad_NT = kTrueMuDirZ_Rad_NT(truth);
148  if(TrueMuDirZ_Rad_NT > TMath::Pi() || TrueMuDirZ_Rad_NT < 0.)
149  return -5.0f;
150  return (float)(TrueMuDirZ_Rad_NT*180./TMath::Pi());
151  });
152 
153  //=========================================================================
154  //True Pi0 Kinetic Vars
155  //=========================================================================
157  ([](const caf::SRNeutrinoProxy* truth)
158  {
159  int nprims = truth->prim.size();
160  int pi0ID = -1;
161  float pi0E = -5.0;
162  for(int i = 0; i < nprims; ++i){
163  if(truth->prim[i].pdg == 111 && pi0E < truth->prim[i].p.E){
164  pi0ID = i;
165  pi0E = truth->prim[i].p.E;
166  }
167  }
168  if(pi0ID < 0) return -5;
169  return pi0ID;
170  });
171 
173  ([](const caf::SRNeutrinoProxy* truth)
174  {
175  int TruePi0_ID = kTruePrimPi0_ID_NT(truth);
176  if(TruePi0_ID < 0) return -5.0f;
177  return (float)truth->prim[TruePi0_ID].p.E;
178  });
179 
181  ([](const caf::SRNeutrinoProxy* truth)
182  {
183  int TruePi0_ID = kTruePrimPi0_ID_NT(truth);
184  if(TruePi0_ID < 0) return -5.0f;
185  return (float)truth->prim[TruePi0_ID].p.Mag();
186  });
187 
189  ([](const caf::SRNeutrinoProxy* truth)
190  {
191  int TruePi0_ID = kTruePrimPi0_ID_NT(truth);
192  if(TruePi0_ID < 0) return -5.0f;
193  TVector3 Pi0P = truth->prim[TruePi0_ID].p.Vect();
195  return (float)Pi0P.Unit().Dot(beamDirND.Unit());
196  });
197 
199  ([](const caf::SRNeutrinoProxy* truth)
200  {
201  int TruePi0_ID = kTruePrimPi0_ID_NT(truth);
202  if(TruePi0_ID < 0) return -5.0f;
203  TVector3 Pi0P = truth->prim[TruePi0_ID].p.Vect();
204  TVector3 beamDirND = NuMIBeamDirection(caf::kNEARDET);
205  return (float)Pi0P.Unit().Angle(beamDirND.Unit());
206  });
207 
209  ([](const caf::SRNeutrinoProxy* truth)
210  {
211  float TruePi0DirZ_Rad_NT = kTruePrimPi0DirZ_Rad_NT(truth);
212  if(TruePi0DirZ_Rad_NT > TMath::Pi() || TruePi0DirZ_Rad_NT < 0.) return -5.0f;
213  return (float)(TruePi0DirZ_Rad_NT*180./TMath::Pi());
214  });
215 
217  ([](const caf::SRNeutrinoProxy* truth)
218  {
219  float nPrimaryPi0 = kNPrimaryPi0_NT(truth);
220  if(nPrimaryPi0 < 1) return -5.f;
221 
222  //---------------------------
223  //NB: For now this only returns 1 value for the event
224  // specifically the energy of the most energetic Pi0 in an event
225  //---------------------------
226 
227  float pi0_t = -5;
228 
229  const int nprims = truth->prim.size();
230 
231  for(int iprim = 0; iprim < nprims; ++iprim){
232  if(truth->prim[iprim].pdg != 111) continue;
233  double E= truth->prim[iprim].p.T();
234  float ke = E-Pi0Mass();
235  if(ke > pi0_t) pi0_t = ke;
236  }
237 
238  return float(pi0_t);
239  });
240 
241  //=========================================================================
242 
243  //True Pi0 Kinetic Vars -- Secondary Pi0
244  //=========================================================================
246  {
247  int nprims = truth->prim.size();
248  int pi0ID1 = -1;
249  float pi0E = -5.0;
250  for(int i = 0; i < nprims; ++i){
251  auto& daughterPDG = truth->prim[i].daughterlist;
252  auto& daughterE = truth->prim[i].daughterEnergies;
253  if(daughterPDG.empty()) continue;
254  for(int j = 0; j < (int)daughterPDG.size(); ++j){
255  if(daughterPDG[j] == 111 && pi0E < daughterE[j]){
256  pi0ID1 = i;
257  pi0E = daughterE[j];
258  }
259  }
260  }
261  if(pi0ID1 < 0) return -5;
262  return pi0ID1;
263  });
264 
266  {
267  int nprims = truth->prim.size();
268  int pi0ID2 = -1;
269  float pi0E = -5.0;
270  for(int i = 0; i < nprims; ++i){
271  auto& daughterPDG = truth->prim[i].daughterlist;
272  auto& daughterE = truth->prim[i].daughterEnergies;
273  if(daughterPDG.empty()) continue;
274  for(int j = 0; j < (int)daughterPDG.size(); ++j){
275  if(daughterPDG[j] == 111 && pi0E < daughterE[j]){
276  pi0ID2 = j;
277  pi0E = daughterE[j];
278  }
279  }
280  }
281  if(pi0ID2 < 0) return -5;
282  return pi0ID2;
283  });
284 
285  const NuTruthVar kTrueSecoPi0E_NT([](const caf::SRNeutrinoProxy* truth)
286  {
287  int TruePi0_ID1 = kTrueSecoPi0_ID1_NT(truth);
288  int TruePi0_ID2 = kTrueSecoPi0_ID2_NT(truth);
289  if(TruePi0_ID1 < 0 || TruePi0_ID2 < 0) return -5.0f;
290  return (float)truth->prim[TruePi0_ID1].daughterEnergies[TruePi0_ID2];
291  });
292 
293  const NuTruthVar kTrueSecoPi0P_NT([](const caf::SRNeutrinoProxy* truth)
294  {
295  int TruePi0_ID1 = kTrueSecoPi0_ID1_NT(truth);
296  int TruePi0_ID2 = kTrueSecoPi0_ID2_NT(truth);
297  double pi0P = -5.0;
298  if(TruePi0_ID1 < 0 || TruePi0_ID2 < 0)
299  return (float)pi0P;
300  double pi0E = truth->prim[TruePi0_ID1].daughterEnergies[TruePi0_ID2];
301  if(pi0E < Pi0Mass()) pi0P = 0.;
302  else pi0P = sqrt(std::pow(pi0E, 2) - std::pow(Pi0Mass(), 2));
303  return (float)pi0P;
304  });
305 
306 
307  //=========================================================================
308  //True Single Transverse Variables - Work in Progress need to test a bit
309  //=========================================================================
310  const TVector3 kNu_Momentum(const caf::SRNeutrinoProxy* truth)
311  {
312  return truth->p.Vect();
313  }
314 
315  const TVector3 kMuon_Momentum(const caf::SRNeutrinoProxy* truth)
316  {
317  int nprims = truth->prim.size();
318 
319  //Lepton is always the first listed primary
320  //This should help with NC and all CC events
321  if(nprims > 0)
322  return truth->prim[0].p.Vect();
323 
324  const TVector3 miss(0,0,0);
325  return miss;
326  }
327 
328  const TVector3 kNucleon_Momentum(const caf::SRNeutrinoProxy* truth)
329  {
330  int nprims = truth->prim.size();
331  for(int iprim = 0; iprim < nprims; iprim++){
332  if( ((truth->pdgorig > 0) && (abs(truth->prim[iprim].pdg) == 2212)) ||
333  ((truth->pdgorig < 0) && (abs(truth->prim[iprim].pdg) == 2112))){
334  return truth->prim[iprim].p.Vect();
335  }
336  }
337 
338  const TVector3 miss(0,0,0);
339  return miss;
340  }
341 
342  const TVector3 kPi0_Momentum(const caf::SRNeutrinoProxy* truth)
343  {
344  const int primPi0ID = kTruePrimPi0_ID_NT(truth);
345 
346  const TVector3 miss(0,0,0);
347  if(primPi0ID < 0) return miss;
348 
349  return truth->prim[primPi0ID].p.Vect();
350  }
351 
352  const TVector3 kH_Momentum(const caf::SRNeutrinoProxy* truth)
353  {
354  TVector3 p_n = kNucleon_Momentum(truth);
355  TVector3 p_pi0 = kPi0_Momentum(truth);
356 
357  TVector3 p_result(0,0,0);
358 
359  p_result += p_n;
360  p_result += p_pi0;
361 
362  return p_result;
363  }
364 
366  ([](const caf::SRNeutrinoProxy* truth)
367  {
368  TVector3 p_h = kH_Momentum(truth);
369 
370  return p_h.Mag();
371  });
372 
373  const TVector3 kDeltaPt(const caf::SRNeutrinoProxy* truth)
374  {
375  TVector3 p_mu = kMuon_Momentum(truth);
376  TVector3 p_h = kH_Momentum(truth);
377  TVector3 p_nu = kNu_Momentum(truth);
378 
379  TVector3 p_mu_T = p_nu.Cross(p_mu);
380  TVector3 p_h_T = p_nu.Cross(p_h);
381 
382  TVector3 delta_p_T(p_mu_T.X(), p_mu_T.Y(), p_mu_T.Z());
383  delta_p_T += p_h_T;
384 
385  return delta_p_T;
386  }
387 
389  ([](const caf::SRNeutrinoProxy* truth)
390  {
391  TVector3 p_mu = kMuon_Momentum(truth);
392  TVector3 p_h = kH_Momentum(truth);
393  TVector3 p_nu = kNu_Momentum(truth);
394 
395  TVector3 p_mu_T = p_nu.Cross(p_mu);
396  TVector3 p_h_T = p_nu.Cross(p_h);
397 
398  TVector3 delta_p_T(p_mu_T.X(), p_mu_T.Y(), p_mu_T.Z());
399  delta_p_T += p_h_T;
400 
401  TVector3 z_TT = p_mu_T.Unit();
402  z_TT.SetX(0.);
403  z_TT.SetY(0.);
404 
405  double delta_p_TT = delta_p_T.Dot(z_TT);
406 
407  if(std::isnan(delta_p_TT)) return -5.0;
408  if(std::isinf(delta_p_TT)) return -5.0;
409  return delta_p_TT;
410  });
411 
412 
414  ([](const caf::SRNeutrinoProxy* truth)
415 
416  {
417  TVector3 p_mu = kMuon_Momentum(truth);
418  TVector3 p_h = kH_Momentum(truth);
419  TVector3 p_nu = kNu_Momentum(truth);
420 
421  TVector3 p_mu_T = p_nu.Cross(p_mu);
422  TVector3 p_h_T = p_nu.Cross(p_h);
423 
424  TVector3 delta_p_T(0,0,0);
425  delta_p_T += p_mu_T;
426  delta_p_T += p_h_T;
427 
428  //Calculating \delta \alpha_{T}
429  double numerator = delta_p_T.Dot(-p_mu_T);
430  double mag_p_mu_T = p_mu_T.Mag();
431  double mag_delta_p_T = delta_p_T.Mag();
432 
433  if(mag_p_mu_T == 0 || mag_delta_p_T == 0)
434  return -5.0;
435 
436  double val = numerator / (mag_p_mu_T * mag_delta_p_T);
437  if(std::isnan(val)) return -5.0;
438  if(std::isinf(val)) return -5.0;
439  double delta_alpha = TMath::ACos(val);
440 
441  if(std::isnan(delta_alpha)) return -5.0;
442  if(std::isinf(delta_alpha)) return -5.0;
443  if(delta_alpha > TMath::Pi() || delta_alpha < 0.) return -5.0;
444  return (double)(delta_alpha*180./TMath::Pi());
445  });
446 
447  //=========================================================================
448  //===================//
449  // Reco vars //
450  //===================//
451 
452  // Number of prongs in one slice
453  const Var kNProngs([](const caf::SRProxy* sr)
454  {
455  if(!sr->vtx.elastic.IsValid) return -5;
456  return (int)sr->vtx.elastic.fuzzyk.npng;
457  });
458 
459  // Best muon track / MuonID
460  //=========================================================================
461  const Var kBestTrack([](const caf::SRProxy* sr)
462  {
463  return (int)sr->trk.kalman.idxmuonid;
464  });
465  const Var kMuonID([](const caf::SRProxy* sr)
466  {
467  int ibesttrk = kBestTrack(sr);
468  if(sr->trk.kalman.ntracks < 1)
469  return -5.0f;
470  if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
471  return -5.0f;
472  return (float)sr->trk.kalman.tracks[ibesttrk].muonid;
473  });
474  //=========================================================================
475 
476  //---------------------------------
477  // Muon Energy for Active region
478  //---------------------------------
479  inline float MuonEAct(double trklenact){
480  double p0 = 1.67012e-01;
481  double p1 = 1.79305e-01;
482  double p2 = 3.74708e-03;
483  double p3 = -1.54232e-04;
484  float MuonE = 0.0;
485  if (trklenact <= 0.0) return 0.0;
486  MuonE = p0
487  + p1 * trklenact
488  + p2 * std::pow(trklenact, 2)
489  + p3 * std::pow(trklenact, 3);
490  return MuonE;
491  }
492  //---------------------------------------
493  // Muon Energy for Active+Catcher region
494  //(for Active region Track length)
495  //---------------------------------------
496  inline float MuonEActandCat(double trklenactandcat){
497  double p0 = 1.21130e-02;
498  double p1 = 1.97903e-01;
499  double p2 = 7.82459e-04;
500  float MuonE = 0.0;
501  if (trklenactandcat <= 0.0) return 0.0;
502  MuonE = p0
503  + p1 * trklenactandcat
504  + p2 * std::pow(trklenactandcat, 2);
505  return MuonE;
506  }
507  //---------------------------------------
508  // Muon Energy for Active+Catcher region
509  //(for catcher region Track length)
510  //---------------------------------------
511  inline float MuonECat(double trklencat){
512  double offset = 1.31325e-01;
513  double slope = 5.35146e-01;
514  float MuonE = 0.0;
515  if (trklencat <= 0.0) return 0.0;
516  MuonE = slope*trklencat + offset;
517  return MuonE;
518  }
519  //---------------------------------------
520  // Haronic Energy
521  //----------------------------------------
522  inline float VisibleHadE(float vishadE) {
523  double p0 = 5.85254e-02;
524  double p1 = 1.27796e+00;
525  double p2 = 3.75457e-01;
526  double p3 = -5.45618e-01;
527  double p4 = 1.65975e-01;
528  float HadE = 0.0;
529  HadE = p0
530  + p1 * vishadE
531  + p2 * std::pow(vishadE, 2)
532  + p3 * std::pow(vishadE, 3)
533  + p4 * std::pow(vishadE, 4);
534  return HadE;
535  }
536 
537  //==============================================================================
538  // Redefine track lengths in active and catcher volumes with the new PID.
539  // Note that at this moment what is called leninact is actually distance from
540  // the start of the track to the transition plane. If the transition plane is
541  // outside of the track, that number is the track length plus the length of the
542  // projection from the last hit along the direction determined by the last two
543  // hits to the transition plane. The lenincat behaves similarly.
544  // In the next production the two variables should be changed to what they are
545  // supposed to mean literally.
546  const Var kTrkLenAct([](const caf::SRProxy* sr)
547  {
548  int ibesttrk = kBestTrack(sr);
549  if(sr->trk.kalman.ntracks < 1)
550  return -1000.f;
551  if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
552  return -1000.f;
553  // If leninact is positive and lenincat is negative,
554  // the transition plane is to the right of the track.
555  if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 &&
556  sr->trk.kalman.tracks[ibesttrk].lenincat < 0)
557  return float((sr->trk.kalman.tracks[ibesttrk].leninact / 100.)
558  +(sr->trk.kalman.tracks[ibesttrk].lenincat / 100.));
559  // If leninact is positive and lenincat is positive,
560  // the transition plane is in the middle of the track.
561  if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 &&
562  sr->trk.kalman.tracks[ibesttrk].lenincat > 0)
563  return float(sr->trk.kalman.tracks[ibesttrk].leninact / 100.);
564  // If leninact is negative and lenincat is positive,
565  // the transition plane is to the left of the track.
566  if(sr->trk.kalman.tracks[ibesttrk].leninact < 0 &&
567  sr->trk.kalman.tracks[ibesttrk].lenincat > 0)
568  return 0.f;
569  return -1000.f;
570  });
571 
572  const Var kTrkLenCat([](const caf::SRProxy* sr)
573  {
574  int ibesttrk = kBestTrack(sr);
575  if(sr->trk.kalman.ntracks < 1)
576  return -1000.f;
577  if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
578  return -1000.f;
579  // If leninact is positive and lenincat is negative,
580  // the transition plane is to the right of the track.
581  if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 &&
582  sr->trk.kalman.tracks[ibesttrk].lenincat < 0)
583  return 0.f;
584  // If leninact is positive and lenincat is positive,
585  // the transition plane is in the middle of the track.
586  if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 &&
587  sr->trk.kalman.tracks[ibesttrk].lenincat > 0)
588  return float(sr->trk.kalman.tracks[ibesttrk].lenincat / 100.);
589  // If leninact is negative and lenincat is positive,
590  // the transition plane is to the left of the track.
591  if(sr->trk.kalman.tracks[ibesttrk].leninact < 0 &&
592  sr->trk.kalman.tracks[ibesttrk].lenincat > 0)
593  return float((sr->trk.kalman.tracks[ibesttrk].leninact / 100.)
594  +(sr->trk.kalman.tracks[ibesttrk].lenincat / 100.));
595  return -1000.f;
596  });
597 
598  // Muon Energy for Inclusive numuCC
599  //===============================================================
600  const Var kIncXsecMuonE([](const caf::SRProxy *sr)
601  {
602  float muonE = 0.0;
603  float muonEact = 0.0;
604  float muonEcat = 0.0;
605  float muonEactandcat = 0.0;
606  float trkLenAct = 0.f;
607  float trkLenCat = 0.f;
608  trkLenAct = kTrkLenAct(sr);
609  trkLenCat = kTrkLenCat(sr);
610  int ibesttrk = kBestTrack(sr);
611  if(sr->trk.kalman.ntracks < 1)
612  return -1000.f;
613  if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
614  return -1000.f;
615  if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 &&
616  sr->trk.kalman.tracks[ibesttrk].lenincat < 0)
617  muonEact = MuonEAct(trkLenAct);
618  else if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 &&
619  sr->trk.kalman.tracks[ibesttrk].lenincat > 0){
620  muonEcat = MuonECat(trkLenCat);
621  muonEactandcat = MuonEActandCat(trkLenAct);
622  muonE = muonEactandcat + muonEcat;
623  }
624  return muonE + muonEact;
625  });
626 
627  // Hadronic Energy for Inclusive numuCC
628  //===============================================================
629  const Var kIncXsecHadE([](const caf::SRProxy* sr)
630  {
631  int ibesttrk = kBestTrack(sr);
632  if(sr->trk.kalman.ntracks < 1)
633  return -1000.f;
634  if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
635  return -1000.f;
636  // Extra energy (hadronic contamination) associated with muon track
637  float ExtraHadE = sr->trk.kalman.tracks[ibesttrk].overlapE;
638  // calorimetric slice energy - Calorimetric Energy of Muon Track
639  float CalhadE = sr->slc.calE - sr->trk.kalman.tracks[ibesttrk].calE;
640  // Add calorimetric hadE and Overlap energy in the muon track
641  float vishadE = CalhadE + ExtraHadE;
642  float hadE = 0.0;
643  hadE = VisibleHadE(vishadE);
644  return hadE;
645  });
646 
647 
648  // Muon-like prong index
649  //==================================================================================
650  const Var kMuonPngIdx([](const caf::SRProxy *sr)
651  {
652  int longest_idx = -5;
653  float longest_len = -5.0;
654  // loop over all verticies, prongs, and tracks to find the longest muon track
655  if(sr->vtx.elastic.IsValid){
656  for(unsigned int png_idx = 0; png_idx < sr->vtx.elastic.fuzzyk.png.size(); png_idx++) {
657  auto &png = sr->vtx.elastic.fuzzyk.png[png_idx];
658  if(png.len > longest_len) {
659  longest_len = png.len;
660  longest_idx = png_idx;
661  }
662  } // end loop over prongs
663  } // end loop over verticies
664  // If there are ANY prongs longer than 500 cm, use the longest prong:
665  if(longest_len > 500.0) return longest_idx;
666  // Otherwise find prong which is most colinear with muon track...
667  int trk_idx = kBestTrack(sr);
668  if(sr->trk.kalman.ntracks < 1 || trk_idx < 0 || trk_idx >= int(sr->trk.kalman.ntracks)) return -5;
669  auto &trk = sr->trk.kalman.tracks[trk_idx];
670  int best_idx = -5;
671  float best_cos = -5.0;
672  if(sr->vtx.elastic.IsValid){
673  for(unsigned int png_idx = 0; png_idx < sr->vtx.elastic.fuzzyk.png.size(); png_idx++) {
674  auto &png = sr->vtx.elastic.fuzzyk.png[png_idx];
675  if(trk.dir.Unit().Dot(png.dir.Unit()) > best_cos) {
676  best_cos = trk.dir.Unit().Dot(png.dir.Unit());
677  best_idx = png_idx;
678  }
679  } // end loop over prongs
680  } // end loop over verticies
681  return best_idx;
682  });
683 
684  const Var kCVNMuonIdx_4view([](const caf::SRProxy *sr)
685  {
686  int longest_idx = -5;
687  float longest_len = -5.0;
688  // loop over all verticies, prongs, and tracks to find the longest muon track
689  if(sr->vtx.elastic.IsValid){
690  for(unsigned int png_idx = 0; png_idx < sr->vtx.elastic.fuzzyk.png.size(); png_idx++) {
691  auto &png = sr->vtx.elastic.fuzzyk.png[png_idx];
692  if(png.len > longest_len) {
693  longest_len = png.len;
694  longest_idx = png_idx;
695  }
696  } // end loop over prongs
697  } // end loop over verticies
698  // If there are ANY prongs longer than 500 cm, use the longest prong:
699  if(longest_len > 500.0) return longest_idx;
700  // Otherwise find the highest scoring prong...
701  int best_idx = -5;
702  float best_score = -5.0;
703  if(sr->vtx.elastic.IsValid){
704  for(unsigned int png_idx = 0; png_idx < sr->vtx.elastic.fuzzyk.png.size(); png_idx++) {
705  auto &png = sr->vtx.elastic.fuzzyk.png[png_idx];
706  if(png.cvnpart.muonid > best_score) {
707  best_score = png.cvnpart.muonid;
708  best_idx = png_idx;
709  }
710  } // end loop over prongs
711  } // end loop over verticies
712  return best_idx;
713  });
714 
715  const Var kCVNMuonIdx_5label([](const caf::SRProxy *sr)
716  {
717  int longest_idx = -5;
718  float longest_len = -5.0;
719  // loop over all verticies, prongs, and tracks to find the longest muon track
720  if(sr->vtx.elastic.IsValid){
721  for(unsigned int png_idx = 0; png_idx < sr->vtx.elastic.fuzzyk.png.size(); png_idx++) {
722  auto &png = sr->vtx.elastic.fuzzyk.png[png_idx];
723  if(png.len > longest_len) {
724  longest_len = png.len;
725  longest_idx = png_idx;
726  }
727  } // end loop over prongs
728  } // end loop over verticies
729  // If there are ANY prongs longer than 500 cm, use the longest prong:
730  if(longest_len > 500.0) return longest_idx;
731  // Otherwise find the highest scoring prong...
732  int best_idx = -5;
733  float best_score = -5.0;
734  if(sr->vtx.elastic.IsValid){
735  for(unsigned int png_idx = 0; png_idx < sr->vtx.elastic.fuzzyk.png.size(); png_idx++) {
736  auto &png = sr->vtx.elastic.fuzzyk.png[png_idx];
737  if(png.spprongcvnpart5label.muonid > best_score) {
738  best_score = png.spprongcvnpart5label.muonid;
739  best_idx = png_idx;
740  }
741  } // end loop over prongs
742  } // end loop over verticies
743  return best_idx;
744  });
745 
746  /// Function to return "adjusted" CVN-prong muon score:
747  double GetCVNProngMuonScore_4view(int ProngIdx, const caf::SRProxy *sr) {
748  // Sanity checks first to prevent segfaults...
749  if(!sr->vtx.elastic.IsValid) return -5.0;
750  if(sr->vtx.elastic.fuzzyk.npng == 0) return -5.0;
751  if(ProngIdx >= (int)sr->vtx.elastic.fuzzyk.npng || ProngIdx < 0) return -5.0;
752  double score = sr->vtx.elastic.fuzzyk.png[ProngIdx].cvnpart.muonid;
753  double len = sr->vtx.elastic.fuzzyk.png[ProngIdx].len;
754  if(len > 500.0) return 1.0;
755  else return score;
756  }
757  double GetCVNProngMuonScore_5label(int ProngIdx, const caf::SRProxy *sr) {
758  // Sanity checks first to prevent segfaults...
759  if(!sr->vtx.elastic.IsValid) return -5.0;
760  if(sr->vtx.elastic.fuzzyk.npng == 0) return -5.0;
761  if(ProngIdx >= (int)sr->vtx.elastic.fuzzyk.npng || ProngIdx < 0) return -5.0;
762  double score = sr->vtx.elastic.fuzzyk.png[ProngIdx].spprongcvnpart5label.muonid;
763  double len = sr->vtx.elastic.fuzzyk.png[ProngIdx].len;
764  if(len > 500.0) return 1.0;
765  else return score;
766  }
767 
769  {
770  return (float)GetCVNProngMuonScore_4view(kMuonPngIdx(sr), sr);
771  });
772  const Var kMuonPng_MuonScore_5label([](const caf::SRProxy *sr)
773  {
774  return (float)GetCVNProngMuonScore_5label(kMuonPngIdx(sr), sr);
775  });
777  {
778  return (float)GetCVNProngMuonScore_4view(kCVNMuonIdx_4view(sr), sr);
779  });
781  {
782  return (float)GetCVNProngMuonScore_5label(kCVNMuonIdx_4view(sr), sr);
783  });
785  {
786  return (float)GetCVNProngMuonScore_4view(kCVNMuonIdx_5label(sr), sr);
787  });
789  {
790  return (float)GetCVNProngMuonScore_5label(kCVNMuonIdx_5label(sr), sr);
791  });
792 
793 
794  // Prong1 & Prong2 index generator
795  //=======================================================================
796  int Prong1ID_Generator (const caf::SRProxy* sr, std::string id_type) {
797  int ProngNum = kNProngs(sr);
798  if(ProngNum < 2) return -5;
799  //int MuonPngID = kCVNMuonIdx(sr);
800  int MuonPngID = kMuonPngIdx(sr);
801  int EMPngID = 0;
802  float EMcvn = -5.0;
803  float tempcvn = -5.0;
804  for(int i = 0; i < ProngNum; ++i) {
805  if(i == MuonPngID) continue;
806  if(id_type == "4view"){
807  tempcvn = sr->vtx.elastic.fuzzyk.png[i].cvnpart.emid;
808  }
809  else if(id_type == "emid"){
810  tempcvn = sr->vtx.elastic.fuzzyk.png[i].spprongcvnpartnumuccemid.emid;
811  }
812  else if(id_type == "5label"){
813  tempcvn = sr->vtx.elastic.fuzzyk.png[i].spprongcvnpart5label.emid;
814  }
815  else {
816  tempcvn = -5.0;
817  }
818  if(tempcvn > EMcvn) {
819  EMPngID = i;
820  EMcvn = tempcvn;
821  }
822  }
823  return EMPngID;
824  }
825 
826  int Prong2ID_Generator (const caf::SRProxy* sr, std::string id_type) {
827  int ProngNum = kNProngs(sr);
828  if(ProngNum < 3) return -5;
829  //int MuonPngID = kCVNMuonIdx(sr);
830  int MuonPngID = kMuonPngIdx(sr);
831  int Prong1ID = Prong1ID_Generator(sr, id_type);
832  int EMPngID = 0;
833  float EMcvn = -5.0;
834  float tempcvn = -5.0;
835  for(int i = 0; i < ProngNum; ++i) {
836  if(i == MuonPngID || i == Prong1ID) continue;
837  if(id_type == "4view"){
838  tempcvn = sr->vtx.elastic.fuzzyk.png[i].cvnpart.emid;
839  }
840  else if(id_type == "emid"){
841  tempcvn = sr->vtx.elastic.fuzzyk.png[i].spprongcvnpartnumuccemid.emid;
842  }
843  else if(id_type == "5label"){
844  tempcvn = sr->vtx.elastic.fuzzyk.png[i].spprongcvnpart5label.emid;
845  }
846  else {
847  tempcvn = -5.0;
848  }
849  if(tempcvn > EMcvn) {
850  EMPngID = i;
851  EMcvn = tempcvn;
852  }
853  }
854  return EMPngID;
855  }
856  //=======================================================================
857 
858  //====================================================================
859  // Prong1 / Prong2 index
860  //====================================================================
861  const Var kProng1ID_4view([](const caf::SRProxy* sr){
862  return Prong1ID_Generator(sr, "4view");
863  });
864  const Var kProng2ID_4view([](const caf::SRProxy* sr){
865  return Prong2ID_Generator(sr, "4view");
866  });
867 
868  const Var kProng1ID_emid([](const caf::SRProxy* sr){
869  return Prong1ID_Generator(sr, "emid");
870  });
871  const Var kProng2ID_emid([](const caf::SRProxy* sr){
872  return Prong2ID_Generator(sr, "emid");
873  });
874 
875  const Var kProng1ID_5label([](const caf::SRProxy* sr){
876  return Prong1ID_Generator(sr, "5label");
877  });
878  const Var kProng2ID_5label([](const caf::SRProxy* sr){
879  return Prong2ID_Generator(sr, "5label");
880  });
881  //=======================================================
882 
883 
884  // Prong NHits (total, x-view, y-view)
885  // === 4view === //
886  const Var kProng1NHit_4view([](const caf::SRProxy* sr)
887  {
888  int png1 = kProng1ID_4view(sr);
889  if(png1 < 0) return -5;
890  return (int)sr->vtx.elastic.fuzzyk.png[png1].nhit;
891  });
892  const Var kProng1NHitX_4view([](const caf::SRProxy* sr)
893  {
894  int png1 = kProng1ID_4view(sr);
895  if(png1 < 0) return -5;
896  return (int)sr->vtx.elastic.fuzzyk.png[png1].nhitx;
897  });
898  const Var kProng1NHitY_4view([](const caf::SRProxy* sr)
899  {
900  int png1 = kProng1ID_4view(sr);
901  if(png1 < 0) return -5;
902  return (int)sr->vtx.elastic.fuzzyk.png[png1].nhity;
903  });
904  const Var kProng2NHit_4view([](const caf::SRProxy* sr)
905  {
906  int png2 = kProng2ID_4view(sr);
907  if(png2 < 0) return -5;
908  return (int)sr->vtx.elastic.fuzzyk.png[png2].nhit;
909  });
910  const Var kProng2NHitX_4view([](const caf::SRProxy* sr)
911  {
912  int png2 = kProng2ID_4view(sr);
913  if(png2 < 0) return -5;
914  return (int)sr->vtx.elastic.fuzzyk.png[png2].nhitx;
915  });
916  const Var kProng2NHitY_4view([](const caf::SRProxy* sr)
917  {
918  int png2 = kProng2ID_4view(sr);
919  if(png2 < 0) return -5;
920  return (int)sr->vtx.elastic.fuzzyk.png[png2].nhity;
921  });
922 
923  // === emid === //
924  const Var kProng1NHit_emid([](const caf::SRProxy* sr)
925  {
926  int png1 = kProng1ID_emid(sr);
927  if(png1 < 0) return -5;
928  return (int)sr->vtx.elastic.fuzzyk.png[png1].nhit;
929  });
930  const Var kProng1NHitX_emid([](const caf::SRProxy* sr)
931  {
932  int png1 = kProng1ID_emid(sr);
933  if(png1 < 0) return -5;
934  return (int)sr->vtx.elastic.fuzzyk.png[png1].nhitx;
935  });
936  const Var kProng1NHitY_emid([](const caf::SRProxy* sr)
937  {
938  int png1 = kProng1ID_emid(sr);
939  if(png1 < 0) return -5;
940  return (int)sr->vtx.elastic.fuzzyk.png[png1].nhity;
941  });
942  const Var kProng2NHit_emid([](const caf::SRProxy* sr)
943  {
944  int png2 = kProng2ID_emid(sr);
945  if(png2 < 0) return -5;
946  return (int)sr->vtx.elastic.fuzzyk.png[png2].nhit;
947  });
948  const Var kProng2NHitX_emid([](const caf::SRProxy* sr)
949  {
950  int png2 = kProng2ID_emid(sr);
951  if(png2 < 0) return -5;
952  return (int)sr->vtx.elastic.fuzzyk.png[png2].nhitx;
953  });
954  const Var kProng2NHitY_emid([](const caf::SRProxy* sr)
955  {
956  int png2 = kProng2ID_emid(sr);
957  if(png2 < 0) return -5;
958  return (int)sr->vtx.elastic.fuzzyk.png[png2].nhity;
959  });
960 
961  // === 5label === //
962  const Var kProng1NHit_5label([](const caf::SRProxy* sr)
963  {
964  int png1 = kProng1ID_5label(sr);
965  if(png1 < 0) return -5;
966  return (int)sr->vtx.elastic.fuzzyk.png[png1].nhit;
967  });
968  const Var kProng1NHitX_5label([](const caf::SRProxy* sr)
969  {
970  int png1 = kProng1ID_5label(sr);
971  if(png1 < 0) return -5;
972  return (int)sr->vtx.elastic.fuzzyk.png[png1].nhitx;
973  });
974  const Var kProng1NHitY_5label([](const caf::SRProxy* sr)
975  {
976  int png1 = kProng1ID_5label(sr);
977  if(png1 < 0) return -5;
978  return (int)sr->vtx.elastic.fuzzyk.png[png1].nhity;
979  });
980  const Var kProng2NHit_5label([](const caf::SRProxy* sr)
981  {
982  int png2 = kProng2ID_5label(sr);
983  if(png2 < 0) return -5;
984  return (int)sr->vtx.elastic.fuzzyk.png[png2].nhit;
985  });
986  const Var kProng2NHitX_5label([](const caf::SRProxy* sr)
987  {
988  int png2 = kProng2ID_5label(sr);
989  if(png2 < 0) return -5;
990  return (int)sr->vtx.elastic.fuzzyk.png[png2].nhitx;
991  });
992  const Var kProng2NHitY_5label([](const caf::SRProxy* sr)
993  {
994  int png2 = kProng2ID_5label(sr);
995  if(png2 < 0) return -5;
996  return (int)sr->vtx.elastic.fuzzyk.png[png2].nhity;
997  });
998 
999  //====================================================================
1000  // Prong1/Prong2 CVN Scores
1001  //====================================================================
1002  const Var kProng1Score_4view([](const caf::SRProxy* sr)
1003  {
1004  int png1 = kProng1ID_4view(sr);
1005  if(png1 < 0) return -5.0f;
1006  const auto& p1 = sr->vtx.elastic.fuzzyk.png[png1];
1007  return (float)(p1.cvnpart.emid);
1008  });
1009  const Var kProng2Score_4view([](const caf::SRProxy* sr)
1010  {
1011  int png2 = kProng2ID_4view(sr);
1012  if(png2 < 0) return -5.0f;
1013  const auto& p2 = sr->vtx.elastic.fuzzyk.png[png2];
1014  return (float)(p2.cvnpart.emid);
1015  });
1016 
1017  const Var kProng1Score_emid([](const caf::SRProxy* sr)
1018  {
1019  int png1 = kProng1ID_emid(sr);
1020  if(png1 < 0) return -5.0f;
1021  const auto& p1 = sr->vtx.elastic.fuzzyk.png[png1];
1022  return (float)(p1.spprongcvnpartnumuccemid.emid);
1023  });
1024  const Var kProng2Score_emid([](const caf::SRProxy* sr)
1025  {
1026  int png2 = kProng2ID_emid(sr);
1027  if(png2 < 0) return -5.0f;
1028  const auto& p2 = sr->vtx.elastic.fuzzyk.png[png2];
1029  return (float)(p2.spprongcvnpartnumuccemid.emid);
1030  });
1031 
1032  const Var kProng1Score_5label([](const caf::SRProxy* sr)
1033  {
1034  int png1 = kProng1ID_5label(sr);
1035  if(png1 < 0) return -5.0f;
1036  const auto& p1 = sr->vtx.elastic.fuzzyk.png[png1];
1037  return (float)(p1.spprongcvnpart5label.emid);
1038  });
1039  const Var kProng2Score_5label([](const caf::SRProxy* sr)
1040  {
1041  int png2 = kProng2ID_5label(sr);
1042  if(png2 < 0) return -5.0f;
1043  const auto& p2 = sr->vtx.elastic.fuzzyk.png[png2];
1044  return (float)(p2.spprongcvnpart5label.emid);
1045  });
1046 
1047 
1048  // Prong Truth PDG
1049  // === 4view === //
1050  const Var kProng1TruePDG_4view([](const caf::SRProxy* sr)
1051  {
1052  int png1 = kProng1ID_4view(sr);
1053  if(png1 < 0) return 0;
1054  return (int)sr->vtx.elastic.fuzzyk.png[png1].truth.pdg;
1055  });
1056  const Var kProng2TruePDG_4view([](const caf::SRProxy* sr)
1057  {
1058  int png2 = kProng2ID_4view(sr);
1059  if(png2 < 0) return 0;
1060  return (int)sr->vtx.elastic.fuzzyk.png[png2].truth.pdg;
1061  });
1062 
1063  // === emid === //
1064  const Var kProng1TruePDG_emid([](const caf::SRProxy* sr)
1065  {
1066  int png1 = kProng1ID_emid(sr);
1067  if(png1 < 0) return 0;
1068  return (int)sr->vtx.elastic.fuzzyk.png[png1].truth.pdg;
1069  });
1070  const Var kProng2TruePDG_emid([](const caf::SRProxy* sr)
1071  {
1072  int png2 = kProng2ID_emid(sr);
1073  if(png2 < 0) return 0;
1074  return (int)sr->vtx.elastic.fuzzyk.png[png2].truth.pdg;
1075  });
1076 
1077  //=================================================================
1078  // Prong1/Prong2 5-label Proton/Photon/Electron Scores
1079  //=================================================================
1080  const Var kProng1ProtonScore_emid([](const caf::SRProxy* sr)
1081  {
1082  int png1 = kProng1ID_emid(sr);
1083  if(png1 < 0) return -5.0f;
1084  const auto& p1 = sr->vtx.elastic.fuzzyk.png[png1];
1085  return (float)(p1.spprongcvnpart5label.protonid);
1086  });
1087  const Var kProng2ProtonScore_emid([](const caf::SRProxy* sr)
1088  {
1089  int png2 = kProng2ID_emid(sr);
1090  if(png2 < 0) return -5.0f;
1091  const auto& p2 = sr->vtx.elastic.fuzzyk.png[png2];
1092  return (float)(p2.spprongcvnpart5label.protonid);
1093  });
1094 
1095  const Var kProng1PhotonScore_emid([](const caf::SRProxy* sr)
1096  {
1097  int png1 = kProng1ID_emid(sr);
1098  if(png1 < 0) return -5.0f;
1099  const auto& p1 = sr->vtx.elastic.fuzzyk.png[png1];
1100  return (float)(p1.spprongcvnpart5label.photonid);
1101  });
1102  const Var kProng2PhotonScore_emid([](const caf::SRProxy* sr)
1103  {
1104  int png2 = kProng2ID_emid(sr);
1105  if(png2 < 0) return -5.0f;
1106  const auto& p2 = sr->vtx.elastic.fuzzyk.png[png2];
1107  return (float)(p2.spprongcvnpart5label.photonid);
1108  });
1109 
1110  const Var kProng1ElectronScore_emid([](const caf::SRProxy* sr)
1111  {
1112  int png1 = kProng1ID_emid(sr);
1113  if(png1 < 0) return -5.0f;
1114  const auto& p1 = sr->vtx.elastic.fuzzyk.png[png1];
1115  return (float)(p1.spprongcvnpart5label.electronid);
1116  });
1117  const Var kProng2ElectronScore_emid([](const caf::SRProxy* sr)
1118  {
1119  int png2 = kProng2ID_emid(sr);
1120  if(png2 < 0) return -5.0f;
1121  const auto& p2 = sr->vtx.elastic.fuzzyk.png[png2];
1122  return (float)(p2.spprongcvnpart5label.electronid);
1123  });
1124 
1125 
1126  //=========================//
1127  // Prong Variables //
1128  //=========================//
1129  // === emid === //
1130 
1131  // Prong efficiency / purity
1132  const Var kProng1TruthEff_emid([](const caf::SRProxy* sr)
1133  {
1134  int png1 = kProng1ID_emid(sr);
1135  if(png1 < 0) return -5.0f;
1136  return (float)sr->vtx.elastic.fuzzyk.png[png1].truth.eff;
1137  });
1138  const Var kProng2TruthEff_emid([](const caf::SRProxy* sr)
1139  {
1140  int png2 = kProng2ID_emid(sr);
1141  if(png2 < 0) return -5.0f;
1142  return (float)sr->vtx.elastic.fuzzyk.png[png2].truth.eff;
1143  });
1144 
1145  const Var kProng1TruthPur_emid([](const caf::SRProxy* sr)
1146  {
1147  int png1 = kProng1ID_emid(sr);
1148  if(png1 < 0) return -5.0f;
1149  return (float)sr->vtx.elastic.fuzzyk.png[png1].truth.pur;
1150  });
1151  const Var kProng2TruthPur_emid([](const caf::SRProxy* sr)
1152  {
1153  int png2 = kProng2ID_emid(sr);
1154  if(png2 < 0) return -5.0f;
1155  return (float)sr->vtx.elastic.fuzzyk.png[png2].truth.pur;
1156  });
1157 
1158  // Prong energy
1159  const Var kProng1CalE_emid([](const caf::SRProxy* sr)
1160  {
1161  int png1 = kProng1ID_emid(sr);
1162  if(png1 < 0) return -5.0f;
1163  return (float)sr->vtx.elastic.fuzzyk.png[png1].calE;
1164  });
1165  const Var kProng2CalE_emid([](const caf::SRProxy* sr)
1166  {
1167  int png2 = kProng2ID_emid(sr);
1168  if(png2 < 0) return -5.0f;
1169  return (float)sr->vtx.elastic.fuzzyk.png[png2].calE;
1170  });
1171 
1172  const Var kProng1TruthE_emid([](const caf::SRProxy* sr)
1173  {
1174  int png1 = kProng1ID_emid(sr);
1175  if(png1 < 0) return -5.0f;
1176  return (float)sr->vtx.elastic.fuzzyk.png[png1].truth.p.E;
1177  });
1178  const Var kProng2TruthE_emid([](const caf::SRProxy* sr)
1179  {
1180  int png2 = kProng2ID_emid(sr);
1181  if(png2 < 0) return -5.0f;
1182  return (float)sr->vtx.elastic.fuzzyk.png[png2].truth.p.E;
1183  });
1184 
1185  const Var kProng1E_Reso_emid([](const caf::SRProxy* sr)
1186  {
1187  int png1 = kProng1ID_emid(sr);
1188  if(png1 < 0) return -5.0f;
1189  float calE = sr->vtx.elastic.fuzzyk.png[png1].calE;
1190  float trueE = sr->vtx.elastic.fuzzyk.png[png1].truth.p.E;
1191  if (trueE == 0) return -5.0f;
1192  return (float)(calE - trueE)/trueE;
1193  });
1194  const Var kProng2E_Reso_emid([](const caf::SRProxy* sr)
1195  {
1196  int png2 = kProng2ID_emid(sr);
1197  if(png2 < 0) return -5.0f;
1198  float calE = sr->vtx.elastic.fuzzyk.png[png2].calE;
1199  float trueE = sr->vtx.elastic.fuzzyk.png[png2].truth.p.E;
1200  if (trueE == 0) return -5.0f;
1201  return (float)(calE - trueE)/trueE;
1202  });
1203 
1204  // Prong direction
1205  const Var kProng1DirZ_Cos_emid([](const caf::SRProxy* sr)
1206  {
1207  int png1 = kProng1ID_emid(sr);
1208  if(png1 < 0) return -5.0f;
1209  auto& pngDir = sr->vtx.elastic.fuzzyk.png[png1].dir;
1210  TVector3 beamDirND = NuMIBeamDirection(caf::kNEARDET);
1211  return (float)pngDir.Unit().Dot(beamDirND.Unit());
1212  });
1213  const Var kProng2DirZ_Cos_emid([](const caf::SRProxy* sr)
1214  {
1215  int png2 = kProng2ID_emid(sr);
1216  if(png2 < 0) return -5.0f;
1217  auto& pngDir = sr->vtx.elastic.fuzzyk.png[png2].dir;
1218  TVector3 beamDirND = NuMIBeamDirection(caf::kNEARDET);
1219  return (float)pngDir.Unit().Dot(beamDirND.Unit());
1220  });
1221  const Var kProng1DirZ_Rad_emid([](const caf::SRProxy* sr)
1222  {
1223  int png1 = kProng1ID_emid(sr);
1224  if(png1 < 0) return -5.0f;
1225  auto& pngDir = sr->vtx.elastic.fuzzyk.png[png1].dir;
1226  TVector3 beamDirND = NuMIBeamDirection(caf::kNEARDET);
1227  return (float)pngDir.Unit().Angle(beamDirND.Unit());
1228  });
1229  const Var kProng2DirZ_Rad_emid([](const caf::SRProxy* sr)
1230  {
1231  int png2 = kProng2ID_emid(sr);
1232  if(png2 < 0) return -5.0f;
1233  auto& pngDir = sr->vtx.elastic.fuzzyk.png[png2].dir;
1234  TVector3 beamDirND = NuMIBeamDirection(caf::kNEARDET);
1235  return (float)pngDir.Unit().Angle(beamDirND.Unit());
1236  });
1237 
1238  const Var kProng1TruthDirZ_Cos_emid([](const caf::SRProxy* sr)
1239  {
1240  int png1 = kProng1ID_emid(sr);
1241  if(png1 < 0) return -5.0f;
1242  TVector3 pngDir = sr->vtx.elastic.fuzzyk.png[png1].truth.p.Vect();
1243  TVector3 beamDirND = NuMIBeamDirection(caf::kNEARDET);
1244  return (float)pngDir.Unit().Dot(beamDirND.Unit());
1245  });
1246  const Var kProng2TruthDirZ_Cos_emid([](const caf::SRProxy* sr)
1247  {
1248  int png2 = kProng2ID_emid(sr);
1249  if(png2 < 0) return -5.0f;
1250  TVector3 pngDir = sr->vtx.elastic.fuzzyk.png[png2].truth.p.Vect();
1251  TVector3 beamDirND = NuMIBeamDirection(caf::kNEARDET);
1252  return (float)pngDir.Unit().Dot(beamDirND.Unit());
1253  });
1254  const Var kProng1TruthDirZ_Rad_emid([](const caf::SRProxy* sr)
1255  {
1256  int png1 = kProng1ID_emid(sr);
1257  if(png1 < 0) return -5.0f;
1258  TVector3 pngDir = sr->vtx.elastic.fuzzyk.png[png1].truth.p.Vect();
1259  TVector3 beamDirND = NuMIBeamDirection(caf::kNEARDET);
1260  return (float)pngDir.Unit().Angle(beamDirND.Unit());
1261  });
1262  const Var kProng2TruthDirZ_Rad_emid([](const caf::SRProxy* sr)
1263  {
1264  int png2 = kProng2ID_emid(sr);
1265  if(png2 < 0) return -5.0f;
1266  TVector3 pngDir = sr->vtx.elastic.fuzzyk.png[png2].truth.p.Vect();
1267  TVector3 beamDirND = NuMIBeamDirection(caf::kNEARDET);
1268  return (float)pngDir.Unit().Angle(beamDirND.Unit());
1269  });
1270 
1271  const Var kProng1DirZ_Diff_Cos_emid([](const caf::SRProxy* sr)
1272  {
1273  int png1 = kProng1ID_emid(sr);
1274  if(png1 < 0) return -5.0f;
1275  auto& pngDir = sr->vtx.elastic.fuzzyk.png[png1].dir;
1276  TVector3 truthDir = sr->vtx.elastic.fuzzyk.png[png1].truth.p.Vect();
1277  return (float)pngDir.Unit().Dot(truthDir.Unit());
1278  });
1279  const Var kProng2DirZ_Diff_Cos_emid([](const caf::SRProxy* sr)
1280  {
1281  int png2 = kProng2ID_emid(sr);
1282  if(png2 < 0) return -5.0f;
1283  auto& pngDir = sr->vtx.elastic.fuzzyk.png[png2].dir;
1284  TVector3 truthDir = sr->vtx.elastic.fuzzyk.png[png2].truth.p.Vect();
1285  return (float)pngDir.Unit().Dot(truthDir.Unit());
1286  });
1287  const Var kProng1DirZ_Diff_Rad_emid([](const caf::SRProxy* sr)
1288  {
1289  int png1 = kProng1ID_emid(sr);
1290  if(png1 < 0) return -5.0f;
1291  auto& pngDir = sr->vtx.elastic.fuzzyk.png[png1].dir;
1292  TVector3 truthDir = sr->vtx.elastic.fuzzyk.png[png1].truth.p.Vect();
1293  return (float)pngDir.Unit().Angle(truthDir.Unit());
1294  });
1295  const Var kProng2DirZ_Diff_Rad_emid([](const caf::SRProxy* sr)
1296  {
1297  int png2 = kProng2ID_emid(sr);
1298  if(png2 < 0) return -5.0f;
1299  auto& pngDir = sr->vtx.elastic.fuzzyk.png[png2].dir;
1300  TVector3 truthDir = sr->vtx.elastic.fuzzyk.png[png2].truth.p.Vect();
1301  return (float)pngDir.Unit().Angle(truthDir.Unit());
1302  });
1303  const Var kProng1DirZ_Diff_Deg_emid([](const caf::SRProxy* sr)
1304  {
1305  float Prong1DirZ_Diff_Rad = kProng1DirZ_Diff_Rad_emid(sr);
1306  if(Prong1DirZ_Diff_Rad > TMath::Pi() || Prong1DirZ_Diff_Rad < 0.) return -5.0f;
1307  return (float)(Prong1DirZ_Diff_Rad*180./TMath::Pi());
1308  });
1309  const Var kProng2DirZ_Diff_Deg_emid([](const caf::SRProxy* sr)
1310  {
1311  float Prong2DirZ_Diff_Rad = kProng2DirZ_Diff_Rad_emid(sr);
1312  if(Prong2DirZ_Diff_Rad > TMath::Pi() || Prong2DirZ_Diff_Rad < 0.) return -5.0f;
1313  return (float)(Prong2DirZ_Diff_Rad*180./TMath::Pi());
1314  });
1315 
1316  // Two prong opening angle
1317  const Var kTwoProngOpenAngle_Cos_emid([](const caf::SRProxy* sr)
1318  {
1319  int png1 = kProng1ID_emid(sr);
1320  int png2 = kProng2ID_emid(sr);
1321  if(png1 < 0 || png2 < 0) return -5.0f;
1322  auto& png1Dir = sr->vtx.elastic.fuzzyk.png[png1].dir;
1323  auto& png2Dir = sr->vtx.elastic.fuzzyk.png[png2].dir;
1324  return (float)png1Dir.Unit().Dot(png2Dir.Unit());
1325  });
1326  const Var kTwoProngOpenAngle_Rad_emid([](const caf::SRProxy* sr)
1327  {
1328  int png1 = kProng1ID_emid(sr);
1329  int png2 = kProng2ID_emid(sr);
1330  if(png1 < 0 || png2 < 0) return -5.0f;
1331  auto& png1Dir = sr->vtx.elastic.fuzzyk.png[png1].dir;
1332  auto& png2Dir = sr->vtx.elastic.fuzzyk.png[png2].dir;
1333  return (float)png1Dir.Unit().Angle(png2Dir.Unit());
1334  });
1335  const Var kTwoProngOpenAngle_Deg_emid([](const caf::SRProxy* sr)
1336  {
1337  float TwoProngOpenAngle_Rad = kTwoProngOpenAngle_Rad_emid(sr);
1338  if(TwoProngOpenAngle_Rad > TMath::Pi() || TwoProngOpenAngle_Rad < 0.) return -5.0f;
1339  return (float)(TwoProngOpenAngle_Rad*180./TMath::Pi());
1340  });
1341 
1343  {
1344  int png1 = kProng1ID_emid(sr);
1345  int png2 = kProng2ID_emid(sr);
1346  if(png1 < 0 || png2 < 0) return -5.0f;
1347  TVector3 png1Dir = sr->vtx.elastic.fuzzyk.png[png1].truth.p.Vect();
1348  TVector3 png2Dir = sr->vtx.elastic.fuzzyk.png[png2].truth.p.Vect();
1349  return (float)png1Dir.Unit().Dot(png2Dir.Unit());
1350  });
1352  {
1353  int png1 = kProng1ID_emid(sr);
1354  int png2 = kProng2ID_emid(sr);
1355  if(png1 < 0 || png2 < 0) return -5.0f;
1356  TVector3 png1Dir = sr->vtx.elastic.fuzzyk.png[png1].truth.p.Vect();
1357  TVector3 png2Dir = sr->vtx.elastic.fuzzyk.png[png2].truth.p.Vect();
1358  return (float)png1Dir.Unit().Angle(png2Dir.Unit());
1359  });
1361  {
1362  float TwoProngTruthOpenAngle_Rad = kTwoProngTruthOpenAngle_Rad_emid(sr);
1363  if(TwoProngTruthOpenAngle_Rad > TMath::Pi() || TwoProngTruthOpenAngle_Rad < 0.) return -5.0f;
1364  return (float)(TwoProngTruthOpenAngle_Rad*180./TMath::Pi());
1365  });
1366 
1368  {
1369  float openAngle = kTwoProngOpenAngle_Rad_emid(sr);
1370  float truthOpenAngle = kTwoProngTruthOpenAngle_Rad_emid(sr);
1371  if(openAngle < 0. || openAngle > TMath::Pi() || truthOpenAngle < 0. || truthOpenAngle > TMath::Pi()) return -5.0f;
1372  return (float)(openAngle - truthOpenAngle);
1373  });
1374 
1375 
1376  //=================================================================
1377  //Reco pi0 invariant mass
1378  //=================================================================
1379  const Var kPi0Mass_4view([](const caf::SRProxy* sr)
1380  {
1381  int png1 = kProng1ID_4view(sr);
1382  int png2 = kProng2ID_4view(sr);
1383  if(png1 < 0 || png2 < 0) return -5.0f;
1384  const auto& p1 = sr->vtx.elastic.fuzzyk.png[png1];
1385  const auto& p2 = sr->vtx.elastic.fuzzyk.png[png2];
1386  const float dot = p1.dir.Unit().Dot(p2.dir.Unit());
1387  return 1000*sqrt(2*p1.calE*p2.calE*(1-dot));
1388  });
1389 
1390  const Var kPi0Mass_emid([](const caf::SRProxy* sr)
1391  {
1392  int png1 = kProng1ID_emid(sr);
1393  int png2 = kProng2ID_emid(sr);
1394  if(png1 < 0 || png2 < 0) return -5.0f;
1395  const auto& p1 = sr->vtx.elastic.fuzzyk.png[png1];
1396  const auto& p2 = sr->vtx.elastic.fuzzyk.png[png2];
1397  const float dot = p1.dir.Unit().Dot(p2.dir.Unit());
1398  return 1000*sqrt(2*p1.calE*p2.calE*(1-dot));
1399  });
1400 
1401  const Var kPi0Mass_5label([](const caf::SRProxy* sr)
1402  {
1403  int png1 = kProng1ID_5label(sr);
1404  int png2 = kProng2ID_5label(sr);
1405  if(png1 < 0 || png2 < 0) return -5.0f;
1406  const auto& p1 = sr->vtx.elastic.fuzzyk.png[png1];
1407  const auto& p2 = sr->vtx.elastic.fuzzyk.png[png2];
1408  const float dot = p1.dir.Unit().Dot(p2.dir.Unit());
1409  return 1000*sqrt(2*p1.calE*p2.calE*(1-dot));
1410  });
1411  //=================================================================
1412  //Reco pi0 kinetic variables
1413  //=================================================================
1414  // === 4view === //
1415  const Var kRecoPi0P_4view([](const caf::SRProxy* sr)
1416  {
1417  int png1 = kProng1ID_4view(sr);
1418  int png2 = kProng2ID_4view(sr);
1419  if(png1 < 0 || png2 < 0) return -5.0f;
1420  const auto& p1 = sr->vtx.elastic.fuzzyk.png[png1];
1421  const auto& p2 = sr->vtx.elastic.fuzzyk.png[png2];
1422  TVector3 p = p1.calE*p1.dir.Unit()+p2.calE*p2.dir.Unit();
1423  return (float)p.Mag();
1424  });
1425 
1426  const Var kRecoPi0DirZ_Cos_4view([](const caf::SRProxy* sr)
1427  {
1428  int png1 = kProng1ID_4view(sr);
1429  int png2 = kProng2ID_4view(sr);
1430  if(png1 < 0 || png2 < 0) return -5.0f;
1431  const auto& p1 = sr->vtx.elastic.fuzzyk.png[png1];
1432  const auto& p2 = sr->vtx.elastic.fuzzyk.png[png2];
1433  TVector3 p = p1.calE*p1.dir.Unit() + p2.calE*p2.dir.Unit();
1434  TVector3 beamDirND = NuMIBeamDirection(caf::kNEARDET);
1435  return (float)p.Unit().Dot(beamDirND.Unit());
1436  });
1437 
1438  const Var kRecoPi0DirZ_Rad_4view([](const caf::SRProxy* sr)
1439  {
1440  int png1 = kProng1ID_4view(sr);
1441  int png2 = kProng2ID_4view(sr);
1442  if(png1 < 0 || png2 < 0) return -5.0f;
1443  const auto& p1 = sr->vtx.elastic.fuzzyk.png[png1];
1444  const auto& p2 = sr->vtx.elastic.fuzzyk.png[png2];
1445  TVector3 p = p1.calE*p1.dir.Unit() + p2.calE*p2.dir.Unit();
1446  TVector3 beamDirND = NuMIBeamDirection(caf::kNEARDET);
1447  return (float)p.Unit().Angle(beamDirND.Unit());
1448  });
1449 
1450  const Var kRecoPi0DirZ_Deg_4view([](const caf::SRProxy* sr)
1451  {
1452  float RecoPi0DirZ_Rad = kRecoPi0DirZ_Rad_4view(sr);
1453  if(RecoPi0DirZ_Rad > TMath::Pi() || RecoPi0DirZ_Rad < 0.) return -5.0f;
1454  return (float)(RecoPi0DirZ_Rad*180./TMath::Pi());
1455  });
1456 
1457  // === emid === //
1458  const Var kRecoPi0P_1_emid([](const caf::SRProxy* sr)
1459  {
1460  int png1 = kProng1ID_emid(sr);
1461  int png2 = kProng2ID_emid(sr);
1462  if(png1 < 0 || png2 < 0) return -5.0f;
1463  const auto& p1 = sr->vtx.elastic.fuzzyk.png[png1];
1464  const auto& p2 = sr->vtx.elastic.fuzzyk.png[png2];
1465  TVector3 p = p1.calE*p1.dir.Unit()+p2.calE*p2.dir.Unit();
1466  return (float)p.Mag();
1467  });
1468  const Var kRecoPi0P_2_emid([](const caf::SRProxy* sr)
1469  {
1470  int png1 = kProng1ID_emid(sr);
1471  int png2 = kProng2ID_emid(sr);
1472  if(png1 < 0 || png2 < 0) return -5.0f;
1473  const auto& p1 = sr->vtx.elastic.fuzzyk.png[png1];
1474  const auto& p2 = sr->vtx.elastic.fuzzyk.png[png2];
1475  float pi0E = p1.calE + p2.calE;
1476  if(pi0E < Pi0Mass()) return -1.0f;
1477  return (float)(sqrt(std::pow(pi0E, 2) - std::pow(Pi0Mass(), 2)));
1478  });
1479 
1480  const Var kRecoPi0DirZ_Cos_emid([](const caf::SRProxy* sr)
1481  {
1482  int png1 = kProng1ID_emid(sr);
1483  int png2 = kProng2ID_emid(sr);
1484  if(png1 < 0 || png2 < 0) return -5.0f;
1485  const auto& p1 = sr->vtx.elastic.fuzzyk.png[png1];
1486  const auto& p2 = sr->vtx.elastic.fuzzyk.png[png2];
1487  TVector3 p = p1.calE*p1.dir.Unit() + p2.calE*p2.dir.Unit();
1488  TVector3 beamDirND = NuMIBeamDirection(caf::kNEARDET);
1489  return (float)p.Unit().Dot(beamDirND.Unit());
1490  });
1491 
1492  const Var kRecoPi0DirZ_Rad_emid([](const caf::SRProxy* sr)
1493  {
1494  int png1 = kProng1ID_emid(sr);
1495  int png2 = kProng2ID_emid(sr);
1496  if(png1 < 0 || png2 < 0) return -5.0f;
1497  const auto& p1 = sr->vtx.elastic.fuzzyk.png[png1];
1498  const auto& p2 = sr->vtx.elastic.fuzzyk.png[png2];
1499  TVector3 p = p1.calE*p1.dir.Unit() + p2.calE*p2.dir.Unit();
1500  TVector3 beamDirND = NuMIBeamDirection(caf::kNEARDET);
1501  return (float)p.Unit().Angle(beamDirND.Unit());
1502  });
1503 
1504  const Var kRecoPi0DirZ_Deg_emid([](const caf::SRProxy* sr)
1505  {
1506  float RecoPi0DirZ_Rad = kRecoPi0DirZ_Rad_emid(sr);
1507  if(RecoPi0DirZ_Rad > TMath::Pi() || RecoPi0DirZ_Rad < 0.) return -5.0f;
1508  return (float)(RecoPi0DirZ_Rad*180./TMath::Pi());
1509  });
1510 
1511 
1512 
1513  const Var kRecoMuonP
1514  ([](const caf::SRProxy* sr)
1515  {
1516  int png_mu = kMuonPngIdx(sr);
1517  if(png_mu < 0) return -5.0f;
1518  const auto& png = sr->vtx.elastic.fuzzyk.png[png_mu];
1519  TVector3 p = png.calE*png.dir.Unit();
1520  return (float)p.Mag();
1521  });
1522 
1523  const Var kRecoMuonDirZ_Rad
1524  ([](const caf::SRProxy* sr)
1525  {
1526  int png_mu = kMuonPngIdx(sr);
1527  if(png_mu < 0) return -5.0f;
1528  const auto& png = sr->vtx.elastic.fuzzyk.png[png_mu];
1529  TVector3 p = png.calE*png.dir.Unit();
1530  TVector3 beamDirND = NuMIBeamDirection(caf::kNEARDET);
1531  return (float)p.Unit().Angle(beamDirND.Unit());
1532  });
1533 
1534  const Var kRecoMuonDirZ_Cos
1535  ([](const caf::SRProxy* sr)
1536  {
1537  int png_mu = kMuonPngIdx(sr);
1538  if(png_mu < 0) return -5.0f;
1539  const auto& png = sr->vtx.elastic.fuzzyk.png[png_mu];
1540  TVector3 p = png.calE*png.dir.Unit();
1541  TVector3 beamDirND = NuMIBeamDirection(caf::kNEARDET);
1542  return (float)p.Unit().Dot(beamDirND.Unit());
1543  });
1544 
1545  const Var kRecoMuonDirZ_Deg
1546  ([](const caf::SRProxy* sr)
1547  {
1548  float RecoMuonDirZ_Rad = kRecoMuonDirZ_Rad(sr);
1549  if(RecoMuonDirZ_Rad > TMath::Pi() || RecoMuonDirZ_Rad < 0.) return -5.0f;
1550  return (float)(RecoMuonDirZ_Rad*180./TMath::Pi());
1551  });
1552 
1553  const Var kRecoQ2
1554  ([](const caf::SRProxy* sr)
1555  {
1556  float nuE = kRecoE(sr);
1557  float muE = kIncXsecMuonE(sr);
1558  float pmu = kRecoMuonP(sr);
1559  float cosmu = kRecoMuonDirZ_Cos(sr);
1560 
1561  if(nuE < 0 || muE < 0 || pmu < 0 || cosmu == -5.0 ) return -5.0f;
1562 
1563  return (float)(2 * nuE * (muE - pmu * cosmu));
1564  });
1565 
1566  const Var kRecoW2
1567  ([](const caf::SRProxy* sr)
1568  {
1569  float nuE = kRecoE(sr);
1570  float muE = kIncXsecMuonE(sr);
1571  float q2 = kRecoQ2(sr);
1572 
1573  if(nuE < 0 || muE < 0 || q2 < 0) return -5.0f;
1574 
1575  float mN = (sr->spill.isFHC)? kNeutronMass(sr) : kProtonMass(sr);
1576 
1577  if((pow(mN,2) + (2 * mN * (nuE - muE)) - q2) < 0) return -5.0f;
1578 
1579  return (float)(pow(mN,2) + (2 * mN * (nuE - muE)) - q2);
1580  });
1581 
1582  const Var kRecoW
1583  ([](const caf::SRProxy* sr)
1584  {
1585  float nuE = kRecoE(sr);
1586  float muE = kIncXsecMuonE(sr);
1587  float q2 = kRecoQ2(sr);
1588 
1589  if(nuE < 0 || muE < 0 || q2 < 0) return -5.0f;
1590 
1591  float mN = (sr->spill.isFHC)? kNeutronMass(sr) : kProtonMass(sr);
1592 
1593  float W2 = (pow(mN,2) + (2 * mN * (nuE - muE)) - q2);
1594  if( W2 < 0) return -5.0f;
1595 
1596  return (float)( sqrt(W2) );
1597  });
1598 
1599  ///////////////////////////////////////////
1600  //Resolution Vars
1601  ///////////////////////////////////////////
1602  const Var kNuE_Diff
1603  ([] (const caf::SRProxy* sr)
1604  {
1605  float truth = kTrueNuE(sr);
1606  float reco = kRecoE(sr);
1607 
1608 
1609  if(truth < 0 || reco < 0 ) return (-5.0f);
1610 
1611  return (float)(reco - truth);
1612  });
1613 
1614  const Var kQ2_Diff
1615  ([] (const caf::SRProxy* sr)
1616  {
1617  float truth = kTrueNuQ2(sr);
1618  float reco = kRecoQ2(sr);
1619 
1620  if(truth < 0 || reco < 0 ) return (-5.0f);
1621 
1622 
1623  return (float)(reco - truth);
1624  });
1625 
1626  const Var kW_Diff
1627  ([] (const caf::SRProxy* sr)
1628  {
1629  float truth = kTrueNuW(sr);
1630  float reco = kRecoW(sr);
1631 
1632  if(truth < 0 || reco < 0 ) return (-5.0f);
1633 
1634  return (float)(reco - truth);
1635  });
1636 
1637  const Var kMuP_Diff
1638  ([] (const caf::SRProxy* sr)
1639  {
1640  float truth = kTrueMuP(sr);
1641  float reco = kRecoMuonP(sr);
1642 
1643  if(truth < 0 || reco < 0 ) return (-5.0f);
1644 
1645  return (float)(reco - truth);
1646  });
1647 
1648  const Var kMuDeg_Diff
1649  ([] (const caf::SRProxy* sr)
1650  {
1651  float truth = kTrueMuDirZ_Deg(sr);
1652  float reco = kRecoMuonDirZ_Deg(sr);
1653 
1654  if(truth < 0 || reco < 0 ) return (-5.0f);
1655 
1656  return (float)(reco - truth);
1657  });
1658 
1659  const Var kPi0P_Diff
1660  ([] (const caf::SRProxy* sr)
1661  {
1662  float truth = kTruePrimPi0P(sr);
1663  float reco = kRecoPi0P_1_emid(sr);
1664 
1665  if(truth < 0 || reco < 0 ) return (-5.0f);
1666 
1667  return (float)(reco - truth);
1668  });
1669 
1670  const Var kPi0Deg_Diff
1671  ([] (const caf::SRProxy* sr)
1672  {
1673  float truth = kTruePrimPi0DirZ_Deg(sr);
1674  float reco = kRecoPi0DirZ_Deg_emid(sr);
1675 
1676  if(truth < 0 || reco < 0 ) return (-5.0f);
1677 
1678  return (float)(reco - truth);
1679  });
1680 
1681 
1682  } // end of napespace numubarccpi0
1683 } // end of namespace ana
const XML_Char int len
Definition: expat.h:262
const Var kProng2ID_4view([](const caf::SRProxy *sr){return Prong2ID_Generator(sr,"4view");})
const Var k4viewMuonPng_MuonScore_4view([](const caf::SRProxy *sr){return(float) GetCVNProngMuonScore_4view(kCVNMuonIdx_4view(sr), sr);})
const NuTruthVar kNPrimaryPi0_NT([](const caf::SRNeutrinoProxy *truth){int nprims=truth->prim.size();int countpi0=0;for(int i=0;i< nprims;++i){if(truth->prim[i].pdg==111) countpi0++;}return countpi0;})
caf::Proxy< size_t > npng
Definition: SRProxy.h:2038
caf::Proxy< caf::SRSpill > spill
Definition: SRProxy.h:2143
const Var kProng1NHitX_emid([](const caf::SRProxy *sr){int png1=kProng1ID_emid(sr);if(png1< 0) return-5;return(int) sr->vtx.elastic.fuzzyk.png[png1].nhitx;})
Near Detector underground.
Definition: SREnums.h:10
const Var kProng1Score_emid([](const caf::SRProxy *sr){int png1=kProng1ID_emid(sr);if(png1< 0) return-5.0f;const auto &p1=sr->vtx.elastic.fuzzyk.png[png1];return(float)(p1.spprongcvnpartnumuccemid.emid);})
const Var kNuE_Diff([](const caf::SRProxy *sr){float truth=kTrueNuE(sr);float reco=kRecoE(sr);if(truth< 0||reco< 0) return(-5.0f);return(float)(reco-truth);})
int isinf(const stan::math::var &a)
Definition: std_isinf.hpp:16
const NuTruthVar kTrueMuDirZ_Rad_NT([](const caf::SRNeutrinoProxy *truth){if(abs(truth->pdg)!=14||!truth->iscc) return-5.0;int nprims=truth->prim.size();for(int iprim=0;iprim< nprims;iprim++){if(abs(truth->prim[iprim].pdg)==13){TVector3 mudir=truth->prim[iprim].p.Vect();TVector3 beamdir=NuMIBeamDirection(caf::kNEARDET);return mudir.Unit().Angle(beamdir.Unit());}}return-5.0;})
const NuTruthVar kTrueSecoPi0E_NT([](const caf::SRNeutrinoProxy *truth){int TruePi0_ID1=kTrueSecoPi0_ID1_NT(truth);int TruePi0_ID2=kTrueSecoPi0_ID2_NT(truth);if(TruePi0_ID1< 0||TruePi0_ID2< 0) return-5.0f;return(float) truth->prim[TruePi0_ID1].daughterEnergies[TruePi0_ID2];})
caf::Proxy< caf::SRFuzzyK > fuzzyk
Definition: SRProxy.h:2059
const Var kRecoPi0P_2_emid([](const caf::SRProxy *sr){int png1=kProng1ID_emid(sr);int png2=kProng2ID_emid(sr);if(png1< 0||png2< 0) return-5.0f;const auto &p1=sr->vtx.elastic.fuzzyk.png[png1];const auto &p2=sr->vtx.elastic.fuzzyk.png[png2];float pi0E=p1.calE+p2.calE;if(pi0E< Pi0Mass()) return-1.0f;return(float)(sqrt(std::pow(pi0E, 2)-std::pow(Pi0Mass(), 2)));})
const NuTruthVar kTrueMuP_NT([](const caf::SRNeutrinoProxy *truth){float pmag=-5;if(abs(truth->pdg)!=14||!truth->iscc) return pmag;int nprims=truth->prim.size();for(int iprim=0;iprim< nprims;iprim++){if(abs(truth->prim[iprim].pdg)==13){pmag=truth->prim[iprim].p.Mag();}}return pmag;})
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const Var kRecoPi0DirZ_Deg_4view([](const caf::SRProxy *sr){float RecoPi0DirZ_Rad=kRecoPi0DirZ_Rad_4view(sr);if(RecoPi0DirZ_Rad > TMath::Pi()||RecoPi0DirZ_Rad< 0.) return-5.0f;return(float)(RecoPi0DirZ_Rad *180./TMath::Pi());})
const Var kPi0P_Diff([](const caf::SRProxy *sr){float truth=kTruePrimPi0P(sr);float reco=kRecoPi0P_1_emid(sr);if(truth< 0||reco< 0) return(-5.0f);return(float)(reco-truth);})
const Var kProng2NHitY_5label([](const caf::SRProxy *sr){int png2=kProng2ID_5label(sr);if(png2< 0) return-5;return(int) sr->vtx.elastic.fuzzyk.png[png2].nhity;})
const Var kProng2DirZ_Diff_Rad_emid([](const caf::SRProxy *sr){int png2=kProng2ID_emid(sr);if(png2< 0) return-5.0f;auto &pngDir=sr->vtx.elastic.fuzzyk.png[png2].dir;TVector3 truthDir=sr->vtx.elastic.fuzzyk.png[png2].truth.p.Vect();return(float) pngDir.Unit().Angle(truthDir.Unit());})
const NuTruthCut kNuSig_NT
const Var kProng2TruthPur_emid([](const caf::SRProxy *sr){int png2=kProng2ID_emid(sr);if(png2< 0) return-5.0f;return(float) sr->vtx.elastic.fuzzyk.png[png2].truth.pur;})
const Var kRecoPi0DirZ_Cos_emid([](const caf::SRProxy *sr){int png1=kProng1ID_emid(sr);int png2=kProng2ID_emid(sr);if(png1< 0||png2< 0) return-5.0f;const auto &p1=sr->vtx.elastic.fuzzyk.png[png1];const auto &p2=sr->vtx.elastic.fuzzyk.png[png2];TVector3 p=p1.calE *p1.dir.Unit()+p2.calE *p2.dir.Unit();TVector3 beamDirND=NuMIBeamDirection(caf::kNEARDET);return(float) p.Unit().Dot(beamDirND.Unit());})
const NuTruthVar kTrueMuKE_NT([](const caf::SRNeutrinoProxy *nu){float ke=-5;if(abs(nu->pdg)!=14||!nu->iscc) return ke;int nprims=nu->prim.size();for(int iprim=0;iprim< nprims;iprim++){if(abs(nu->prim[iprim].pdg)==13){double E=nu->prim[iprim].p.T();ke=E-MuonMass();}}return ke;})
const Var kBestTrack([](const caf::SRProxy *sr){return(int) sr->trk.kalman.idxmuonid;})
caf::Proxy< size_t > ntracks
Definition: SRProxy.h:1778
const Var kProng1PhotonScore_emid([](const caf::SRProxy *sr){int png1=kProng1ID_emid(sr);if(png1< 0) return-5.0f;const auto &p1=sr->vtx.elastic.fuzzyk.png[png1];return(float)(p1.spprongcvnpart5label.photonid);})
Proxy for caf::SRNeutrino.
Definition: SRProxy.h:510
const Var kProng1TruthE_emid([](const caf::SRProxy *sr){int png1=kProng1ID_emid(sr);if(png1< 0) return-5.0f;return(float) sr->vtx.elastic.fuzzyk.png[png1].truth.p.E;})
const NuTruthCut kSecoPi0_NT
const Var kProng2TruePDG_4view([](const caf::SRProxy *sr){int png2=kProng2ID_4view(sr);if(png2< 0) return 0;return(int) sr->vtx.elastic.fuzzyk.png[png2].truth.pdg;})
const Var kProng1NHitY_emid([](const caf::SRProxy *sr){int png1=kProng1ID_emid(sr);if(png1< 0) return-5;return(int) sr->vtx.elastic.fuzzyk.png[png1].nhity;})
const Var kProng2Score_emid([](const caf::SRProxy *sr){int png2=kProng2ID_emid(sr);if(png2< 0) return-5.0f;const auto &p2=sr->vtx.elastic.fuzzyk.png[png2];return(float)(p2.spprongcvnpartnumuccemid.emid);})
const Var kPi0Mass_5label([](const caf::SRProxy *sr){int png1=kProng1ID_5label(sr);int png2=kProng2ID_5label(sr);if(png1< 0||png2< 0) return-5.0f;const auto &p1=sr->vtx.elastic.fuzzyk.png[png1];const auto &p2=sr->vtx.elastic.fuzzyk.png[png2];const float dot=p1.dir.Unit().Dot(p2.dir.Unit());return 1000 *sqrt(2 *p1.calE *p2.calE *(1-dot));})
const char * p
Definition: xmltok.h:285
const Var kProng2TruthDirZ_Cos_emid([](const caf::SRProxy *sr){int png2=kProng2ID_emid(sr);if(png2< 0) return-5.0f;TVector3 pngDir=sr->vtx.elastic.fuzzyk.png[png2].truth.p.Vect();TVector3 beamDirND=NuMIBeamDirection(caf::kNEARDET);return(float) pngDir.Unit().Dot(beamDirND.Unit());})
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
T sqrt(T number)
Definition: d0nt_math.hpp:156
const NuTruthVar kTruePrimPi0P_NT([](const caf::SRNeutrinoProxy *truth){int TruePi0_ID=kTruePrimPi0_ID_NT(truth);if(TruePi0_ID< 0) return-5.0f;return(float) truth->prim[TruePi0_ID].p.Mag();})
constexpr T pow(T x)
Definition: pow.h:75
const Var kProng1DirZ_Diff_Cos_emid([](const caf::SRProxy *sr){int png1=kProng1ID_emid(sr);if(png1< 0) return-5.0f;auto &pngDir=sr->vtx.elastic.fuzzyk.png[png1].dir;TVector3 truthDir=sr->vtx.elastic.fuzzyk.png[png1].truth.p.Vect();return(float) pngDir.Unit().Dot(truthDir.Unit());})
const NuTruthVar kTruePrimPi0DirZ_Cos_NT([](const caf::SRNeutrinoProxy *truth){int TruePi0_ID=kTruePrimPi0_ID_NT(truth);if(TruePi0_ID< 0) return-5.0f;TVector3 Pi0P=truth->prim[TruePi0_ID].p.Vect();TVector3 beamDirND=NuMIBeamDirection(caf::kNEARDET);return(float) Pi0P.Unit().Dot(beamDirND.Unit());})
const Var kMuP_Diff([](const caf::SRProxy *sr){float truth=kTrueMuP(sr);float reco=kRecoMuonP(sr);if(truth< 0||reco< 0) return(-5.0f);return(float)(reco-truth);})
double GetCVNProngMuonScore_4view(int ProngIdx, const caf::SRProxy *sr)
Function to return "adjusted" CVN-prong muon score:
const Var kProng2NHitX_4view([](const caf::SRProxy *sr){int png2=kProng2ID_4view(sr);if(png2< 0) return-5;return(int) sr->vtx.elastic.fuzzyk.png[png2].nhitx;})
void abs(TH1 *hist)
const Var kProng2Score_5label([](const caf::SRProxy *sr){int png2=kProng2ID_5label(sr);if(png2< 0) return-5.0f;const auto &p2=sr->vtx.elastic.fuzzyk.png[png2];return(float)(p2.spprongcvnpart5label.emid);})
const Var kProng1ID_emid([](const caf::SRProxy *sr){return Prong1ID_Generator(sr,"emid");})
const Var kProng2NHit_5label([](const caf::SRProxy *sr){int png2=kProng2ID_5label(sr);if(png2< 0) return-5;return(int) sr->vtx.elastic.fuzzyk.png[png2].nhit;})
const Var kProng1DirZ_Cos_emid([](const caf::SRProxy *sr){int png1=kProng1ID_emid(sr);if(png1< 0) return-5.0f;auto &pngDir=sr->vtx.elastic.fuzzyk.png[png1].dir;TVector3 beamDirND=NuMIBeamDirection(caf::kNEARDET);return(float) pngDir.Unit().Dot(beamDirND.Unit());})
const NuTruthVar kTrueNuQ2_NT([](const caf::SRNeutrinoProxy *nu){return nu->q2;})
caf::Proxy< float > E
Definition: SRProxy.h:520
const Var kMuonPng_MuonScore_5label([](const caf::SRProxy *sr){return(float) GetCVNProngMuonScore_5label(kMuonPngIdx(sr), sr);})
const Var kProng2NHit_emid([](const caf::SRProxy *sr){int png2=kProng2ID_emid(sr);if(png2< 0) return-5;return(int) sr->vtx.elastic.fuzzyk.png[png2].nhit;})
const Var kRecoW2([](const caf::SRProxy *sr){float nuE=kRecoE(sr);float muE=kIncXsecMuonE(sr);float q2=kRecoQ2(sr);if(nuE< 0||muE< 0||q2< 0) return-5.0f;float mN=(sr->spill.isFHC)?kNeutronMass(sr):kProtonMass(sr);if((pow(mN, 2)+(2 *mN *(nuE-muE))-q2)< 0) return-5.0f;return(float)(pow(mN, 2)+(2 *mN *(nuE-muE))-q2);})
const Var kTwoProngTruthOpenAngle_Deg_emid([](const caf::SRProxy *sr){float TwoProngTruthOpenAngle_Rad=kTwoProngTruthOpenAngle_Rad_emid(sr);if(TwoProngTruthOpenAngle_Rad > TMath::Pi()||TwoProngTruthOpenAngle_Rad< 0.) return-5.0f;return(float)(TwoProngTruthOpenAngle_Rad *180./TMath::Pi());})
const Var kRecoPi0DirZ_Deg_emid([](const caf::SRProxy *sr){float RecoPi0DirZ_Rad=kRecoPi0DirZ_Rad_emid(sr);if(RecoPi0DirZ_Rad > TMath::Pi()||RecoPi0DirZ_Rad< 0.) return-5.0f;return(float)(RecoPi0DirZ_Rad *180./TMath::Pi());})
const Var kRecoMuonDirZ_Rad([](const caf::SRProxy *sr){int png_mu=kMuonPngIdx(sr);if(png_mu< 0) return-5.0f;const auto &png=sr->vtx.elastic.fuzzyk.png[png_mu];TVector3 p=png.calE *png.dir.Unit();TVector3 beamDirND=NuMIBeamDirection(caf::kNEARDET);return(float) p.Unit().Angle(beamDirND.Unit());})
const Var kProng1ElectronScore_emid([](const caf::SRProxy *sr){int png1=kProng1ID_emid(sr);if(png1< 0) return-5.0f;const auto &p1=sr->vtx.elastic.fuzzyk.png[png1];return(float)(p1.spprongcvnpart5label.electronid);})
int Prong1ID_Generator(const caf::SRProxy *sr, std::string id_type)
const NuTruthVar kTrueDeltaAlpha_NT([](const caf::SRNeutrinoProxy *truth){TVector3 p_mu=kMuon_Momentum(truth);TVector3 p_h=kH_Momentum(truth);TVector3 p_nu=kNu_Momentum(truth);TVector3 p_mu_T=p_nu.Cross(p_mu);TVector3 p_h_T=p_nu.Cross(p_h);TVector3 delta_p_T(0, 0, 0);delta_p_T+=p_mu_T;delta_p_T+=p_h_T;double numerator=delta_p_T.Dot(-p_mu_T);double mag_p_mu_T=p_mu_T.Mag();double mag_delta_p_T=delta_p_T.Mag();if(mag_p_mu_T==0||mag_delta_p_T==0) return-5.0;double val=numerator/(mag_p_mu_T *mag_delta_p_T);if(std::isnan(val)) return-5.0;if(std::isinf(val)) return-5.0;double delta_alpha=TMath::ACos(val);if(std::isnan(delta_alpha)) return-5.0;if(std::isinf(delta_alpha)) return-5.0;if(delta_alpha > TMath::Pi()||delta_alpha< 0.) return-5.0;return(double)(delta_alpha *180./TMath::Pi());})
const Var kMuonPngIdx([](const caf::SRProxy *sr){int longest_idx=-5;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;int trk_idx=kBestTrack(sr);if(sr->trk.kalman.ntracks< 1||trk_idx< 0||trk_idx >=int(sr->trk.kalman.ntracks)) return-5;auto &trk=sr->trk.kalman.tracks[trk_idx];int best_idx=-5;float best_cos=-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(trk.dir.Unit().Dot(png.dir.Unit()) > best_cos){best_cos=trk.dir.Unit().Dot(png.dir.Unit());best_idx=png_idx;}}}return best_idx;})
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
caf::Proxy< caf::SRElastic > elastic
Definition: SRProxy.h:2118
Double_t q2[12][num]
Definition: f2_nu.C:137
const Var kProng2TruthEff_emid([](const caf::SRProxy *sr){int png2=kProng2ID_emid(sr);if(png2< 0) return-5.0f;return(float) sr->vtx.elastic.fuzzyk.png[png2].truth.eff;})
float MuonEActandCat(double trklenactandcat)
const Var kProng2NHitY_emid([](const caf::SRProxy *sr){int png2=kProng2ID_emid(sr);if(png2< 0) return-5;return(int) sr->vtx.elastic.fuzzyk.png[png2].nhity;})
const Var kProng2NHitY_4view([](const caf::SRProxy *sr){int png2=kProng2ID_4view(sr);if(png2< 0) return-5;return(int) sr->vtx.elastic.fuzzyk.png[png2].nhity;})
const Var kProng2Score_4view([](const caf::SRProxy *sr){int png2=kProng2ID_4view(sr);if(png2< 0) return-5.0f;const auto &p2=sr->vtx.elastic.fuzzyk.png[png2];return(float)(p2.cvnpart.emid);})
const Var kProng1TruePDG_emid([](const caf::SRProxy *sr){int png1=kProng1ID_emid(sr);if(png1< 0) return 0;return(int) sr->vtx.elastic.fuzzyk.png[png1].truth.pdg;})
const NuTruthVar kTruePrimPi0DirZ_Deg_NT([](const caf::SRNeutrinoProxy *truth){float TruePi0DirZ_Rad_NT=kTruePrimPi0DirZ_Rad_NT(truth);if(TruePi0DirZ_Rad_NT > TMath::Pi()||TruePi0DirZ_Rad_NT< 0.) return-5.0f;return(float)(TruePi0DirZ_Rad_NT *180./TMath::Pi());})
caf::Proxy< bool > isFHC
Definition: SRProxy.h:1375
Track finder for cosmic rays.
const Var kRecoMuonDirZ_Cos([](const caf::SRProxy *sr){int png_mu=kMuonPngIdx(sr);if(png_mu< 0) return-5.0f;const auto &png=sr->vtx.elastic.fuzzyk.png[png_mu];TVector3 p=png.calE *png.dir.Unit();TVector3 beamDirND=NuMIBeamDirection(caf::kNEARDET);return(float) p.Unit().Dot(beamDirND.Unit());})
caf::Proxy< std::vector< caf::SRFuzzyKProng > > png
Definition: SRProxy.h:2043
const Var kCVNMuonIdx_5label([](const caf::SRProxy *sr){int longest_idx=-5;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;int best_idx=-5;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.spprongcvnpart5label.muonid > best_score){best_score=png.spprongcvnpart5label.muonid;best_idx=png_idx;}}}return best_idx;})
const Var kIncXsecHadE([](const caf::SRProxy *sr){int ibesttrk=kBestTrack(sr);if(sr->trk.kalman.ntracks< 1) return-1000.f;if(ibesttrk< 0||ibesttrk >=int(sr->trk.kalman.ntracks)) return-1000.f;float ExtraHadE=sr->trk.kalman.tracks[ibesttrk].overlapE;float CalhadE=sr->slc.calE-sr->trk.kalman.tracks[ibesttrk].calE;float vishadE=CalhadE+ExtraHadE;float hadE=0.0;hadE=VisibleHadE(vishadE);return hadE;})
const Var kProng2NHitX_5label([](const caf::SRProxy *sr){int png2=kProng2ID_5label(sr);if(png2< 0) return-5;return(int) sr->vtx.elastic.fuzzyk.png[png2].nhitx;})
const Var kProng1DirZ_Diff_Deg_emid([](const caf::SRProxy *sr){float Prong1DirZ_Diff_Rad=kProng1DirZ_Diff_Rad_emid(sr);if(Prong1DirZ_Diff_Rad > TMath::Pi()||Prong1DirZ_Diff_Rad< 0.) return-5.0f;return(float)(Prong1DirZ_Diff_Rad *180./TMath::Pi());})
const Var kIncXsecMuonE([](const caf::SRProxy *sr){float muonE=0.0;float muonEact=0.0;float muonEcat=0.0;float muonEactandcat=0.0;float trkLenAct=0.f;float trkLenCat=0.f;trkLenAct=kTrkLenAct(sr);trkLenCat=kTrkLenCat(sr);int ibesttrk=kBestTrack(sr);if(sr->trk.kalman.ntracks< 1) return-1000.f;if(ibesttrk< 0||ibesttrk >=int(sr->trk.kalman.ntracks)) return-1000.f;if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 && sr->trk.kalman.tracks[ibesttrk].lenincat< 0) muonEact=MuonEAct(trkLenAct);else if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 && sr->trk.kalman.tracks[ibesttrk].lenincat > 0){muonEcat=MuonECat(trkLenCat);muonEactandcat=MuonEActandCat(trkLenAct);muonE=muonEactandcat+muonEcat;}return muonE+muonEact;})
const TVector3 kNu_Momentum(const caf::SRNeutrinoProxy *truth)
const NuTruthVar kTrueNuW_NT([](const caf::SRNeutrinoProxy *nu){float val=nu->W2;if(val< 0) return(double)-5.0;return(double) sqrt(nu->W2);})
const Var kNeutronMass([](const caf::SRProxy *sr){return NeutronMass();})
const Var kTruePrimPi0DirZ_Deg
const NuTruthVar kTruthTag_NT([](const caf::SRNeutrinoProxy *truth){bool AntiNuSig=kAntiNuSig_NT(truth);bool NuSig=kNuSig_NT(truth);bool SecoPi0=kNumuCC_NT(truth)&&kSecoPi0_NT(truth);bool NoPi0=kNumuCC_NT(truth)&&kNoPi0_NT(truth);bool NC=kNC_NT(truth);bool Other=kIsOther_NT(truth);if(AntiNuSig) return 1;else if(NuSig) return 2;else if(SecoPi0) return 3;else if(NoPi0) return 4;else if(NC) return 5;else if(Other) return 6;else return 0;})
float MuonECat(double trklencat)
const Var kRecoW([](const caf::SRProxy *sr){float nuE=kRecoE(sr);float muE=kIncXsecMuonE(sr);float q2=kRecoQ2(sr);if(nuE< 0||muE< 0||q2< 0) return-5.0f;float mN=(sr->spill.isFHC)?kNeutronMass(sr):kProtonMass(sr);float W2=(pow(mN, 2)+(2 *mN *(nuE-muE))-q2);if(W2< 0) return-5.0f;return(float)(sqrt(W2));})
const TVector3 kH_Momentum(const caf::SRNeutrinoProxy *truth)
float VisibleHadE(float vishadE)
const Var kMuDeg_Diff([](const caf::SRProxy *sr){float truth=kTrueMuDirZ_Deg(sr);float reco=kRecoMuonDirZ_Deg(sr);if(truth< 0||reco< 0) return(-5.0f);return(float)(reco-truth);})
const Var kPi0Mass_4view([](const caf::SRProxy *sr){int png1=kProng1ID_4view(sr);int png2=kProng2ID_4view(sr);if(png1< 0||png2< 0) return-5.0f;const auto &p1=sr->vtx.elastic.fuzzyk.png[png1];const auto &p2=sr->vtx.elastic.fuzzyk.png[png2];const float dot=p1.dir.Unit().Dot(p2.dir.Unit());return 1000 *sqrt(2 *p1.calE *p2.calE *(1-dot));})
const Var kProng1NHitY_4view([](const caf::SRProxy *sr){int png1=kProng1ID_4view(sr);if(png1< 0) return-5;return(int) sr->vtx.elastic.fuzzyk.png[png1].nhity;})
caf::Proxy< caf::SRTrackBranch > trk
Definition: SRProxy.h:2145
const Var k5labelMuonPng_MuonScore_5label([](const caf::SRProxy *sr){return(float) GetCVNProngMuonScore_5label(kCVNMuonIdx_5label(sr), sr);})
const Var kRecoPi0P_4view([](const caf::SRProxy *sr){int png1=kProng1ID_4view(sr);int png2=kProng2ID_4view(sr);if(png1< 0||png2< 0) return-5.0f;const auto &p1=sr->vtx.elastic.fuzzyk.png[png1];const auto &p2=sr->vtx.elastic.fuzzyk.png[png2];TVector3 p=p1.calE *p1.dir.Unit()+p2.calE *p2.dir.Unit();return(float) p.Mag();})
Float_t E
Definition: plot.C:20
const Var kRecoMuonDirZ_Deg([](const caf::SRProxy *sr){float RecoMuonDirZ_Rad=kRecoMuonDirZ_Rad(sr);if(RecoMuonDirZ_Rad > TMath::Pi()||RecoMuonDirZ_Rad< 0.) return-5.0f;return(float)(RecoMuonDirZ_Rad *180./TMath::Pi());})
int Prong2ID_Generator(const caf::SRProxy *sr, std::string id_type)
const Var kProng2DirZ_Diff_Cos_emid([](const caf::SRProxy *sr){int png2=kProng2ID_emid(sr);if(png2< 0) return-5.0f;auto &pngDir=sr->vtx.elastic.fuzzyk.png[png2].dir;TVector3 truthDir=sr->vtx.elastic.fuzzyk.png[png2].truth.p.Vect();return(float) pngDir.Unit().Dot(truthDir.Unit());})
const Var kProng2PhotonScore_emid([](const caf::SRProxy *sr){int png2=kProng2ID_emid(sr);if(png2< 0) return-5.0f;const auto &p2=sr->vtx.elastic.fuzzyk.png[png2];return(float)(p2.spprongcvnpart5label.photonid);})
const Var kProng1TruthDirZ_Rad_emid([](const caf::SRProxy *sr){int png1=kProng1ID_emid(sr);if(png1< 0) return-5.0f;TVector3 pngDir=sr->vtx.elastic.fuzzyk.png[png1].truth.p.Vect();TVector3 beamDirND=NuMIBeamDirection(caf::kNEARDET);return(float) pngDir.Unit().Angle(beamDirND.Unit());})
const NuTruthCut kAntiNuSig_NT
const Var kW_Diff([](const caf::SRProxy *sr){float truth=kTrueNuW(sr);float reco=kRecoW(sr);if(truth< 0||reco< 0) return(-5.0f);return(float)(reco-truth);})
const Var kTwoProngOpenAngle_Reso_Rad_emid([](const caf::SRProxy *sr){float openAngle=kTwoProngOpenAngle_Rad_emid(sr);float truthOpenAngle=kTwoProngTruthOpenAngle_Rad_emid(sr);if(openAngle< 0.||openAngle > TMath::Pi()||truthOpenAngle< 0.||truthOpenAngle > TMath::Pi()) return-5.0f;return(float)(openAngle-truthOpenAngle);})
caf::Proxy< bool > iscc
Definition: SRProxy.h:538
const NuTruthVar kTrueMuDirZ_Cos_NT([](const caf::SRNeutrinoProxy *truth){if(abs(truth->pdg)!=14||!truth->iscc) return-5.0;int nprims=truth->prim.size();for(int iprim=0;iprim< nprims;iprim++){if(abs(truth->prim[iprim].pdg)==13){TVector3 mudir=truth->prim[iprim].p.Vect();TVector3 beamdir=NuMIBeamDirection(caf::kNEARDET);return mudir.Unit().Dot(beamdir.Unit());}}return-5.0;})
_Var< caf::SRNeutrinoProxy > NuTruthVar
Var designed to be used over the nuTree, ie all neutrinos, not just those that got slices...
Definition: Var.h:82
const Var kProng1CalE_emid([](const caf::SRProxy *sr){int png1=kProng1ID_emid(sr);if(png1< 0) return-5.0f;return(float) sr->vtx.elastic.fuzzyk.png[png1].calE;})
const Var kProng1ID_5label([](const caf::SRProxy *sr){return Prong1ID_Generator(sr,"5label");})
const NuTruthCut kNoPi0_NT
const NuTruthVar kTrueDeltaPtt_NT([](const caf::SRNeutrinoProxy *truth){TVector3 p_mu=kMuon_Momentum(truth);TVector3 p_h=kH_Momentum(truth);TVector3 p_nu=kNu_Momentum(truth);TVector3 p_mu_T=p_nu.Cross(p_mu);TVector3 p_h_T=p_nu.Cross(p_h);TVector3 delta_p_T(p_mu_T.X(), p_mu_T.Y(), p_mu_T.Z());delta_p_T+=p_h_T;TVector3 z_TT=p_mu_T.Unit();z_TT.SetX(0.);z_TT.SetY(0.);double delta_p_TT=delta_p_T.Dot(z_TT);if(std::isnan(delta_p_TT)) return-5.0;if(std::isinf(delta_p_TT)) return-5.0;return delta_p_TT;})
const TVector3 kMuon_Momentum(const caf::SRNeutrinoProxy *truth)
const Var kProng1Score_4view([](const caf::SRProxy *sr){int png1=kProng1ID_4view(sr);if(png1< 0) return-5.0f;const auto &p1=sr->vtx.elastic.fuzzyk.png[png1];return(float)(p1.cvnpart.emid);})
const Var k4viewMuonPng_MuonScore_5label([](const caf::SRProxy *sr){return(float) GetCVNProngMuonScore_5label(kCVNMuonIdx_4view(sr), sr);})
caf::StandardRecord * sr
const Var kProng2DirZ_Diff_Deg_emid([](const caf::SRProxy *sr){float Prong2DirZ_Diff_Rad=kProng2DirZ_Diff_Rad_emid(sr);if(Prong2DirZ_Diff_Rad > TMath::Pi()||Prong2DirZ_Diff_Rad< 0.) return-5.0f;return(float)(Prong2DirZ_Diff_Rad *180./TMath::Pi());})
const double j
Definition: BetheBloch.cxx:29
const Var kProng1TruePDG_4view([](const caf::SRProxy *sr){int png1=kProng1ID_4view(sr);if(png1< 0) return 0;return(int) sr->vtx.elastic.fuzzyk.png[png1].truth.pdg;})
const Var kMuonPng_MuonScore_4view([](const caf::SRProxy *sr){return(float) GetCVNProngMuonScore_4view(kMuonPngIdx(sr), sr);})
const NuTruthVar kTrueSecoPi0_ID1_NT([](const caf::SRNeutrinoProxy *truth){int nprims=truth->prim.size();int pi0ID1=-1;float pi0E=-5.0;for(int i=0;i< nprims;++i){auto &daughterPDG=truth->prim[i].daughterlist;auto &daughterE=truth->prim[i].daughterEnergies;if(daughterPDG.empty()) continue;for(int j=0;j< (int) daughterPDG.size();++j){if(daughterPDG[j]==111 &&pi0E< daughterE[j]){pi0ID1=i;pi0E=daughterE[j];}}}if(pi0ID1< 0) return-5;return pi0ID1;})
const NuTruthVar kTruePrimPi0E_NT([](const caf::SRNeutrinoProxy *truth){int TruePi0_ID=kTruePrimPi0_ID_NT(truth);if(TruePi0_ID< 0) return-5.0f;return(float) truth->prim[TruePi0_ID].p.E;})
const NuTruthVar kTrueMuDirZ_Deg_NT([](const caf::SRNeutrinoProxy *truth){float TrueMuDirZ_Rad_NT=kTrueMuDirZ_Rad_NT(truth);if(TrueMuDirZ_Rad_NT > TMath::Pi()||TrueMuDirZ_Rad_NT< 0.) return-5.0f;return(float)(TrueMuDirZ_Rad_NT *180./TMath::Pi());})
const Var kMuonID([](const caf::SRProxy *sr){int ibesttrk=kBestTrack(sr);if(sr->trk.kalman.ntracks< 1) return-5.0f;if(ibesttrk< 0||ibesttrk >=int(sr->trk.kalman.ntracks)) return-5.0f;return(float) sr->trk.kalman.tracks[ibesttrk].muonid;})
const Var kTwoProngTruthOpenAngle_Cos_emid([](const caf::SRProxy *sr){int png1=kProng1ID_emid(sr);int png2=kProng2ID_emid(sr);if(png1< 0||png2< 0) return-5.0f;TVector3 png1Dir=sr->vtx.elastic.fuzzyk.png[png1].truth.p.Vect();TVector3 png2Dir=sr->vtx.elastic.fuzzyk.png[png2].truth.p.Vect();return(float) png1Dir.Unit().Dot(png2Dir.Unit());})
const Var kProng2E_Reso_emid([](const caf::SRProxy *sr){int png2=kProng2ID_emid(sr);if(png2< 0) return-5.0f;float calE=sr->vtx.elastic.fuzzyk.png[png2].calE;float trueE=sr->vtx.elastic.fuzzyk.png[png2].truth.p.E;if(trueE==0) return-5.0f;return(float)(calE-trueE)/trueE;})
caf::Proxy< unsigned int > idxmuonid
Definition: SRProxy.h:1776
const Var kCVNMuonIdx_4view([](const caf::SRProxy *sr){int longest_idx=-5;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;int best_idx=-5;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;})
const TVector3 kPi0_Momentum(const caf::SRNeutrinoProxy *truth)
const Var kProng1NHit_4view([](const caf::SRProxy *sr){int png1=kProng1ID_4view(sr);if(png1< 0) return-5;return(int) sr->vtx.elastic.fuzzyk.png[png1].nhit;})
const TVector3 beamDirND
const Var kRecoMuonP([](const caf::SRProxy *sr){int png_mu=kMuonPngIdx(sr);if(png_mu< 0) return-5.0f;const auto &png=sr->vtx.elastic.fuzzyk.png[png_mu];TVector3 p=png.calE *png.dir.Unit();return(float) p.Mag();})
const Var kTwoProngOpenAngle_Cos_emid([](const caf::SRProxy *sr){int png1=kProng1ID_emid(sr);int png2=kProng2ID_emid(sr);if(png1< 0||png2< 0) return-5.0f;auto &png1Dir=sr->vtx.elastic.fuzzyk.png[png1].dir;auto &png2Dir=sr->vtx.elastic.fuzzyk.png[png2].dir;return(float) png1Dir.Unit().Dot(png2Dir.Unit());})
const Var kProng1TruthPur_emid([](const caf::SRProxy *sr){int png1=kProng1ID_emid(sr);if(png1< 0) return-5.0f;return(float) sr->vtx.elastic.fuzzyk.png[png1].truth.pur;})
const Var kTrkLenAct([](const caf::SRProxy *sr){int ibesttrk=kBestTrack(sr);if(sr->trk.kalman.ntracks< 1) return-1000.f;if(ibesttrk< 0||ibesttrk >=int(sr->trk.kalman.ntracks)) return-1000.f; if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 && sr->trk.kalman.tracks[ibesttrk].lenincat< 0) return float((sr->trk.kalman.tracks[ibesttrk].leninact/100.) +(sr->trk.kalman.tracks[ibesttrk].lenincat/100.)); if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 && sr->trk.kalman.tracks[ibesttrk].lenincat > 0) return float(sr->trk.kalman.tracks[ibesttrk].leninact/100.); if(sr->trk.kalman.tracks[ibesttrk].leninact< 0 && sr->trk.kalman.tracks[ibesttrk].lenincat > 0) return 0.f;return-1000.f;})
const Var kQ2_Diff([](const caf::SRProxy *sr){float truth=kTrueNuQ2(sr);float reco=kRecoQ2(sr);if(truth< 0||reco< 0) return(-5.0f);return(float)(reco-truth);})
const NuTruthVar kTrueHadronicSystemP_NT([](const caf::SRNeutrinoProxy *truth){TVector3 p_h=kH_Momentum(truth);return p_h.Mag();})
double dot(const std::vector< double > &x, const std::vector< double > &y)
Definition: dot.hpp:10
const Var kProng1E_Reso_emid([](const caf::SRProxy *sr){int png1=kProng1ID_emid(sr);if(png1< 0) return-5.0f;float calE=sr->vtx.elastic.fuzzyk.png[png1].calE;float trueE=sr->vtx.elastic.fuzzyk.png[png1].truth.p.E;if(trueE==0) return-5.0f;return(float)(calE-trueE)/trueE;})
const NuTruthCut kIsOther_NT
const Var kTwoProngOpenAngle_Rad_emid([](const caf::SRProxy *sr){int png1=kProng1ID_emid(sr);int png2=kProng2ID_emid(sr);if(png1< 0||png2< 0) return-5.0f;auto &png1Dir=sr->vtx.elastic.fuzzyk.png[png1].dir;auto &png2Dir=sr->vtx.elastic.fuzzyk.png[png2].dir;return(float) png1Dir.Unit().Angle(png2Dir.Unit());})
caf::Proxy< caf::SRKalman > kalman
Definition: SRProxy.h:1797
const int NC
const Var kTwoProngOpenAngle_Deg_emid([](const caf::SRProxy *sr){float TwoProngOpenAngle_Rad=kTwoProngOpenAngle_Rad_emid(sr);if(TwoProngOpenAngle_Rad > TMath::Pi()||TwoProngOpenAngle_Rad< 0.) return-5.0f;return(float)(TwoProngOpenAngle_Rad *180./TMath::Pi());})
const Var kProng1ID_4view([](const caf::SRProxy *sr){return Prong1ID_Generator(sr,"4view");})
caf::Proxy< bool > IsValid
Definition: SRProxy.h:2058
caf::Proxy< short int > pdgorig
Definition: SRProxy.h:553
const Var kProng2NHitX_emid([](const caf::SRProxy *sr){int png2=kProng2ID_emid(sr);if(png2< 0) return-5;return(int) sr->vtx.elastic.fuzzyk.png[png2].nhitx;})
const Var kProng1ProtonScore_emid([](const caf::SRProxy *sr){int png1=kProng1ID_emid(sr);if(png1< 0) return-5.0f;const auto &p1=sr->vtx.elastic.fuzzyk.png[png1];return(float)(p1.spprongcvnpart5label.protonid);})
caf::Proxy< caf::SRSlice > slc
Definition: SRProxy.h:2142
const Var kProng1NHitX_5label([](const caf::SRProxy *sr){int png1=kProng1ID_5label(sr);if(png1< 0) return-5;return(int) sr->vtx.elastic.fuzzyk.png[png1].nhitx;})
const Var kNProngs([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-5;return(int) sr->vtx.elastic.fuzzyk.npng;})
const Var kProng1DirZ_Rad_emid([](const caf::SRProxy *sr){int png1=kProng1ID_emid(sr);if(png1< 0) return-5.0f;auto &pngDir=sr->vtx.elastic.fuzzyk.png[png1].dir;TVector3 beamDirND=NuMIBeamDirection(caf::kNEARDET);return(float) pngDir.Unit().Angle(beamDirND.Unit());})
caf::Proxy< float > calE
Definition: SRProxy.h:1292
const NuTruthCut kNumuCC_NT
const Var kTrkLenCat([](const caf::SRProxy *sr){int ibesttrk=kBestTrack(sr);if(sr->trk.kalman.ntracks< 1) return-1000.f;if(ibesttrk< 0||ibesttrk >=int(sr->trk.kalman.ntracks)) return-1000.f; if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 && sr->trk.kalman.tracks[ibesttrk].lenincat< 0) return 0.f; if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 && sr->trk.kalman.tracks[ibesttrk].lenincat > 0) return float(sr->trk.kalman.tracks[ibesttrk].lenincat/100.); if(sr->trk.kalman.tracks[ibesttrk].leninact< 0 && sr->trk.kalman.tracks[ibesttrk].lenincat > 0) return float((sr->trk.kalman.tracks[ibesttrk].leninact/100.) +(sr->trk.kalman.tracks[ibesttrk].lenincat/100.));return-1000.f;})
const Var kTwoProngTruthOpenAngle_Rad_emid([](const caf::SRProxy *sr){int png1=kProng1ID_emid(sr);int png2=kProng2ID_emid(sr);if(png1< 0||png2< 0) return-5.0f;TVector3 png1Dir=sr->vtx.elastic.fuzzyk.png[png1].truth.p.Vect();TVector3 png2Dir=sr->vtx.elastic.fuzzyk.png[png2].truth.p.Vect();return(float) png1Dir.Unit().Angle(png2Dir.Unit());})
const Var kProng2ID_5label([](const caf::SRProxy *sr){return Prong2ID_Generator(sr,"5label");})
const Var kProng1Score_5label([](const caf::SRProxy *sr){int png1=kProng1ID_5label(sr);if(png1< 0) return-5.0f;const auto &p1=sr->vtx.elastic.fuzzyk.png[png1];return(float)(p1.spprongcvnpart5label.emid);})
const NuTruthVar kTruePrimPi0DirZ_Rad_NT([](const caf::SRNeutrinoProxy *truth){int TruePi0_ID=kTruePrimPi0_ID_NT(truth);if(TruePi0_ID< 0) return-5.0f;TVector3 Pi0P=truth->prim[TruePi0_ID].p.Vect();TVector3 beamDirND=NuMIBeamDirection(caf::kNEARDET);return(float) Pi0P.Unit().Angle(beamDirND.Unit());})
const NuTruthVar kTruePrimPi0KE_NT([](const caf::SRNeutrinoProxy *truth){float nPrimaryPi0=kNPrimaryPi0_NT(truth);if(nPrimaryPi0< 1) return-5.f; float pi0_t=-5;const int nprims=truth->prim.size();for(int iprim=0;iprim< nprims;++iprim){if(truth->prim[iprim].pdg!=111) continue;double E=truth->prim[iprim].p.T();float ke=E-Pi0Mass();if(ke > pi0_t) pi0_t=ke;}return float(pi0_t);})
const Var kProng2NHit_4view([](const caf::SRProxy *sr){int png2=kProng2ID_4view(sr);if(png2< 0) return-5;return(int) sr->vtx.elastic.fuzzyk.png[png2].nhit;})
double GetCVNProngMuonScore_5label(int ProngIdx, const caf::SRProxy *sr)
const NuTruthVar kTrueSecoPi0P_NT([](const caf::SRNeutrinoProxy *truth){int TruePi0_ID1=kTrueSecoPi0_ID1_NT(truth);int TruePi0_ID2=kTrueSecoPi0_ID2_NT(truth);double pi0P=-5.0;if(TruePi0_ID1< 0||TruePi0_ID2< 0) return(float) pi0P;double pi0E=truth->prim[TruePi0_ID1].daughterEnergies[TruePi0_ID2];if(pi0E< Pi0Mass()) pi0P=0.;else pi0P=sqrt(std::pow(pi0E, 2)-std::pow(Pi0Mass(), 2));return(float) pi0P;})
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:329
const Var kProng1NHit_5label([](const caf::SRProxy *sr){int png1=kProng1ID_5label(sr);if(png1< 0) return-5;return(int) sr->vtx.elastic.fuzzyk.png[png1].nhit;})
caf::Proxy< caf::SRLorentzVector > p
Definition: SRProxy.h:551
const TVector3 kNucleon_Momentum(const caf::SRNeutrinoProxy *truth)
const NuTruthVar kNSecondaryPi0_NT([](const caf::SRNeutrinoProxy *truth){int nprims=truth->prim.size();int countpi0=0;for(int i=0;i< nprims;++i){auto &daughters=truth->prim[i].daughterlist;if(daughters.empty()) continue;for(const auto &pdg:daughters){if(pdg==111) countpi0++;}}return countpi0;})
const Var k5labelMuonPng_MuonScore_4view([](const caf::SRProxy *sr){return(float) GetCVNProngMuonScore_4view(kCVNMuonIdx_5label(sr), sr);})
const Var kProng2TruthDirZ_Rad_emid([](const caf::SRProxy *sr){int png2=kProng2ID_emid(sr);if(png2< 0) return-5.0f;TVector3 pngDir=sr->vtx.elastic.fuzzyk.png[png2].truth.p.Vect();TVector3 beamDirND=NuMIBeamDirection(caf::kNEARDET);return(float) pngDir.Unit().Angle(beamDirND.Unit());})
caf::Proxy< float > q2
Definition: SRProxy.h:557
caf::Proxy< short int > pdg
Definition: SRProxy.h:552
const Var kRecoPi0P_1_emid([](const caf::SRProxy *sr){int png1=kProng1ID_emid(sr);int png2=kProng2ID_emid(sr);if(png1< 0||png2< 0) return-5.0f;const auto &p1=sr->vtx.elastic.fuzzyk.png[png1];const auto &p2=sr->vtx.elastic.fuzzyk.png[png2];TVector3 p=p1.calE *p1.dir.Unit()+p2.calE *p2.dir.Unit();return(float) p.Mag();})
caf::Proxy< float > W2
Definition: SRProxy.h:522
caf::Proxy< std::vector< caf::SRKalmanTrack > > tracks
Definition: SRProxy.h:1780
const NuTruthVar kTruePrimPi0_ID_NT([](const caf::SRNeutrinoProxy *truth){int nprims=truth->prim.size();int pi0ID=-1;float pi0E=-5.0;for(int i=0;i< nprims;++i){if(truth->prim[i].pdg==111 &&pi0E< truth->prim[i].p.E){pi0ID=i;pi0E=truth->prim[i].p.E;}}if(pi0ID< 0) return-5;return pi0ID;})
const Var kProng1NHitY_5label([](const caf::SRProxy *sr){int png1=kProng1ID_5label(sr);if(png1< 0) return-5;return(int) sr->vtx.elastic.fuzzyk.png[png1].nhity;})
caf::Proxy< caf::SRVertexBranch > vtx
Definition: SRProxy.h:2146
const Var kProng1TruthEff_emid([](const caf::SRProxy *sr){int png1=kProng1ID_emid(sr);if(png1< 0) return-5.0f;return(float) sr->vtx.elastic.fuzzyk.png[png1].truth.eff;})
Template for Var and SpillVar.
const NuTruthCut kNC_NT
const TVector3 kDeltaPt(const caf::SRNeutrinoProxy *truth)
const Var kRecoPi0DirZ_Cos_4view([](const caf::SRProxy *sr){int png1=kProng1ID_4view(sr);int png2=kProng2ID_4view(sr);if(png1< 0||png2< 0) return-5.0f;const auto &p1=sr->vtx.elastic.fuzzyk.png[png1];const auto &p2=sr->vtx.elastic.fuzzyk.png[png2];TVector3 p=p1.calE *p1.dir.Unit()+p2.calE *p2.dir.Unit();TVector3 beamDirND=NuMIBeamDirection(caf::kNEARDET);return(float) p.Unit().Dot(beamDirND.Unit());})
const Var kRecoPi0DirZ_Rad_4view([](const caf::SRProxy *sr){int png1=kProng1ID_4view(sr);int png2=kProng2ID_4view(sr);if(png1< 0||png2< 0) return-5.0f;const auto &p1=sr->vtx.elastic.fuzzyk.png[png1];const auto &p2=sr->vtx.elastic.fuzzyk.png[png2];TVector3 p=p1.calE *p1.dir.Unit()+p2.calE *p2.dir.Unit();TVector3 beamDirND=NuMIBeamDirection(caf::kNEARDET);return(float) p.Unit().Angle(beamDirND.Unit());})
const NuTruthVar kTrueSecoPi0_ID2_NT([](const caf::SRNeutrinoProxy *truth){int nprims=truth->prim.size();int pi0ID2=-1;float pi0E=-5.0;for(int i=0;i< nprims;++i){auto &daughterPDG=truth->prim[i].daughterlist;auto &daughterE=truth->prim[i].daughterEnergies;if(daughterPDG.empty()) continue;for(int j=0;j< (int) daughterPDG.size();++j){if(daughterPDG[j]==111 &&pi0E< daughterE[j]){pi0ID2=j;pi0E=daughterE[j];}}}if(pi0ID2< 0) return-5;return pi0ID2;})
const Var kProng1DirZ_Diff_Rad_emid([](const caf::SRProxy *sr){int png1=kProng1ID_emid(sr);if(png1< 0) return-5.0f;auto &pngDir=sr->vtx.elastic.fuzzyk.png[png1].dir;TVector3 truthDir=sr->vtx.elastic.fuzzyk.png[png1].truth.p.Vect();return(float) pngDir.Unit().Angle(truthDir.Unit());})
const Var kProng1TruthDirZ_Cos_emid([](const caf::SRProxy *sr){int png1=kProng1ID_emid(sr);if(png1< 0) return-5.0f;TVector3 pngDir=sr->vtx.elastic.fuzzyk.png[png1].truth.p.Vect();TVector3 beamDirND=NuMIBeamDirection(caf::kNEARDET);return(float) pngDir.Unit().Dot(beamDirND.Unit());})
float MuonEAct(double trklenact)
const Var kProng2TruePDG_emid([](const caf::SRProxy *sr){int png2=kProng2ID_emid(sr);if(png2< 0) return 0;return(int) sr->vtx.elastic.fuzzyk.png[png2].truth.pdg;})
const Var kProtonMass([](const caf::SRProxy *sr){return ProtonMass();})
const Var kProng2CalE_emid([](const caf::SRProxy *sr){int png2=kProng2ID_emid(sr);if(png2< 0) return-5.0f;return(float) sr->vtx.elastic.fuzzyk.png[png2].calE;})
const Var kPi0Mass_emid([](const caf::SRProxy *sr){int png1=kProng1ID_emid(sr);int png2=kProng2ID_emid(sr);if(png1< 0||png2< 0) return-5.0f;const auto &p1=sr->vtx.elastic.fuzzyk.png[png1];const auto &p2=sr->vtx.elastic.fuzzyk.png[png2];const float dot=p1.dir.Unit().Dot(p2.dir.Unit());return 1000 *sqrt(2 *p1.calE *p2.calE *(1-dot));})
const Var kRecoPi0DirZ_Rad_emid([](const caf::SRProxy *sr){int png1=kProng1ID_emid(sr);int png2=kProng2ID_emid(sr);if(png1< 0||png2< 0) return-5.0f;const auto &p1=sr->vtx.elastic.fuzzyk.png[png1];const auto &p2=sr->vtx.elastic.fuzzyk.png[png2];TVector3 p=p1.calE *p1.dir.Unit()+p2.calE *p2.dir.Unit();TVector3 beamDirND=NuMIBeamDirection(caf::kNEARDET);return(float) p.Unit().Angle(beamDirND.Unit());})
const Var kProng2ProtonScore_emid([](const caf::SRProxy *sr){int png2=kProng2ID_emid(sr);if(png2< 0) return-5.0f;const auto &p2=sr->vtx.elastic.fuzzyk.png[png2];return(float)(p2.spprongcvnpart5label.protonid);})
const NuTruthVar kTrueNuE_NT([](const caf::SRNeutrinoProxy *nu){return nu->E;})
const Var kProng1NHit_emid([](const caf::SRProxy *sr){int png1=kProng1ID_emid(sr);if(png1< 0) return-5;return(int) sr->vtx.elastic.fuzzyk.png[png1].nhit;})
const Var kProng2TruthE_emid([](const caf::SRProxy *sr){int png2=kProng2ID_emid(sr);if(png2< 0) return-5.0f;return(float) sr->vtx.elastic.fuzzyk.png[png2].truth.p.E;})
const Var kProng2ID_emid([](const caf::SRProxy *sr){return Prong2ID_Generator(sr,"emid");})
caf::Proxy< std::vector< caf::SRTrueParticle > > prim
Definition: SRProxy.h:555
const Var kProng2DirZ_Rad_emid([](const caf::SRProxy *sr){int png2=kProng2ID_emid(sr);if(png2< 0) return-5.0f;auto &pngDir=sr->vtx.elastic.fuzzyk.png[png2].dir;TVector3 beamDirND=NuMIBeamDirection(caf::kNEARDET);return(float) pngDir.Unit().Angle(beamDirND.Unit());})
const Var kPi0Deg_Diff([](const caf::SRProxy *sr){float truth=kTruePrimPi0DirZ_Deg(sr);float reco=kRecoPi0DirZ_Deg_emid(sr);if(truth< 0||reco< 0) return(-5.0f);return(float)(reco-truth);})
const Var kProng2DirZ_Cos_emid([](const caf::SRProxy *sr){int png2=kProng2ID_emid(sr);if(png2< 0) return-5.0f;auto &pngDir=sr->vtx.elastic.fuzzyk.png[png2].dir;TVector3 beamDirND=NuMIBeamDirection(caf::kNEARDET);return(float) pngDir.Unit().Dot(beamDirND.Unit());})
const Var kRecoQ2([](const caf::SRProxy *sr){float nuE=kRecoE(sr);float muE=kIncXsecMuonE(sr);float pmu=kRecoMuonP(sr);float cosmu=kRecoMuonDirZ_Cos(sr);if(nuE< 0||muE< 0||pmu< 0||cosmu==-5.0) return-5.0f;return(float)(2 *nuE *(muE-pmu *cosmu));})
const Var kProng2ElectronScore_emid([](const caf::SRProxy *sr){int png2=kProng2ID_emid(sr);if(png2< 0) return-5.0f;const auto &p2=sr->vtx.elastic.fuzzyk.png[png2];return(float)(p2.spprongcvnpart5label.electronid);})
const Var kProng1NHitX_4view([](const caf::SRProxy *sr){int png1=kProng1ID_4view(sr);if(png1< 0) return-5;return(int) sr->vtx.elastic.fuzzyk.png[png1].nhitx;})
enum BeamMode string