NumuEFSpecialxs.cxx
Go to the documentation of this file.
2 
3 /*
4  * TODO: properly refactor code to have a header with required declarations
5  */
7 
9 
10 #define NumuMuE2017FDW(f) \
11  [] (const caf::SRProxy *sr) -> double \
12  { \
13  if (sr->trk.kalman.ntracks == 0) { \
14  return 0; \
15  } \
16  \
17  double trkLen = sr->trk.kalman.tracks[0].len; \
18  bool ismc = sr->hdr.ismc; \
19  \
20  return f(trkLen, ismc); \
21  }
22 
23 #define NumuHadE2017FDW(f) \
24  [] (const caf::SRProxy *sr) -> double \
25  { \
26  if (sr->trk.kalman.ntracks == 0) { \
27  return 0; \
28  } \
29  \
30  bool ismc = sr->hdr.ismc; \
31  double hadVisE = sr->energy.numu.hadtrkE + sr->energy.numu.hadcalE; \
32  \
33  return f(hadVisE, ismc); \
34  }
35 
36 #define NumuHadE2017NDW(f) \
37  [] (const caf::SRProxy *sr) -> double \
38  { \
39  if (sr->trk.kalman.ntracks == 0) { \
40  return 0; \
41  } \
42  \
43  double hadVisE = sr->energy.numu.hadtrkE + sr->energy.numu.hadcalE; \
44  \
45  return f(hadVisE); \
46  }
47 
48 #define NumuMuE2017NDW(f_act, f_cat) \
49  [] (const caf::SRProxy *sr) -> double \
50  { \
51  if (sr->trk.kalman.ntracks == 0) { \
52  return 0; \
53  } \
54  \
55  double trkLenAct = sr->energy.numu.ndtrklenact; \
56  double trkLenCat = sr->energy.numu.ndtrklencat; \
57  \
58  if ( (sr->energy.numu.ndtrkcalactE == 0) \
59  && (sr->energy.numu.ndtrkcaltranE == 0)) \
60  { \
61  return 0; \
62  } \
63  \
64  return f_act(trkLenAct) + f_cat(trkLenCat); \
65  }
66 
67 namespace ana {
68 
69  /*
70  * special energy estimators
71  */
72 
73  double predict_special_fd_p1_muon_energy(double trkLen)
74  {
75 
76  const double offset = 1.613851970352435661e-01;
77  const double slope0 = 1.884065471495919020e-01;
78  const double stitch1 = 2.466666646099010762e+00;
79  const double slope1 = 1.989826976703109596e-01;
80  const double stitch2 = 4.422209871371248546e+00;
81  const double slope2 = 2.056268485963292092e-01;
82  const double stitch3 = 9.933336098068100029e+00;
83  const double slope3 = 2.187431512281476642e-01;
84 
85  if (trkLen <= 0) {
86  return 0;
87  }
88 
89  return spline3(trkLen / 100,
90  offset, slope0,
91  stitch1, slope1,
92  stitch2, slope2,
93  stitch3, slope3);
94  }
95 
96 
97 
98  double predict_special_fd_p2_muon_energy(double trkLen)
99  {
100 
101  const double offset = 1.640715025504196900e-01;
102  const double slope0 = 1.884091893782043614e-01;
103  const double stitch1 = 2.376534942967446717e+00;
104  const double slope1 = 1.986844669481659387e-01;
105  const double stitch2 = 4.704987401094306243e+00;
106  const double slope2 = 2.065885016956024722e-01;
107  const double stitch3 = 1.000449201373624142e+01;
108  const double slope3 = 2.181253627897896397e-01;
109 
110  if (trkLen <= 0) {
111  return 0;
112  }
113 
114  return spline3(trkLen / 100,
115  offset, slope0,
116  stitch1, slope1,
117  stitch2, slope2,
118  stitch3, slope3);
119  }
120 
121 
122 
123  double predict_special_fd_p3_muon_energy(double trkLen)
124  {
125 
126  const double offset = 1.467148337502355293e-01;
127  const double slope0 = 1.960474766711914618e-01;
128  const double stitch1 = 3.996862510511446587e+00;
129  const double slope1 = 2.053499719645651977e-01;
130  const double stitch2 = 9.703458476274496647e+00;
131  const double slope2 = 2.180968259059623271e-01;
132 
133  if (trkLen <= 0) {
134  return 0;
135  }
136 
137  return spline2(trkLen / 100,
138  offset, slope0,
139  stitch1, slope1,
140  stitch2, slope2);
141  }
142 
143 
144  double predict_special_fd_p5_muon_energy(double trkLen)
145  {
146 
147  const double offset = 1.472489248766962744e-01;
148  const double slope0 = 1.959403715530378398e-01;
149  const double stitch1 = 4.131691877836759730e+00;
150  const double slope1 = 2.058580582558999228e-01;
151  const double stitch2 = 9.609108710793506702e+00;
152  const double slope2 = 2.169615956961473968e-01;
153 
154  if (trkLen <= 0) {
155  return 0;
156  }
157 
158  return spline2(trkLen / 100,
159  offset, slope0,
160  stitch1, slope1,
161  stitch2, slope2);
162  }
163 
164 
165 
166  double predict_special_fd_p1_had_energy(double hadVisE)
167  {
168 
169  const double offset = -1.115503071653284328e-02;
170  const double slope0 = 1.981675201010266285e+00;
171  const double stitch1 = 6.250000000537712097e-01;
172  const double slope1 = 1.608178619859576219e+00;
173  const double stitch2 = 9.019538685837620307e-01;
174  const double slope2 = 2.247024194052140711e+00;
175 
176  if (hadVisE <= 0) {
177  return 0;
178  }
179 
180  return spline2(hadVisE,
181  offset, slope0,
182  stitch1, slope1,
183  stitch2, slope2);
184  }
185 
186 
187 
188  double predict_special_fd_p2_had_energy(double hadVisE)
189  {
190 
191  const double offset = 4.099301758354712000e-02;
192  const double slope0 = 1.938155325867301659e+00;
193  const double stitch1 = 2.801477052019014091e-01;
194  const double slope1 = 2.236099595431565668e+00;
195  const double stitch2 = 4.385370235138835171e-01;
196  const double slope2 = 1.750514963254879985e+00;
197  const double stitch3 = 8.316561158197010029e-01;
198  const double slope3 = 2.063851716405838310e+00;
199 
200  if (hadVisE <= 0) {
201  return 0;
202  }
203 
204  return spline3(hadVisE,
205  offset, slope0,
206  stitch1, slope1,
207  stitch2, slope2,
208  stitch3, slope3);
209  }
210 
211 
212 
213  double predict_special_fd_p3_had_energy(double hadVisE)
214  {
215 
216  const double offset = 3.745449190173599785e-02;
217  const double slope0 = 1.892075873818970910e+00;
218  const double stitch1 = 2.974351615574292729e-01;
219  const double slope1 = 2.156697626694576098e+00;
220  const double stitch2 = 4.574523715735526186e-01;
221  const double slope2 = 1.699958562404741302e+00;
222  const double stitch3 = 8.487244850111413941e-01;
223  const double slope3 = 1.973030303034408517e+00;
224 
225  if (hadVisE <= 0) {
226  return 0;
227  }
228 
229  return spline3(hadVisE,
230  offset, slope0,
231  stitch1, slope1,
232  stitch2, slope2,
233  stitch3, slope3);
234  }
235 
236 
237 
238  double predict_special_fd_p5_had_energy(double hadVisE)
239  {
240 
241  const double offset = 3.684482597641203228e-02;
242  const double slope0 = 1.890974360981876323e+00;
243  const double stitch1 = 2.772971163373517678e-01;
244  const double slope1 = 2.145786130402019154e+00;
245  const double stitch2 = 4.600878655820809238e-01;
246  const double slope2 = 1.689944140975524700e+00;
247  const double stitch3 = 8.545467147032163036e-01;
248  const double slope3 = 2.026605695978441091e+00;
249 
250  if (hadVisE <= 0) {
251  return 0;
252  }
253 
254  return spline3(hadVisE,
255  offset, slope0,
256  stitch1, slope1,
257  stitch2, slope2,
258  stitch3, slope3);
259  }
260 
261 
262 
263  double predict_special_nd_p3_act_energy(double ndtrklenact)
264  {
265 
266  const double offset = 1.563378889497577529e-01;
267  const double slope0 = 1.912184900372692065e-01;
268  const double stitch1 = 2.620527402376078285e+00;
269  const double slope1 = 1.999735087857159310e-01;
270  const double stitch2 = 5.027393643667430467e+00;
271  const double slope2 = 2.078465326350594222e-01;
272 
273  if (ndtrklenact <= 0) {
274  return 0;
275  }
276 
277  return spline2(ndtrklenact / 100,
278  offset, slope0,
279  stitch1, slope1,
280  stitch2, slope2);
281  }
282 
283 
284 
285  double predict_special_nd_p3_cat_energy(double ndtrklencat)
286  {
287 
288  const double offset = 7.184563085776096703e-02;
289  const double slope0 = 7.125191883683058836e-02;
290  const double stitch1 = 2.141333257577366922e-01;
291  const double slope1 = 7.735382807819963791e-01;
292  const double stitch2 = 4.158036082639919861e-01;
293  const double slope2 = 1.441826598112122548e-01;
294  const double stitch3 = 4.466999142739985773e-01;
295  const double slope3 = 5.585469746098317145e-01;
296 
297  if (ndtrklencat <= 0) {
298  return 0;
299  }
300 
301  return spline3(ndtrklencat / 100,
302  offset, slope0,
303  stitch1, slope1,
304  stitch2, slope2,
305  stitch3, slope3);
306  }
307 
308 
309 
310  double predict_special_nd_p3_had_energy(double hadVisE)
311  {
312 
313  const double offset = 4.270943510939317900e-02;
314  const double slope0 = 2.053687935268237119e+00;
315  const double stitch1 = 5.551117447250514259e-01;
316  const double slope1 = 1.795794073953269843e+00;
317 
318  if (hadVisE <= 0) {
319  return 0;
320  }
321 
322  return spline1(hadVisE,
323  offset, slope0,
324  stitch1, slope1);
325  }
326 
327 
328 
329  double predict_special_fd_p1_muon_energy_corr(double trkLen, bool ismc)
330  {
331  trkLen = correct_prod3_fd_muon_length(trkLen, ismc);
332  return predict_special_fd_p1_muon_energy(trkLen);
333  }
334 
335 
336 
337  double predict_special_fd_p2_muon_energy_corr(double trkLen, bool ismc)
338  {
339  trkLen = correct_prod3_fd_muon_length(trkLen, ismc);
340  return predict_special_fd_p2_muon_energy(trkLen);
341  }
342 
343 
344 
345  double predict_special_fd_p3_muon_energy_corr(double trkLen, bool ismc)
346  {
347  trkLen = correct_prod3_fd_muon_length(trkLen, ismc);
348  return predict_special_fd_p3_muon_energy(trkLen);
349  }
350 
351 
352  double predict_special_fd_p5_muon_energy_corr(double trkLen, bool ismc)
353  {
354  trkLen = correct_prod3_fd_muon_length(trkLen, ismc);
355  return predict_special_fd_p5_muon_energy(trkLen);
356  }
357 
358 
359 
360  double predict_special_fd_p1_had_energy_corr(double hadVisE, bool ismc)
361  {
362  hadVisE = correct_prod3_fd_calE_low_gain(hadVisE, ismc);
363  return predict_special_fd_p1_had_energy(hadVisE);
364  }
365 
366 
367 
368  double predict_special_fd_p2_had_energy_corr(double hadVisE, bool ismc)
369  {
370  hadVisE = correct_prod3_fd_calE_low_gain(hadVisE, ismc);
371  return predict_special_fd_p2_had_energy(hadVisE);
372  }
373 
374 
375 
376  double predict_special_fd_p3_had_energy_corr(double hadVisE, bool ismc)
377  {
378  hadVisE = correct_prod3_fd_calE_high_gain(hadVisE, ismc);
379  return predict_special_fd_p3_had_energy(hadVisE);
380  }
381 
382 
383 
384  double predict_special_fd_p5_had_energy_corr(double hadVisE, bool ismc)
385  {
386  hadVisE = correct_prod3_fd_calE_high_gain(hadVisE, ismc);
387  return predict_special_fd_p5_had_energy(hadVisE);
388  }
389 
390 
391 
392  double predict_special_fd_muon_energy(double trkLen, int run, bool ismc)
393  {
394  trkLen = correct_prod3_fd_muon_length(trkLen, ismc);
395 
396  if (run <= fd_period_1_run_end) {
397  return predict_special_fd_p1_muon_energy(trkLen);
398  }
399 
400  if (run <= fd_period_2_run_end) {
401  return predict_special_fd_p2_muon_energy(trkLen);
402  }
403 
404  if (run <= fd_period_3_run_end) {
405  return predict_special_fd_p3_muon_energy(trkLen);
406  }
407 
408  return predict_special_fd_p5_muon_energy(trkLen);
409  }
410 
411 
412 
413  double predict_special_fd_had_energy(double hadVisE, int run, bool ismc)
414  {
415  hadVisE = correct_prod3_fd_calE(hadVisE, run, ismc);
416 
417  if (run <= fd_period_1_run_end) {
418  return predict_special_fd_p1_had_energy(hadVisE);
419  }
420 
421  if (run <= fd_period_2_run_end) {
422  return predict_special_fd_p2_had_energy(hadVisE);
423  }
424 
425  if (run <= fd_period_3_run_end) {
426  return predict_special_fd_p3_had_energy(hadVisE);
427  }
428 
429  return predict_special_fd_p5_had_energy(hadVisE);
430  }
431 
432 
433  double predict_special_nd_act_energy(double ndtrklenact)
434  {
435  return predict_special_nd_p3_act_energy(ndtrklenact);
436  }
437 
438 
439  double predict_special_nd_cat_energy(double ndtrklencat)
440  {
441  return predict_special_nd_p3_cat_energy(ndtrklencat);
442  }
443 
444 
445  double predict_special_nd_had_energy(double hadVisE)
446  {
447  return predict_special_nd_p3_had_energy(hadVisE);
448  }
449 
452  );
455  );
458  );
461  );
462 
465  );
468  );
471  );
474  );
475 
479  );
480 
483  );
484 
485  const Var kNumuMuE2017Special(
486  [] (const caf::SRProxy *sr) -> double
487  {
488  if (sr->trk.kalman.ntracks == 0) {
489  return 0;
490  }
491 
492  caf::Det_t det = sr->hdr.det;
493 
494  if (det == caf::kFARDET)
495  {
496  int run = sr->hdr.run;
497  bool ismc = sr->hdr.ismc;
498  double trkLen = sr->trk.kalman.tracks[0].len;
499 
500  return predict_special_fd_muon_energy(trkLen, run, ismc);
501  }
502 
503  if (det == caf::kNEARDET)
504  {
505  double trkLenAct = sr->energy.numu.ndtrklenact;
506  double trkLenCat = sr->energy.numu.ndtrklencat;
507 
508  if ( (sr->energy.numu.ndtrkcalactE == 0)
509  && (sr->energy.numu.ndtrkcaltranE == 0))
510  {
511  return 0;
512  }
513 
514  return predict_special_nd_act_energy(trkLenAct) +
516  }
517 
518  /* neither near nor far detector */
519  return 0;
520  }
521  );
522 
523 
525  [] (const caf::SRProxy *sr) -> double
526  {
527  if (sr->trk.kalman.ntracks == 0) {
528  return 0;
529  }
530 
531  caf::Det_t det = sr->hdr.det;
532  double hadVisE = sr->energy.numu.hadtrkE + sr->energy.numu.hadcalE;
533 
534  if (det == caf::kFARDET) {
535  int run = sr->hdr.run;
536  bool ismc = sr->hdr.ismc;
537  return predict_special_fd_had_energy(hadVisE, run, ismc);
538  }
539 
540  if (det == caf::kNEARDET) {
541  return predict_special_nd_had_energy(hadVisE);
542  }
543 
544  /* neither near nor far detector */
545  return 0;
546  }
547  );
548 
549 }
550 
Near Detector underground.
Definition: SREnums.h:10
double predict_special_fd_p3_had_energy_corr(double hadVisE, bool ismc)
Det_t
Which NOvA detector?
Definition: SREnums.h:7
double predict_special_fd_p2_had_energy(double hadVisE)
double correct_prod3_fd_muon_length(double trkLen, bool ismc)
const Var kNumuHadE2017SpecialFDp3(NumuHadE2017FDW(predict_special_fd_p3_had_energy_corr))
Far Detector at Ash River.
Definition: SREnums.h:11
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
#define NumuMuE2017FDW(f)
const Var kNumuMuE2017SpecialFDp1(NumuMuE2017FDW(predict_special_fd_p1_muon_energy_corr))
caf::Proxy< size_t > ntracks
Definition: SRProxy.h:1778
const Var kNumuMuE2017Special([](const caf::SRProxy *sr) -> double{if(sr->trk.kalman.ntracks==0){return 0;}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_special_fd_muon_energy(trkLen, run, ismc);}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_special_nd_act_energy(trkLenAct)+predict_special_nd_cat_energy(trkLenCat);}return 0;})
const Var kNumuHadE2017SpecialFDp5(NumuHadE2017FDW(predict_special_fd_p5_had_energy_corr))
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 kNumuHadE2017Special([](const caf::SRProxy *sr) -> double{if(sr->trk.kalman.ntracks==0){return 0;}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_special_fd_had_energy(hadVisE, run, ismc);}if(det==caf::kNEARDET){return predict_special_nd_had_energy(hadVisE);}return 0;})
double predict_special_fd_p3_muon_energy(double trkLen)
double predict_special_fd_p1_had_energy(double hadVisE)
caf::Proxy< float > hadtrkE
Definition: SRProxy.h:171
double predict_special_fd_p1_muon_energy_corr(double trkLen, bool ismc)
double predict_special_nd_p3_act_energy(double ndtrklenact)
double predict_special_nd_had_energy(double hadVisE)
#define NumuHadE2017NDW(f)
const Var kNumuMuE2017SpecialFDp5(NumuMuE2017FDW(predict_special_fd_p5_muon_energy_corr))
double predict_special_nd_p3_had_energy(double hadVisE)
caf::Proxy< caf::SREnergyBranch > energy
Definition: SRProxy.h:2136
double predict_special_fd_p5_muon_energy_corr(double trkLen, bool ismc)
caf::Proxy< unsigned int > run
Definition: SRProxy.h:248
const Var kNumuHadE2017SpecialFDp2(NumuHadE2017FDW(predict_special_fd_p2_had_energy_corr))
const int fd_period_2_run_end
Definition: Numu2017Fits.cxx:5
static double spline1(double x, double offset, double slope0, double stitch1, double slope1)
const int fd_period_1_run_end
Definition: Numu2017Fits.cxx:4
const int fd_period_3_run_end
Definition: Numu2017Fits.cxx:6
if(dump)
double predict_special_nd_act_energy(double ndtrklenact)
caf::Proxy< caf::SRTrackBranch > trk
Definition: SRProxy.h:2145
double predict_special_fd_p1_had_energy_corr(double hadVisE, bool ismc)
double predict_special_nd_p3_cat_energy(double ndtrklencat)
double predict_special_fd_had_energy(double hadVisE, int run, bool ismc)
caf::Proxy< float > ndtrkcaltranE
Definition: SRProxy.h:185
caf::StandardRecord * sr
caf::Proxy< float > ndtrkcalactE
Definition: SRProxy.h:183
double predict_special_fd_p5_had_energy(double hadVisE)
caf::Proxy< float > hadcalE
Definition: SRProxy.h:169
double correct_prod3_fd_calE(double calE, int run, bool ismc)
double correct_prod3_fd_calE_high_gain(double calE, bool ismc)
Definition: run.py:1
double predict_special_fd_p2_muon_energy_corr(double trkLen, bool ismc)
static double spline3(double x, double offset, double slope0, double stitch1, double slope1, double stitch2, double slope2, double stitch3, double slope3)
#define NumuMuE2017NDW(f_act, f_cat)
double predict_special_fd_muon_energy(double trkLen, int run, bool ismc)
double predict_special_fd_p2_muon_energy(double trkLen)
double predict_special_fd_p5_had_energy_corr(double hadVisE, bool ismc)
double predict_special_fd_p3_had_energy(double hadVisE)
caf::Proxy< caf::SRKalman > kalman
Definition: SRProxy.h:1797
const Var kNumuHadE2017SpecialND(NumuHadE2017NDW(predict_special_nd_p3_had_energy))
#define NumuHadE2017FDW(f)
double predict_special_fd_p3_muon_energy_corr(double trkLen, bool ismc)
double predict_special_fd_p1_muon_energy(double trkLen)
caf::Proxy< float > ndtrklenact
Definition: SRProxy.h:186
caf::Proxy< bool > ismc
Definition: SRProxy.h:242
double predict_special_fd_p5_muon_energy(double trkLen)
double predict_special_nd_cat_energy(double ndtrklencat)
double correct_prod3_fd_calE_low_gain(double calE, bool ismc)
caf::Proxy< std::vector< caf::SRKalmanTrack > > tracks
Definition: SRProxy.h:1780
const Var kNumuMuE2017SpecialFDp2(NumuMuE2017FDW(predict_special_fd_p2_muon_energy_corr))
const Var kNumuMuE2017SpecialND(NumuMuE2017NDW(predict_special_nd_p3_act_energy, predict_special_nd_p3_cat_energy))
double predict_special_fd_p2_had_energy_corr(double hadVisE, bool ismc)
const Var kNumuHadE2017SpecialFDp1(NumuHadE2017FDW(predict_special_fd_p1_had_energy_corr))
caf::Proxy< float > ndtrklencat
Definition: SRProxy.h:187
static double spline2(double x, double offset, double slope0, double stitch1, double slope1, double stitch2, double slope2)
const Var kNumuMuE2017SpecialFDp3(NumuMuE2017FDW(predict_special_fd_p3_muon_energy_corr))
caf::Proxy< caf::Det_t > det
Definition: SRProxy.h:232