NDXSecTrackPID.cxx
Go to the documentation of this file.
1 #include "NDAna/numucc_ccPi/NDXSecTrackPID.h"
2 
3 #include "NDAna/numucc_inc/NDXSecMuonPID.h"
4 #include "3FlavorAna/Vars/NumuVars.h"
5 
6 #include "CAFAna/Core/Utilities.h"
7 
8 #include "StandardRecord/Proxy/SRProxy.h"
9 
10 #include "TMVA/Reader.h"
11 
12 namespace{
13 const TVector3 beamDirND = ana::NuMIBeamDirection(caf::kNEARDET);
14 }
15 
16 namespace ana
17 {
18  static float TMVAvars[6] = {-999.0};
19  static TMVA::Reader* fReaderBDTPion = 0;
20  static TMVA::Reader* fReaderBDTProton = 0;
21  static TMVA::Reader* fReaderBDTGamma = 0;
22 
23  const TVector3 beamDirND = ana::NuMIBeamDirection(caf::kNEARDET);
24 
25 
26  //-------------------------------------------------------------
27  float GetPionID::operator() ( const caf::SRProxy* sr) const
28  {
29 
30  float highestpid = -999.0;
31  float pionbdtvalue = -999.0;
32  float protonbdtvalue = -999.0;
33  float gammabdtvalue = -999.0;
34 
35  if(!fReaderBDTPion) InitTMVAPion();
36  if(!fReaderBDTProton) InitTMVAProton();
37  if(!fReaderBDTGamma) InitTMVAGamma();
38  float dedxll = -999.0, scatll = -999.0;
39  float fAvededxlast10cm = -999.0;
40  float fAvededxlast40cm = -999.0;
41  float fTrackStartGap = -999.0;
42  float fGapDensity = -999.0;
43 
44  int nkals = sr->trk.kalman.ntracks;
45 
46  for (int itrk = 0; itrk < nkals; ++itrk){
47  if (sr->trk.kalman.tracks[itrk].rempid == -1) continue;
48  if (itrk==GetBestTrack()(sr)) continue;
49 
50  scatll = sr->trk.kalman.tracks[itrk].scatllh2017;
51  dedxll = sr->trk.kalman.tracks[itrk].dedxllh2017;
52 
53  if (sr->trk.kalman.tracks[itrk].avedEdxlast10cm == -5) continue;
54  if (sr->trk.kalman.tracks[itrk].avedEdxlast30cm == -5) continue;
55  if (sr->trk.kalman.tracks[itrk].avedEdxlast20cm == -5) continue;
56  if (sr->trk.kalman.tracks[itrk].avedEdxlast40cm == -5) continue;
57 
58  if (sr->trk.kalman.tracks[itrk].avedEdxlast10cm > 30){
59  fAvededxlast10cm = 30;
60  }else{
61  fAvededxlast10cm = sr->trk.kalman.tracks[itrk].avedEdxlast10cm;
62  }
63 
64  if (sr->trk.kalman.tracks[itrk].avedEdxlast40cm > 30){
65  fAvededxlast40cm = 30;
66  }else{
67  fAvededxlast40cm = sr->trk.kalman.tracks[itrk].avedEdxlast40cm;
68  }
69 
70  TVector3 vtxpos(sr->trk.kalman.tracks[sr->trk.kalman.idxremid].start.X(),sr->trk.kalman.tracks[sr->trk.kalman.idxremid].start.Y(),sr->trk.kalman.tracks[sr->trk.kalman.idxremid].start.Z());
71  TVector3 trkpos(sr->trk.kalman.tracks[itrk].start.X(),sr->trk.kalman.tracks[itrk].start.Y(),sr->trk.kalman.tracks[itrk].start.Z());
72  float trackstartgap = (vtxpos-trkpos).Mag();
73  TVector3 trkvec = sr->trk.kalman.tracks[itrk].dir;
74  float angle = trkvec.Angle(beamDirND);
75  float length = sr->trk.kalman.tracks[itrk].len;
76  float coslength = length*cos(angle);
77  float ngaps = sr->trk.kalman.tracks[itrk].nplanegap;
78 
79  fTrackStartGap=trackstartgap;
80  fGapDensity=ngaps/coslength;
81 
82  TMVAvars[0] = dedxll;
83  TMVAvars[1] = scatll;
84  TMVAvars[2] = fAvededxlast10cm;
85  TMVAvars[3] = fAvededxlast40cm;
86  TMVAvars[4] = fTrackStartGap;
87  TMVAvars[5] = fGapDensity;
88 
89  if(!fReaderBDTPion) InitTMVAPion();
90  if(!fReaderBDTProton) InitTMVAProton();
91  if(!fReaderBDTGamma) InitTMVAGamma();
92 
93  pionbdtvalue = fReaderBDTPion->EvaluateMVA("BDTGPion");
94  protonbdtvalue = fReaderBDTProton->EvaluateMVA("BDTGProton");
95  gammabdtvalue = fReaderBDTGamma->EvaluateMVA("BDTGGamma");
96 
97  if (pionbdtvalue > protonbdtvalue && pionbdtvalue > gammabdtvalue && pionbdtvalue > highestpid){
98  highestpid = pionbdtvalue;
99  }else continue;
100 
101  }
102  return highestpid;
103  }
104 
105 
106  float GetPionTrack::operator() (const caf::SRProxy* sr) const
107  {
108 
109  float highestpid = -999.0;
110  int imax = -1; //this can be used as your track index with highest PID score.
111  float pionbdtvalue = -999.0;
112  float protonbdtvalue = -999.0;
113  float gammabdtvalue = -999.0;
114 
115  if(!fReaderBDTPion) InitTMVAPion();
116  if(!fReaderBDTProton) InitTMVAProton();
117  if(!fReaderBDTGamma) InitTMVAGamma();
118  float dedxll = -999.0, scatll = -999.0;
119  float fAvededxlast10cm = -999.0;
120  float fAvededxlast40cm = -999.0;
121  float fTrackStartGap = -999.0;
122  float fGapDensity = -999.0;
123 
124  int nkals = sr->trk.kalman.ntracks;
125 
126  for (int itrk = 0; itrk < nkals; ++itrk){
127  if (sr->trk.kalman.tracks[itrk].rempid == -1) continue;
128  if (itrk==GetBestTrack()(sr)) continue;
129 
130  scatll = sr->trk.kalman.tracks[itrk].scatllh2017;
131  dedxll = sr->trk.kalman.tracks[itrk].dedxllh2017;
132 
133  if (sr->trk.kalman.tracks[itrk].avedEdxlast10cm == -5) continue;
134  if (sr->trk.kalman.tracks[itrk].avedEdxlast30cm == -5) continue;
135  if (sr->trk.kalman.tracks[itrk].avedEdxlast20cm == -5) continue;
136  if (sr->trk.kalman.tracks[itrk].avedEdxlast40cm == -5) continue;
137 
138  if (sr->trk.kalman.tracks[itrk].avedEdxlast10cm > 30){
139  fAvededxlast10cm = 30;
140  }else{
141  fAvededxlast10cm = sr->trk.kalman.tracks[itrk].avedEdxlast10cm;
142  }
143 
144  if (sr->trk.kalman.tracks[itrk].avedEdxlast40cm > 30){
145  fAvededxlast40cm = 30;
146  }else{
147  fAvededxlast40cm = sr->trk.kalman.tracks[itrk].avedEdxlast40cm;
148  }
149 
150  TVector3 vtxpos(sr->trk.kalman.tracks[sr->trk.kalman.idxremid].start.X(),sr->trk.kalman.tracks[sr->trk.kalman.idxremid].start.Y(),sr->trk.kalman.tracks[sr->trk.kalman.idxremid].start.Z());
151  TVector3 trkpos(sr->trk.kalman.tracks[itrk].start.X(),sr->trk.kalman.tracks[itrk].start.Y(),sr->trk.kalman.tracks[itrk].start.Z());
152  float trackstartgap = (vtxpos-trkpos).Mag();
153  TVector3 trkvec = sr->trk.kalman.tracks[itrk].dir;
154  float angle = trkvec.Angle(beamDirND);
155  float length = sr->trk.kalman.tracks[itrk].len;
156  float coslength = length*cos(angle);
157  float ngaps = sr->trk.kalman.tracks[itrk].nplanegap;
158 
159  fTrackStartGap=trackstartgap;
160  fGapDensity=ngaps/coslength;
161 
162  TMVAvars[0] = dedxll;
163  TMVAvars[1] = scatll;
164  TMVAvars[2] = fAvededxlast10cm;
165  TMVAvars[3] = fAvededxlast40cm;
166  TMVAvars[4] = fTrackStartGap;
167  TMVAvars[5] = fGapDensity;
168 
169  if(!fReaderBDTPion) InitTMVAPion();
170  if(!fReaderBDTProton) InitTMVAProton();
171  if(!fReaderBDTGamma) InitTMVAGamma();
172 
173  pionbdtvalue = fReaderBDTPion->EvaluateMVA("BDTGPion");
174  protonbdtvalue = fReaderBDTProton->EvaluateMVA("BDTGProton");
175  gammabdtvalue = fReaderBDTGamma->EvaluateMVA("BDTGGamma");
176 
177  if (pionbdtvalue > protonbdtvalue && pionbdtvalue > gammabdtvalue && pionbdtvalue > highestpid){
178  imax = itrk;
179  highestpid = pionbdtvalue;
180  }else continue;
181  }
182  return imax;
183  }
184 
185 
186  //-------------------------------------------------------------
187  float GetBestPionID::operator() ( const caf::SRProxy* sr) const
188  {
189 
190  float highestpid = -999.0;
191  float pionbdtvalue = -999.0;
192 
193  if(!fReaderBDTPion) InitTMVAPion();
194  float dedxll = -999.0, scatll = -999.0;
195  float fAvededxlast10cm = -999.0;
196  float fAvededxlast40cm = -999.0;
197  float fTrackStartGap = -999.0;
198  float fGapDensity = -999.0;
199 
200  int nkals = sr->trk.kalman.ntracks;
201 
202  for (int itrk = 0; itrk < nkals; ++itrk){
203  if (sr->trk.kalman.tracks[itrk].rempid == -1) continue;
204  if (itrk==GetBestTrack()(sr)) continue;
205 
206  scatll = sr->trk.kalman.tracks[itrk].scatllh2017;
207  dedxll = sr->trk.kalman.tracks[itrk].dedxllh2017;
208 
209  if (sr->trk.kalman.tracks[itrk].avedEdxlast10cm == -5) continue;
210  if (sr->trk.kalman.tracks[itrk].avedEdxlast30cm == -5) continue;
211  if (sr->trk.kalman.tracks[itrk].avedEdxlast20cm == -5) continue;
212  if (sr->trk.kalman.tracks[itrk].avedEdxlast40cm == -5) continue;
213 
214  if (sr->trk.kalman.tracks[itrk].avedEdxlast10cm > 30){
215  fAvededxlast10cm = 30;
216  }else{
217  fAvededxlast10cm = sr->trk.kalman.tracks[itrk].avedEdxlast10cm;
218  }
219 
220  if (sr->trk.kalman.tracks[itrk].avedEdxlast40cm > 30){
221  fAvededxlast40cm = 30;
222  }else{
223  fAvededxlast40cm = sr->trk.kalman.tracks[itrk].avedEdxlast40cm;
224  }
225 
226  TVector3 vtxpos(sr->trk.kalman.tracks[sr->trk.kalman.idxremid].start.X(),sr->trk.kalman.tracks[sr->trk.kalman.idxremid].start.Y(),sr->trk.kalman.tracks[sr->trk.kalman.idxremid].start.Z());
227  TVector3 trkpos(sr->trk.kalman.tracks[itrk].start.X(),sr->trk.kalman.tracks[itrk].start.Y(),sr->trk.kalman.tracks[itrk].start.Z());
228  float trackstartgap = (vtxpos-trkpos).Mag();
229  TVector3 trkvec = sr->trk.kalman.tracks[itrk].dir;
230  float angle = trkvec.Angle(beamDirND);
231  float length = sr->trk.kalman.tracks[itrk].len;
232  float coslength = length*cos(angle);
233  float ngaps = sr->trk.kalman.tracks[itrk].nplanegap;
234 
235  fTrackStartGap=trackstartgap;
236  fGapDensity=ngaps/coslength;
237 
238  TMVAvars[0] = dedxll;
239  TMVAvars[1] = scatll;
240  TMVAvars[2] = fAvededxlast10cm;
241  TMVAvars[3] = fAvededxlast40cm;
242  TMVAvars[4] = fTrackStartGap;
243  TMVAvars[5] = fGapDensity;
244 
245  if(!fReaderBDTPion) InitTMVAPion();
246 
247  pionbdtvalue = fReaderBDTPion->EvaluateMVA("BDTGPion");
248 
249  if (pionbdtvalue > highestpid){
250  highestpid = pionbdtvalue;
251  }else continue;
252 
253  }
254  return highestpid;
255  }
256 
257 
258  float GetBestPionTrack::operator() (const caf::SRProxy* sr) const
259  {
260 
261  float highestpid = -999.0;
262  int imax = -1; //this can be used as your track index with highest PID score.
263  float pionbdtvalue = -999.0;
264 
265  if(!fReaderBDTPion) InitTMVAPion();
266  float dedxll = -999.0, scatll = -999.0;
267  float fAvededxlast10cm = -999.0;
268  float fAvededxlast40cm = -999.0;
269  float fTrackStartGap = -999.0;
270  float fGapDensity = -999.0;
271 
272  int nkals = sr->trk.kalman.ntracks;
273 
274  for (int itrk = 0; itrk < nkals; ++itrk){
275  if (sr->trk.kalman.tracks[itrk].rempid == -1) continue;
276  if (itrk==GetBestTrack()(sr)) continue;
277 
278  scatll = sr->trk.kalman.tracks[itrk].scatllh2017;
279  dedxll = sr->trk.kalman.tracks[itrk].dedxllh2017;
280 
281  if (sr->trk.kalman.tracks[itrk].avedEdxlast10cm == -5) continue;
282  if (sr->trk.kalman.tracks[itrk].avedEdxlast30cm == -5) continue;
283  if (sr->trk.kalman.tracks[itrk].avedEdxlast20cm == -5) continue;
284  if (sr->trk.kalman.tracks[itrk].avedEdxlast40cm == -5) continue;
285 
286  if (sr->trk.kalman.tracks[itrk].avedEdxlast10cm > 30){
287  fAvededxlast10cm = 30;
288  }else{
289  fAvededxlast10cm = sr->trk.kalman.tracks[itrk].avedEdxlast10cm;
290  }
291 
292  if (sr->trk.kalman.tracks[itrk].avedEdxlast40cm > 30){
293  fAvededxlast40cm = 30;
294  }else{
295  fAvededxlast40cm = sr->trk.kalman.tracks[itrk].avedEdxlast40cm;
296  }
297 
298  TVector3 vtxpos(sr->trk.kalman.tracks[sr->trk.kalman.idxremid].start.X(),sr->trk.kalman.tracks[sr->trk.kalman.idxremid].start.Y(),sr->trk.kalman.tracks[sr->trk.kalman.idxremid].start.Z());
299  TVector3 trkpos(sr->trk.kalman.tracks[itrk].start.X(),sr->trk.kalman.tracks[itrk].start.Y(),sr->trk.kalman.tracks[itrk].start.Z());
300  float trackstartgap = (vtxpos-trkpos).Mag();
301  TVector3 trkvec = sr->trk.kalman.tracks[itrk].dir;
302  float angle = trkvec.Angle(beamDirND);
303  float length = sr->trk.kalman.tracks[itrk].len;
304  float coslength = length*cos(angle);
305  float ngaps = sr->trk.kalman.tracks[itrk].nplanegap;
306 
307  fTrackStartGap=trackstartgap;
308  fGapDensity=ngaps/coslength;
309 
310  TMVAvars[0] = dedxll;
311  TMVAvars[1] = scatll;
312  TMVAvars[2] = fAvededxlast10cm;
313  TMVAvars[3] = fAvededxlast40cm;
314  TMVAvars[4] = fTrackStartGap;
315  TMVAvars[5] = fGapDensity;
316 
317  if(!fReaderBDTPion) InitTMVAPion();
318 
319  pionbdtvalue = fReaderBDTPion->EvaluateMVA("BDTGPion");
320 
321  if (pionbdtvalue > highestpid){
322  imax = itrk;
323  highestpid = pionbdtvalue;
324  }else continue;
325  }
326  return imax;
327  }
328 
329  float GetPionIDVal::operator() ( float dedxll, float scatll, float avededxlast10cm, float avededxlast40cm, float trackstartgap, float gapdensity) const
330  {
331  float pionbdtvalue = -999.0;
332 
333  if(!fReaderBDTPion) InitTMVAPion();
334  /*
335  float dedxll = -999.0, scatll = -999.0;
336  float fAvededxlast10cm = -999.0;
337  float fAvededxlast40cm = -999.0;
338  float fTrackStartGap = -999.0;
339  float fGapDensity = -999.0;
340 
341  //int nkals = sr->trk.kalman.ntracks;
342 
343 
344  if (sr->trk.kalman.tracks[idx].rempid == -1) return -2.0f;
345  //if (idx==GetBestTrack()(sr)) return -2.0f;
346 
347  scatll = sr->trk.kalman.tracks[idx].scatllh2017;
348  dedxll = sr->trk.kalman.tracks[idx].dedxllh2017;
349 
350  if (sr->trk.kalman.tracks[idx].avedEdxlast10cm == -5) return -2.0f;
351  if (sr->trk.kalman.tracks[idx].avedEdxlast30cm == -5) return -2.0f;
352  if (sr->trk.kalman.tracks[idx].avedEdxlast20cm == -5) return -2.0f;
353  if (sr->trk.kalman.tracks[idx].avedEdxlast40cm == -5) return -2.0f;
354 
355  if (sr->trk.kalman.tracks[idx].avedEdxlast10cm > 30){
356  fAvededxlast10cm = 30;
357  }else{
358  fAvededxlast10cm = sr->trk.kalman.tracks[idx].avedEdxlast10cm;
359  }
360 
361  if (sr->trk.kalman.tracks[idx].avedEdxlast40cm > 30){
362  fAvededxlast40cm = 30;
363  }else{
364  fAvededxlast40cm = sr->trk.kalman.tracks[idx].avedEdxlast40cm;
365  }
366 
367  TVector3 vtxpos(sr->trk.kalman.tracks[sr->trk.kalman.idxremid].start.X(),sr->trk.kalman.tracks[sr->trk.kalman.idxremid].start.Y(),sr->trk.kalman.tracks[sr->trk.kalman.idxremid].start.Z());
368  TVector3 trkpos(sr->trk.kalman.tracks[idx].start.X(),sr->trk.kalman.tracks[idx].start.Y(),sr->trk.kalman.tracks[idx].start.Z());
369  float trackstartgap = (vtxpos-trkpos).Mag();
370  TVector3 trkvec = sr->trk.kalman.tracks[idx].dir;
371  float angle = trkvec.Angle(beamDirND);
372  float length = sr->trk.kalman.tracks[idx].len;
373  float coslength = length*cos(angle);
374  float ngaps = sr->trk.kalman.tracks[idx].nplanegap;
375 
376  fTrackStartGap=trackstartgap;
377  fGapDensity=ngaps/coslength;*/
378 
379  TMVAvars[0] = dedxll;
380  TMVAvars[1] = scatll;
381  TMVAvars[2] = avededxlast10cm;
382  TMVAvars[3] = avededxlast40cm;
383  TMVAvars[4] = trackstartgap;
384  TMVAvars[5] = gapdensity;
385 
386  if(!fReaderBDTPion) InitTMVAPion();
387 
388  pionbdtvalue = fReaderBDTPion->EvaluateMVA("BDTGPion");
389 
390  return pionbdtvalue;
391  }
392 
393  float GetProtonIDVal::operator() ( float dedxll, float scatll, float avededxlast10cm, float avededxlast40cm, float trackstartgap, float gapdensity) const
394  {
395 
396  float protonbdtvalue = -999.0;
397 
398  if(!fReaderBDTProton) InitTMVAProton();
399  /*
400  float dedxll = -999.0, scatll = -999.0;
401  float fAvededxlast10cm = -999.0;
402  float fAvededxlast40cm = -999.0;
403  float fTrackStartGap = -999.0;
404  float fGapDensity = -999.0;
405 
406  //int nkals = sr->trk.kalman.ntracks;
407 
408 
409  if (sr->trk.kalman.tracks[idx].rempid == -1) return -2.0f;
410  //if (idx==GetBestTrack()(sr)) return -2.0f;
411 
412  scatll = sr->trk.kalman.tracks[idx].scatllh2017;
413  dedxll = sr->trk.kalman.tracks[idx].dedxllh2017;
414 
415  if (sr->trk.kalman.tracks[idx].avedEdxlast10cm == -5) return -2.0f;
416  if (sr->trk.kalman.tracks[idx].avedEdxlast30cm == -5) return -2.0f;
417  if (sr->trk.kalman.tracks[idx].avedEdxlast20cm == -5) return -2.0f;
418  if (sr->trk.kalman.tracks[idx].avedEdxlast40cm == -5) return -2.0f;
419 
420  if (sr->trk.kalman.tracks[idx].avedEdxlast10cm > 30){
421  fAvededxlast10cm = 30;
422  }else{
423  fAvededxlast10cm = sr->trk.kalman.tracks[idx].avedEdxlast10cm;
424  }
425 
426  if (sr->trk.kalman.tracks[idx].avedEdxlast40cm > 30){
427  fAvededxlast40cm = 30;
428  }else{
429  fAvededxlast40cm = sr->trk.kalman.tracks[idx].avedEdxlast40cm;
430  }
431 
432  TVector3 vtxpos(sr->trk.kalman.tracks[sr->trk.kalman.idxremid].start.X(),sr->trk.kalman.tracks[sr->trk.kalman.idxremid].start.Y(),sr->trk.kalman.tracks[sr->trk.kalman.idxremid].start.Z());
433  TVector3 trkpos(sr->trk.kalman.tracks[idx].start.X(),sr->trk.kalman.tracks[idx].start.Y(),sr->trk.kalman.tracks[idx].start.Z());
434  float trackstartgap = (vtxpos-trkpos).Mag();
435  TVector3 trkvec = sr->trk.kalman.tracks[idx].dir;
436  float angle = trkvec.Angle(beamDirND);
437  float length = sr->trk.kalman.tracks[idx].len;
438  float coslength = length*cos(angle);
439  float ngaps = sr->trk.kalman.tracks[idx].nplanegap;
440 
441  fTrackStartGap=trackstartgap;
442  fGapDensity=ngaps/coslength;*/
443 
444  TMVAvars[0] = dedxll;
445  TMVAvars[1] = scatll;
446  TMVAvars[2] = avededxlast10cm;
447  TMVAvars[3] = avededxlast40cm;
448  TMVAvars[4] = trackstartgap;
449  TMVAvars[5] = gapdensity;
450 
451  if(!fReaderBDTProton) InitTMVAProton();
452 
453  protonbdtvalue = fReaderBDTProton->EvaluateMVA("BDTGProton");
454 
455  return protonbdtvalue;
456  }
457 
458  float GetGammaIDVal::operator() ( float dedxll, float scatll, float avededxlast10cm, float avededxlast40cm, float trackstartgap, float gapdensity) const
459  {
460 
461  float gammabdtvalue = -999.0;
462 
463  if(!fReaderBDTGamma) InitTMVAGamma();
464  /*
465  float dedxll = -999.0, scatll = -999.0;
466  float fAvededxlast10cm = -999.0;
467  float fAvededxlast40cm = -999.0;
468  float fTrackStartGap = -999.0;
469  float fGapDensity = -999.0;
470 
471  //int nkals = sr->trk.kalman.ntracks;
472 
473 
474  if (sr->trk.kalman.tracks[idx].rempid == -1) return -2.0f;
475  //if (idx==GetBestTrack()(sr)) return -2.0f;
476 
477  scatll = sr->trk.kalman.tracks[idx].scatllh2017;
478  dedxll = sr->trk.kalman.tracks[idx].dedxllh2017;
479 
480  if (sr->trk.kalman.tracks[idx].avedEdxlast10cm == -5) return -2.0f;
481  if (sr->trk.kalman.tracks[idx].avedEdxlast30cm == -5) return -2.0f;
482  if (sr->trk.kalman.tracks[idx].avedEdxlast20cm == -5) return -2.0f;
483  if (sr->trk.kalman.tracks[idx].avedEdxlast40cm == -5) return -2.0f;
484 
485  if (sr->trk.kalman.tracks[idx].avedEdxlast10cm > 30){
486  fAvededxlast10cm = 30;
487  }else{
488  fAvededxlast10cm = sr->trk.kalman.tracks[idx].avedEdxlast10cm;
489  }
490 
491  if (sr->trk.kalman.tracks[idx].avedEdxlast40cm > 30){
492  fAvededxlast40cm = 30;
493  }else{
494  fAvededxlast40cm = sr->trk.kalman.tracks[idx].avedEdxlast40cm;
495  }
496 
497  TVector3 vtxpos(sr->trk.kalman.tracks[sr->trk.kalman.idxremid].start.X(),sr->trk.kalman.tracks[sr->trk.kalman.idxremid].start.Y(),sr->trk.kalman.tracks[sr->trk.kalman.idxremid].start.Z());
498  TVector3 trkpos(sr->trk.kalman.tracks[idx].start.X(),sr->trk.kalman.tracks[idx].start.Y(),sr->trk.kalman.tracks[idx].start.Z());
499  float trackstartgap = (vtxpos-trkpos).Mag();
500  TVector3 trkvec = sr->trk.kalman.tracks[idx].dir;
501  float angle = trkvec.Angle(beamDirND);
502  float length = sr->trk.kalman.tracks[idx].len;
503  float coslength = length*cos(angle);
504  float ngaps = sr->trk.kalman.tracks[idx].nplanegap;
505 
506  fTrackStartGap=trackstartgap;
507  fGapDensity=ngaps/coslength;*/
508 
509  TMVAvars[0] = dedxll;
510  TMVAvars[1] = scatll;
511  TMVAvars[2] = avededxlast10cm;
512  TMVAvars[3] = avededxlast40cm;
513  TMVAvars[4] = trackstartgap;
514  TMVAvars[5] = gapdensity;
515 
516  if(!fReaderBDTGamma) InitTMVAGamma();
517 
518  gammabdtvalue = fReaderBDTGamma->EvaluateMVA("BDTGGamma");
519 
520  return gammabdtvalue;
521  }
522 
523  //---------------------------------------------------
524  float GetPionID::InitTMVAPion() const
525  {
526  if(!fReaderBDTPion) fReaderBDTPion = new TMVA::Reader( "!Color:!Silent" );
527  const char* path = getenv("SRT_PUBLIC_CONTEXT");
528  std::string pidlib = std::string(path) + "/NDAna/numucc_ccPi/PionFeatureID_BDT.weights.xml";
529  fReaderBDTPion->AddVariable("TrackStartGap", &TMVAvars[4]);
530  fReaderBDTPion->AddVariable("nGaps/ProjLength", &TMVAvars[5]);
531  fReaderBDTPion->AddVariable("ScatLL", &TMVAvars[1]);
532  fReaderBDTPion->AddVariable("DedxLL", &TMVAvars[0]);
533  fReaderBDTPion->AddVariable("Avededxlast10cm", &TMVAvars[2]);
534  fReaderBDTPion->AddVariable("Avededxlast40cm", &TMVAvars[3]);
535  fReaderBDTPion->BookMVA("BDTGPion", pidlib);
536  return -5.0f;
537  }
538 
539  float GetPionID::InitTMVAProton() const
540  {
541  const char* path = getenv("SRT_PUBLIC_CONTEXT");
542  if(!fReaderBDTProton) fReaderBDTProton = new TMVA::Reader( "!Color:!Silent" );
543  std::string pidlib1 = std::string(path) + "/NDAna/numucc_ccPi/ProtonFeatureID_BDT.weights.xml";
544  fReaderBDTProton->AddVariable("TrackStartGap", &TMVAvars[4]);
545  fReaderBDTProton->AddVariable("ScatLL", &TMVAvars[1]);
546  fReaderBDTProton->AddVariable("DedxLL", &TMVAvars[0]);
547  fReaderBDTProton->AddVariable("Avededxlast10cm", &TMVAvars[2]);
548  fReaderBDTProton->AddVariable("Avededxlast40cm", &TMVAvars[3]);
549  fReaderBDTProton->AddVariable("nGaps/ProjLength", &TMVAvars[5]);
550  fReaderBDTProton->BookMVA("BDTGProton", pidlib1);
551  return -5.0f;
552  }
553 
554  float GetPionID::InitTMVAGamma() const
555  {
556  const char* path = getenv("SRT_PUBLIC_CONTEXT");
557  if(!fReaderBDTGamma) fReaderBDTGamma = new TMVA::Reader( "!Color:!Silent" );
558  std::string pidlib2 = std::string(path) + "/NDAna/numucc_ccPi/GammaFeatureID_BDT.weights.xml";
559  fReaderBDTGamma->AddVariable("TrackStartGap", &TMVAvars[4]);
560  fReaderBDTGamma->AddVariable("ScatLL", &TMVAvars[1]);
561  fReaderBDTGamma->AddVariable("DedxLL", &TMVAvars[0]);
562  fReaderBDTGamma->AddVariable("Avededxlast10cm", &TMVAvars[2]);
563  fReaderBDTGamma->AddVariable("Avededxlast40cm", &TMVAvars[3]);
564  fReaderBDTGamma->AddVariable("nGaps/ProjLength", &TMVAvars[5]);
565  fReaderBDTGamma->BookMVA("BDTGGamma", pidlib2);
566  return -5.0f;
567  }
568 
569  float GetPionTrack::InitTMVAPion() const
570  {
571  if(!fReaderBDTPion) fReaderBDTPion = new TMVA::Reader( "!Color:!Silent" );
572  const char* path = getenv("SRT_PUBLIC_CONTEXT");
573  std::string pidlib = std::string(path) + "/NDAna/numucc_ccPi/PionFeatureID_BDT.weights.xml";
574  fReaderBDTPion->AddVariable("TrackStartGap", &TMVAvars[4]);
575  fReaderBDTPion->AddVariable("nGaps/ProjLength", &TMVAvars[5]);
576  fReaderBDTPion->AddVariable("ScatLL", &TMVAvars[1]);
577  fReaderBDTPion->AddVariable("DedxLL", &TMVAvars[0]);
578  fReaderBDTPion->AddVariable("Avededxlast10cm", &TMVAvars[2]);
579  fReaderBDTPion->AddVariable("Avededxlast40cm", &TMVAvars[3]);
580  fReaderBDTPion->BookMVA("BDTGPion", pidlib);
581  return -5.0f;
582  }
583 
584  float GetPionTrack::InitTMVAProton() const
585  {
586  const char* path = getenv("SRT_PUBLIC_CONTEXT");
587  if(!fReaderBDTProton) fReaderBDTProton = new TMVA::Reader( "!Color:!Silent" );
588  std::string pidlib1 = std::string(path) + "/NDAna/numucc_ccPi/ProtonFeatureID_BDT.weights.xml";
589  fReaderBDTProton->AddVariable("TrackStartGap", &TMVAvars[4]);
590  fReaderBDTProton->AddVariable("ScatLL", &TMVAvars[1]);
591  fReaderBDTProton->AddVariable("DedxLL", &TMVAvars[0]);
592  fReaderBDTProton->AddVariable("Avededxlast10cm", &TMVAvars[2]);
593  fReaderBDTProton->AddVariable("Avededxlast40cm", &TMVAvars[3]);
594  fReaderBDTProton->AddVariable("nGaps/ProjLength", &TMVAvars[5]);
595  fReaderBDTProton->BookMVA("BDTGProton", pidlib1);
596  return -5.0f;
597  }
598  float GetPionTrack::InitTMVAGamma() const
599  {
600  const char* path = getenv("SRT_PUBLIC_CONTEXT");
601  if(!fReaderBDTGamma) fReaderBDTGamma = new TMVA::Reader( "!Color:!Silent" );
602  std::string pidlib2 = std::string(path) + "/NDAna/numucc_ccPi/GammaFeatureID_BDT.weights.xml";
603  fReaderBDTGamma->AddVariable("TrackStartGap", &TMVAvars[4]);
604  fReaderBDTGamma->AddVariable("ScatLL", &TMVAvars[1]);
605  fReaderBDTGamma->AddVariable("DedxLL", &TMVAvars[0]);
606  fReaderBDTGamma->AddVariable("Avededxlast10cm", &TMVAvars[2]);
607  fReaderBDTGamma->AddVariable("Avededxlast40cm", &TMVAvars[3]);
608  fReaderBDTGamma->AddVariable("nGaps/ProjLength", &TMVAvars[5]);
609  fReaderBDTGamma->BookMVA("BDTGGamma", pidlib2);
610  return -5.0f;
611  }
612 
613  float GetBestPionID::InitTMVAPion() const
614  {
615  if(!fReaderBDTPion) fReaderBDTPion = new TMVA::Reader( "!Color:!Silent" );
616  const char* path = getenv("SRT_PUBLIC_CONTEXT");
617  std::string pidlib = std::string(path) + "/NDAna/numucc_ccPi/PionFeatureID_BDT.weights.xml";
618  fReaderBDTPion->AddVariable("TrackStartGap", &TMVAvars[4]);
619  fReaderBDTPion->AddVariable("nGaps/ProjLength", &TMVAvars[5]);
620  fReaderBDTPion->AddVariable("ScatLL", &TMVAvars[1]);
621  fReaderBDTPion->AddVariable("DedxLL", &TMVAvars[0]);
622  fReaderBDTPion->AddVariable("Avededxlast10cm", &TMVAvars[2]);
623  fReaderBDTPion->AddVariable("Avededxlast40cm", &TMVAvars[3]);
624  fReaderBDTPion->BookMVA("BDTGPion", pidlib);
625  return -5.0f;
626  }
627 
628  float GetBestPionTrack::InitTMVAPion() const
629  {
630  if(!fReaderBDTPion) fReaderBDTPion = new TMVA::Reader( "!Color:!Silent" );
631  const char* path = getenv("SRT_PUBLIC_CONTEXT");
632  std::string pidlib = std::string(path) + "/NDAna/numucc_ccPi/PionFeatureID_BDT.weights.xml";
633  fReaderBDTPion->AddVariable("TrackStartGap", &TMVAvars[4]);
634  fReaderBDTPion->AddVariable("nGaps/ProjLength", &TMVAvars[5]);
635  fReaderBDTPion->AddVariable("ScatLL", &TMVAvars[1]);
636  fReaderBDTPion->AddVariable("DedxLL", &TMVAvars[0]);
637  fReaderBDTPion->AddVariable("Avededxlast10cm", &TMVAvars[2]);
638  fReaderBDTPion->AddVariable("Avededxlast40cm", &TMVAvars[3]);
639  fReaderBDTPion->BookMVA("BDTGPion", pidlib);
640  return -5.0f;
641  }
642 
643  float GetPionIDVal::InitTMVAPion() const
644  {
645  if(!fReaderBDTPion) fReaderBDTPion = new TMVA::Reader( "!Color:!Silent" );
646  const char* path = getenv("SRT_PUBLIC_CONTEXT");
647  std::string pidlib = std::string(path) + "/NDAna/numucc_ccPi/PionFeatureID_BDT.weights.xml";
648  fReaderBDTPion->AddVariable("TrackStartGap", &TMVAvars[4]);
649  fReaderBDTPion->AddVariable("nGaps/ProjLength", &TMVAvars[5]);
650  fReaderBDTPion->AddVariable("ScatLL", &TMVAvars[1]);
651  fReaderBDTPion->AddVariable("DedxLL", &TMVAvars[0]);
652  fReaderBDTPion->AddVariable("Avededxlast10cm", &TMVAvars[2]);
653  fReaderBDTPion->AddVariable("Avededxlast40cm", &TMVAvars[3]);
654  fReaderBDTPion->BookMVA("BDTGPion", pidlib);
655  return -5.0f;
656  }
657 
658  float GetProtonIDVal::InitTMVAProton() const
659  {
660  const char* path = getenv("SRT_PUBLIC_CONTEXT");
661  if(!fReaderBDTProton) fReaderBDTProton = new TMVA::Reader( "!Color:!Silent" );
662  std::string pidlib1 = std::string(path) + "/NDAna/numucc_ccPi/ProtonFeatureID_BDT.weights.xml";
663  fReaderBDTProton->AddVariable("TrackStartGap", &TMVAvars[4]);
664  fReaderBDTProton->AddVariable("ScatLL", &TMVAvars[1]);
665  fReaderBDTProton->AddVariable("DedxLL", &TMVAvars[0]);
666  fReaderBDTProton->AddVariable("Avededxlast10cm", &TMVAvars[2]);
667  fReaderBDTProton->AddVariable("Avededxlast40cm", &TMVAvars[3]);
668  fReaderBDTProton->AddVariable("nGaps/ProjLength", &TMVAvars[5]);
669  fReaderBDTProton->BookMVA("BDTGProton", pidlib1);
670  return -5.0f;
671  }
672 
673  float GetGammaIDVal::InitTMVAGamma() const
674  {
675  const char* path = getenv("SRT_PUBLIC_CONTEXT");
676  if(!fReaderBDTGamma) fReaderBDTGamma = new TMVA::Reader( "!Color:!Silent" );
677  std::string pidlib2 = std::string(path) + "/NDAna/numucc_ccPi/GammaFeatureID_BDT.weights.xml";
678  fReaderBDTGamma->AddVariable("TrackStartGap", &TMVAvars[4]);
679  fReaderBDTGamma->AddVariable("ScatLL", &TMVAvars[1]);
680  fReaderBDTGamma->AddVariable("DedxLL", &TMVAvars[0]);
681  fReaderBDTGamma->AddVariable("Avededxlast10cm", &TMVAvars[2]);
682  fReaderBDTGamma->AddVariable("Avededxlast40cm", &TMVAvars[3]);
683  fReaderBDTGamma->AddVariable("nGaps/ProjLength", &TMVAvars[5]);
684  fReaderBDTGamma->BookMVA("BDTGGamma", pidlib2);
685  return -5.0f;
686  }
687 
688 } // name space ana
689 
caf::kNEARDET
Near Detector underground.
Definition: SREnums.h:10
ana::GetPionTrack::InitTMVAGamma
float InitTMVAGamma() const
Definition: NDXSecTrackPID.cxx:598
ana::fReaderBDTGamma
static TMVA::Reader * fReaderBDTGamma
Definition: NDXSecTrackPID.cxx:21
angle
Double_t angle
Definition: plot.C:86
ana
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
ana::GetPionID::InitTMVAPion
float InitTMVAPion() const
Definition: NDXSecTrackPID.cxx:524
caf::Proxy< caf::SRKalman >::idxremid
caf::Proxy< unsigned int > idxremid
Definition: SRProxy.h:1777
caf::Proxy< caf::SRKalman >::ntracks
caf::Proxy< size_t > ntracks
Definition: SRProxy.h:1778
ana::GetPionTrack::InitTMVAPion
float InitTMVAPion() const
Definition: NDXSecTrackPID.cxx:569
beamDirND
const TVector3 beamDirND
Definition: caf_numu_validation.C:93
caf::Proxy< caf::StandardRecord >
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
make_root_from_grid_output.scatll
scatll
Definition: make_root_from_grid_output.py:67
ana::fReaderBDTProton
static TMVA::Reader * fReaderBDTProton
Definition: NDXSecTrackPID.cxx:20
NDXSecMuonPID.h
ana::GetPionIDVal::operator()
float operator()(float dedxll, float scatll, float avededxlast10cm, float avededxlast40cm, float trackstartgap, float gapdensity) const
Definition: NDXSecTrackPID.cxx:329
ana::GetBestPionTrack::InitTMVAPion
float InitTMVAPion() const
Definition: NDXSecTrackPID.cxx:628
demo0.length
length
Definition: demo0.py:21
caf::Proxy< caf::StandardRecord >::trk
caf::Proxy< caf::SRTrackBranch > trk
Definition: SRProxy.h:2145
make_root_from_grid_output.dedxll
dedxll
Definition: make_root_from_grid_output.py:66
ana::GetPionTrack::InitTMVAProton
float InitTMVAProton() const
Definition: NDXSecTrackPID.cxx:584
cet::getenv
std::string getenv(std::string const &name)
ana::TMVAvars
static float TMVAvars[4]
Definition: NDNCPi0Xsec.cxx:10
ana::GetPionIDVal::InitTMVAPion
float InitTMVAPion() const
Definition: NDXSecTrackPID.cxx:643
sr
caf::StandardRecord * sr
Definition: FillTrainingTrees.h:18
ana::GetPionID::InitTMVAGamma
float InitTMVAGamma() const
Definition: NDXSecTrackPID.cxx:554
ana::GetPionTrack::operator()
float operator()(const caf::SRProxy *sr) const
Definition: NDXSecTrackPID.cxx:106
ana::fReaderBDTPion
static TMVA::Reader * fReaderBDTPion
Definition: NDXSecTrackPID.cxx:19
ana::GetPionID::operator()
float operator()(const caf::SRProxy *sr) const
Definition: NDXSecTrackPID.cxx:27
ana::beamDirND
const TVector3 beamDirND
Definition: NDXSecTrackPID.cxx:23
path
const std::string path
Definition: plot_BEN.C:43
ana::GetGammaIDVal::InitTMVAGamma
float InitTMVAGamma() const
Definition: NDXSecTrackPID.cxx:673
ana::GetBestPionID::InitTMVAPion
float InitTMVAPion() const
Definition: NDXSecTrackPID.cxx:613
caf::Proxy< caf::SRTrackBranch >::kalman
caf::Proxy< caf::SRKalman > kalman
Definition: SRProxy.h:1797
Mag
float Mag() const
Definition: SRLorentzVectorAddons.h:10
ana::GetBestPionTrack::operator()
float operator()(const caf::SRProxy *sr) const
Definition: NDXSecTrackPID.cxx:258
ana::GetPionID::InitTMVAProton
float InitTMVAProton() const
Definition: NDXSecTrackPID.cxx:539
std::cos
T cos(T number)
Definition: d0nt_math.hpp:78
ana::GetGammaIDVal::operator()
float operator()(float dedxll, float scatll, float avededxlast10cm, float avededxlast40cm, float trackstartgap, float gapdensity) const
Definition: NDXSecTrackPID.cxx:458
ana::NuMIBeamDirection
TVector3 NuMIBeamDirection(caf::Det_t det)
Average direction of NuMI neutrinos in a given detector This function is a copy of geo::GeometryBase:...
Definition: Utilities.cxx:329
ana::GetProtonIDVal::InitTMVAProton
float InitTMVAProton() const
Definition: NDXSecTrackPID.cxx:658
caf::Proxy< caf::SRKalman >::tracks
caf::Proxy< std::vector< caf::SRKalmanTrack > > tracks
Definition: SRProxy.h:1780
ana::GetProtonIDVal::operator()
float operator()(float dedxll, float scatll, float avededxlast10cm, float avededxlast40cm, float trackstartgap, float gapdensity) const
Definition: NDXSecTrackPID.cxx:393
ana::GetBestPionID::operator()
float operator()(const caf::SRProxy *sr) const
Definition: NDXSecTrackPID.cxx:187
string
enum BeamMode string
Generated on Sun Mar 7 2021 14:13:09 for NOvA by   doxygen 1.8.11