NumuEFxs.cxx
Go to the documentation of this file.
2 
4 
5 #define NumuMuE2017FDW(f) \
6  [] (const caf::SRProxy *sr) -> double \
7  { \
8  if (sr->trk.kalman.ntracks == 0) { \
9  return 0; \
10  } \
11  \
12  double trkLen = sr->trk.kalman.tracks[0].len; \
13  bool ismc = sr->hdr.ismc; \
14  \
15  return f(trkLen, ismc); \
16  }
17 
18 #define NumuHadE2017FDW(f) \
19  [] (const caf::SRProxy *sr) -> double \
20  { \
21  if (sr->trk.kalman.ntracks == 0) { \
22  return 0; \
23  } \
24  \
25  bool ismc = sr->hdr.ismc; \
26  double hadVisE = sr->energy.numu.hadtrkE + sr->energy.numu.hadcalE; \
27  \
28  return f(hadVisE, ismc); \
29  }
30 
31 #define NumuHadE2017NDW(f) \
32  [] (const caf::SRProxy *sr) -> double \
33  { \
34  if (sr->trk.kalman.ntracks == 0) { \
35  return 0; \
36  } \
37  \
38  double hadVisE = sr->energy.numu.hadtrkE + sr->energy.numu.hadcalE; \
39  \
40  return f(hadVisE); \
41  }
42 
43 #define NumuMuE2017NDW(f_act, f_cat) \
44  [] (const caf::SRProxy *sr) -> double \
45  { \
46  if (sr->trk.kalman.ntracks == 0) { \
47  return 0; \
48  } \
49  \
50  double trkLenAct = sr->energy.numu.ndtrklenact; \
51  double trkLenCat = sr->energy.numu.ndtrklencat; \
52  \
53  if ( (sr->energy.numu.ndtrkcalactE == 0) \
54  && (sr->energy.numu.ndtrkcaltranE == 0)) \
55  { \
56  return 0; \
57  } \
58  \
59  return f_act(trkLenAct) + f_cat(trkLenCat); \
60  }
61 
62 
63 namespace ana
64 {
65 
66  // Energy Estimator trained on prod5 MC
67  Var initNumuMuE2020FDpXVar( const std::function<NumuEnergyFunc::prod5_fd_energy_p_func_t> &f )
68  {
69  return Var
70  ([f] (const caf::SRProxy *sr) -> double
71  {
72  if (sr->trk.kalman.ntracks == 0)
73  {
74  return -5;
75  }
76 
77  double trkLen = sr->trk.kalman.tracks[0].len;
78  bool ismc = sr->hdr.ismc;
79 
80  return f(trkLen, ismc);
81  }
82  );
83  }
84 
85  Var initNumuHadE2020FDpXVar( const std::function<NumuEnergyFunc::prod5_fd_energy_p_func_t> &f )
86  {
87  return Var
88  ([f] (const caf::SRProxy *sr) -> double
89  {
90  if (sr->trk.kalman.ntracks == 0)
91  {
92  return -5;
93  }
94 
95  double hadVisE = sr->energy.numu.hadtrkE + sr->energy.numu.hadcalE;
96  bool ismc = sr->hdr.ismc;
97 
98  return f(hadVisE, ismc);
99  }
100  );
101  }
102 
103  Var initNumuMuE2020NDpXVar( const std::function<NumuEnergyFunc::prod5_nd_energy_p_func_t> &f_act, const std::function<NumuEnergyFunc::prod5_nd_energy_p_func_t> &f_cat )
104  {
105  return Var
106  ([f_act,f_cat] (const caf::SRProxy *sr) -> double
107  {
108  if (sr->trk.kalman.ntracks == 0)
109  {
110  return -5;
111  }
112 
113  double trkLenAct = sr->energy.numu.ndtrklenact;
114  double trkLenCat = sr->energy.numu.ndtrklencat;
115 
116  if ( (sr->energy.numu.ndtrkcalactE == 0)
117  && (sr->energy.numu.ndtrkcaltranE == 0))
118  {
119  return -5;
120  }
121 
122  return f_act(trkLenAct) + f_cat(trkLenCat);
123  }
124  );
125  }
126 
127  Var initNumuHadE2020NDpXVar( const std::function<NumuEnergyFunc::prod5_nd_energy_p_func_t> &f )
128  {
129  return Var
130  ([f] (const caf::SRProxy *sr) -> double
131  {
132  if (sr->trk.kalman.ntracks == 0)
133  {
134  return -5;
135  }
136 
137  double hadVisE = sr->energy.numu.hadtrkE + sr->energy.numu.hadcalE;
138 
139  return f(hadVisE);
140  }
141  );
142  }
143 
144 
145  Var initNumuMuE2020Var( const std::function<NumuEnergyFunc::prod5_fd_energy_func_t> &f_fd, const std::function<NumuEnergyFunc::prod5_nd_energy_func_t> &f_nd_act, const std::function<NumuEnergyFunc::prod5_nd_energy_func_t> &f_nd_cat)
146  {
147  return Var
148  ([f_fd,f_nd_act,f_nd_cat] (const caf::SRProxy *sr) -> double
149  {
150  if (sr->trk.kalman.ntracks == 0)
151  {
152  return -5;
153  }
154 
155  bool isRHC = sr->spill.isRHC;
156  caf::Det_t det = sr->hdr.det;
157 
158  if (det == caf::kFARDET)
159  {
160  int run = sr->hdr.run;
161  bool ismc = sr->hdr.ismc;
162  double trkLen = sr->trk.kalman.tracks[0].len;
163 
164  return f_fd(trkLen, run, ismc, isRHC);
165  }
166 
167  if (det == caf::kNEARDET)
168  {
169  double trkLenAct = sr->energy.numu.ndtrklenact;
170  double trkLenCat = sr->energy.numu.ndtrklencat;
171 
172  if ( (sr->energy.numu.ndtrkcalactE == 0)
173  && (sr->energy.numu.ndtrkcaltranE == 0))
174  {
175  return -5;
176  }
177 
178  return f_nd_act(trkLenAct, isRHC) + f_nd_cat(trkLenCat, isRHC);
179  }
180 
181  /* neither near nor far detector */
182  return -5;
183  }
184  );
185  }
186 
187 
188  Var initNumuHadE2020Var( const std::function<NumuEnergyFunc::prod5_fd_energy_func_t> &f_fd, const std::function<NumuEnergyFunc::prod5_nd_energy_func_t> &f_nd)
189  {
190  return Var
191  ([f_fd,f_nd] (const caf::SRProxy *sr) -> double
192  {
193  if (sr->trk.kalman.ntracks == 0)
194  {
195  return -5;
196  }
197 
198  bool isRHC = sr->spill.isRHC;
199  caf::Det_t det = sr->hdr.det;
200  double hadVisE = sr->energy.numu.hadtrkE + sr->energy.numu.hadcalE;
201 
202  if (det == caf::kFARDET)
203  {
204  int run = sr->hdr.run;
205  bool ismc = sr->hdr.ismc;
206  return f_fd(hadVisE, run, ismc, isRHC);
207  }
208 
209  if (det == caf::kNEARDET)
210  {
211  return f_nd(hadVisE, isRHC);
212  }
213 
214  /* neither near nor far detector */
215  return -5;
216  }
217  );
218  }
219 
220 
221 
222 /*...........................................................................
223  Third Analysis Energy Estimator For ND
224  ...........................................................................*/
226  const std::function<NumuEnergyFunc::prod4_fd_energy_p_func_t> &f
227  )
228  {
229  return Var([f] (const caf::SRProxy *sr) -> double
230  {
231  if (sr->trk.kalman.ntracks == 0) {
232  return -5;
233  }
234 
235  double trkLen = sr->trk.kalman.tracks[0].len;
236  bool ismc = sr->hdr.ismc;
237 
238  return f(trkLen, ismc);
239  });
240  }
241 
243  const std::function<NumuEnergyFunc::prod4_fd_energy_p_func_t> &f
244  )
245  {
246  return Var([f] (const caf::SRProxy *sr) -> double
247  {
248  if (sr->trk.kalman.ntracks == 0) {
249  return -5;
250  }
251 
252  double hadVisE = sr->energy.numu.hadtrkE + sr->energy.numu.hadcalE;
253  bool ismc = sr->hdr.ismc;
254 
255  return f(hadVisE, ismc);
256  });
257  }
258 
260  const std::function<NumuEnergyFunc::prod4_nd_energy_p_func_t> &f_act,
261  const std::function<NumuEnergyFunc::prod4_nd_energy_p_func_t> &f_cat
262  )
263  {
264  return Var([f_act,f_cat] (const caf::SRProxy *sr) -> double
265  {
266  if (sr->trk.kalman.ntracks == 0) {
267  return -5;
268  }
269 
270  double trkLenAct = sr->energy.numu.ndtrklenact;
271  double trkLenCat = sr->energy.numu.ndtrklencat;
272 
273  if ( (sr->energy.numu.ndtrkcalactE == 0)
274  && (sr->energy.numu.ndtrkcaltranE == 0))
275  {
276  return -5;
277  }
278 
279  return f_act(trkLenAct) + f_cat(trkLenCat);
280  });
281  }
282 
284  const std::function<NumuEnergyFunc::prod4_nd_energy_p_func_t> &f
285  )
286  {
287  return Var([f] (const caf::SRProxy *sr) -> double
288  {
289  if (sr->trk.kalman.ntracks == 0) {
290  return -5;
291  }
292 
293  double hadVisE = sr->energy.numu.hadtrkE + sr->energy.numu.hadcalE;
294 
295  return f(hadVisE);
296  });
297  }
298 
299 
301  const std::function<NumuEnergyFunc::prod4_fd_energy_func_t> &f_fd,
302  const std::function<NumuEnergyFunc::prod4_nd_energy_func_t> &f_nd_act,
303  const std::function<NumuEnergyFunc::prod4_nd_energy_func_t> &f_nd_cat
304  )
305  {
306  return Var([f_fd,f_nd_act,f_nd_cat] (const caf::SRProxy *sr) -> double
307  {
308  if (sr->trk.kalman.ntracks == 0) {
309  return -5;
310  }
311 
312  bool isRHC = sr->spill.isRHC;
313  caf::Det_t det = sr->hdr.det;
314 
315  if (det == caf::kFARDET)
316  {
317  int run = sr->hdr.run;
318  bool ismc = sr->hdr.ismc;
319  double trkLen = sr->trk.kalman.tracks[0].len;
320 
321  return f_fd(trkLen, run, ismc, isRHC);
322  }
323 
324  if (det == caf::kNEARDET)
325  {
326  double trkLenAct = sr->energy.numu.ndtrklenact;
327  double trkLenCat = sr->energy.numu.ndtrklencat;
328 
329  if ( (sr->energy.numu.ndtrkcalactE == 0)
330  && (sr->energy.numu.ndtrkcaltranE == 0))
331  {
332  return -5;
333  }
334 
335  return f_nd_act(trkLenAct, isRHC) + f_nd_cat(trkLenCat, isRHC);
336  }
337 
338  /* neither near nor far detector */
339  return -5;
340  });
341  }
342 
343 
345  const std::function<NumuEnergyFunc::prod4_fd_energy_func_t> &f_fd,
346  const std::function<NumuEnergyFunc::prod4_nd_energy_func_t> &f_nd
347  )
348  {
349  return Var([f_fd,f_nd] (const caf::SRProxy *sr) -> double
350  {
351  if (sr->trk.kalman.ntracks == 0) {
352  return -5;
353  }
354 
355  bool isRHC = sr->spill.isRHC;
356  caf::Det_t det = sr->hdr.det;
357  double hadVisE = sr->energy.numu.hadtrkE + sr->energy.numu.hadcalE;
358 
359  if (det == caf::kFARDET)
360  {
361  int run = sr->hdr.run;
362  bool ismc = sr->hdr.ismc;
363  return f_fd(hadVisE, run, ismc, isRHC);
364  }
365 
366  if (det == caf::kNEARDET) {
367  return f_nd(hadVisE, isRHC);
368  }
369 
370  /* neither near nor far detector */
371  return -5;
372  });
373  }
374 
375 
376 
377  float TAMuE(double trklen) {
378  double stitch1 = 143; // cm
379  double stitch2 = 335; // cm
380  double stitch3 = 542; // cm
381  double offset = 0.1642; // GeV
382  double slope1 = 0.001831; // GeV/cm
383  double slope2 = 0.001963; // GeV/cm
384  double slope3 = 0.002025; // GeV/cm
385  double slope4 = 0.002078; // GeV/cm
386 
387  float MuonE = 0.0;
388 
389  if (trklen <= 0.0) return 0.0;
390 
391  if (trklen < stitch1){
392  MuonE = slope1*trklen + offset;
393  }
394  else if (trklen < stitch2){
395  MuonE = slope2*trklen + ((slope1-slope2)*stitch1 + offset);
396  }
397  else if (trklen <stitch3){
398  MuonE = slope3*trklen + ((slope2-slope3)*stitch2 +(slope1-slope2)*stitch1 + offset);
399  }
400  else{
401  MuonE = slope4*trklen + ((slope3-slope4)*stitch3 +(slope2-slope3)*stitch2 +(slope1-slope2)*stitch1 + offset);
402  }
403  return MuonE;
404  }
405 
406 
407  float TAGetNDActLen(double trkLenact, double trkLencat) {
408  double stitch1 = 163; // cm
409  double stitch2 = 535; // cm
410  double stitch3 = 884; // cm
411  double offset = 0.0429; // GeV
412  double slope1 = 0.001714; // GeV/cm
413  double slope2 = 0.002043; // GeV/cm
414  double slope3 = 0.002075; // GeV/cm
415  double slope4 = 0.002141; // GeV/cm
416 
417  double offsetMC = 0.139; // GeV
418  double slope1MC = 0.00540; // GeV/cm
419 
420  double MuonEMC = 0.0;
421  MuonEMC = slope1MC*trkLencat + offsetMC;
422 
423  float muE = 0.0;
424 
425  if (trkLenact <= 0.0) return 0.0;
426 
427  if (trkLenact < stitch1){
428  muE = slope1*trkLenact + offset;
429  }
430  else if (trkLenact < stitch2){
431  muE = slope2*trkLenact + ((slope1-slope2)*stitch1 + offset);
432  }
433  else if (trkLenact <stitch3){
434  muE = slope3*trkLenact + ((slope2-slope3)*stitch2 +(slope1-slope2)*stitch1 + offset);
435  }
436  else{
437  muE = slope4*trkLenact + ((slope3-slope4)*stitch3 +(slope2-slope3)*stitch2 +(slope1-slope2)*stitch1 + offset);
438  }
439  return (muE + MuonEMC);
440  }
441 
442  float TAHadEND(double HadVisE) {
443 
444  float HadE = 0.0;
445  double stitchH1 = 0.045; // GeV
446  double stitchH2 = 0.54; // GeV
447  double stitchH3 = 1.05; // GeV
448  double offsetH = 0.0616; // GeV
449  double slopeH1 = 1.66; // unitless
450  double slopeH2 = 2.182; // unitless
451  double slopeH3 = 1.738; // unitless
452  double slopeH4 = 2.226; // unitless
453 
454  if (HadVisE < 0.0){
455  HadE = offsetH;
456  }
457  else if (HadVisE < stitchH1){
458  HadE = slopeH1*HadVisE + offsetH;
459  }
460  else if (HadVisE < stitchH2){
461  HadE = slopeH2*HadVisE + ((slopeH1-slopeH2)*stitchH1 + offsetH);
462  }
463  else if (HadVisE <stitchH3){
464  HadE = slopeH3*HadVisE + ((slopeH2-slopeH3)*stitchH2 +(slopeH1-slopeH2)*stitchH1 + offsetH);
465  }
466  else{
467  HadE = slopeH4*HadVisE + ((slopeH3-slopeH4)*stitchH3 +(slopeH2-slopeH3)*stitchH2 +(slopeH1-slopeH2)*stitchH1 + offsetH);
468  }
469  return HadE;
470  }
471 
472 
473 
474  const Var kNumuCCXsecMuonE(
475  [](const caf::SRProxy *sr)
476  {
477  if(sr->trk.kalman.ntracks < 1) return 0.0f;
478  if(sr->energy.numu.E <= 0) return 0.0f;
479 
480  caf::Det_t detector = sr->hdr.det;
481  // bool isFD = (detector == caf::kFARDET);
482  bool isND = (detector == caf::kNEARDET);
483  if ((detector == caf::kNDOS) || (detector == caf::kUNKNOWN)) return 0.0f;
484  if (isND == false) return 0.0f;
485 
486  float muonE = 0.0;
487  if (isND){
488  double trkLenAct = sr->energy.numu.ndtrklenact;
489  double trkLenCat = sr->energy.numu.ndtrklencat;
490  if ((trkLenAct <= 0.0)&&(trkLenCat <= 0.0)) return 0.0f;
491 
492  bool actOnly = (((sr->energy.numu.ndtrkcalactE > 0.0)||(sr->energy.numu.ndtrkcaltranE > 0.0))&&(sr->energy.numu.ndtrkcalcatE == 0.0));
493  bool actAndCat = (((sr->energy.numu.ndtrkcalactE > 0.0)||(sr->energy.numu.ndtrkcaltranE > 0.0))&&(sr->energy.numu.ndtrkcalcatE > 0.0));
494 
495  if((actOnly == false)&&(actAndCat == false)) return 0.0f;
496  if((actOnly == true)&&(actAndCat == true)) return 0.0f;
497 
498  double trklen = 0.0f;
499 
500  if (actOnly){
501  trklen = trkLenAct;
502  muonE = TAMuE(trklen);
503  }
504 
505  if (actAndCat){
506  muonE = TAGetNDActLen(trkLenAct,trkLenCat);
507  }//End of actAndCat
508  }//ND energy
509  return muonE;
510  });
511 
512  const Var kNumuCCXsecHadE(
513  [](const caf::SRProxy *sr)
514  {
515  if(sr->trk.kalman.ntracks < 1) return 0.0f;
516  if(sr->energy.numu.E <= 0) return 0.0f;
517 
518  caf::Det_t detector = sr->hdr.det;
519  bool isND = (detector == caf::kNEARDET);
520  if ((detector == caf::kNDOS) || (detector == caf::kUNKNOWN)) return 0.0f;
521  if (isND == false) return 0.0f;
522 
523  float hadE = 0.0;
524  if (isND){
525  double hadVisE = sr->energy.numu.hadtrkE + sr->energy.numu.hadcalE;
526  hadE = TAHadEND(hadVisE);
527  }//ND energy
528  return hadE;
529  });
530 
531 
532 
533  const Var kNumuMuE2017FDp1(
535  );
536  const Var kNumuMuE2017FDp2(
538  );
539  const Var kNumuMuE2017FDp3(
541  );
542  const Var kNumuMuE2017FDp4(
544  );
545  const Var kNumuMuE2017FDp5(
547  );
548 
549  const Var kNumuHadE2017FDp1(
551  );
552  const Var kNumuHadE2017FDp2(
554  );
555  const Var kNumuHadE2017FDp3(
557  );
558  const Var kNumuHadE2017FDp4(
560  );
561  const Var kNumuHadE2017FDp5(
563  );
564 
565 
568 
570 
573 
576 
577 
578  const Var kNumuMuE2017(
579  [] (const caf::SRProxy *sr) -> double
580  {
581  if (sr->trk.kalman.ntracks == 0) {
582  return 0;
583  }
584 
585  bool isRHC = sr->spill.isRHC;
586  caf::Det_t det = sr->hdr.det;
587 
588  if (det == caf::kFARDET)
589  {
590  int run = sr->hdr.run;
591  bool ismc = sr->hdr.ismc;
592  double trkLen = sr->trk.kalman.tracks[0].len;
593 
594  return predict_prod3_fd_muon_energy(trkLen, run, ismc, isRHC);
595  }
596 
597  if (det == caf::kNEARDET)
598  {
599  double trkLenAct = sr->energy.numu.ndtrklenact;
600  double trkLenCat = sr->energy.numu.ndtrklencat;
601 
602  if ( (sr->energy.numu.ndtrkcalactE == 0)
603  && (sr->energy.numu.ndtrkcaltranE == 0))
604  {
605  return 0;
606  }
607 
608  return predict_prod3_nd_act_energy(trkLenAct, isRHC) +
610  }
611 
612  /* neither near nor far detector */
613  return 0;
614  }
615  );
616 
617 
618  const Var kNumuHadE2017(
619  [] (const caf::SRProxy *sr) -> double
620  {
621  if (sr->trk.kalman.ntracks == 0) {
622  return 0;
623  }
624 
625  bool isRHC = sr->spill.isRHC;
626  caf::Det_t det = sr->hdr.det;
627  double hadVisE = sr->energy.numu.hadtrkE + sr->energy.numu.hadcalE;
628 
629  if (det == caf::kFARDET) {
630  int run = sr->hdr.run;
631  bool ismc = sr->hdr.ismc;
632  return predict_prod3_fd_had_energy(hadVisE, run, ismc, isRHC);
633  }
634 
635  if (det == caf::kNEARDET) {
636  return predict_prod3_nd_had_energy(hadVisE, isRHC);
637  }
638 
639  /* neither near nor far detector */
640  return 0;
641  }
642  );
643 
644 
645  // old SA tunings
646 
647  float SAMuE(double trkLen) {
648  //New muon tuning values
649  double stitch1 = 342; // cm
650  double stitch2 = 520; // cm
651  double stitch3 = 1090; // cm
652  double offset = 0.1491; // GeV
653  double slope1 = 0.001951; // GeV/cm
654  double slope2 = 0.002022; // GeV/cm
655  double slope3 = 0.002067; // GeV/cm
656  double slope4 = 0.002196; // GeV/cm
657 
658  float muonE = 0.0;
659 
660  if (trkLen <= 0.0) return 0.0;
661 
662  if (trkLen < stitch1){
663  muonE = slope1*trkLen + offset;
664  }
665  else if (trkLen < stitch2){
666  muonE = slope2*trkLen + ((slope1-slope2)*stitch1 + offset);
667  }
668  else if (trkLen <stitch3){
669  muonE = slope3*trkLen + ((slope2-slope3)*stitch2 +(slope1-slope2)*stitch1 + offset);
670  }
671  else{
672  muonE = slope4*trkLen + ((slope3-slope4)*stitch3 +(slope2-slope3)*stitch2 +(slope1-slope2)*stitch1 + offset);
673  }
674 
675  return muonE;
676  }
677 
678  float SAGetNDCatEffLen(double trkLenCat) {
679 
680  double stitch1 = 342; // cm
681  double stitch2 = 520; // cm
682  double stitch3 = 1090; // cm
683  double offset = 0.1491; // GeV
684  double slope1 = 0.001951; // GeV/cm
685  double slope2 = 0.002022; // GeV/cm
686  double slope3 = 0.002067; // GeV/cm
687  double slope4 = 0.002196; // GeV/cm
688 
689  //New tuning values rounded for errors
690  double offsetMC = 0.139; // GeV
691  double slope1MC = 0.00537; // GeV/cm
692 
693  double muonEMC = slope1MC*trkLenCat + offsetMC;
694 
695  float effActTrkLen = 0.0;
696 
697  if (muonEMC <= 0.0) return 0;
698 
699  if (muonEMC <= offset) {
700  effActTrkLen = 0.0;
701  }
702  else if (muonEMC < (slope1*stitch1 + offset)){
703  effActTrkLen = (muonEMC - offset)/slope1;
704  }
705  else if (muonEMC < (slope2*stitch2 + ((slope1-slope2)*stitch1 + offset))){
706  effActTrkLen = (muonEMC - ((slope1-slope2)*stitch1 + offset))/slope2;
707  }
708  else if (muonEMC < (slope3*stitch3 + ((slope2-slope3)*stitch2 + (slope1-slope2)*stitch1 + offset))){
709  effActTrkLen = (muonEMC - ((slope2-slope3)*stitch2 +(slope1-slope2)*stitch1 + offset))/slope3;
710  }
711  else {
712  effActTrkLen = (muonEMC - ((slope3-slope4)*stitch3 + (slope2-slope3)*stitch2 + (slope1-slope2)*stitch1 + offset))/slope4;
713  }
714 
715  return effActTrkLen;
716  }
717 
718  float SAHadEFD(double hadVisE, int run) {
719 
720  double hadE = 0.0;
721 
722  //New Tune Values with Error Rounding - Period 1
723  double stitchH1 = 0.0799; // GeV
724  double stitchH2 = 0.373; // GeV
725  double stitchH3 = 0.925; // GeV
726  double offsetH = 0.0419; // GeV
727  double slopeH1 = 0.816; // unitless
728  double slopeH2 = 2.04; // unitless
729  double slopeH3 = 1.68; // unitless
730  double slopeH4 = 2.45; // unitless
731 
732  if ((run >= 17140)&&(run < 20753)){
733  //New Tune Values with Error Rounding - Period 2
734  stitchH1 = 0.0694; // GeV
735  stitchH2 = 0.288; // GeV
736  stitchH3 = 0.930; // GeV
737  offsetH = 0.0392; // GeV
738  slopeH1 = 1.45; // unitless
739  slopeH2 = 2.13; // unitless
740  slopeH3 = 1.76; // unitless
741  slopeH4 = 2.17; // unitless
742  }
743 
744  if (run >= 20753){
745  //New Tune Values with Error Rounding - Period 3
746  stitchH1 = 0.0702; // GeV
747  stitchH2 = 0.321; // GeV
748  stitchH3 = 1.06; // GeV
749  offsetH = 0.0365; // GeV
750  slopeH1 = 1.34; // unitless
751  slopeH2 = 2.05; // unitless
752  slopeH3 = 1.68; // unitless
753  slopeH4 = 2.28; // unitless
754  }
755  if (hadVisE < 0.0){
756  hadE = offsetH;
757  }
758  else if (hadVisE < stitchH1){
759  hadE = slopeH1*hadVisE + offsetH;
760  }
761  else if (hadVisE < stitchH2){
762  hadE = slopeH2*hadVisE + ((slopeH1-slopeH2)*stitchH1 + offsetH);
763  }
764  else if (hadVisE <stitchH3){
765  hadE = slopeH3*hadVisE + ((slopeH2-slopeH3)*stitchH2 +(slopeH1-slopeH2)*stitchH1 + offsetH);
766  }
767  else{
768  hadE = slopeH4*hadVisE + ((slopeH3-slopeH4)*stitchH3 +(slopeH2-slopeH3)*stitchH2 +(slopeH1-slopeH2)*stitchH1 + offsetH);
769  }
770  return hadE;
771  }
772 
773 
774  float SAHadEND(double hadVisE) {
775 
776  float hadE = 0.0;
777 
778  //New Tune Values with Error Rounding
779  double stitchH1 = 0.0394; // GeV
780  double stitchH2 = 0.362; // GeV
781  double stitchH3 = 0.955; // GeV
782  double offsetH = 0.0534; // GeV
783  double slopeH1 = 0.956; // unitless
784  double slopeH2 = 2.14; // unitless
785  double slopeH3 = 1.78; // unitless
786  double slopeH4 = 2.00; // unitless
787 
788  if (hadVisE < 0.0){
789  hadE = offsetH;
790  }
791  else if (hadVisE < stitchH1){
792  hadE = slopeH1*hadVisE + offsetH;
793  }
794  else if (hadVisE < stitchH2){
795  hadE = slopeH2*hadVisE + ((slopeH1-slopeH2)*stitchH1 + offsetH);
796  }
797  else if (hadVisE <stitchH3){
798  hadE = slopeH3*hadVisE + ((slopeH2-slopeH3)*stitchH2 +(slopeH1-slopeH2)*stitchH1 + offsetH);
799  }
800  else{
801  hadE = slopeH4*hadVisE + ((slopeH3-slopeH4)*stitchH3 +(slopeH2-slopeH3)*stitchH2 +(slopeH1-slopeH2)*stitchH1 + offsetH);
802  }
803  return hadE;
804  }
805 
806  double SAMuEND(double lenact, double lencat,
807  double actE, double tranE,
808  double catE)
809  {
810 
811  if ((lenact <= 0.0)&&(lencat <= 0.0)) return 0.0;
812 
813  bool actOnly = (((actE > 0.0)||(tranE > 0.0))&&(catE == 0.0));
814  bool actAndCat = (((actE > 0.0)||(tranE > 0.0))&&(catE > 0.0));
815 
816  if((actOnly == false)&&(actAndCat == false)) return 0.0;
817  if((actOnly == true)&&(actAndCat == true)) return 0.0;
818 
819  double trkLen = 0.0;
820 
821  if (actOnly){
822  trkLen = lenact;
823  }
824 
825  if (actAndCat){
826  trkLen = lenact + SAGetNDCatEffLen(lencat);
827  }//End of actAndCat
828 
829  return SAMuE(trkLen);
830 
831  }
832 }
double predict_prod3_nd_p4_cat_energy(double ndtrklencat)
caf::Proxy< caf::SRSpill > spill
Definition: SRProxy.h:2143
#define NumuHadE2017NDW(f)
Definition: NumuEFxs.cxx:31
Near Detector underground.
Definition: SREnums.h:10
Unknown detector.
Definition: SREnums.h:9
Det_t
Which NOvA detector?
Definition: SREnums.h:7
const Var kNumuHadE2017RHCND(NumuHadE2017NDW( predict_prod3_nd_p4_had_energy))
Definition: NumuEFxs.h:58
Far Detector at Ash River.
Definition: SREnums.h:11
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
double predict_prod3_fd_p4_muon_energy_corr(double trkLen, bool ismc)
double predict_prod3_fd_p3_muon_energy_corr(double trkLen, bool ismc)
double predict_prod3_fd_p5_had_energy_corr(double hadVisE, bool ismc)
double predict_prod3_fd_p2_muon_energy_corr(double trkLen, bool ismc)
caf::Proxy< float > E
Definition: SRProxy.h:164
Var initNumuMuE2020NDpXVar(const std::function< NumuEnergyFunc::prod5_nd_energy_p_func_t > &f_act, const std::function< NumuEnergyFunc::prod5_nd_energy_p_func_t > &f_cat)
Definition: NumuEFxs.cxx:103
caf::Proxy< size_t > ntracks
Definition: SRProxy.h:1778
double predict_prod3_fd_p4_had_energy_corr(double hadVisE, bool ismc)
caf::Proxy< caf::SRHeader > hdr
Definition: SRProxy.h:2137
caf::Proxy< caf::SRNumuEnergy > numu
Definition: SRProxy.h:214
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
const Var kNumuMuE2017([](const caf::SRProxy *sr) -> double{if(sr->trk.kalman.ntracks==0){return 0;}bool isRHC=sr->spill.isRHC;caf::Det_t det=sr->hdr.det;if(det==caf::kFARDET){int run=sr->hdr.run;bool ismc=sr->hdr.ismc;double trkLen=sr->trk.kalman.tracks[0].len;return predict_prod3_fd_muon_energy(trkLen, run, ismc, isRHC);}if(det==caf::kNEARDET){double trkLenAct=sr->energy.numu.ndtrklenact;double trkLenCat=sr->energy.numu.ndtrklencat;if( (sr->energy.numu.ndtrkcalactE==0) &&(sr->energy.numu.ndtrkcaltranE==0)){return 0;}return predict_prod3_nd_act_energy(trkLenAct, isRHC)+predict_prod3_nd_cat_energy(trkLenCat, isRHC);}return 0;})
Definition: NumuEFxs.h:63
double predict_prod3_fd_p5_muon_energy_corr(double trkLen, bool ismc)
#define NumuMuE2017NDW(f_act, f_cat)
Definition: NumuEFxs.cxx:43
const Var kNumuCCXsecHadE([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return 0.0f;if(sr->energy.numu.E<=0) return 0.0f;caf::Det_t detector=sr->hdr.det;bool isND=(detector==caf::kNEARDET);if((detector==caf::kNDOS)||(detector==caf::kUNKNOWN)) return 0.0f;if(isND==false) return 0.0f;float hadE=0.0;if(isND){double hadVisE=sr->energy.numu.hadtrkE+sr->energy.numu.hadcalE;hadE=TAHadEND(hadVisE);}return hadE;})
Definition: NumuEFxs.h:29
caf::Proxy< float > hadtrkE
Definition: SRProxy.h:171
Var initNumuHadE2020NDpXVar(const std::function< NumuEnergyFunc::prod5_nd_energy_p_func_t > &f)
Definition: NumuEFxs.cxx:127
double predict_prod3_fd_muon_energy(double trkLen, int run, bool ismc, bool isRHC)
double predict_prod3_nd_p4_had_energy(double hadVisE)
double predict_prod3_fd_p2_had_energy_corr(double hadVisE, bool ismc)
double predict_prod3_nd_p3_act_energy(double ndtrklenact)
caf::Proxy< caf::SREnergyBranch > energy
Definition: SRProxy.h:2136
Var initNumuMuE2020FDpXVar(const std::function< NumuEnergyFunc::prod5_fd_energy_p_func_t > &f)
Definition: NumuEFxs.cxx:67
caf::Proxy< unsigned int > run
Definition: SRProxy.h:248
_Var< caf::SRProxy > Var
Representation of a variable to be retrieved from a caf::StandardRecord object.
Definition: Var.h:74
float SAGetNDCatEffLen(double trkLenCat)
Definition: NumuEFxs.cxx:678
#define NumuMuE2017FDW(f)
Definition: NumuEFxs.cxx:5
const Var kNumuHadE2017([](const caf::SRProxy *sr) -> double{if(sr->trk.kalman.ntracks==0){return 0;}bool isRHC=sr->spill.isRHC;caf::Det_t det=sr->hdr.det;double hadVisE=sr->energy.numu.hadtrkE+sr->energy.numu.hadcalE;if(det==caf::kFARDET){int run=sr->hdr.run;bool ismc=sr->hdr.ismc;return predict_prod3_fd_had_energy(hadVisE, run, ismc, isRHC);}if(det==caf::kNEARDET){return predict_prod3_nd_had_energy(hadVisE, isRHC);}return 0;})
Definition: NumuEFxs.h:64
const Var kNumuMuE2017FDp1(NumuMuE2017FDW(predict_prod3_fd_p1_muon_energy_corr))
Definition: NumuEFxs.h:34
float SAHadEND(double hadVisE)
Definition: NumuEFxs.cxx:774
Var initNumuHadE2018NDpXVar(const std::function< NumuEnergyFunc::prod4_nd_energy_p_func_t > &f)
Definition: NumuEFxs.cxx:283
Var initNumuHadE2018Var(const std::function< NumuEnergyFunc::prod4_fd_energy_func_t > &f_fd, const std::function< NumuEnergyFunc::prod4_nd_energy_func_t > &f_nd)
Definition: NumuEFxs.cxx:344
double predict_prod3_nd_p3_had_energy(double hadVisE)
Var initNumuMuE2018FDpXVar(const std::function< NumuEnergyFunc::prod4_fd_energy_p_func_t > &f)
Definition: NumuEFxs.cxx:225
const Var kNumuCCXsecMuonE([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return 0.0f;if(sr->energy.numu.E<=0) return 0.0f;caf::Det_t detector=sr->hdr.det;bool isND=(detector==caf::kNEARDET);if((detector==caf::kNDOS)||(detector==caf::kUNKNOWN)) return 0.0f;if(isND==false) return 0.0f;float muonE=0.0;if(isND){double trkLenAct=sr->energy.numu.ndtrklenact;double trkLenCat=sr->energy.numu.ndtrklencat;if((trkLenAct<=0.0)&&(trkLenCat<=0.0)) return 0.0f;bool actOnly=(((sr->energy.numu.ndtrkcalactE > 0.0)||(sr->energy.numu.ndtrkcaltranE > 0.0))&&(sr->energy.numu.ndtrkcalcatE==0.0));bool actAndCat=(((sr->energy.numu.ndtrkcalactE > 0.0)||(sr->energy.numu.ndtrkcaltranE > 0.0))&&(sr->energy.numu.ndtrkcalcatE > 0.0));if((actOnly==false)&&(actAndCat==false)) return 0.0f;if((actOnly==true)&&(actAndCat==true)) return 0.0f;double trklen=0.0f;if(actOnly){trklen=trkLenAct;muonE=TAMuE(trklen);}if(actAndCat){muonE=TAGetNDActLen(trkLenAct, trkLenCat);}}return muonE;})
Definition: NumuEFxs.h:28
if(dump)
double SAMuEND(double lenact, double lencat, double actE, double tranE, double catE)
Definition: NumuEFxs.cxx:806
const Var kNumuMuE2017RHCND(NumuMuE2017NDW(predict_prod3_nd_p4_act_energy, predict_prod3_nd_p4_cat_energy))
Definition: NumuEFxs.h:57
#define NumuHadE2017FDW(f)
Definition: NumuEFxs.cxx:18
caf::Proxy< caf::SRTrackBranch > trk
Definition: SRProxy.h:2145
Prototype Near Detector on the Surface.
Definition: SREnums.h:12
caf::Proxy< float > ndtrkcalcatE
Definition: SRProxy.h:184
const Var kNumuMuE2017FDp2(NumuMuE2017FDW(predict_prod3_fd_p2_muon_energy_corr))
Definition: NumuEFxs.h:35
Var initNumuMuE2018NDpXVar(const std::function< NumuEnergyFunc::prod4_nd_energy_p_func_t > &f_act, const std::function< NumuEnergyFunc::prod4_nd_energy_p_func_t > &f_cat)
Definition: NumuEFxs.cxx:259
caf::Proxy< float > ndtrkcaltranE
Definition: SRProxy.h:185
Var initNumuMuE2020Var(const std::function< NumuEnergyFunc::prod5_fd_energy_func_t > &f_fd, const std::function< NumuEnergyFunc::prod5_nd_energy_func_t > &f_nd_act, const std::function< NumuEnergyFunc::prod5_nd_energy_func_t > &f_nd_cat)
Definition: NumuEFxs.cxx:145
const Var kNumuHadE2017FDp3(NumuHadE2017FDW(predict_prod3_fd_p3_had_energy_corr))
Definition: NumuEFxs.h:42
double predict_prod3_fd_p3_had_energy_corr(double hadVisE, bool ismc)
caf::StandardRecord * sr
caf::Proxy< float > ndtrkcalactE
Definition: SRProxy.h:183
const Var kNumuMuE2017FDp4(NumuMuE2017FDW(predict_prod3_fd_p4_muon_energy_corr))
Definition: NumuEFxs.h:37
caf::Proxy< float > hadcalE
Definition: SRProxy.h:169
const Var kNumuHadE2017FDp5(NumuHadE2017FDW(predict_prod3_fd_p5_had_energy_corr))
Definition: NumuEFxs.h:44
float TAMuE(double trklen)
Definition: NumuEFxs.cxx:377
double predict_prod3_nd_p4_act_energy(double ndtrklenact)
const Var kNumuHadE2017FDp1(NumuHadE2017FDW(predict_prod3_fd_p1_had_energy_corr))
Definition: NumuEFxs.h:40
Var initNumuHadE2020Var(const std::function< NumuEnergyFunc::prod5_fd_energy_func_t > &f_fd, const std::function< NumuEnergyFunc::prod5_nd_energy_func_t > &f_nd)
Definition: NumuEFxs.cxx:188
Definition: run.py:1
double predict_prod3_fd_p1_had_energy_corr(double hadVisE, bool ismc)
float SAMuE(double trkLen)
Definition: NumuEFxs.cxx:647
const Var kNumuMuE2017FDp3(NumuMuE2017FDW(predict_prod3_fd_p3_muon_energy_corr))
Definition: NumuEFxs.h:36
Var initNumuHadE2020FDpXVar(const std::function< NumuEnergyFunc::prod5_fd_energy_p_func_t > &f)
Definition: NumuEFxs.cxx:85
caf::Proxy< caf::SRKalman > kalman
Definition: SRProxy.h:1797
double predict_prod3_nd_p3_cat_energy(double ndtrklencat)
float TAGetNDActLen(double trkLenact, double trkLencat)
Definition: NumuEFxs.cxx:407
double predict_prod3_nd_cat_energy(double ndtrklencat, bool isRHC)
double predict_prod3_nd_had_energy(double hadVisE, bool isRHC)
float SAHadEFD(double hadVisE, int run)
Definition: NumuEFxs.cxx:718
caf::Proxy< float > ndtrklenact
Definition: SRProxy.h:186
Var initNumuMuE2018Var(const std::function< NumuEnergyFunc::prod4_fd_energy_func_t > &f_fd, const std::function< NumuEnergyFunc::prod4_nd_energy_func_t > &f_nd_act, const std::function< NumuEnergyFunc::prod4_nd_energy_func_t > &f_nd_cat)
Definition: NumuEFxs.cxx:300
Var initNumuHadE2018FDpXVar(const std::function< NumuEnergyFunc::prod4_fd_energy_p_func_t > &f)
Definition: NumuEFxs.cxx:242
caf::Proxy< bool > ismc
Definition: SRProxy.h:242
const Var kNumuHadE2017FDp4(NumuHadE2017FDW(predict_prod3_fd_p4_had_energy_corr))
Definition: NumuEFxs.h:43
const Var kNumuHadE2017ND(NumuHadE2017NDW(predict_prod3_nd_p3_had_energy))
Definition: NumuEFxs.h:53
float TAHadEND(double HadVisE)
Definition: NumuEFxs.cxx:442
caf::Proxy< bool > isRHC
Definition: SRProxy.h:1376
caf::Proxy< std::vector< caf::SRKalmanTrack > > tracks
Definition: SRProxy.h:1780
const Var kNumuMuE2017FDp5(NumuMuE2017FDW(predict_prod3_fd_p5_muon_energy_corr))
Definition: NumuEFxs.h:38
double predict_prod3_fd_p1_muon_energy_corr(double trkLen, bool ismc)
caf::Proxy< float > ndtrklencat
Definition: SRProxy.h:187
const Var kNumuHadE2017FDp2(NumuHadE2017FDW(predict_prod3_fd_p2_had_energy_corr))
Definition: NumuEFxs.h:41
const Var kNumuMuE2017ND(NumuMuE2017NDW(predict_prod3_nd_p3_act_energy, predict_prod3_nd_p3_cat_energy))
Definition: NumuEFxs.h:52
caf::Proxy< caf::Det_t > det
Definition: SRProxy.h:232
double predict_prod3_fd_had_energy(double hadVisE, int run, bool ismc, bool isRHC)
double predict_prod3_nd_act_energy(double ndtrklenact, bool isRHC)