NDXSecMuonPID.cxx
Go to the documentation of this file.
2 
4 
5 #include "TFile.h"
6 
7 #include "TMVA/Reader.h"
8 
9 #include <cassert>
10 #include <math.h>
11 
12 namespace ana
13 {
14  static float TMVAvars[4] = {-999.0};
15  static TMVA::Reader* fReaderBDT = 0;
16 
17 
18  //-------------------------------------------------------------
20  {
21 
22  float heighstpid = -999.0;
23  float bdtvalue = -999.0;
24 
25  if(!fReaderBDT) InitTMVA();
26  float dedxll = -999.0, scatll = -999.0;
27  float fAvededxlast10cm = -999.0;
28  float fAvededxlast40cm = -999.0;
29 
30  int nkals = sr->trk.kalman.ntracks;
31 
32  for (int itrk = 0; itrk < nkals; ++itrk){
33  if (sr->trk.kalman.tracks[itrk].rempid == -1) continue;
34 
35  scatll = sr->trk.kalman.tracks[itrk].scatllh2017;
36  dedxll = sr->trk.kalman.tracks[itrk].dedxllh2017;
37 
38  if (sr->trk.kalman.tracks[itrk].avedEdxlast10cm == -5) continue;
39  if (sr->trk.kalman.tracks[itrk].avedEdxlast30cm == -5) continue;
40  if (sr->trk.kalman.tracks[itrk].avedEdxlast20cm == -5) continue;
41  if (sr->trk.kalman.tracks[itrk].avedEdxlast40cm == -5) continue;
42 
43  if (sr->trk.kalman.tracks[itrk].avedEdxlast10cm > 30){
44  fAvededxlast10cm = 30;
45  }else{
46  fAvededxlast10cm = sr->trk.kalman.tracks[itrk].avedEdxlast10cm;
47  }
48 
49  if (sr->trk.kalman.tracks[itrk].avedEdxlast40cm > 30){
50  fAvededxlast40cm = 30;
51  }else{
52  fAvededxlast40cm = sr->trk.kalman.tracks[itrk].avedEdxlast40cm;
53  }
54 
55  TMVAvars[0] = dedxll;
56  TMVAvars[1] = scatll;
57  TMVAvars[2] = fAvededxlast10cm;
58  TMVAvars[3] = fAvededxlast40cm;
59 
60  if(!fReaderBDT) InitTMVA();
61 
62  bdtvalue = fReaderBDT->EvaluateMVA("BDTG_1k");
63 
64  if (bdtvalue > heighstpid){
65  heighstpid = bdtvalue;
66  }else continue;
67 
68  }
69  return heighstpid;
70  }
71 
72 
74  {
75 
76  float heighstpid = -999.0;
77  int imax = -1; //this can be used as your track index with highest PID score.
78  float bdtvalue = -999.0;
79 
80  if(!fReaderBDT) InitTMVA();
81  float dedxll = -999.0, scatll = -999.0;
82  float fAvededxlast10cm = -999.0;
83  float fAvededxlast40cm = -999.0;
84 
85  int nkals = sr->trk.kalman.ntracks;
86 
87  for (int itrk = 0; itrk < nkals; ++itrk){
88  if (sr->trk.kalman.tracks[itrk].rempid == -1) continue;
89 
90  scatll = sr->trk.kalman.tracks[itrk].scatllh2017;
91  dedxll = sr->trk.kalman.tracks[itrk].dedxllh2017;
92 
93  if (sr->trk.kalman.tracks[itrk].avedEdxlast10cm == -5) continue;
94  if (sr->trk.kalman.tracks[itrk].avedEdxlast30cm == -5) continue;
95  if (sr->trk.kalman.tracks[itrk].avedEdxlast20cm == -5) continue;
96  if (sr->trk.kalman.tracks[itrk].avedEdxlast40cm == -5) continue;
97 
98  if (sr->trk.kalman.tracks[itrk].avedEdxlast10cm > 30){
99  fAvededxlast10cm = 30;
100  }else{
101  fAvededxlast10cm = sr->trk.kalman.tracks[itrk].avedEdxlast10cm;
102  }
103 
104  if (sr->trk.kalman.tracks[itrk].avedEdxlast40cm > 30){
105  fAvededxlast40cm = 30;
106  }else{
107  fAvededxlast40cm = sr->trk.kalman.tracks[itrk].avedEdxlast40cm;
108  }
109 
110  TMVAvars[0] = dedxll;
111  TMVAvars[1] = scatll;
112  TMVAvars[2] = fAvededxlast10cm;
113  TMVAvars[3] = fAvededxlast40cm;
114 
115  if(!fReaderBDT) InitTMVA();
116 
117  bdtvalue = fReaderBDT->EvaluateMVA("BDTG_1k");
118  if (bdtvalue > heighstpid){
119  imax = itrk;
120  heighstpid = bdtvalue;
121  }else continue;
122  }
123  return imax;
124  }
125 
126  //---------------------------------------------------
127  float GetMuonID::InitTMVA() const
128  {
129  if(!fReaderBDT) fReaderBDT = new TMVA::Reader( "!Color:!Silent" );
130  const char* path = getenv("SRT_PUBLIC_CONTEXT");
131  std::string pidlib = std::string(path) + "/NDAna/numucc_inc/MuonID_BDTG.weights.xml";
132 
133  fReaderBDT->AddVariable("DedxLL", &TMVAvars[0]);
134  fReaderBDT->AddVariable("ScatLL", &TMVAvars[1]);
135  fReaderBDT->AddVariable("Avededxlast10cm", &TMVAvars[2]);
136  fReaderBDT->AddVariable("Avededxlast40cm", &TMVAvars[3]);
137  fReaderBDT->BookMVA("BDTG_1k", pidlib);
138  return -5.0f;
139  }
140 
142  {
143  if(!fReaderBDT) fReaderBDT = new TMVA::Reader( "!Color:!Silent" );
144  const char* path = getenv("SRT_PUBLIC_CONTEXT");
145  std::string pidlib = std::string(path) + "/NDAna/numucc_inc/MuonID_BDTG.weights.xml";
146 
147  fReaderBDT->AddVariable("DedxLL", &TMVAvars[0]);
148  fReaderBDT->AddVariable("ScatLL", &TMVAvars[1]);
149  fReaderBDT->AddVariable("Avededxlast10cm", &TMVAvars[2]);
150  fReaderBDT->AddVariable("Avededxlast40cm", &TMVAvars[3]);
151  fReaderBDT->BookMVA("BDTG_1k", pidlib);
152  return -5.0f;
153  }
154 } // name space ana
155 
static TMVA::Reader * fReaderBDT
Definition: NDNCPi0Xsec.cxx:11
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
float InitTMVA() const
caf::Proxy< size_t > ntracks
Definition: SRProxy.h:1777
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2125
float InitTMVA() const
caf::Proxy< caf::SRTrackBranch > trk
Definition: SRProxy.h:2144
std::string getenv(std::string const &name)
static float TMVAvars[4]
Definition: NDNCPi0Xsec.cxx:10
caf::StandardRecord * sr
float operator()(const caf::SRProxy *sr) const
const std::string path
Definition: plot_BEN.C:43
caf::Proxy< caf::SRKalman > kalman
Definition: SRProxy.h:1796
float operator()(const caf::SRProxy *sr) const
caf::Proxy< std::vector< caf::SRKalmanTrack > > tracks
Definition: SRProxy.h:1779
enum BeamMode string