CLAApplication.C
Go to the documentation of this file.
1 /**********************************************************************************
2  * Project : TMVA - a Root-integrated toolkit for multivariate data analysis *
3  * Package : TMVA *
4  * Root Macro: TMVAMulticlassApplication *
5  * *
6  * This macro provides a simple example on how to use the trained multiclass *
7  * classifiers within an analysis module *
8  **********************************************************************************/
9 
10 #include <cstdlib>
11 #include <iostream>
12 #include <map>
13 #include <string>
14 #include <vector>
15 
16 #include "TFile.h"
17 #include "TTree.h"
18 #include "TString.h"
19 #include "TSystem.h"
20 #include "TROOT.h"
21 #include "TStopwatch.h"
22 #include "TH1F.h"
23 
24 #if not defined(__CINT__) || defined(__MAKECINT__)
25 #include "TMVA/Tools.h"
26 #include "TMVA/Reader.h"
27 #endif
28 
29 using namespace TMVA;
30 
31 void CLAApplication( TString myMethodList = "" )
32 {
33 #ifdef __CINT__
34  gROOT->ProcessLine( ".O0" ); // turn off optimization in CINT
35 #endif
36 
37  TMVA::Tools::Instance();
38 
39  //---------------------------------------------------------------
40  // default MVA methods to be trained + tested
41  std::map<std::string,int> Use;
42 
43  Use["BDTG"] = 0;
44  Use["BDTA"] = 1;
45  Use["BDTB"] = 0;
46  Use["BDTD"] = 0;
47  Use["BDTMitFisher"] = 0;
48  Use["MLP"] = 0;
49  Use["MLPBFGS"] = 0;
50  //---------------------------------------------------------------
51 
53  std::cout << "==> Start CLAApplication" << std::endl;
54  if (myMethodList != "") {
55  for (std::map<std::string,int>::iterator it = Use.begin(); it != Use.end(); it++) it->second = 0;
56 
57  std::vector<TString> mlist = gTools().SplitString( myMethodList, ',' );
58  for (UInt_t i=0; i<mlist.size(); i++) {
59  std::string regMethod(mlist[i]);
60  if (Use.find(regMethod) == Use.end()) {
61  std::cout << "Method \"" << regMethod << "\" not known in TMVA under this name. Choose among the following:" << std::endl;
62  for (std::map<std::string,int>::iterator it = Use.begin(); it != Use.end(); it++) std::cout << it->first << " " << std::endl;
63  std::cout << std::endl;
64  return;
65  }
66  Use[regMethod] = 1;
67  }
68  }
69 
70  // create the Reader object
71  //TMVA::Reader *reader = new TMVA::Reader( "!Color:!Silent" );
72  TMVA::Reader* reader = new TMVA::Reader( "" );
73 
74  // create a set of variables and declare them to the reader
75  // the variable names must corresponds in name and type to
76  // those given in the weight file(s) that you use
77  //Float_t anglekal, dirFY, boxmaxFY, kalnhit, kallen, kalfwdcell, kalbakcell, scatt, nhit, energy, boxminFY, nkal, ncid, nueid, numuid, nutauid;
79  reader->AddVariable( "ncid", &ncid );
80  reader->AddVariable( "partptp", &partptp );
81  reader->AddVariable( "shwnhit", &shwnhit );
82  reader->AddVariable( "shwxminusy", &shwxminusy );
83  reader->AddVariable( "shwxplusy", &shwxplusy );
84  reader->AddVariable( "shwxovery", &shwxovery );
85  reader->AddVariable( "shwcalE", &shwcalE );
86  reader->AddVariable( "shwdirY", &shwdirY );
87  reader->AddVariable( "shwlen", &shwlen );
88  reader->AddVariable( "shwwwidth", &shwwwidth );
89  reader->AddVariable( "shwGap", &shwGap );
90  reader->AddVariable( "nshwlid", &nshwlid );
91  reader->AddVariable( "nmiphit", &nmiphit );
92  //reader->AddVariable( "nueid", &nueid );
93  //reader->AddVariable( "numuid", &numuid );
94  //reader->AddVariable( "nutauid", &nutauid );
95  //reader->AddVariable( "shwnplane", &shwnplane );
96  //reader->AddVariable( "shwmaxplanecont", &shwmaxplanecont );
97  //reader->AddVariable( "shwmaxplanegap", &shwmaxplanegap );
98  //reader->AddVariable( "shwstartX", &shwstartX );
99  //reader->AddVariable( "shwstartY", &shwstartY );
100  //reader->AddVariable( "shwstartZ", &shwstartZ );
101  //reader->AddVariable( "shwstopX", &shwstopX );
102  //reader->AddVariable( "shwstopY", &shwstopY );
103  //reader->AddVariable( "shwstopZ", &shwstopZ );
104  //reader->AddVariable( "vtxX", &vtxX );
105  //reader->AddVariable( "vtxY", &vtxY );
106  //reader->AddVariable( "vtxZ", &vtxZ );
107  //reader->AddVariable( "cosdang", &cosdang );
108  //reader->AddVariable( "sparsenessasymm", &sparsenessasymm );
109  //reader->AddVariable( "hitsperplaneasymm", &hitsperplaneasymm );
110  //reader->AddVariable( "energy", &energy );
111  //reader->AddVariable( "xAll", &xAll );
112  //reader->AddVariable( "ehit", &ehit );
113  /*
114  theTree->SetBranchAddress( "cosdang", &cosdang);
115  theTree->SetBranchAddress( "sparsenessasymm", &sparsenessasymm);
116  theTree->SetBranchAddress( "hitsperplaneasymm", &hitsperplaneasymm);
117  theTree->SetBranchAddress( "partptp", &partptp);
118  theTree->SetBranchAddress( "nhit", &nhit);
119  theTree->SetBranchAddress( "energy", &energy);
120  reader->AddVariable( "anglekal", &anglekal );
121  reader->AddVariable( "dirFY", &dirFY );
122  reader->AddVariable( "boxmaxFY", &boxmaxFY );
123  reader->AddVariable( "kalnhit", &kalnhit );
124  reader->AddVariable( "kallen", &kallen );
125  reader->AddVariable( "kalfwdcell",&kalfwdcell);
126  reader->AddVariable( "kalbakcell",&kalbakcell );
127  reader->AddVariable( "scatt", &scatt );
128  reader->AddVariable( "nhit", &nhit );
129  reader->AddVariable( "energy", &energy );
130  reader->AddVariable( "boxminFY", &boxminFY );
131  reader->AddVariable( "nkal", &nkal );
132  */
133  // book the MVA methods
134  TString dir = "weights/";
135  TString prefix = "CLAShower";
136 
137  for (std::map<std::string,int>::iterator it = Use.begin(); it != Use.end(); it++) {
138  if (it->second) {
139  TString methodName = TString(it->first) + TString(" method");
140  TString weightfile = dir + prefix + TString("_") + TString(it->first) + TString(".weights.xml");
141  reader->BookMVA( methodName, weightfile );
142  }
143  }
144 
145  // book output histograms
146  UInt_t nbin = 100;
147  TH1F *BDTA_NC;
148  TH1F *BDTA_CS;
149  //TH1F *histMLP_signal(0), *histBDTG_signal(0), *histFDAGA_signal(0), *histPDEFoam_signal(0);
150  if (Use["BDTA"])
151  //histBDTA_NC = new TH1F( "MVA_BDTA_NC", "MVA_BDT_NC", nbin, 0., 1.1 );
152  BDTA_NC = new TH1F( "MVA_BDTA_NC", "MVA_BDT_NC", nbin, 0.0, 1.1 );
153  BDTA_CS = new TH1F( "MVA_BDTA_CS", "MVA_BDT_CS", nbin, 0.0, 1.1 );
154  /*
155  if (Use["MLP"])
156  histMLP_signal = new TH1F( "MVA_MLP_signal", "MVA_MLP_signal", nbin, 0., 1.1 );
157  if (Use["BDTG"])
158  histBDTG_signal = new TH1F( "MVA_BDTG_signal", "MVA_BDTG_signal", nbin, 0., 1.1 );
159  if (Use["MLPBFGS"])
160  histMLPBFGS_signal = new TH1F( "MVA_MLPBFGS_signal", "MVA_MLPBFGS_signal", nbin, 0., 1.1 );
161  if (Use["BDTB"])
162  histBDTB_signal = new TH1F( "MVA_BDTB_signal", "MVA_BDTB_signal", nbin, 0., 1.1 );
163  if (Use["BDTD"])
164  histBDTD_signal = new TH1F( "MVA_BDTD_signal", "MVA_BDTD_signal", nbin, 0., 1.1 );
165  if (Use["BDTMitFisher"])
166  histBDTMitFisher_signal = new TH1F( "MVA_BDTMitFisher_signal", "MVA_BDTMitFisher_signal", nbin, 0., 1.1 );
167  */
168 
169  TFile *input(0);
170  //TString fname = "TEST_Track.root";
171  TString fname = "APP_Track.root";
172 
173  if (!gSystem->AccessPathName( fname )) {
174  input = TFile::Open( fname ); // check if file in local directory exists
175  }
176  if (!input) {
177  std::cout << "ERROR: could not open data file, please generate example data first!" << std::endl;
178  exit(1);
179  }
180  std::cout << "--- TMVAMulticlassApp : Using input file: " << input->GetName() << std::endl;
181 
182  // prepare the tree
183  // - here the variable names have to corresponds to your tree
184  // - you can use the same variables as above which is slightly faster,
185  // but of course you can use different ones and copy the values inside the event loop
186 
187  TTree* theTree = (TTree*)input->Get("ncTree");
188  //TTree* theTree = (TTree*)input->Get("csTree");
189  std::cout << "--- Select signal sample" << std::endl;
190  theTree->SetBranchAddress( "ncid", &ncid);
191  theTree->SetBranchAddress( "partptp", &partptp);
192  theTree->SetBranchAddress( "shwnhit", &shwnhit);
193  theTree->SetBranchAddress( "shwxminusy", &shwxminusy);
194  theTree->SetBranchAddress( "shwxplusy", &shwxplusy);
195  theTree->SetBranchAddress( "shwxovery", &shwxovery);
196  theTree->SetBranchAddress( "shwcalE", &shwcalE);
197  theTree->SetBranchAddress( "shwdirY", &shwdirY);
198  theTree->SetBranchAddress( "shwlen", &shwlen);
199  theTree->SetBranchAddress( "shwwwidth", &shwwwidth);
200  theTree->SetBranchAddress( "shwGap", &shwGap);
201  theTree->SetBranchAddress( "nshwlid", &nshwlid);
202  theTree->SetBranchAddress( "nmiphit", &nmiphit);
203 
204  //The following variables are from Kirk's BDT
205  //theTree->SetBranchAddress( "anglekal", &anglekal );
206  //theTree->SetBranchAddress( "dirFY", &dirFY );
207  //theTree->SetBranchAddress( "boxmaxFY", &boxmaxFY );
208  //theTree->SetBranchAddress( "kalnhit", &kalnhit);
209  //theTree->SetBranchAddress( "kallen", &kallen);
210  //theTree->SetBranchAddress( "kalfwdcell",&kalfwdcell);
211  //theTree->SetBranchAddress( "kalbakcell",&kalbakcell);
212  //theTree->SetBranchAddress( "scatt", &scatt);
213  //theTree->SetBranchAddress( "boxminFY", &boxminFY);
214  //theTree->SetBranchAddress( "nkal", &nkal);
215  //theTree->SetBranchAddress( "nueid", &nueid);
216  //theTree->SetBranchAddress( "numuid", &numuid);
217  //theTree->SetBranchAddress( "nutauid", &nutauid);
218  //theTree->SetBranchAddress( "cosdang", &cosdang);
219  //theTree->SetBranchAddress( "sparsenessasymm", &sparsenessasymm);
220  //theTree->SetBranchAddress( "hitsperplaneasymm", &hitsperplaneasymm);
221  //theTree->SetBranchAddress( "energy", &energy);
222  //theTree->SetBranchAddress( "xAll := shwnhitx/shwnhit", &xAll);
223  //theTree->SetBranchAddress( "ehit := shwcalE/shwnhit",&ehit);
224  //theTree->SetBranchAddress( "shwcalE", &shwcalE );
225  //theTree->SetBranchAddress( "shwnplane", &shwnplane);
226  //theTree->SetBranchAddress( "shwmaxplanecont", &shwmaxplanecont);
227  //theTree->SetBranchAddress( "shwmaxplanegap", &shwmaxplanegap);
228  //theTree->SetBranchAddress( "shwstartX", &shwstartX);
229  //theTree->SetBranchAddress( "shwstartY", &shwstartY);
230  //theTree->SetBranchAddress( "shwstartZ", &shwstartZ);
231  //theTree->SetBranchAddress( "shwstopX", &shwstopX);
232  //theTree->SetBranchAddress( "shwstopY", &shwstopY);
233  //theTree->SetBranchAddress( "shwstopZ", &shwstopZ);
234  //theTree->SetBranchAddress( "vtxX", &vtxX);
235  //theTree->SetBranchAddress( "vtxY", &vtxY);
236  //theTree->SetBranchAddress( "vtxZ", &vtxZ);
237 
238  std::cout << "--- Processing: " << theTree->GetEntries() << "events" << std::endl;
239  TStopwatch sw;
240  sw.Start();
241 
242  for (Long64_t ievt=0; ievt<theTree->GetEntries();ievt++) {
243  if (ievt%1000 == 0){
244  std::cout << "--- ... Processing event: " << ievt << std::endl;
245  }
246  //if (ievt == 3600) break;
247  if (ievt == 133400) break;
248  theTree->GetEntry(ievt);
249 
250  if (Use["BDTA"])
251  //histBDTA_NC->Fill((reader->EvaluateMulticlass( "BDTA method" ))[0]);
252  BDTA_NC->Fill((reader->EvaluateMVA( "BDTA method" ))[0]);
253  //BDTA_CS->Fill((reader->EvaluateMVA( "BDTA method" ))[0]);
254  //histBDTA_NC-Fill(reader->EvaluateMVA("BDTA"));
255  /* if (Use["MLP"])
256  histMLP_signal->Fill((reader->EvaluateMulticlass( "MLP method" ))[0]);
257  if (Use["BDTG"])
258  histBDTG_signal->Fill((reader->EvaluateMulticlass( "BDTG method" ))[0]);
259  if (Use["MLPBFGS"])
260  histMLPBFGS_signal->Fill((reader->EvaluateMulticlass( "MLPBFGS method" ))[0]);
261  if (Use["BDTB"])
262  histBDTB_signal->Fill((reader->EvaluateMulticlass( "BDTB method" ))[0]);
263  if (Use["BDTD"])
264  histBDTD_signal->Fill((reader->EvaluateMulticlass( "BDTD method" ))[0]);
265  if (Use["BDTMitFisher"])
266  histBDTMitFisher_signal->Fill((reader->EvaluateMulticlass( "BDTMitFisher method" ))[0]);*/
267  }
268 
269  // get elapsed time
270  sw.Stop();
271  std::cout << "--- End of event loop: "; sw.Print();
272 
273  TFile *target = new TFile( "RealBDTNC.root","RECREATE" );
274  //TFile *target = new TFile( "RealBDTCS.root","RECREATE" );
275  if (Use["BDTA"])
276  BDTA_NC->Write();
277  //BDTA_CS->Write();
278 
279  /* if (Use["MLP"])
280  histMLP_signal->Write();
281  if (Use["BDTG"])
282  histBDTG_signal->Write();
283  if (Use["MLPBFGS"])
284  histMLPBFGS_signal->Write();
285  if (Use["BDTB"])
286  histBDTB_signal->Write();
287  if (Use["BDTD"])
288  histBDTD_signal->Write();
289  if (Use["BDTMitFisher"])
290  histBDTMitFisher_signal->Write();*/
291 
292  target->Close();
293  std::cout << "--- Created root file: \"RealBDTNC.root\" containing the MVA output histograms" << std::endl;
294 
295  delete reader;
296 
297  std::cout << "==> TMVAApplication is done!" << std::endl << std::endl;
298 }
float nshwlid
const XML_Char * target
Definition: expat.h:268
set< int >::iterator it
float shwGap
float shwxplusy
float shwcalE
float shwlen
float nmiphit
float partptp
void CLAApplication(TString myMethodList="")
const XML_Char * prefix
Definition: expat.h:380
float shwdirY
OStream cout
Definition: OStream.cxx:6
float ncid
Float_t sw
Definition: plot.C:20
float shwxovery
TDirectory * dir
Definition: macro.C:5
float shwxminusy
exit(0)
Definition: tmvaglob.h:28
float shwnhit
enum BeamMode string