NusVarsTemp.cxx
Go to the documentation of this file.
2 
4 
5 #include <cstdlib>
6 #include "TString.h"
7 #include "TROOT.h"
8 #include "TSystem.h"
9 #include "TObjString.h"
10 #include "TMVA/Reader.h"
11 #include "TMVA/IMethod.h"
12 #include "TMVA/Factory.h"
13 #include "TMVA/DataLoader.h"
14 #include "TMVA/Tools.h"
15 
16 #pragma GCC diagnostic push
17 #pragma GCC diagnostic ignored "-Wunused-but-set-variable"
18 #include "NCID/func/KerasModel.h"
19 #pragma GCC diagnostic pop
20 
21 namespace ana
22 {
23 
24  // new for 2020
25  static TMVA::Reader* fReader2020BDTG = 0;
26 
27  float TMVA2020vars[21] = {};
28  float TMVA2020specs[4] = {};
29 
30  float vtxx = -5.f;
31  float vtxy = -5.f;
32  float vtxz = -5.f;
33  float closestslicemindist = -5.f;
34  float ncontplanes = -5.f;
35  float shwhitx = -5.f;
36  float shwhity = -5.f;
37  float shwhittot = -5.f;
38  float shwhitasymm = -5.f;
39  float shwhitratio = -5.f;
40  float nshowers = -5.f;
41  float showerwidth = -5.f;
42  float showerlength = -5.f;
43  float showerdirycosine = -5.f;
44  float showergap = -5.f;
45  float showercale = -5.f;
46  float nhitsperplane = -5.f;
47  float nmiphits = -5.f;
48  float nhitsperslice = -5.f;
49  float closestslicetime = -5.f;
50  float partptp = -5.f;
51  int run = -5;
52  int subrun = -5;
53  int evt = -5;
54  int subevt = -5;
55 
57  {
58  if (!fReader2020BDTG) InitTMVA();
59 
60  if(!sr->vtx.elastic.IsValid) return -5.f;
61  if(sr->vtx.elastic.fuzzyk.npng == 0) return -5.f;
62  if(sr->vtx.elastic.fuzzyk.nshwlid == 0) return -5.f;
63  if(std::isnan(1.*sr->sel.nuecosrej.partptp) ||
64  std::isinf(1.*sr->sel.nuecosrej.partptp)) return -5.f;
65 
66  // vetex variables
67  vtxx = sr->vtx.elastic.vtx.x;
68  vtxy = sr->vtx.elastic.vtx.y;
69  vtxz = sr->vtx.elastic.vtx.z;
70 
71  // shower hit variables
72  shwhitx = sr->vtx.elastic.fuzzyk.png[0].shwlid.nhitx;
73  shwhity = sr->vtx.elastic.fuzzyk.png[0].shwlid.nhity;
74  shwhittot = shwhitx + shwhity;
75  shwhitasymm = shwhitx - shwhity;
76  shwhitratio = (float)shwhitx/(float)shwhity;
77 
78  // other shower variables
79  nshowers = sr->vtx.elastic.fuzzyk.nshwlid;
80  showerwidth = sr->vtx.elastic.fuzzyk.png[0].shwlid.width;
81  showerlength = sr->vtx.elastic.fuzzyk.png[0].shwlid.len;
82  showerdirycosine = sr->vtx.elastic.fuzzyk.png[0].shwlid.dir.y;
83  showergap = sr->vtx.elastic.fuzzyk.png[0].shwlid.gap;
84  showercale = sr->vtx.elastic.fuzzyk.png[0].shwlid.calE;
85 
86  // other variables
87  nhitsperplane = sr->sel.nuecosrej.hitsperplane;
88  partptp = sr->sel.nuecosrej.partptp;
89  nmiphits = sr->slc.nmiphit;
90  nhitsperslice = sr->slc.nhit;
91  ncontplanes = sr->slc.ncontplanes;
92  closestslicetime = sr->slc.closestslicetime;
93  closestslicemindist = sr->slc.closestslicemindist;
94 
95  // spectators
96  run = sr->hdr.run;
97  subrun = sr->hdr.subrun;
98  evt = sr->hdr.evt;
99  subevt = sr->hdr.subevt;
100 
101  TMVA2020vars[0] = vtxx ;
102  TMVA2020vars[1] = vtxy ;
103  TMVA2020vars[2] = vtxz ;
104  TMVA2020vars[3] = closestslicemindist ;
105  TMVA2020vars[4] = ncontplanes ;
106  TMVA2020vars[5] = shwhitx ;
107  TMVA2020vars[6] = shwhity ;
108  TMVA2020vars[7] = shwhittot ;
109  TMVA2020vars[8] = shwhitasymm ;
110  TMVA2020vars[9] = shwhitratio ;
111  TMVA2020vars[10] = nshowers ;
112  TMVA2020vars[11] = showerwidth ;
113  TMVA2020vars[12] = showerlength ;
114  TMVA2020vars[13] = showerdirycosine ;
115  TMVA2020vars[14] = showergap ;
116  TMVA2020vars[15] = showercale ;
117  TMVA2020vars[16] = nhitsperplane ;
118  TMVA2020vars[17] = nmiphits ;
119  TMVA2020vars[18] = nhitsperslice ;
120  TMVA2020vars[19] = partptp ;
121  TMVA2020vars[20] = closestslicetime ;
122  TMVA2020specs[0] = run ;
123  TMVA2020specs[1] = subrun ;
124  TMVA2020specs[2] = evt ;
125  TMVA2020specs[3] = subevt ;
126 
127  float score = fReader2020BDTG->EvaluateMVA("BDTAGrad");
128  memset(TMVA2020vars, 0, sizeof(TMVA2020vars));
129  memset(TMVA2020specs, 0, sizeof(TMVA2020specs));
130 
131  return score;
132 
133  }
134 
136  {
137  if(!fReader2020BDTG) fReader2020BDTG = new TMVA::Reader;
138 
139  const char* libpath = getenv("NCID_LIB_PATH");
140  assert(libpath);
141  std::string pidlib = std::string(libpath)+"/TMVAClassifier_BDTAGrad.weights.xml";
142 
143  fReader2020BDTG->AddVariable("vtxx" , &TMVA2020vars[0]);
144  fReader2020BDTG->AddVariable("vtxy" , &TMVA2020vars[1]);
145  fReader2020BDTG->AddVariable("vtxz" , &TMVA2020vars[2]);
146  fReader2020BDTG->AddVariable("closestslicemindist" , &TMVA2020vars[3]);
147  fReader2020BDTG->AddVariable("ncontplanes" , &TMVA2020vars[4]);
148  fReader2020BDTG->AddVariable("shwhitx" , &TMVA2020vars[5]);
149  fReader2020BDTG->AddVariable("shwhity" , &TMVA2020vars[6]);
150  fReader2020BDTG->AddVariable("shwhittot" , &TMVA2020vars[7]);
151  fReader2020BDTG->AddVariable("shwhitasymm" , &TMVA2020vars[8]);
152  fReader2020BDTG->AddVariable("shwhitratio" , &TMVA2020vars[9]);
153  fReader2020BDTG->AddVariable("nshowers" , &TMVA2020vars[10]);
154  fReader2020BDTG->AddVariable("showerwidth" , &TMVA2020vars[11]);
155  fReader2020BDTG->AddVariable("showerlength" , &TMVA2020vars[12]);
156  fReader2020BDTG->AddVariable("showerdirycosine" , &TMVA2020vars[13]);
157  fReader2020BDTG->AddVariable("showergap" , &TMVA2020vars[14]);
158  fReader2020BDTG->AddVariable("showercale" , &TMVA2020vars[15]);
159  fReader2020BDTG->AddVariable("nhitsperplane" , &TMVA2020vars[16]);
160  fReader2020BDTG->AddVariable("nmiphits" , &TMVA2020vars[17]);
161  fReader2020BDTG->AddVariable("nhitsperslice" , &TMVA2020vars[18]);
162  fReader2020BDTG->AddVariable("partptp" , &TMVA2020vars[19]);
163  fReader2020BDTG->AddVariable("closestslicetime" , &TMVA2020vars[20]);
164  fReader2020BDTG->AddSpectator("run" , &TMVA2020specs[0]);
165  fReader2020BDTG->AddSpectator("subrun", &TMVA2020specs[1]);
166  fReader2020BDTG->AddSpectator("evt" , &TMVA2020specs[2]);
167  fReader2020BDTG->AddSpectator("subevt", &TMVA2020specs[3]);
168 
169  fReader2020BDTG->BookMVA("BDTAGrad", pidlib);
170 
171  }
172 
173 
174  //-------------------------------
175  // Pre-2020 analyses below here
176  //-------------------------------
177 
178  static TMVA::Reader* fReaderGBDT = 0;
179  static TMVA::Reader* fReaderBDT = 0;
180  static TMVA::Reader* fReaderBDT1 = 0;
181  static TMVA::Reader* fReaderBDT2 = 0;
182  static TMVA::Reader* fReaderBDT3 = 0;
183  static TMVA::Reader* fReaderBDT4 = 0;
184  float TMVAGvars[13] = {0}; // for GBDT-variant
185  float TMVAvars[13] = {0}; // for BDT
186 
187  float cosmicid= -5.f, shwnhit= -5.f, shwnhitx = -5.f, shwnhity= -5.f, shwxminusy= -5.f, shwxplusy= -5.f, shwxovery = -5.f, shwcalE= -5.f, shwlen = -5.f, shwwidth = -5.f, shwGap= -5.f, nshwlid = -5.f, nmiphit= -5.f, shwdirY= -5.f;
188  float ncid = -5.f;
189 
191 
193  {
194  if(!fReaderGBDT) InitTMVA();
195 
196  float ncid;
197  if(!sr->vtx.elastic.IsValid) return -5.f;
198  if(sr->vtx.elastic.fuzzyk.npng == 0) return -5.f;
199  if(sr->vtx.elastic.fuzzyk.nshwlid == 0) return -5.f;
200  std::cout << "ERROR::GetNCCosRejG. Looking for cvn2017. Branch no longer exists." << std::endl;
201  abort();
202  //if(sr->sel.cvn2017.noutput == 0) return -5.f; // Karl::Deleted branch
203  if(std::isnan(1.*sr->sel.nuecosrej.partptp) ||
204  std::isinf(1.*sr->sel.nuecosrej.partptp)) return -5.f;
205  //cosmicid = sr->sel.cvn2017.output[14]; // Karl::Deleted branch
206  partptp = sr->sel.nuecosrej.partptp;
207  shwnhit = (float)sr->vtx.elastic.fuzzyk.png[0].shwlid.nhit;
208  shwnhitx = (float)sr->vtx.elastic.fuzzyk.png[0].shwlid.nhitx;
209  shwnhity = (float)sr->vtx.elastic.fuzzyk.png[0].shwlid.nhity;
213  shwcalE = sr->vtx.elastic.fuzzyk.png[0].shwlid.calE;
214  shwdirY = sr->vtx.elastic.fuzzyk.png[0].shwlid.dir.y;
215  shwlen = sr->vtx.elastic.fuzzyk.png[0].shwlid.len;
216  shwwidth = sr->vtx.elastic.fuzzyk.png[0].shwlid.width;
217  shwGap = sr->vtx.elastic.fuzzyk.png[0].shwlid.gap;
219  nmiphit = sr->slc.nmiphit;
220 
221  TMVAGvars[0] = cosmicid; // ncid
222  TMVAGvars[1] = partptp; // partptp
223  TMVAGvars[2] = shwnhit; // shwnhit;
224  TMVAGvars[3] = shwxminusy; // shwnhitx - shwnhity
225  TMVAGvars[4] = shwxplusy; // shnhitx + shwnhity
226  TMVAGvars[5] = shwxovery; // shwxminusy/shwxplusy
227  TMVAGvars[6] = shwcalE; // shwcalE
228  TMVAGvars[7] = shwdirY; // shwdirY
229  TMVAGvars[8] = shwlen; // shwlen
230  TMVAGvars[9] = shwwidth; // shwwidth
231  TMVAGvars[10] = shwGap; // shwgap
232  TMVAGvars[11] = nshwlid; // nshwlid
233  TMVAGvars[12] = nmiphit; // nmiphit;
234 
235  ncid = fReaderGBDT->EvaluateMVA("BDTG");
236 
237  return ncid;
238 
239  }
240 
242  {
243  if(!fReaderBDT) InitTMVA();
244 
245  if(!sr->vtx.elastic.IsValid) return -5.f;
246  if(sr->vtx.elastic.fuzzyk.npng == 0) return -5.f;
247  if(sr->vtx.elastic.fuzzyk.nshwlid == 0) return -5.f;
248  std::cout << "ERROR::GetNCCosRej. Looking for cvn2017. Branch no longer exists." << std::endl;
249  abort();
250  //if(sr->sel.cvn2017.noutput == 0) return -5.f; // Karl::Deleted branch
251  if(std::isnan(1.*sr->sel.nuecosrej.partptp) ||
252  std::isinf(1.*sr->sel.nuecosrej.partptp)) return -5.f;
253  //cosmicid = sr->sel.cvn2017.output[14]; // Karl::Deleted branch
254  partptp = sr->sel.nuecosrej.partptp;
255  shwnhit = (float)sr->vtx.elastic.fuzzyk.png[0].shwlid.nhit;
256  shwnhitx = (float)sr->vtx.elastic.fuzzyk.png[0].shwlid.nhitx;
257  shwnhity = (float)sr->vtx.elastic.fuzzyk.png[0].shwlid.nhity;
261  shwcalE = sr->vtx.elastic.fuzzyk.png[0].shwlid.calE;
262  shwdirY = sr->vtx.elastic.fuzzyk.png[0].shwlid.dir.y;
263  shwlen = sr->vtx.elastic.fuzzyk.png[0].shwlid.len;
264  shwwidth = sr->vtx.elastic.fuzzyk.png[0].shwlid.width;
265  shwGap = sr->vtx.elastic.fuzzyk.png[0].shwlid.gap;
267  nmiphit = sr->slc.nmiphit;
268 
269  TMVAvars[0] = cosmicid; // ncid
270  TMVAvars[1] = partptp; // partptp
271  TMVAvars[2] = shwnhit; // shwnhit;
272  TMVAvars[3] = shwxminusy; // shwnhitx - shwnhity
273  TMVAvars[4] = shwxplusy; // shnhitx + shwnhity
274  TMVAvars[5] = shwxovery; // shwxminusy/shwxplusy
275  TMVAvars[6] = shwcalE; // shwcalE
276  TMVAvars[7] = shwdirY; // shwdirY
277  TMVAvars[8] = shwlen; // shwlen
278  TMVAvars[9] = shwwidth; // shwwidth
279  TMVAvars[10] = shwGap; // shwgap
280  TMVAvars[11] = nshwlid; // nshwlid
281  TMVAvars[12] = nmiphit; // nmiphit;
282 
283  ncid = fReaderBDT->EvaluateMVA("BDT");
284  memset(TMVAvars, 0, sizeof(TMVAvars));
285  return ncid;
286 
287  }
288 
290  {
291  if(!fReaderBDT1) InitTMVA();
292 
293  if(!sr->vtx.elastic.IsValid) return -5.f;
294  if(sr->vtx.elastic.fuzzyk.npng == 0) return -5.f;
295  if(sr->vtx.elastic.fuzzyk.nshwlid == 0) return -5.f;
296  std::cout << "ERROR::GetNCCosRejP1. Looking for cvn2017. Branch no longer exists." << std::endl;
297  abort();
298  //if(sr->sel.cvn2017.noutput == 0) return -5.f; // Karl::Deleted branch
299  if(std::isnan(1.*sr->sel.nuecosrej.partptp) ||
300  std::isinf(1.*sr->sel.nuecosrej.partptp)) return -5.f;
301  //cosmicid = sr->sel.cvn2017.output[14]; // Karl::Deleted branch
302  partptp = sr->sel.nuecosrej.partptp;
303  shwnhit = (float)sr->vtx.elastic.fuzzyk.png[0].shwlid.nhit;
304  shwnhitx = (float)sr->vtx.elastic.fuzzyk.png[0].shwlid.nhitx;
305  shwnhity = (float)sr->vtx.elastic.fuzzyk.png[0].shwlid.nhity;
309  shwcalE = sr->vtx.elastic.fuzzyk.png[0].shwlid.calE;
310  shwdirY = sr->vtx.elastic.fuzzyk.png[0].shwlid.dir.y;
311  shwlen = sr->vtx.elastic.fuzzyk.png[0].shwlid.len;
312  shwwidth = sr->vtx.elastic.fuzzyk.png[0].shwlid.width;
313  shwGap = sr->vtx.elastic.fuzzyk.png[0].shwlid.gap;
315  nmiphit = sr->slc.nmiphit;
316 
317  TMVAvars[0] = cosmicid; // ncid
318  TMVAvars[1] = partptp; // partptp
319  TMVAvars[2] = shwnhit; // shwnhit;
320  TMVAvars[3] = shwxminusy; // shwnhitx - shwnhity
321  TMVAvars[4] = shwxplusy; // shnhitx + shwnhity
322  TMVAvars[5] = shwxovery; // shwxminusy/shwxplusy
323  TMVAvars[6] = shwcalE; // shwcalE
324  TMVAvars[7] = shwdirY; // shwdirY
325  TMVAvars[8] = shwlen; // shwlen
326  TMVAvars[9] = shwwidth; // shwwidth
327  TMVAvars[10] = shwGap; // shwgap
328  TMVAvars[11] = nshwlid; // nshwlid
329  TMVAvars[12] = nmiphit; // nmiphit;
330 
331  ncid = fReaderBDT1->EvaluateMVA("BDT");
332  memset(TMVAvars, 0, sizeof(TMVAvars));
333  return ncid;
334 
335  }
336 
338  {
339  if(!fReaderBDT2) InitTMVA();
340 
341  if(!sr->vtx.elastic.IsValid) return -5.f;
342  if(sr->vtx.elastic.fuzzyk.npng == 0) return -5.f;
343  if(sr->vtx.elastic.fuzzyk.nshwlid == 0) return -5.f;
344  std::cout << "ERROR::GetNCCosRejP2. Looking for cvn2017. Branch no longer exists." << std::endl;
345  abort();
346  //if(sr->sel.cvn2017.noutput == 0) return -5.f; // Karl::Deleted branch
347  if(std::isnan(1.*sr->sel.nuecosrej.partptp) ||
348  std::isinf(1.*sr->sel.nuecosrej.partptp)) return -5.f;
349  //cosmicid = sr->sel.cvn2017.output[14]; // Karl::Deleted branch
350  partptp = sr->sel.nuecosrej.partptp;
351  shwnhit = (float)sr->vtx.elastic.fuzzyk.png[0].shwlid.nhit;
352  shwnhitx = (float)sr->vtx.elastic.fuzzyk.png[0].shwlid.nhitx;
353  shwnhity = (float)sr->vtx.elastic.fuzzyk.png[0].shwlid.nhity;
357  shwcalE = sr->vtx.elastic.fuzzyk.png[0].shwlid.calE;
358  shwdirY = sr->vtx.elastic.fuzzyk.png[0].shwlid.dir.y;
359  shwlen = sr->vtx.elastic.fuzzyk.png[0].shwlid.len;
360  shwwidth = sr->vtx.elastic.fuzzyk.png[0].shwlid.width;
361  shwGap = sr->vtx.elastic.fuzzyk.png[0].shwlid.gap;
363  nmiphit = sr->slc.nmiphit;
364 
365  TMVAvars[0] = cosmicid; // ncid
366  TMVAvars[1] = partptp; // partptp
367  TMVAvars[2] = shwnhit; // shwnhit;
368  TMVAvars[3] = shwxminusy; // shwnhitx - shwnhity
369  TMVAvars[4] = shwxplusy; // shnhitx + shwnhity
370  TMVAvars[5] = shwxovery; // shwxminusy/shwxplusy
371  TMVAvars[6] = shwcalE; // shwcalE
372  TMVAvars[7] = shwdirY; // shwdirY
373  TMVAvars[8] = shwlen; // shwlen
374  TMVAvars[9] = shwwidth; // shwwidth
375  TMVAvars[10] = shwGap; // shwgap
376  TMVAvars[11] = nshwlid; // nshwlid
377  TMVAvars[12] = nmiphit; // nmiphit;
378 
379  ncid = fReaderBDT2->EvaluateMVA("BDT");
380  memset(TMVAvars, 0, sizeof(TMVAvars));
381  return ncid;
382 
383  }
384 
385 
387  {
388  if(!fReaderBDT3) InitTMVA();
389 
390  if(!sr->vtx.elastic.IsValid) return -5.f;
391  if(sr->vtx.elastic.fuzzyk.npng == 0) return -5.f;
392  if(sr->vtx.elastic.fuzzyk.nshwlid == 0) return -5.f;
393  std::cout << "ERROR::GetNCCosRejP3_5. Looking for cvn2017. Branch no longer exists." << std::endl;
394  abort();
395  //if(sr->sel.cvn2017.noutput == 0) return -5.f; // Karl::Deleted branch
396  if(std::isnan(1.*sr->sel.nuecosrej.partptp) ||
397  std::isinf(1.*sr->sel.nuecosrej.partptp)) return -5.f;
398  //cosmicid = sr->sel.cvn2017.output[14]; // Karl::Deleted branch
399  partptp = sr->sel.nuecosrej.partptp;
400  shwnhit = (float)sr->vtx.elastic.fuzzyk.png[0].shwlid.nhit;
401  shwnhitx = (float)sr->vtx.elastic.fuzzyk.png[0].shwlid.nhitx;
402  shwnhity = (float)sr->vtx.elastic.fuzzyk.png[0].shwlid.nhity;
406  shwcalE = sr->vtx.elastic.fuzzyk.png[0].shwlid.calE;
407  shwdirY = sr->vtx.elastic.fuzzyk.png[0].shwlid.dir.y;
408  shwlen = sr->vtx.elastic.fuzzyk.png[0].shwlid.len;
409  shwwidth = sr->vtx.elastic.fuzzyk.png[0].shwlid.width;
410  shwGap = sr->vtx.elastic.fuzzyk.png[0].shwlid.gap;
412  nmiphit = sr->slc.nmiphit;
413 
414  TMVAvars[0] = cosmicid; // ncid
415  TMVAvars[1] = partptp; // partptp
416  TMVAvars[2] = shwnhit; // shwnhit;
417  TMVAvars[3] = shwxminusy; // shwnhitx - shwnhity
418  TMVAvars[4] = shwxplusy; // shnhitx + shwnhity
419  TMVAvars[5] = shwxovery; // shwxminusy/shwxplusy
420  TMVAvars[6] = shwcalE; // shwcalE
421  TMVAvars[7] = shwdirY; // shwdirY
422  TMVAvars[8] = shwlen; // shwlen
423  TMVAvars[9] = shwwidth; // shwwidth
424  TMVAvars[10] = shwGap; // shwgap
425  TMVAvars[11] = nshwlid; // nshwlid
426  TMVAvars[12] = nmiphit; // nmiphit;
427 
428  ncid = fReaderBDT3->EvaluateMVA("BDT");
429  memset(TMVAvars, 0, sizeof(TMVAvars));
430  return ncid;
431 
432  }
433 
435  {
436  if(!fReaderBDT4) InitTMVA();
437 
438  if(!sr->vtx.elastic.IsValid) return -5.f;
439  if(sr->vtx.elastic.fuzzyk.npng == 0) return -5.f;
440  if(sr->vtx.elastic.fuzzyk.nshwlid == 0) return -5.f;
441  std::cout << "ERROR::GetNCCosRejP4_6. Looking for cvn2017. Branch no longer exists." << std::endl;
442  abort();
443  //if(sr->sel.cvn2017.noutput == 0) return -5.f; // Karl::Deleted branch
444  if(std::isnan(1.*sr->sel.nuecosrej.partptp) ||
445  std::isinf(1.*sr->sel.nuecosrej.partptp)) return -5.f;
446  //cosmicid = sr->sel.cvn2017.output[14]; // Karl::Deleted branch
447  partptp = sr->sel.nuecosrej.partptp;
448  shwnhit = (float)sr->vtx.elastic.fuzzyk.png[0].shwlid.nhit;
449  shwnhitx = (float)sr->vtx.elastic.fuzzyk.png[0].shwlid.nhitx;
450  shwnhity = (float)sr->vtx.elastic.fuzzyk.png[0].shwlid.nhity;
454  shwcalE = sr->vtx.elastic.fuzzyk.png[0].shwlid.calE;
455  shwdirY = sr->vtx.elastic.fuzzyk.png[0].shwlid.dir.y;
456  shwlen = sr->vtx.elastic.fuzzyk.png[0].shwlid.len;
457  shwwidth = sr->vtx.elastic.fuzzyk.png[0].shwlid.width;
458  shwGap = sr->vtx.elastic.fuzzyk.png[0].shwlid.gap;
460  nmiphit = sr->slc.nmiphit;
461 
462  TMVAvars[0] = cosmicid; // ncid
463  TMVAvars[1] = partptp; // partptp
464  TMVAvars[2] = shwnhit; // shwnhit;
465  TMVAvars[3] = shwxminusy; // shwnhitx - shwnhity
466  TMVAvars[4] = shwxplusy; // shnhitx + shwnhity
467  TMVAvars[5] = shwxovery; // shwxminusy/shwxplusy
468  TMVAvars[6] = shwcalE; // shwcalE
469  TMVAvars[7] = shwdirY; // shwdirY
470  TMVAvars[8] = shwlen; // shwlen
471  TMVAvars[9] = shwwidth; // shwwidth
472  TMVAvars[10] = shwGap; // shwgap
473  TMVAvars[11] = nshwlid; // nshwlid
474  TMVAvars[12] = nmiphit; // nmiphit;
475 
476  ncid = fReaderBDT4->EvaluateMVA("BDT");
477  memset(TMVAvars, 0, sizeof(TMVAvars));
478  return ncid;
479 
480  }
481 
483  {
484  keras::DataChunk *sample = new keras::DataChunkFlat();
485 
486  if(!fModel) InitKeras();
487 
488  if(!sr->vtx.elastic.IsValid) return -5.f;
489  if(sr->vtx.elastic.fuzzyk.npng == 0) return -5.f;
490  if(sr->vtx.elastic.fuzzyk.nshwlid == 0) return -5.f;
491  std::cout << "ERROR::GetNCCosRejKeras. Looking for cvn2017. Branch no longer exists." << std::endl;
492  abort();
493  //if(sr->sel.cvn2017.noutput == 0) return -5.f; // Karl::Deleted branch
494  if(std::isnan(1.*sr->sel.nuecosrej.partptp) ||
495  std::isinf(1.*sr->sel.nuecosrej.partptp)) return -5.f;
496  //cosmicid = sr->sel.cvn2017.output[14]; // Karl::Deleted branch
497  partptp = sr->sel.nuecosrej.partptp;
498  shwnhit = (float)sr->vtx.elastic.fuzzyk.png[0].shwlid.nhit;
499  shwnhitx = (float)sr->vtx.elastic.fuzzyk.png[0].shwlid.nhitx;
500  shwnhity = (float)sr->vtx.elastic.fuzzyk.png[0].shwlid.nhity;
504  shwcalE = sr->vtx.elastic.fuzzyk.png[0].shwlid.calE;
505  shwdirY = sr->vtx.elastic.fuzzyk.png[0].shwlid.dir.y;
506  shwlen = sr->vtx.elastic.fuzzyk.png[0].shwlid.len;
507  shwwidth = sr->vtx.elastic.fuzzyk.png[0].shwlid.width;
508  shwGap = sr->vtx.elastic.fuzzyk.png[0].shwlid.gap;
510  nmiphit = sr->slc.nmiphit;
511 
512  float myfloats[] = {partptp,
513  cosmicid,
514  nmiphit,
515  shwcalE,
516  shwnhit,
517  nshwlid,
518  shwnhitx,
519  shwnhity,
520  shwxminusy,
521  shwxplusy,
522  shwxovery,
523  shwdirY,
524  shwlen,
525  shwwidth,
526  shwGap};
527 
528  std::vector<float> vec_myfloats(myfloats, myfloats + sizeof(myfloats) / sizeof(float));
529  sample->set_data(vec_myfloats);
530  ncid = fModel->compute_output(sample)[0];
531  if(std::isnan(ncid) || std::isinf(ncid)){
532  ncid = -5.;
533  }
534  delete sample;
535 
536  return ncid;
537 
538  }
539 
540 
542  {
543  if(!fReaderGBDT) fReaderGBDT = new TMVA::Reader;
544 
545  const char* libpath = getenv("NCID_LIB_PATH");
546  assert(libpath);
547  std::string pidlib = std::string(libpath)+"/NCCosRejPID_BDTG_nus.weights.xml";
548 
549  fReaderGBDT->AddVariable("cosmicid", &TMVAGvars[0]); //cosmicid
550  fReaderGBDT->AddVariable("partptp", &TMVAGvars[1]); //partptp
551  fReaderGBDT->AddVariable("shwnhit", &TMVAGvars[2]); //shwnhit
552  fReaderGBDT->AddVariable("shwxminusy", &TMVAGvars[3]); //shwminusy
553  fReaderGBDT->AddVariable("shwxplusy", &TMVAGvars[4]); //shwxplusy
554  fReaderGBDT->AddVariable("shwxovery", &TMVAGvars[5]); //shwxpvery
555  fReaderGBDT->AddVariable("shwcalE", &TMVAGvars[6]); //shwcalE
556  fReaderGBDT->AddVariable("shwdirY", &TMVAGvars[7]); //shwdirY
557  fReaderGBDT->AddVariable("shwlen", &TMVAGvars[8]); //shwlen
558  fReaderGBDT->AddVariable("shwwwidth", &TMVAGvars[9]); //shwwwidth
559  fReaderGBDT->AddVariable("shwGap", &TMVAGvars[10]); //shwGap
560  fReaderGBDT->AddVariable("nshwlid", &TMVAGvars[11]); //nshwlid
561  fReaderGBDT->AddVariable("nmiphit", &TMVAGvars[12]); //nmiphit
562  fReaderGBDT->BookMVA("BDTG", pidlib);
563  }
564 
566  {
567  if(!fReaderBDT) fReaderBDT = new TMVA::Reader;
568 
569  const char* libpath = getenv("NCID_LIB_PATH");
570  assert(libpath);
571  std::string pidlib = std::string(libpath)+"/NCCosRejPID_BDTA_nus.weights.xml";
572 
573  fReaderBDT->AddVariable("cosmicid", &TMVAvars[0]); //cosmicid
574  fReaderBDT->AddVariable("partptp", &TMVAvars[1]); //partptp
575  fReaderBDT->AddVariable("shwnhit", &TMVAvars[2]); //shwnhit
576  fReaderBDT->AddVariable("shwxminusy", &TMVAvars[3]); //shwminusy
577  fReaderBDT->AddVariable("shwxplusy", &TMVAvars[4]); //shwxplusy
578  fReaderBDT->AddVariable("shwxovery", &TMVAvars[5]); //shwxpvery
579  fReaderBDT->AddVariable("shwcalE", &TMVAvars[6]); //shwcalE
580  fReaderBDT->AddVariable("shwdirY", &TMVAvars[7]); //shwdirY
581  fReaderBDT->AddVariable("shwlen", &TMVAvars[8]); //shwlen
582  fReaderBDT->AddVariable("shwwwidth", &TMVAvars[9]); //shwwwidth
583  fReaderBDT->AddVariable("shwGap", &TMVAvars[10]); //shwGap
584  fReaderBDT->AddVariable("nshwlid", &TMVAvars[11]); //nshwlid
585  fReaderBDT->AddVariable("nmiphit", &TMVAvars[12]); //nmiphit
586  fReaderBDT->BookMVA("BDT", pidlib);
587  }
588 
589 
590  //BDT Trained with Period-1 cosmic data. "O_BDTA.weights.xml" is the weight file corresponding to that
592  {
593  if(!fReaderBDT1) fReaderBDT1 = new TMVA::Reader;
594  const char* libpath = getenv("NCID_LIB_PATH");
595  assert(libpath);
596  std::string pidlib = std::string(libpath)+"/O_BDTA.weights.xml";
597 
598  fReaderBDT1->AddVariable("cosmicid", &TMVAvars[0]); //cosmicid
599  fReaderBDT1->AddVariable("partptp", &TMVAvars[1]); //partptp
600  fReaderBDT1->AddVariable("shwnhit", &TMVAvars[2]); //shwnhit
601  fReaderBDT1->AddVariable("shwxminusy", &TMVAvars[3]); //shwminusy
602  fReaderBDT1->AddVariable("shwxplusy", &TMVAvars[4]); //shwxplusy
603  fReaderBDT1->AddVariable("shwxovery", &TMVAvars[5]); //shwxpvery
604  fReaderBDT1->AddVariable("shwcalE", &TMVAvars[6]); //shwcalE
605  fReaderBDT1->AddVariable("shwdirY", &TMVAvars[7]); //shwdirY
606  fReaderBDT1->AddVariable("shwlen", &TMVAvars[8]); //shwlen
607  fReaderBDT1->AddVariable("shwwwidth", &TMVAvars[9]); //shwwwidth
608  fReaderBDT1->AddVariable("shwGap", &TMVAvars[10]); //shwGap
609  fReaderBDT1->AddVariable("nshwlid", &TMVAvars[11]); //nshwlid
610  fReaderBDT1->AddVariable("nmiphit", &TMVAvars[12]); //nmiphit
611  fReaderBDT1->BookMVA("BDT", pidlib);
612  }
613 
614 
615  //BDT Trained with Period-2 cosmic data. "T_BDTA.weights.xml" is the weight file corresponding to that
616  //
618  {
619  if(!fReaderBDT2) fReaderBDT2 = new TMVA::Reader;
620  const char* libpath = getenv("NCID_LIB_PATH");
621  assert(libpath);
622  std::string pidlib = std::string(libpath)+"/T_BDTA.weights.xml";
623  fReaderBDT2->AddVariable("cosmicid", &TMVAvars[0]); //cosmicid
624  fReaderBDT2->AddVariable("partptp", &TMVAvars[1]); //partptp
625  fReaderBDT2->AddVariable("shwnhit", &TMVAvars[2]); //shwnhit
626  fReaderBDT2->AddVariable("shwxminusy", &TMVAvars[3]); //shwminusy
627  fReaderBDT2->AddVariable("shwxplusy", &TMVAvars[4]); //shwxplusy
628  fReaderBDT2->AddVariable("shwxovery", &TMVAvars[5]); //shwxpvery
629  fReaderBDT2->AddVariable("shwcalE", &TMVAvars[6]); //shwcalE
630  fReaderBDT2->AddVariable("shwdirY", &TMVAvars[7]); //shwdirY
631  fReaderBDT2->AddVariable("shwlen", &TMVAvars[8]); //shwlen
632  fReaderBDT2->AddVariable("shwwwidth", &TMVAvars[9]); //shwwwidth
633  fReaderBDT2->AddVariable("shwGap", &TMVAvars[10]); //shwGap
634  fReaderBDT2->AddVariable("nshwlid", &TMVAvars[11]); //nshwlid
635  fReaderBDT2->AddVariable("nmiphit", &TMVAvars[12]); //nmiphit
636  fReaderBDT2->BookMVA("BDT", pidlib);
637  }
638 
639 
640  //BDT Trained with Period-3 and period-5 cosmic data. "TF_BDTA.weights.xml" is the weight file corresponding to that
642  {
643  if(!fReaderBDT3) fReaderBDT3 = new TMVA::Reader;
644  const char* libpath = getenv("NCID_LIB_PATH");
645  assert(libpath);
646  std::string pidlib = std::string(libpath)+"/TF_BDTA.weights.xml";
647  fReaderBDT3->AddVariable("cosmicid", &TMVAvars[0]); //cosmicid
648  fReaderBDT3->AddVariable("partptp", &TMVAvars[1]); //partptp
649  fReaderBDT3->AddVariable("shwnhit", &TMVAvars[2]); //shwnhit
650  fReaderBDT3->AddVariable("shwxminusy", &TMVAvars[3]); //shwminusy
651  fReaderBDT3->AddVariable("shwxplusy", &TMVAvars[4]); //shwxplusy
652  fReaderBDT3->AddVariable("shwxovery", &TMVAvars[5]); //shwxpvery
653  fReaderBDT3->AddVariable("shwcalE", &TMVAvars[6]); //shwcalE
654  fReaderBDT3->AddVariable("shwdirY", &TMVAvars[7]); //shwdirY
655  fReaderBDT3->AddVariable("shwlen", &TMVAvars[8]); //shwlen
656  fReaderBDT3->AddVariable("shwwwidth", &TMVAvars[9]); //shwwwidth
657  fReaderBDT3->AddVariable("shwGap", &TMVAvars[10]); //shwGap
658  fReaderBDT3->AddVariable("nshwlid", &TMVAvars[11]); //nshwlid
659  fReaderBDT3->AddVariable("nmiphit", &TMVAvars[12]); //nmiphit
660  fReaderBDT3->BookMVA("BDT", pidlib);
661  }
662 
663  //BDT Trained with Period-4 and period-6 cosmic data. "SF_BDTA.weights.xml" is the weight file corresponding to that (RHC Specific)
664 
666  {
667  if(!fReaderBDT4) fReaderBDT4 = new TMVA::Reader;
668 
669  const char* libpath = getenv("NCID_LIB_PATH");
670  assert(libpath);
671  std::string pidlib = std::string(libpath)+"/SF_BDTA.weights.xml";
672  fReaderBDT4->AddVariable("cosmicid", &TMVAvars[0]); //cosmicid
673  fReaderBDT4->AddVariable("partptp", &TMVAvars[1]); //partptp
674  fReaderBDT4->AddVariable("shwnhit", &TMVAvars[2]); //shwnhit
675  fReaderBDT4->AddVariable("shwxminusy", &TMVAvars[3]); //shwminusy
676  fReaderBDT4->AddVariable("shwxplusy", &TMVAvars[4]); //shwxplusy
677  fReaderBDT4->AddVariable("shwxovery", &TMVAvars[5]); //shwxpvery
678  fReaderBDT4->AddVariable("shwcalE", &TMVAvars[6]); //shwcalE
679  fReaderBDT4->AddVariable("shwdirY", &TMVAvars[7]); //shwdirY
680  fReaderBDT4->AddVariable("shwlen", &TMVAvars[8]); //shwlen
681  fReaderBDT4->AddVariable("shwwwidth", &TMVAvars[9]); //shwwwidth
682  fReaderBDT4->AddVariable("shwGap", &TMVAvars[10]); //shwGap
683  fReaderBDT4->AddVariable("nshwlid", &TMVAvars[11]); //nshwlid
684  fReaderBDT4->AddVariable("nmiphit", &TMVAvars[12]); //nmiphit
685  fReaderBDT4->BookMVA("BDT", pidlib);
686  }
688  {
689 
690  if(!fModel){
691  const char* libpath = getenv("NCID_LIB_PATH");
692  assert(libpath);
693  std::string keras_model = std::string(libpath)+"/keras_nn_5layer_v2_weight.nnet";
694  fModel = new keras::KerasModel(keras_model,true);
695  }
696 
697  }
698 }
699 
float nshwlid
float shwhitratio
Definition: NusVarsTemp.cxx:39
caf::Proxy< size_t > npng
Definition: SRProxy.h:2021
caf::Proxy< unsigned int > nshwlid
Definition: SRProxy.h:2023
void InitTMVA() const
Definition: FillPIDs.h:18
int isinf(const stan::math::var &a)
Definition: std_isinf.hpp:16
static TMVA::Reader * fReaderBDT
Definition: NDNCPi0Xsec.cxx:11
caf::Proxy< caf::SRFuzzyK > fuzzyk
Definition: SRProxy.h:2042
float partptp
Definition: NusVarsTemp.cxx:50
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
static TMVA::Reader * fReaderGBDT
float shwwidth
float operator()(const caf::SRProxy *sr) const
float nshowers
Definition: NusVarsTemp.cxx:40
float closestslicetime
Definition: NusVarsTemp.cxx:49
float shwnhity
caf::Proxy< unsigned int > subrun
Definition: SRProxy.h:253
float shwGap
float shwhittot
Definition: NusVarsTemp.cxx:37
caf::Proxy< caf::SRHeader > hdr
Definition: SRProxy.h:2120
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2109
float operator()(const caf::SRProxy *sr) const
static TMVA::Reader * fReaderBDT3
caf::Proxy< unsigned int > ncontplanes
Definition: SRProxy.h:1293
float shwhitasymm
Definition: NusVarsTemp.cxx:38
float showergap
Definition: NusVarsTemp.cxx:44
float shwxplusy
float vtxz
Definition: NusVarsTemp.cxx:32
float shwnhitx
float operator()(const caf::SRProxy *sr) const
Definition: NusVarsTemp.cxx:56
std::vector< float > compute_output(keras::DataChunk *dc)
Definition: KerasModel.cxx:354
float shwcalE
caf::Proxy< unsigned int > run
Definition: SRProxy.h:247
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
static TMVA::Reader * fReader2020BDTG
Definition: NusVarsTemp.cxx:25
float cosmicid
caf::Proxy< float > closestslicetime
Definition: SRProxy.h:1281
caf::Proxy< caf::SRElastic > elastic
Definition: SRProxy.h:2101
caf::Proxy< caf::SRNueCosRej > nuecosrej
Definition: SRProxy.h:1244
caf::Proxy< std::vector< caf::SRFuzzyKProng > > png
Definition: SRProxy.h:2026
float TMVA2020specs[4]
Definition: NusVarsTemp.cxx:28
caf::Proxy< unsigned int > nhit
Definition: SRProxy.h:1294
void InitTMVA() const
caf::Proxy< float > z
Definition: SRProxy.h:107
float nhitsperslice
Definition: NusVarsTemp.cxx:48
float shwlen
float nmiphit
caf::Proxy< short unsigned int > subevt
Definition: SRProxy.h:249
float nhitsperplane
Definition: NusVarsTemp.cxx:46
caf::Proxy< float > x
Definition: SRProxy.h:105
static keras::KerasModel * fModel
float TMVA2020vars[21]
Definition: NusVarsTemp.cxx:27
caf::Proxy< unsigned int > nmiphit
Definition: SRProxy.h:1295
std::string getenv(std::string const &name)
static TMVA::Reader * fReaderBDT2
float shwhitx
Definition: NusVarsTemp.cxx:35
static float TMVAvars[4]
Definition: NDNCPi0Xsec.cxx:10
caf::Proxy< float > hitsperplane
Definition: SRProxy.h:1025
void InitTMVA() const
void InitTMVA() const
caf::StandardRecord * sr
caf::Proxy< float > partptp
Definition: SRProxy.h:1036
int subevt
Definition: NusVarsTemp.cxx:54
int evt
Definition: NusVarsTemp.cxx:53
float operator()(const caf::SRProxy *sr) const
void InitTMVA() const
float operator()(const caf::SRProxy *sr) const
float showerwidth
Definition: NusVarsTemp.cxx:41
float shwdirY
Definition: run.py:1
static TMVA::Reader * fReaderBDT1
caf::Proxy< unsigned int > evt
Definition: SRProxy.h:236
OStream cout
Definition: OStream.cxx:6
float ncid
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void InitTMVA() const
int run
Definition: NusVarsTemp.cxx:51
caf::Proxy< bool > IsValid
Definition: SRProxy.h:2041
caf::Proxy< caf::SRSlice > slc
Definition: SRProxy.h:2125
float vtxy
Definition: NusVarsTemp.cxx:31
float shwxovery
caf::Proxy< float > y
Definition: SRProxy.h:106
void InitKeras() const
caf::Proxy< float > closestslicemindist
Definition: SRProxy.h:1273
float shwhity
Definition: NusVarsTemp.cxx:36
float nmiphits
Definition: NusVarsTemp.cxx:47
float shwxminusy
float operator()(const caf::SRProxy *sr) const
float ncontplanes
Definition: NusVarsTemp.cxx:34
caf::Proxy< caf::SRVector3D > vtx
Definition: SRProxy.h:2056
assert(nhit_max >=nhit_nbins)
float vtxx
Definition: NusVarsTemp.cxx:30
float shwnhit
caf::Proxy< caf::SRIDBranch > sel
Definition: SRProxy.h:2124
float closestslicemindist
Definition: NusVarsTemp.cxx:33
caf::Proxy< caf::SRVertexBranch > vtx
Definition: SRProxy.h:2129
float operator()(const caf::SRProxy *sr) const
static TMVA::Reader * fReaderBDT4
float TMVAGvars[13]
float showercale
Definition: NusVarsTemp.cxx:45
float operator()(const caf::SRProxy *sr) const
int subrun
Definition: NusVarsTemp.cxx:52
float showerlength
Definition: NusVarsTemp.cxx:42
float showerdirycosine
Definition: NusVarsTemp.cxx:43