NumuCCIncVars.cxx
Go to the documentation of this file.
2 
4 
5 #include "sys/stat.h"
6 #include "TCanvas.h"
7 #include "TFile.h"
8 
9 namespace ana{
10 
13 
14  float MuonLLR(const float dedxll,
15  const float scatll)
16  {
17 
18  if(!hprob_dedx_vs_scat_sig || !hprob_dedx_vs_scat_bg){
19  const char* path = getenv("SRT_PUBLIC_CONTEXT");
20  std::string fname = std::string(path) + "/NDAna/numucc_inc/muon_pdf.root";
21 
22  TFile *muon_prob_file = new TFile( fname.c_str(), "r");
23 
24  hprob_dedx_vs_scat_sig = (TH2F*)muon_prob_file->Get("prob_sig");
25  hprob_dedx_vs_scat_bg = (TH2F*)muon_prob_file->Get("prob_bg");
26 
27  }
28 
29  const int dedx_vs_scat_nxbins = hprob_dedx_vs_scat_bg->GetXaxis()->GetNbins();
30  const int dedx_vs_scat_nybins = hprob_dedx_vs_scat_bg->GetYaxis()->GetNbins();
31 
32  int xbin_sig = hprob_dedx_vs_scat_sig->GetXaxis()->FindBin(dedxll);
33  int ybin_sig = hprob_dedx_vs_scat_sig->GetYaxis()->FindBin(scatll);
34 
35  int xbin_bg = hprob_dedx_vs_scat_bg->GetXaxis()->FindBin(dedxll);
36  int ybin_bg = hprob_dedx_vs_scat_bg->GetYaxis()->FindBin(scatll);
37 
38  float llr = -1e12;
39  // if the bins are underflow or overflow, it's beyond the range of our model
40  // skip over to the next track.
41  if (xbin_sig == 0 || xbin_sig > dedx_vs_scat_nxbins ||
42  ybin_sig == 0 || ybin_sig > dedx_vs_scat_nybins ||
43  xbin_bg == 0 || xbin_bg > dedx_vs_scat_nxbins ||
44  ybin_bg == 0 || ybin_bg > dedx_vs_scat_nybins)
45  return llr;
46 
47  float psig = std::log(hprob_dedx_vs_scat_sig->GetBinContent(xbin_sig, ybin_sig));
48  float pbg = std::log(hprob_dedx_vs_scat_bg->GetBinContent(xbin_bg, ybin_bg));
49 
50  llr = psig - pbg;
51  return llr;
52  }
53 
54 
55  //---------------------------------
56  // Muon Energy for Active region
57  //---------------------------------
58 
59  float MuonEAct(double trklenact){
60 
61  double p0 = 1.67012e-01;
62  double p1 = 1.79305e-01;
63  double p2 = 3.74708e-03;
64  double p3 = -1.54232e-04;
65 
66  float MuonE = 0.0;
67 
68  if (trklenact <= 0.0) return 0.0;
69 
70  MuonE = p0
71  + p1 * trklenact
72  + p2 * std::pow(trklenact, 2)
73  + p3 * std::pow(trklenact, 3);
74 
75  return MuonE;
76  }
77 
78 
79  //---------------------------------------
80  // Muon Energy for Active+Catcher region
81  //(for Active region Track length)
82  //---------------------------------------
83 
84  float MuonEActandCat(double trklenactandcat){
85 
86  double p0 = 1.21130e-02;
87  double p1 = 1.97903e-01;
88  double p2 = 7.82459e-04;
89 
90  float MuonE = 0.0;
91 
92  if (trklenactandcat <= 0.0) return 0.0;
93 
94  MuonE = p0
95  + p1 * trklenactandcat
96  + p2 * std::pow(trklenactandcat, 2);
97 
98  return MuonE;
99  }
100 
101 
102  //---------------------------------------
103  // Muon Energy for Active+Catcher region
104  //(for catcher region Track length)
105  //---------------------------------------
106 
107  float MuonECat(double trklencat){
108 
109  double offset = 1.31325e-01;
110  double slope = 5.35146e-01;
111 
112  float MuonE = 0.0;
113 
114  if (trklencat <= 0.0) return 0.0;
115 
116  MuonE = slope*trklencat + offset;
117 
118  return MuonE;
119  }
120 
121  //---------------------------------------
122  // Haronic Energy
123  //----------------------------------------
124 
125  float VisibleHadE(float vishadE) {
126 
127  double p0 = 5.85254e-02;
128  double p1 = 1.27796e+00;
129  double p2 = 3.75457e-01;
130  double p3 = -5.45618e-01;
131  double p4 = 1.65975e-01;
132 
133  float HadE = 0.0;
134 
135  HadE = p0
136  + p1 * vishadE
137  + p2 * std::pow(vishadE, 2)
138  + p3 * std::pow(vishadE, 3)
139  + p4 * std::pow(vishadE, 4);
140 
141  return HadE;
142  }
143 
144 
145  //-------------------------------------------------------------------------------------------------
146  // PID related Vars
147 
148  /// Return the muon to non-muon likelihod ration, given
149  /// the dedx and scattering likelihood ratios for the track
150  float MuonLLR(const float dedxll, const float scatll);
151 
152  const Var kMostMuonLikeTrackIdx([](const caf::SRProxy* sr)
153  {
154  float maxllr = -1e12;
155  int bestidx = -5;
156  // find the track with max LLR for being a muon
157  for(unsigned int itrk = 0; itrk < sr->trk.kalman.ntracks; itrk++){
158  float llr = MuonLLR(sr->trk.kalman.tracks[itrk].dedxllh,
159  sr->trk.kalman.tracks[itrk].scatllh);
160  if(llr > maxllr){
161  maxllr = llr;
162  bestidx = itrk;
163  }
164  } // end loop over tracks
165  return bestidx;
166  });
167 
168  const Var kMuonLLRPID([](const caf::SRProxy* sr)
169  {
170  int bestidx = kMostMuonLikeTrackIdx(sr);
171  float llr = -1e12;
172 
173  if(bestidx < 0 || bestidx >= int(sr->trk.kalman.ntracks))
174  return llr;
175 
176  // if we get here, the bestidx is a valid index into the tracks vector
177  llr = MuonLLR(sr->trk.kalman.tracks[bestidx].dedxllh,
178  sr->trk.kalman.tracks[bestidx].scatllh);
179  return llr;
180  });
181 
182  //----------------------------
183  // Energy Estimation
184  //----------------------------
185 
186  float MuonEActandCat(double trklenactandcat);
187  float MuonEAct(double trklenact);
188  float MuonECat(double trklencat);
189  float VisibleHadE(float vishadE);
190 
191  const Var kMuStartX([](const caf::SRProxy* sr)
192  {
193  int ibesttrk = kBestTrack(sr);
194 
195  if(sr->trk.kalman.ntracks < 1)
196  return -1000.f;
197 
198  if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
199  return -1000.f;
200 
201  return sr->trk.kalman.tracks[ibesttrk].start.X();
202  });
203 
204 
205  const Var kMuStartY([](const caf::SRProxy* sr)
206  {
207  int ibesttrk = kBestTrack(sr);
208 
209  if(sr->trk.kalman.ntracks < 1)
210  return -1000.f;
211 
212  if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
213  return -1000.f;
214 
215  return sr->trk.kalman.tracks[ibesttrk].start.Y();
216  });
217 
218 
219  const Var kMuStartZ([](const caf::SRProxy* sr)
220  {
221  int ibesttrk = kBestTrack(sr);
222 
223  if(sr->trk.kalman.ntracks < 1)
224  return -1000.f;
225 
226  if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
227  return -1000.f;
228 
229  return sr->trk.kalman.tracks[ibesttrk].start.Z();
230  });
231 
232  const Var kMuStopX([](const caf::SRProxy* sr)
233  {
234  int ibesttrk = kBestTrack(sr);
235 
236  if(sr->trk.kalman.ntracks < 1)
237  return -1000.f;
238 
239  if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
240  return -1000.f;
241 
242  return sr->trk.kalman.tracks[ibesttrk].stop.X();
243  });
244 
245 
246  const Var kMuStopY([](const caf::SRProxy* sr)
247  {
248  int ibesttrk = kBestTrack(sr);
249 
250  if(sr->trk.kalman.ntracks < 1)
251  return -1000.f;
252 
253  if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
254  return -1000.f;
255 
256  return sr->trk.kalman.tracks[ibesttrk].stop.Y();
257  });
258 
259 
260  const Var kMuStopZ([](const caf::SRProxy* sr)
261  {
262  int ibesttrk = kBestTrack(sr);
263 
264  if(sr->trk.kalman.ntracks < 1)
265  return -1000.f;
266 
267  if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
268  return -1000.f;
269 
270  return sr->trk.kalman.tracks[ibesttrk].stop.Z();
271  });
272 
273 
274  const Var kMuonLength([](const caf::SRProxy* sr)
275  {
276  int ibesttrk = kBestTrack(sr);
277 
278  if(sr->trk.kalman.ntracks < 1)
279  return -1000.f;
280 
281  if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
282  return -1000.f;
283 
284  return float(sr->trk.kalman.tracks[ibesttrk].len/100.);
285  });
286 
287  const Var kVisHadE([](const caf::SRProxy* sr)
288  {
289  int ibesttrk = kBestTrack(sr);
290 
291  if(sr->trk.kalman.ntracks < 1)
292  return -1000.f;
293 
294  if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
295  return -1000.f;
296 
297  // Extra energy (hadronic contamination) associated with muon track
298  float ExtraHadE = sr->trk.kalman.tracks[ibesttrk].overlapE;
299 
300  // calorimetric slice energy - Calorimetric Energy of Muon Track
301  float CalhadE = sr->slc.calE - sr->trk.kalman.tracks[ibesttrk].calE;
302 
303  // Add calorimetric hadE and Overlap energy in the muon track
304  float vishadE = CalhadE + ExtraHadE;
305 
306  return vishadE;
307  });
308 
309  const Var kHadEonTrack([](const caf::SRProxy* sr)
310  {
311  int ibesttrk = kBestTrack(sr);
312 
313  if(sr->trk.kalman.ntracks < 1)
314  return -1000.f;
315 
316  if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
317  return -1000.f;
318 
319  // Extra energy (hadronic contamination) associated with muon track
320  return (float)sr->trk.kalman.tracks[ibesttrk].overlapE;
321  });
322 
323 
324  const Var kAveragededx40([](const caf::SRProxy* sr)
325  {
326  int ibesttrk = kBestTrack(sr);
327 
328  if(sr->trk.kalman.ntracks < 1)
329  return -1000.f;
330 
331  if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
332  return -1000.f;
333 
334  if (sr->trk.kalman.tracks[ibesttrk].avedEdxlast40cm == -5)
335  return -1000.f;
336 
337  return float(sr->trk.kalman.tracks[ibesttrk].avedEdxlast40cm);
338  });
339 
340 //==============================================================================
341 // Redefine track lengths in active and catcher volumes with the new PID.
342 // Note that at this moment what is called leninact is actually distance from
343 // the start of the track to the transition plane. If the transition plane is
344 // outside of the track, that number is the track length plus the length of the
345 // projection from the last hit along the direction determined by the last two
346 // hits to the transition plane. The lenincat behaves similarly.
347 // In the next production the two variables should be changed to what they are
348 // supposed to mean literally.
349 
350  const Var kTrkLenAct([](const caf::SRProxy* sr)
351  {
352  int ibesttrk = kBestTrack(sr);
353  if(sr->trk.kalman.ntracks < 1)
354  return -1000.f;
355 
356  if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
357  return -1000.f;
358 
359  // If leninact is positive and lenincat is negative,
360  // the transition plane is to the right of the track.
361  if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 &&
362  sr->trk.kalman.tracks[ibesttrk].lenincat < 0)
363  return float((sr->trk.kalman.tracks[ibesttrk].leninact / 100.)
364  +(sr->trk.kalman.tracks[ibesttrk].lenincat / 100.));
365 
366  // If leninact is positive and lenincat is positive,
367  // the transition plane is in the middle of the track.
368  if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 &&
369  sr->trk.kalman.tracks[ibesttrk].lenincat > 0)
370  return float(sr->trk.kalman.tracks[ibesttrk].leninact / 100.);
371 
372  // If leninact is negative and lenincat is positive,
373  // the transition plane is to the left of the track.
374  if(sr->trk.kalman.tracks[ibesttrk].leninact < 0 &&
375  sr->trk.kalman.tracks[ibesttrk].lenincat > 0)
376  return 0.f;
377 
378  return -1000.f;
379  });
380 
381  const Var kTrkLenCat([](const caf::SRProxy* sr)
382  {
383  int ibesttrk = kBestTrack(sr);
384  if(sr->trk.kalman.ntracks < 1)
385  return -1000.f;
386 
387  if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
388  return -1000.f;
389 
390  // If leninact is positive and lenincat is negative,
391  // the transition plane is to the right of the track.
392  if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 &&
393  sr->trk.kalman.tracks[ibesttrk].lenincat < 0)
394  return 0.f;
395 
396  // If leninact is positive and lenincat is positive,
397  // the transition plane is in the middle of the track.
398  if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 &&
399  sr->trk.kalman.tracks[ibesttrk].lenincat > 0)
400  return float(sr->trk.kalman.tracks[ibesttrk].lenincat / 100.);
401 
402  // If leninact is negative and lenincat is positive,
403  // the transition plane is to the left of the track.
404  if(sr->trk.kalman.tracks[ibesttrk].leninact < 0 &&
405  sr->trk.kalman.tracks[ibesttrk].lenincat > 0)
406  return float((sr->trk.kalman.tracks[ibesttrk].leninact / 100.)
407  +(sr->trk.kalman.tracks[ibesttrk].lenincat / 100.));
408 
409  return -1000.f;
410  });
411 
412 
413  // Muon Energy for Inclusive numuCC
414  //------------------------------------
415  const Var kIncXsecMuonE([](const caf::SRProxy *sr)
416  {
417  float muonE = 0.0;
418  float muonEact = 0.0;
419  float muonEcat = 0.0;
420  float muonEactandcat = 0.0;
421 
422  float trkLenAct = 0.f;
423  float trkLenCat = 0.f;
424 
425  trkLenAct = kTrkLenAct(sr);
426  trkLenCat = kTrkLenCat(sr);
427 
428  int ibesttrk = kBestTrack(sr);
429  if(sr->trk.kalman.ntracks < 1)
430  return -1000.f;
431 
432  if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
433  return -1000.f;
434 
435  if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 &&
436  sr->trk.kalman.tracks[ibesttrk].lenincat < 0)
437  muonEact = MuonEAct(trkLenAct);
438  else if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 &&
439  sr->trk.kalman.tracks[ibesttrk].lenincat > 0){
440  muonEcat = MuonECat(trkLenCat);
441  muonEactandcat = MuonEActandCat(trkLenAct);
442  muonE = muonEactandcat + muonEcat;
443  }
444 
445  return muonE + muonEact;
446  });
447 
448 
449  // Hadronic Energy for Inclusive numuCC
450  //---------------------------------------
451 
452  const Var kIncXsecHadE([](const caf::SRProxy* sr)
453  {
454  int ibesttrk = kBestTrack(sr);
455 
456  if(sr->trk.kalman.ntracks < 1)
457  return -1000.f;
458 
459  if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
460  return -1000.f;
461 
462  // Extra energy (hadronic contamination) associated with muon track
463  float ExtraHadE = sr->trk.kalman.tracks[ibesttrk].overlapE;
464 
465  // calorimetric slice energy - Calorimetric Energy of Muon Track
466  float CalhadE = sr->slc.calE - sr->trk.kalman.tracks[ibesttrk].calE;
467 
468  // Add calorimetric hadE and Overlap energy in the muon track
469  float vishadE = CalhadE + ExtraHadE;
470 
471  float hadE = 0.0;
472  hadE = VisibleHadE(vishadE);
473 
474  return hadE;
475  });
476 
477  //-------------------------------------------------------------------------------------------------
478  // Standard analysis vars and axes
479 
480  // Reconstructed candidate muon angle wrt beam
481  const Var kRecoMuCostheta([](const caf::SRProxy* sr)
482  {
483  int ibesttrk = kBestTrack(sr);
484 
485  if(sr->trk.kalman.ntracks < 1)
486  return -1000.f;
487 
488  if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
489  return -1000.f;
490 
491  TVector3 dir = sr->trk.kalman.tracks[ibesttrk].dir;
492  TVector3 beamdir = NuMIBeamDirection(sr->hdr.det);
493  return (float)dir.Dot(beamdir);
494  });
495 
496  // Angle systematics (+/- 0.0025 mrad wrt beam)
497  const Var kRecoMuTheta([](const caf::SRProxy* sr)
498  {
499  int ibesttrk = kBestTrack(sr);
500 
501  if (sr->trk.kalman.ntracks < 1)
502  return -1000.f;
503 
504  if (ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
505  return -1000.f;
506 
507  TVector3 dir = sr->trk.kalman.tracks[ibesttrk].dir;
508  TVector3 beamdir = NuMIBeamDirection(sr->hdr.det);
509  return (float)dir.Angle(beamdir);
510  });
511 
512  const Var kRecoMuCosthetaUp([](const caf::SRProxy* sr)
513  {
514  return (float)std::cos(kRecoMuThetaUp(sr));
515  });
516 
517  const Var kRecoMuCosthetaDw([](const caf::SRProxy* sr)
518  {
519  return (float)std::cos(kRecoMuThetaDw(sr));
520  });
521 
522  const NuTruthVar kTrueEST([](const caf::SRNeutrinoProxy* nu)
523  {
524  return nu->E;
525  });
526 
527  const NuTruthVar kTrueMuEST([](const caf::SRNeutrinoProxy* nu)
528  {
529  float MuE = -5.;
530 
531  if(abs(nu->pdg) != 14 || !nu->iscc)
532  return MuE;
533 
534  int nprims = nu->prim.size();
535  for(int iprim = 0; iprim < nprims; iprim++){
536  if(abs(nu->prim[iprim].pdg) == 13){
537  MuE = nu->prim[iprim].p.T();
538  }
539  }//end loop over primaries
540 
541  return MuE;
542  });
543 
544  const NuTruthVar kTrueMuKEST([](const caf::SRNeutrinoProxy* nu)
545  {
546 
547  float ke = -5;
548 
549  if(abs(nu->pdg) != 14 || !nu->iscc)
550  return ke;
551 
552  int nprims = nu->prim.size();
553  for(int iprim = 0; iprim < nprims; iprim++){
554  if(abs(nu->prim[iprim].pdg) == 13){
555  double E= nu->prim[iprim].p.T();
556  ke = E-MuonMass();
557  }
558  }//end loop over primaries
559 
560  return ke;
561  });
562 
563  // True muon cos theta wrt beam direction
564 
566  {
567  if(abs(nu->pdg) != 14 || !nu->iscc)
568  return -5.0;
569 
570  int nprims = nu->prim.size();
571  for(int iprim = 0; iprim < nprims; iprim++){
572  if(abs(nu->prim[iprim].pdg) == 13){
573  TVector3 mudir = nu->prim[iprim].p.Vect();
574  TVector3 beamdir = NuMIBeamDirection(caf::kNEARDET);
575 
576  return mudir.Unit().Dot(beamdir.Unit());
577  }
578  }//end loop over primaries
579  return -5.0;
580  });
581 
582  // True hadronic y
583  const NuTruthVar kTrueYST([](const caf::SRNeutrinoProxy* sr)
584  {
585  return float(sr->y);
586  });
587 
588 
589  // True hadronic energy
590  const NuTruthVar kTrueHadEST([](const caf::SRNeutrinoProxy* sr)
591  {
592  return float(sr->y * sr->E);
593  });
594 
595  // Reconstructed available energy (Determined by Travis Olson)
596  const Var kRecoEavail([] (const caf::SRProxy *sr)
597  {
598  return 1.68*kNumuHadVisE(sr)+0.02*pow(kNumuHadVisE(sr),2);
599  });
600 
601 
602 
606  kRecoMuKE,
607  mukebins);
608 
612  kTrueMuKEST,
613  mukebins);
614 
618  kTrueMuKEST,
619  mukebins,
620  kTrueEST,
621  enubins);
622 
626  kRecoMuKE,
627  mukebins,
628  kRecoE,
629  enubins);
630 
634  kTrueMuKEST,
635  mukebins,
637  eavailbins);
638 
642  kRecoMuKE,
643  mukebins,
644  kRecoEavail,
645  eavailbins);
646 
647 /////////////////////////////////////////
648 /// \brief Systematic Shift Vars and HistAxes
649 /////////////////////////////////////////
650 
651  // Vars related to muon energy scale systematics:
652  const Var kIncXsecMuonEUp([](const caf::SRProxy *sr)
653  {
654  //factors to make the muon energy up
655  float mue_act_up = 1.0079;
656  float mue_cat_up = 1.0120;
657 
658  float muonE = 0.0;
659  float muonEact = 0.0;
660  float muonEcat = 0.0;
661  float muonEactandcat = 0.0;
662 
663  float trkLenAct = 0.f;
664  float trkLenCat = 0.f;
665 
666  trkLenAct = kTrkLenAct(sr);
667  trkLenCat = kTrkLenCat(sr);
668 
669  int ibesttrk = kBestTrack(sr);
670  if(sr->trk.kalman.ntracks < 1)
671  return -1000.f;
672 
673  if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
674  return -1000.f;
675 
676  if( sr->trk.kalman.tracks[ibesttrk].leninact > 0 &&
677  sr->trk.kalman.tracks[ibesttrk].lenincat < 0)
678  muonEact = mue_act_up * MuonEAct(trkLenAct);
679  else if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 &&
680  sr->trk.kalman.tracks[ibesttrk].lenincat > 0){
681  muonEcat = MuonECat(trkLenCat);
682  muonEactandcat = MuonEActandCat(trkLenAct);
683  muonE = mue_act_up*muonEactandcat + mue_cat_up*muonEcat;
684  }
685 
686  return muonE + muonEact;
687  });
688 
689  const Var kIncXsecMuonEDw([](const caf::SRProxy *sr)
690  {
691  //factors to make the muon energy down
692  float mue_act_dw = 0.9921;
693  float mue_cat_dw = 0.9880;
694 
695  float muonE = 0.0;
696  float muonEact = 0.0;
697  float muonEcat = 0.0;
698  float muonEactandcat = 0.0;
699 
700  float trkLenAct = 0.f;
701  float trkLenCat = 0.f;
702 
703  trkLenAct = kTrkLenAct(sr);
704  trkLenCat = kTrkLenCat(sr);
705 
706  int ibesttrk = kBestTrack(sr);
707  if(sr->trk.kalman.ntracks < 1)
708  return -1000.f;
709 
710  if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
711  return -1000.f;
712 
713  if( sr->trk.kalman.tracks[ibesttrk].leninact > 0 &&
714  sr->trk.kalman.tracks[ibesttrk].lenincat < 0)
715  muonEact = mue_act_dw * MuonEAct(trkLenAct);
716  else if(sr->trk.kalman.tracks[ibesttrk].leninact > 0 &&
717  sr->trk.kalman.tracks[ibesttrk].lenincat > 0){
718  muonEcat = MuonECat(trkLenCat);
719  muonEactandcat = MuonEActandCat(trkLenAct);
720  muonE = mue_act_dw*muonEactandcat + mue_cat_dw*muonEcat;
721  }
722 
723  return muonE + muonEact;
724  });
725 
726  // 3D vars for muon energy shift (Enu as 3rd axis)
730  kRecoMuKEUp,
731  mukebins,
732  kRecoE,
733  enubins);
737  kRecoMuKEDw,
738  mukebins,
739  kRecoE,
740  enubins);
744  kRecoMuKEUp,
745  mukebins,
746  kRecoEUp,
747  enubins);
751  kRecoMuKEDw,
752  mukebins,
753  kRecoEDw,
754  enubins);
755 
756  // 3D vars for muon energy shift (Eavail as 3rd axis)
760  kRecoMuKEUp,
761  mukebins,
762  kRecoEavail,
763  eavailbins);
767  kRecoMuKEDw,
768  mukebins,
769  kRecoEavail,
770  eavailbins);
771 
772  // 3D vars for muon angle shift (Eavail as 3rd axis)
776  kRecoMuKE,
777  mukebins,
778  kRecoEavail,
779  eavailbins);
780 
784  kRecoMuKE,
785  mukebins,
786  kRecoEavail,
787  eavailbins);
788 
789  //-------------------------------------------------------------------------------------------------
790 
791 
792  const Var kDummy([](const caf::SRProxy* sr)
793  {
794  return 1;
795  });
796 
797 
798 
799  const Var ksigResdown([](const caf::SRProxy* sr)
800  {
801  if(sr->mc.nnu == 0)
802  return 0.;
803 
804  assert(sr->mc.nnu == 1);
805 
806  if( (sr->mc.nu[0].iscc) && sr->mc.nu[0].mode==caf::kRes)
807  return 0.909090;
808 
809  else return 1.;
810  });
811 
812  const Var ksigDISup([](const caf::SRProxy* sr)
813  {
814  if(sr->mc.nnu == 0)
815  return 0.;
816 
817  assert(sr->mc.nnu == 1);
818 
819  if( (sr->mc.nu[0].iscc) && sr->mc.nu[0].mode==caf::kDIS)
820  return 1.1;
821 
822  else return 1.;
823  });
824 
825  const Var kbkgDISdown([](const caf::SRProxy* sr)
826  {
827  if(sr->mc.nnu == 0)
828  return 0.;
829 
830  assert(sr->mc.nnu == 1);
831 
832  if( (!sr->mc.nu[0].iscc) && sr->mc.nu[0].mode==caf::kDIS)
833  return 0.909090;
834 
835  else return 1.;
836  });
837 
838 
839  //Number of prongs
840  const Var kNProngs([](const caf::SRProxy* sr)
841  {
842  return (sr->vtx.elastic[0].fuzzyk.npng);
843  });
844 
845  // True Q2
846 
847  const Var kQsqr([](const caf::SRProxy* sr)
848  {
849  return (sr->mc.nnu == 0) ? 0. : float(sr->mc.nu[0].q2);
850  });
851 
852  // True W2
853 
854  const Var kWsqr([](const caf::SRProxy* sr)
855  {
856  return (sr->mc.nnu == 0) ? 0. : float(sr->mc.nu[0].W2);
857  });
858 
859  const Var kweightedten([](const caf::SRProxy*){return 1.1;});
860 
861  // Primary muon momentum
862 
863  const Var kTrueMuP([](const caf::SRProxy* sr)
864  {
865  if(sr->mc.nnu < 1)
866  return -5.0;
867 
868  if(abs(sr->mc.nu[0].pdg) != 14 || !sr->mc.nu[0].iscc)
869  return -5.0;
870  // unsigned int ipartmu = 0;
871  int nprims = sr->mc.nu[0].prim.size();
872  for(int iprim = 0; iprim < nprims; iprim++){
873  if(abs(sr->mc.nu[0].prim[iprim].pdg) == 13){
874  // ipartmu = iprim;
875  return sr->mc.nu[0].prim[iprim].p.Vect().Mag();
876  }
877  }//end loop over primaries
878  return -5.0;
879  });
880 
881 
882  // Primary muon transverse momentum wrt beam direction
883 
884  const Var kTrueMuPt([](const caf::SRProxy* sr)
885  {
886  if(sr->mc.nnu < 1)
887  return -5.0;
888 
889  if(abs(sr->mc.nu[0].pdg) != 14 || !sr->mc.nu[0].iscc)
890  return -5.0;
891 
892  TVector3 beamdir = NuMIBeamDirection(sr->hdr.det);
893  int nprims = sr->mc.nu[0].prim.size();
894 
895  for(int iprim = 0; iprim < nprims; iprim++){
896  if(abs(sr->mc.nu[0].prim[iprim].pdg) == 13){
897  TVector3 p = sr->mc.nu[0].prim[iprim].p.Vect();
898  return (p-p.Dot(beamdir)*beamdir).Mag();
899  }
900  }//end loop over primaries
901  return -5.0;
902  });
903 
904  // Primary muon momentum along the beam direction
905 
906  const Var kTrueMuPz([](const caf::SRProxy* sr)
907  {
908  if(sr->mc.nnu < 1)
909  return -5.0;
910 
911  if(abs(sr->mc.nu[0].pdg) != 14 || !sr->mc.nu[0].iscc)
912  return -5.0;
913 
914  TVector3 beamdir = NuMIBeamDirection(sr->hdr.det);
915  int nprims = sr->mc.nu[0].prim.size();
916 
917  for(int iprim = 0; iprim < nprims; iprim++){
918  if(abs(sr->mc.nu[0].prim[iprim].pdg) == 13){
919  TVector3 p = sr->mc.nu[0].prim[iprim].p.Vect();
920  return p.Dot(beamdir);
921  }
922  }//end loop over primaries
923  return -5.0;
924  });
925 
926  // Reconstructed muon momentum
927  const Var kRecoMuPsqr([](const caf::SRProxy * sr){
928  float muonpsqr = std::pow(kRecoMuE(sr), 2) - kMuonMass2(sr);
929  return (std::max(muonpsqr, (float)0.0)); // Return muonpsqr if > 0, otherwise 0.
930  });
931 
932  const Var kRecoMuPsqrUp([](const caf::SRProxy * sr){
933  float muonpsqr = std::pow(kRecoMuEUp(sr), 2) - kMuonMass2(sr);
934  return (std::max(muonpsqr, (float)0.0)); // Return muonpsqr if > 0, otherwise 0.
935  });
936 
937  const Var kRecoMuPsqrDw([](const caf::SRProxy * sr){
938  float muonpsqr = std::pow(kRecoMuEDw(sr), 2) - kMuonMass2(sr);
939  return (std::max(muonpsqr, (float)0.0)); // Return muonpsqr if > 0, otherwise 0.
940  });
941 
942 
943  // Reconstructed muon tranvserse momentum wrt beam direction
944 
945  const Var kRecoMuPt([](const caf::SRProxy* sr)
946  {
947  int bestidx = sr->trk.kalman.idxremid;
948  int nkals = sr->trk.kalman.ntracks;
949 
950  if( nkals < 1 ||
951  nkals < bestidx ||
952  bestidx < 0 )
953  return -5.0f;
954 
955  double E = TAMuE(sr->trk.kalman.tracks[0].len);
956 
957  TVector3 dir = sr->trk.kalman.tracks[0].dir;
958  TVector3 pvec = TMath::Sqrt( E*E - util::sqr(MuonMass()))*dir;
959  TVector3 beamdir = NuMIBeamDirection(sr->hdr.det);
960 
961  return float((pvec - pvec.Dot(beamdir)*beamdir).Mag());
962  });
963 
964 
965  // Reconstructed muon momentum along beam direction
966 
967  const Var kRecoMuPz([](const caf::SRProxy* sr)
968  {
969 
970  double E = TAMuE(sr->trk.kalman.tracks[0].len);
971 
972  TVector3 dir = sr->trk.kalman.tracks[0].dir;
973  TVector3 pvec = TMath::Sqrt( E*E - util::sqr(MuonMass()))*dir;
974  TVector3 beamdir = NuMIBeamDirection(sr->hdr.det);
975 
976  return pvec.Dot(beamdir);
977  });
978 
979  // Fractional muon momentum resolution
980  const Var kResMuP([](const caf::SRProxy* sr)
981  {
982  if(sr->mc.nnu < 1)
983  return -5.0f;
984 
985  if( abs(sr->mc.nu[0].pdg) != 14 ||
986  !sr->mc.nu[0].iscc)
987  return -5.0f;
988  int bestidx = sr->trk.kalman.idxremid;
989  int nkals = sr->trk.kalman.ntracks;
990  //int bestidx = sr->sel.remid.bestidx;
991  //int nmups = sr->energy.mutrkE.size();
992  if( nkals < 1 ||
993  nkals < bestidx ||
994  bestidx < 0 )
995  return -5.0f;
996 
997  float mup = -5, muRecoE;
998  int nprims = sr->mc.nu[0].prim.size();
999  for(int iprim = 0; iprim < nprims; iprim++){
1000  if(abs(sr->mc.nu[0].prim[iprim].pdg) == 13){
1001  return float(sr->mc.nu[0].prim[iprim].p.Vect().Mag());
1002  }
1003  // mup = sr->mc.nu[0].prim[iprim].p.P();
1004  // break;
1005  }//end loop over primaries
1006 
1007  if(mup == -5)
1008  return -5.0f;
1009 
1010  muRecoE = TAMuE(sr->trk.kalman.tracks[0].len);
1011 
1012  double recop = TMath::Sqrt( muRecoE*muRecoE - util::sqr(MuonMass()));
1013 
1014  return float((recop-mup)/mup);
1015  });
1016 
1017 
1018  // Fractional muon transverse momentum resolution
1019  const Var kResMuPt([](const caf::SRProxy* sr)
1020  {
1021  if(sr->mc.nnu < 1)
1022  return -5.0f;
1023 
1024  if( abs(sr->mc.nu[0].pdg) != 14 ||
1025  !sr->mc.nu[0].iscc)
1026  return -5.0f;
1027  int bestidx = sr->trk.kalman.idxremid;
1028  int nkals = sr->trk.kalman.ntracks;
1029  //int bestidx = sr->sel.remid.bestidx;
1030  //int nmups = sr->energy.mutrkE.size();
1031  TVector3 beamdir = NuMIBeamDirection(sr->hdr.det);
1032 
1033  if( nkals < 1 ||
1034  nkals < bestidx ||
1035  bestidx < 0 )
1036  return -5.0f;
1037 
1038  float mupt = -5, muRecoE;
1039  int nprims = sr->mc.nu[0].prim.size();
1040  for(int iprim = 0; iprim < nprims; iprim++){
1041  if(abs(sr->mc.nu[0].prim[iprim].pdg) == 13){
1042  TVector3 p = sr->mc.nu[0].prim[iprim].p.Vect();
1043  mupt = (p-p.Dot(beamdir)*beamdir).Mag();
1044  }
1045  break;
1046  }//end loop over primaries
1047 
1048  if(mupt == -5)
1049  return -5.0f;
1050 
1051  muRecoE = TAMuE(sr->trk.kalman.tracks[0].len);
1052 
1053  TVector3 dir = sr->trk.kalman.tracks[0].dir;
1054  TVector3 recop = TMath::Sqrt( muRecoE*muRecoE - util::sqr(MuonMass()))*dir;
1055  double recopt = (recop-recop.Dot(beamdir)*beamdir).Mag();
1056  return float((recopt-mupt)/mupt);
1057  });
1058 
1059 
1060  // Fractional resolution of muon momentum along beam direction
1061  const Var kResMuPz([](const caf::SRProxy* sr)
1062  {
1063  if(sr->mc.nnu < 1)
1064  return -5.0f;
1065 
1066  if( abs(sr->mc.nu[0].pdg) != 14 ||
1067  !sr->mc.nu[0].iscc)
1068  return -5.0f;
1069 
1070  TVector3 beamdir = NuMIBeamDirection(sr->hdr.det);
1071 
1072  float mupz = -5, muRecoE;
1073  int nprims = sr->mc.nu[0].prim.size();
1074  for(int iprim = 0; iprim < nprims; iprim++){
1075  if(abs(sr->mc.nu[0].prim[iprim].pdg) == 13){
1076  TVector3 p = sr->mc.nu[0].prim[iprim].p.Vect();
1077  mupz = p.Dot(beamdir);
1078  }
1079  break;
1080  }//end loop over primaries
1081 
1082  if(mupz == -5)
1083  return -5.0f;
1084 
1085  muRecoE = TAMuE(sr->trk.kalman.tracks[0].len);
1086 
1087  TVector3 dir = sr->trk.kalman.tracks[0].dir;
1088  TVector3 recop = TMath::Sqrt( muRecoE*muRecoE - util::sqr(MuonMass()))*dir;
1089 
1090  double recopz = recop.Dot(beamdir);
1091  return float((recopz-mupz)/mupz);
1092  });
1093 
1094 
1095  // Remid
1096  const Var kRemid([](const caf::SRProxy* sr){
1097  return (sr->sel.remid.pid);
1098  });
1099 
1100  // Remid variables
1101  const Var kScatLL = SIMPLEVAR(sel.remid.scatllh);
1102  const Var kDedxLL = SIMPLEVAR(sel.remid.dedxllh);
1103  const Var kMeasfrac = SIMPLEVAR(sel.remid.measfrac);
1104 
1105  // Length of candidate muon track
1106  const Var kMuLength([](const caf::SRProxy* sr)
1107  {
1108  // return sr->trk.kalman[bestidx].len;
1109  return sr->trk.kalman.tracks[0].len/100;
1110  });
1111 
1112  // Energy per length of the candidate muon track
1113  const Var kEPerLength([](const caf::SRProxy* sr)
1114  {
1115  if (sr->trk.kalman.ntracks > 0 &&
1116  sr->trk.kalman.idxremid != 999)
1117  return TAMuE(sr->trk.kalman.tracks[0].len)/sr->trk.kalman.tracks[0].len;
1118 
1119  return -5.f;
1120  });
1121 
1122 
1123  // Resolution of candidate muon angle wrt beam
1124  const Var kResMuCostheta([](const caf::SRProxy* sr)
1125  {
1126  if(sr->mc.nnu < 1)
1127  return -5.0;
1128 
1129  if( abs(sr->mc.nu[0].pdg) != 14 ||
1130  !sr->mc.nu[0].iscc)
1131  return -5.0;
1132 
1133  int bestidx = sr->trk.kalman.idxremid;
1134  int nkals = sr->trk.kalman.ntracks;
1135  if( nkals < 1 ||
1136  nkals < bestidx ||
1137  bestidx < 0 )
1138  return -5.0;
1139 
1140  double mu = 0, muReco = 0;
1141  TVector3 beamdir = NuMIBeamDirection(sr->hdr.det);
1142  int nprims = sr->mc.nu[0].prim.size();
1143  for(int iprim = 0; iprim < nprims; iprim++){
1144  if(abs(sr->mc.nu[0].prim[iprim].pdg) == 13){
1145  TVector3 mudir = sr->mc.nu[0].prim[iprim].p.Vect();
1146  mu = beamdir.Dot(mudir);
1147  mu = beamdir.Dot(sr->mc.nu[0].prim[iprim].p.Vect().Unit());
1148  }
1149  }//end loop over primaries
1150 
1151  if(mu == 0)
1152  return -5.0;
1153 
1154  TVector3 dir = sr->trk.kalman.tracks[0].dir;
1155 
1156  muReco = beamdir.Dot(dir);
1157 
1158  return (muReco-mu)/mu;
1159  });
1160 
1161 
1162 
1163 
1164 
1165  // True visible energy fraction of the muon
1166  const Var kTrueEFrac([](const caf::SRProxy* sr)
1167  {
1168  if(sr->mc.nnu < 1)
1169  return -5.0f;
1170 
1171  return float(sr->mc.nu[0].visE/sr->mc.nu[0].E);
1172  });
1173 
1174  // 2D var : energy resolution vs true energy
1176  kResE,
1177  resbins,
1178  kTrueE,
1179  ebins);
1180 
1181  // 2D var : energy resolution vs reco energy
1183  kResE,
1184  resbins,
1185  kRecoE,
1186  ebins);
1187 
1188  // 2D var : muon momentum resolution vs true muon momentum
1190  kResMuP,
1191  resbins,
1192  kTrueMuP,
1193  ebins);
1194 
1195  // 2D var : muon momentum resolution vs reco muon momentum
1197  kResMuP,
1198  resbins,
1199  kRecoMuP,
1200  ebins);
1201 
1202  // 2D var : muon transverse momentum resolution vs true muon transverse momentum
1204  kResMuPt,
1205  resbins,
1206  kTrueMuPt,
1207  ebins);
1208 
1209  // 2D var : muon transverse momentum resolution vs reco muon transverse momentum
1211  kResMuPt,
1212  resbins,
1213  kRecoMuPt,
1214  ebins);
1215 
1216  // 2D var : muon longitudinal momentum resolution vs true muon longitudinal momentum
1218  kResMuPz,
1219  resbins,
1220  kTrueMuPz,
1221  ebins);
1222 
1223  // 2D var : muon longitudinal momentum resolution vs reco muon longitudinal momentum
1225  kResMuPz,
1226  resbins,
1227  kRecoMuPz,
1228  ebins);
1229 
1230 
1233  resbins,
1235  angbins);
1236 
1239  resbins,
1241  angbins);
1242 
1244  kRecoMuP,
1245  ebins,
1247  angbins);
1248 
1250  kRecoMuPz,
1251  ebins,
1252  kRecoMuPt,
1253  ebins);
1254 
1256  kTrueMuP,
1257  ebins,
1259  angbins);
1260 
1262  kTrueMuPz,
1263  ebins,
1264  kTrueMuPt,
1265  ebins);
1266 
1267  // background weights
1268  const Var kUncontainedBGWeight([](const caf::SRProxy* sr)
1269  {
1270  if(sr->trk.kalman.ntracks < 1)
1271  return -1000.0;
1272 
1273  return -0.00032*sr->trk.kalman.tracks[0].start.Y() + 1.14;
1274  });
1275 
1276 
1277  // Var For Energy estimation
1278  const Var kTrueCatcherE = SIMPLEVAR(energy.numu.mc.truemuoncatcherE);
1279 
1280  const Var kReMIdTrkLenAct ([](const caf::SRProxy* sr)
1281  {
1282  return (sr->energy.numu.ndtrklenact / 100); // in m
1283  });
1284 
1285  const Var kReMIdTrkLenCat ([](const caf::SRProxy* sr)
1286  {
1287  return (sr->energy.numu.ndtrklencat / 100); // in m
1288  });
1289 
1290 }// end of namespace
const Var kRecoMuKEVsCosVsEavailUp
caf::Proxy< size_t > npng
Definition: SRProxy.h:2020
Near Detector underground.
Definition: SREnums.h:10
T max(const caf::Proxy< T > &a, T b)
const NuTruthVar kTrueMuEST([](const caf::SRNeutrinoProxy *nu){float MuE=-5.;if(abs(nu->pdg)!=14||!nu->iscc) return MuE;int nprims=nu->prim.size();for(int iprim=0;iprim< nprims;iprim++){if(abs(nu->prim[iprim].pdg)==13){MuE=nu->prim[iprim].p.T();}}return MuE;})
const Var kRecoMuKEVsCosVsEavail
const Var kTrueMuPt([](const caf::SRProxy *sr){if(sr->mc.nnu< 1) return-5.0;if(abs(sr->mc.nu[0].pdg)!=14||!sr->mc.nu[0].iscc) return-5.0;TVector3 beamdir=NuMIBeamDirection(sr->hdr.det);int nprims=sr->mc.nu[0].prim.size();for(int iprim=0;iprim< nprims;iprim++){if(abs(sr->mc.nu[0].prim[iprim].pdg)==13){TVector3 p=sr->mc.nu[0].prim[iprim].p.Vect();return(p-p.Dot(beamdir)*beamdir).Mag();}}return-5.0;})
caf::Proxy< caf::SRFuzzyK > fuzzyk
Definition: SRProxy.h:2041
const Var kReMIdTrkLenAct([](const caf::SRProxy *sr){return(sr->energy.numu.ndtrklenact/100);})
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
TH2F * hprob_dedx_vs_scat_bg
const Var kResMuCostheta([](const caf::SRProxy *sr){if(sr->mc.nnu< 1) return-5.0;if(abs(sr->mc.nu[0].pdg)!=14|| !sr->mc.nu[0].iscc) return-5.0;int bestidx=sr->trk.kalman.idxremid;int nkals=sr->trk.kalman.ntracks;if(nkals< 1|| nkals< bestidx|| bestidx< 0) return-5.0;double mu=0, muReco=0;TVector3 beamdir=NuMIBeamDirection(sr->hdr.det);int nprims=sr->mc.nu[0].prim.size();for(int iprim=0;iprim< nprims;iprim++){if(abs(sr->mc.nu[0].prim[iprim].pdg)==13){TVector3 mudir=sr->mc.nu[0].prim[iprim].p.Vect();mu=beamdir.Dot(mudir);mu=beamdir.Dot(sr->mc.nu[0].prim[iprim].p.Vect().Unit());}}if(mu==0) return-5.0;TVector3 dir=sr->trk.kalman.tracks[0].dir;muReco=beamdir.Dot(dir);return(muReco-mu)/mu;})
const Var kIncXsecMuonEUp([](const caf::SRProxy *sr){float mue_act_up=1.0079;float mue_cat_up=1.0120;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=mue_act_up *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=mue_act_up *muonEactandcat+mue_cat_up *muonEcat;}return muonE+muonEact;})
Systematic Shift Vars and HistAxes.
const NuTruthVar kTrueMuCosthetaST([](const caf::SRNeutrinoProxy *nu){if(abs(nu->pdg)!=14||!nu->iscc) return-5.0;int nprims=nu->prim.size();for(int iprim=0;iprim< nprims;iprim++){if(abs(nu->prim[iprim].pdg)==13){TVector3 mudir=nu->prim[iprim].p.Vect();TVector3 beamdir=NuMIBeamDirection(caf::kNEARDET);return mudir.Unit().Dot(beamdir.Unit());}}return-5.0;})
const Var kRecoMuKE
const Var kRecoE
Definition: NumuCCIncVars.h:97
const Var kTrueMuCostheta
const NuTruthVar kTrueHadEST([](const caf::SRNeutrinoProxy *sr){return float(sr->y *sr->E);})
const Var ksigDISup([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return 0.;assert(sr->mc.nnu==1);if((sr->mc.nu[0].iscc)&&sr->mc.nu[0].mode==caf::kDIS) return 1.1;else return 1.;})
const Var kQsqr([](const caf::SRProxy *sr){return(sr->mc.nnu==0)?0.:float(sr->mc.nu[0].q2);})
const Var kVisHadE([](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;return vishadE;})
Definition: NumuCCIncVars.h:57
const NuTruthVar kTrueMuKEVsCosST
const Var kResMuPt([](const caf::SRProxy *sr){if(sr->mc.nnu< 1) return-5.0f;if(abs(sr->mc.nu[0].pdg)!=14|| !sr->mc.nu[0].iscc) return-5.0f;int bestidx=sr->trk.kalman.idxremid;int nkals=sr->trk.kalman.ntracks; TVector3 beamdir=NuMIBeamDirection(sr->hdr.det);if(nkals< 1|| nkals< bestidx|| bestidx< 0) return-5.0f;float mupt=-5, muRecoE;int nprims=sr->mc.nu[0].prim.size();for(int iprim=0;iprim< nprims;iprim++){if(abs(sr->mc.nu[0].prim[iprim].pdg)==13){TVector3 p=sr->mc.nu[0].prim[iprim].p.Vect();mupt=(p-p.Dot(beamdir)*beamdir).Mag();}break;}if(mupt==-5) return-5.0f;muRecoE=TAMuE(sr->trk.kalman.tracks[0].len);TVector3 dir=sr->trk.kalman.tracks[0].dir;TVector3 recop=TMath::Sqrt(muRecoE *muRecoE-util::sqr(MuonMass()))*dir;double recopt=(recop-recop.Dot(beamdir)*beamdir).Mag();return float((recopt-mupt)/mupt);})
caf::Proxy< unsigned int > idxremid
Definition: SRProxy.h:1755
const Var kRecoMuKEVsCosVsEnu_onlyMuKEDw
TH2F * hprob_dedx_vs_scat_sig
caf::Proxy< size_t > ntracks
Definition: SRProxy.h:1756
caf::Proxy< caf::SRHeader > hdr
Definition: SRProxy.h:2119
const Binning resbins
Definition: NumuCCIncBins.h:31
const char * p
Definition: xmltok.h:285
caf::Proxy< caf::SRNumuEnergy > numu
Definition: SRProxy.h:213
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2108
_Var< T > Var2D(const _Var< T > &a, const Binning &binsa, const _Var< T > &b, const Binning &binsb)
Variable formed from two input variables.
Definition: Var.cxx:247
const Var kRecoMuEDw
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:597
constexpr T pow(T x)
Definition: pow.h:75
const Var kRecoEavail([](const caf::SRProxy *sr){return 1.68 *kNumuHadVisE(sr)+0.02 *pow(kNumuHadVisE(sr), 2);})
const Var kResMuP([](const caf::SRProxy *sr){if(sr->mc.nnu< 1) return-5.0f;if(abs(sr->mc.nu[0].pdg)!=14|| !sr->mc.nu[0].iscc) return-5.0f;int bestidx=sr->trk.kalman.idxremid;int nkals=sr->trk.kalman.ntracks; if(nkals< 1|| nkals< bestidx|| bestidx< 0) return-5.0f;float mup=-5, muRecoE;int nprims=sr->mc.nu[0].prim.size();for(int iprim=0;iprim< nprims;iprim++){if(abs(sr->mc.nu[0].prim[iprim].pdg)==13){return float(sr->mc.nu[0].prim[iprim].p.Vect().Mag());} }if(mup==-5) return-5.0f;muRecoE=TAMuE(sr->trk.kalman.tracks[0].len);double recop=TMath::Sqrt(muRecoE *muRecoE-util::sqr(MuonMass()));return float((recop-mup)/mup);})
caf::Proxy< float > pid
Definition: SRProxy.h:1115
const Var kRecoMuPsqrUp([](const caf::SRProxy *sr){float muonpsqr=std::pow(kRecoMuEUp(sr), 2)-kMuonMass2(sr);return(std::max(muonpsqr,(float) 0.0));})
const Binning ebins
Definition: NumuCCIncBins.h:24
const Var kRecoMuPz([](const caf::SRProxy *sr){double E=TAMuE(sr->trk.kalman.tracks[0].len);TVector3 dir=sr->trk.kalman.tracks[0].dir;TVector3 pvec=TMath::Sqrt(E *E-util::sqr(MuonMass()))*dir;TVector3 beamdir=NuMIBeamDirection(sr->hdr.det);return pvec.Dot(beamdir);})
const Var kResMuPz([](const caf::SRProxy *sr){if(sr->mc.nnu< 1) return-5.0f;if(abs(sr->mc.nu[0].pdg)!=14|| !sr->mc.nu[0].iscc) return-5.0f;TVector3 beamdir=NuMIBeamDirection(sr->hdr.det);float mupz=-5, muRecoE;int nprims=sr->mc.nu[0].prim.size();for(int iprim=0;iprim< nprims;iprim++){if(abs(sr->mc.nu[0].prim[iprim].pdg)==13){TVector3 p=sr->mc.nu[0].prim[iprim].p.Vect();mupz=p.Dot(beamdir);}break;}if(mupz==-5) return-5.0f;muRecoE=TAMuE(sr->trk.kalman.tracks[0].len);TVector3 dir=sr->trk.kalman.tracks[0].dir;TVector3 recop=TMath::Sqrt(muRecoE *muRecoE-util::sqr(MuonMass()))*dir;double recopz=recop.Dot(beamdir);return float((recopz-mupz)/mupz);})
const NuTruthVar kTrueMuKEST([](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 kMuStopY([](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;return sr->trk.kalman.tracks[ibesttrk].stop.Y();})
Definition: NumuCCIncVars.h:52
const Var kResVsTrueE
const Var kUncontainedBGWeight([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-1000.0;return-0.00032 *sr->trk.kalman.tracks[0].start.Y()+1.14;})
void abs(TH1 *hist)
const Var kResVsRecoMuP
caf::Proxy< short int > nnu
Definition: SRProxy.h:596
const Var kTrueMuP([](const caf::SRProxy *sr){if(sr->mc.nnu< 1) return-5.0;if(abs(sr->mc.nu[0].pdg)!=14||!sr->mc.nu[0].iscc) return-5.0;int nprims=sr->mc.nu[0].prim.size();for(int iprim=0;iprim< nprims;iprim++){if(abs(sr->mc.nu[0].prim[iprim].pdg)==13){return sr->mc.nu[0].prim[iprim].p.Vect().Mag();}}return-5.0;})
const Binning eavailbins
const Var kRecoMuE
const NuTruthVar kTrueMuKEVsCosVsEnuST
const Var kResE
const Var kTruePtVsPz
const Var kRecoPtVsPz
caf::Proxy< caf::SREnergyBranch > energy
Definition: SRProxy.h:2118
const Var kRecoMuCosthetaUp([](const caf::SRProxy *sr){return(float) std::cos(kRecoMuThetaUp(sr));})
const Var kMuonLLRPID([](const caf::SRProxy *sr){int bestidx=kMostMuonLikeTrackIdx(sr);float llr=-1e12;if(bestidx< 0||bestidx >=int(sr->trk.kalman.ntracks)) return llr;llr=MuonLLR(sr->trk.kalman.tracks[bestidx].dedxllh, sr->trk.kalman.tracks[bestidx].scatllh);return llr;})
Definition: NumuCCIncVars.h:36
const Var kRecoMuP
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
const Var kNumuHadVisE([](const caf::SRProxy *sr){return kNumuHadCalE(sr)+kNumuHadTrkE(sr);})
Definition: NumuVars.h:124
const Var kReMIdTrkLenCat([](const caf::SRProxy *sr){return(sr->energy.numu.ndtrklencat/100);})
_Var< T > Var3D(const _Var< T > &a, const Binning &binsa, const _Var< T > &b, const Binning &binsb, const _Var< T > &c, const Binning &binsc)
This is just like a Var2D, but useful for 3D Spectra.
Definition: Var.cxx:273
_Var< caf::SRProxy > Var
Representation of a variable to be retrieved from a caf::StandardRecord object.
Definition: Var.h:74
const Var kRecoMuCostheta([](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;TVector3 dir=sr->trk.kalman.tracks[ibesttrk].dir;TVector3 beamdir=NuMIBeamDirection(sr->hdr.det);return(float) dir.Dot(beamdir);})
const Var kTrueMuPz([](const caf::SRProxy *sr){if(sr->mc.nnu< 1) return-5.0;if(abs(sr->mc.nu[0].pdg)!=14||!sr->mc.nu[0].iscc) return-5.0;TVector3 beamdir=NuMIBeamDirection(sr->hdr.det);int nprims=sr->mc.nu[0].prim.size();for(int iprim=0;iprim< nprims;iprim++){if(abs(sr->mc.nu[0].prim[iprim].pdg)==13){TVector3 p=sr->mc.nu[0].prim[iprim].p.Vect();return p.Dot(beamdir);}}return-5.0;})
const NuTruthVar kTrueEST([](const caf::SRNeutrinoProxy *nu){return nu->E;})
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;})
Definition: NumuCCIncVars.h:75
caf::Proxy< caf::SRElastic > elastic
Definition: SRProxy.h:2100
const Var kIncXsecMuonEDw([](const caf::SRProxy *sr){float mue_act_dw=0.9921;float mue_cat_dw=0.9880;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=mue_act_dw *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=mue_act_dw *muonEactandcat+mue_cat_dw *muonEcat;}return muonE+muonEact;})
const Var kRecoMuKEVsCos
const Var kMuonLength([](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;return float(sr->trk.kalman.tracks[ibesttrk].len/100.);})
Definition: NumuCCIncVars.h:55
const Binning angbins
Definition: NumuCCIncBins.h:37
float VisibleHadE(float vishadE)
const Var kRecoMuKEVsCosVsEnu_onlyMuKEUp
#define SIMPLEVAR(CAFNAME)
For Vars where literally all you need is a single CAF variable.
Definition: Var.h:88
const Binning enubins
const Var kTrueE([](const caf::SRProxy *sr){return(sr->mc.nnu==0)?0.:float(sr->mc.nu[0].E);})
Definition: Vars.cxx:85
const Var kRecoMuKEVsCosVsEavailAngleUp
const NuTruthVar kTrueMuKEVsCosVsEavailST
const Var kMuStartX([](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;return sr->trk.kalman.tracks[ibesttrk].start.X();})
Definition: NumuCCIncVars.h:47
const Var kAveragededx40([](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].avedEdxlast40cm==-5) return-1000.f;return float(sr->trk.kalman.tracks[ibesttrk].avedEdxlast40cm);})
Definition: NumuCCIncVars.h:61
const Var kMeasfrac
const Var kRemid([](const caf::SRProxy *sr){return(sr->sel.remid.pid);})
caf::Proxy< caf::SRTrackBranch > trk
Definition: SRProxy.h:2127
const Var kResVsRecoMuPt
Float_t E
Definition: plot.C:20
float MuonEActandCat(double trklenactandcat)
const Var kRecoMuKEVsCosVsEnuDw
std::string getenv(std::string const &name)
const Var kweightedten([](const caf::SRProxy *){return 1.1;})
const Var kbkgDISdown([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return 0.;assert(sr->mc.nnu==1);if((!sr->mc.nu[0].iscc)&&sr->mc.nu[0].mode==caf::kDIS) return 0.909090;else return 1.;})
const Var kScatLL
const Var kMostMuonLikeTrackIdx([](const caf::SRProxy *sr){float maxllr=-1e12;int bestidx=-5;for(unsigned int itrk=0;itrk< sr->trk.kalman.ntracks;itrk++){float llr=MuonLLR(sr->trk.kalman.tracks[itrk].dedxllh, sr->trk.kalman.tracks[itrk].scatllh);if(llr > maxllr){maxllr=llr;bestidx=itrk;}}return bestidx;})
Definition: NumuCCIncVars.h:34
double energy
Definition: plottest35.C:25
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;})
Definition: NumuCCIncVars.h:73
const Var kDedxLL
const Var kHadEonTrack([](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;return(float) sr->trk.kalman.tracks[ibesttrk].overlapE;})
Definition: NumuCCIncVars.h:59
const Var kTruePVsCos
double MuonMass()
caf::StandardRecord * sr
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;})
Definition: NumuCCIncVars.h:87
const Var kRecoMuKEVsCosVsEnuUp
const Var kRecoEDw
caf::Proxy< caf::SRRemid > remid
Definition: SRProxy.h:1248
float TAMuE(double trklen)
Definition: NumuEFxs.cxx:377
float MuonEAct(double trklenact)
const NuTruthVar kTrueYST([](const caf::SRNeutrinoProxy *sr){return float(sr->y);})
const Var kRecoMuPt([](const caf::SRProxy *sr){int bestidx=sr->trk.kalman.idxremid;int nkals=sr->trk.kalman.ntracks;if(nkals< 1|| nkals< bestidx|| bestidx< 0) return-5.0f;double E=TAMuE(sr->trk.kalman.tracks[0].len);TVector3 dir=sr->trk.kalman.tracks[0].dir;TVector3 pvec=TMath::Sqrt(E *E-util::sqr(MuonMass()))*dir;TVector3 beamdir=NuMIBeamDirection(sr->hdr.det);return float((pvec-pvec.Dot(beamdir)*beamdir).Mag());})
const Var kRecoPVsCos
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;})
Definition: NumuCCIncVars.h:82
const Var kRecoMuKEDw
const Var kNProngs([](const caf::SRProxy *sr){return(sr->vtx.elastic[0].fuzzyk.npng);})
const Var kRecoMuThetaDw
const Var kResVsTrueMuPz
const Var kEPerLength([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks > 0 && sr->trk.kalman.idxremid!=999) return TAMuE(sr->trk.kalman.tracks[0].len)/sr->trk.kalman.tracks[0].len;return-5.f;})
const std::string path
Definition: plot_BEN.C:43
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2120
const Var kRecoMuKEVsCosVsEavailAngleDw
const Var kRecoMuKEVsCosVsEavailDw
caf::Proxy< caf::SRKalman > kalman
Definition: SRProxy.h:1775
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const Var kRecoMuThetaUp
const Binning mukebins
Definition: NumuCCIncBins.h:90
caf::Proxy< caf::SRSlice > slc
Definition: SRProxy.h:2124
const Var kRecoMuCosthetaDw([](const caf::SRProxy *sr){return(float) std::cos(kRecoMuThetaDw(sr));})
TDirectory * dir
Definition: macro.C:5
caf::Proxy< float > calE
Definition: SRProxy.h:1271
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
const Var kResVsTrueMuP
const Var kRecoEUp
const Var kMuStopZ([](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;return sr->trk.kalman.tracks[ibesttrk].stop.Z();})
Definition: NumuCCIncVars.h:53
caf::Proxy< float > ndtrklenact
Definition: SRProxy.h:185
const Var kRecoMuTheta([](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;TVector3 dir=sr->trk.kalman.tracks[ibesttrk].dir;TVector3 beamdir=NuMIBeamDirection(sr->hdr.det);return(float) dir.Angle(beamdir);})
float MuonECat(double trklencat)
const Var kRecoMuPsqr([](const caf::SRProxy *sr){float muonpsqr=std::pow(kRecoMuE(sr), 2)-kMuonMass2(sr);return(std::max(muonpsqr,(float) 0.0));})
float MuonLLR(const float dedxll, const float scatll)
float Mag() const
const Var kResVsRecoE
const Var kRecoMuPsqrDw([](const caf::SRProxy *sr){float muonpsqr=std::pow(kRecoMuEDw(sr), 2)-kMuonMass2(sr);return(std::max(muonpsqr,(float) 0.0));})
T cos(T number)
Definition: d0nt_math.hpp:78
assert(nhit_max >=nhit_nbins)
const Var kTrueCatcherE
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:333
const Var kDummy([](const caf::SRProxy *sr){return 1;})
const Var kMuStopX([](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;return sr->trk.kalman.tracks[ibesttrk].stop.X();})
Definition: NumuCCIncVars.h:51
const Var kMuLength([](const caf::SRProxy *sr){return sr->trk.kalman.tracks[0].len/100;})
caf::Proxy< std::vector< caf::SRKalmanTrack > > tracks
Definition: SRProxy.h:1758
caf::Proxy< caf::SRIDBranch > sel
Definition: SRProxy.h:2123
const Var kResVsTrueMuPt
const Var kMuStartZ([](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;return sr->trk.kalman.tracks[ibesttrk].start.Z();})
Definition: NumuCCIncVars.h:49
const NuTruthVar kTrueEavailST
const Var kMuStartY([](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;return sr->trk.kalman.tracks[ibesttrk].start.Y();})
Definition: NumuCCIncVars.h:48
caf::Proxy< caf::SRVertexBranch > vtx
Definition: SRProxy.h:2128
const Binning angbinsCustom
Definition: NumuCCIncBins.h:72
const Var ksigResdown([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return 0.;assert(sr->mc.nnu==1);if((sr->mc.nu[0].iscc)&&sr->mc.nu[0].mode==caf::kRes) return 0.909090;else return 1.;})
const Var kWsqr([](const caf::SRProxy *sr){return(sr->mc.nnu==0)?0.:float(sr->mc.nu[0].W2);})
const Var kResVsRecoMuPz
caf::Proxy< float > ndtrklencat
Definition: SRProxy.h:186
const Var kRecoMuKEUp
const Var kRecoMuEUp
const Var kTrueEFrac([](const caf::SRProxy *sr){if(sr->mc.nnu< 1) return-5.0f;return float(sr->mc.nu[0].visE/sr->mc.nu[0].E);})
const Var kBestTrack
Definition: NDXSecMuonPID.h:33
const Var kResVsTrueCostheta
caf::Proxy< caf::Det_t > det
Definition: SRProxy.h:231
const Var kRecoMuKEVsCosVsEnu
const Var kResVsRecoCostheta
const Var kMuonMass2([](const caf::SRProxy *sr){return util::sqr(MuonMass());})