TMVAEIDTraining.cxx
Go to the documentation of this file.
1 
2 /** TMVA-Based Electron ID classification training and diagnostics.
3  *
4  * This class is used to train a specified classifier to recognize charged current nue events and their
5  * backgrounds. This code is based on JMShower EID and is intended to take advantage of the diagnostics and
6  * ability to use different classifiers in TMVA.
7  */
8 #include <TString.h>
9 #include <vector>
10 
11 #include "TMVAEIDTraining.h"
12 
13 
14 #if not defined(__CINT__) || defined(__MAKECINT__)
15 // needs to be included when makecint runs (ACLIC)
16 #endif
17 
18 TMVAEIDTraining::TMVAEIDTraining(TString jobNameIn, TString outPathIn,
19  TString myMethodListIn, TString useENuIn) {
20 
21  //---------------------------------------------------------------
22  // This loads the library
23  TMVA::Tools::Instance();
24  // systemPtr = new TUnixSystem();
25  fOutPath = outPathIn;
26  methodsBooked = false;
27  TMVAEIDTraining::fJobName = jobNameIn;
28  TMVAEIDTraining::myMethodList = myMethodListIn;
30  if (useENuIn == "useENu") {
32  std::cout << "###Using ENu in training" << std::endl;
33  }
34  if (myMethodList == "") {
35  myMethodList = "TMlpANN";
36  }
37  TMVAEIDTraining::weightDir = fOutPath + TString("/output/") + fJobName + TString("/weights/");
38  TMVAEIDTraining::plotDir = fOutPath + TString("/output/") + fJobName + TString("/plots/");
39  TMVAEIDTraining::histDir = fOutPath + TString("/output/") + fJobName + TString("/hist/");
40 
41  gSystem->mkdir(weightDir.Data(), true);
42  gSystem->mkdir(plotDir.Data(), true);
43  gSystem->mkdir(histDir.Data(), true);
44 
45  TMVAEIDTraining::outfileName = histDir + TString("TMVA") + TString(".root");
46 
47  //
48  // Set TMVA options
49  //
50  (TMVA::gConfig().GetIONames()).fWeightFileDir = weightDir;
51 
52  // to get access to the GUI and all tmva macros
53  TString tmva_dir(TString(gRootDir) + "/tmva");
54  if (gSystem->Getenv("TMVASYS"))
55  tmva_dir = TString(gSystem->Getenv("TMVASYS"));
56  std::cout << tmva_dir << std::endl;
57  gROOT->SetMacroPath("$ROOTSYS/root/tmva/test/:.:" + tmva_dir + "/test/:" + gROOT->GetMacroPath());
58  gROOT->ProcessLine(".L TMVAGui.C");
59 
60 }
61 
63 }
64 
66 
67  delete factory;
68  delete signalIn;
69  delete backgroundIn;
70  delete trainingSample;
71 }
72 
73 void TMVAEIDTraining::runTraining(TString cSVFileNameIn) {
74 
75  cSVFileName = cSVFileNameIn;
79  if (!isTransformedInput) {
81  if (cSVFileName != NULL && cSVFileName != "") {
82  return;
83  }
85  }
90 }
91 
93 
95 }
96 
97 void TMVAEIDTraining::setSignalFilesUntransformed(std::vector<TString> inputFileNameIn) {
98  std::cout << "IN setSignalFilesUntransformed" << std::endl;
99  isTransformedInput = false;
100  signalIn = new TChain("showerrec/fEvent");
101  for (unsigned int ii = 0; ii < inputFileNameIn.size(); ii++) {
102  signalIn->Add(inputFileNameIn[ii]);
103  }
104 
105  // if (inputFileNameIn == "FULL") {
106  // // temp hack
107  // signalIn->Add("/export/data/local/novadata/nuevt/mc/S12-10-07/FD/posfkcaf_swap_fd_nhc_[1-9].root");
108  // signalIn->Add("/export/data/local/novadata/nuevt/mc/S12-10-07/FD/posfkcaf_swap_fd_nhc_[1-9][0-9].root");
109  // signalIn->Add("/export/data/local/novadata/nuevt/mc/S12-10-07/FD/posfkcaf_swap_fd_nhc_[1-1][0-4][0-9].root");
110  // signalIn->Add("/export/data/local/novadata/nuevt/mc/S12-10-07/FD/posfkcaf_swap_fd_nhc_150.root");
111  // std::cout << "Added FULL set of 150 signal files" << std::endl;
112  // } else {
113  // signalIn->Add(inputFileNameIn);
114  // std::cout << "Added signal file " << inputFileNameIn << std::endl;
115  // }
116 }
117 
118 void TMVAEIDTraining::setBackgroundFilesUntransformed(std::vector<TString> inputFileNameIn) {
119  std::cout << "IN setBackgroundFilesUntransformed" << std::endl;
120  isTransformedInput = false;
121  backgroundIn = new TChain("showerrec/fEvent");
122  for (unsigned int ii = 0; ii < inputFileNameIn.size(); ii++) {
123  backgroundIn->Add(inputFileNameIn[ii]);
124  }
125  // if (inputFileNameIn == "FULL") {
126  // //temp hack
127  // backgroundIn->Add("/export/data/local/novadata/nuevt/mc/S12-10-07/FD/posfkcaf_fd_nhc_[1-9].root");
128  // backgroundIn->Add("/export/data/local/novadata/nuevt/mc/S12-10-07/FD/posfkcaf_fd_nhc_[1-9][0-9].root");
129  // backgroundIn->Add("/export/data/local/novadata/nuevt/mc/S12-10-07/FD/posfkcaf_fd_nhc_[1-2][0-9][0-9].root");
130  // backgroundIn->Add("/export/data/local/novadata/nuevt/mc/S12-10-07/FD/posfkcaf_fd_nhc_300.root");
131  // std::cout << "Added FULL set of 300 background files" << std::endl;
132  // } else {
133  // backgroundIn->Add(inputFileNameIn);
134  // std::cout << "Added background file " << inputFileNameIn << std::endl;
135  // }
136 }
137 
139  std::cout << "IN setTMVAAlgorithms" << std::endl;
140 
141  // Default MVA methods to be trained + tested
142 
143  // --- Cut optimisation
144  Use["Cuts"] = 0;
145  Use["CutsD"] = 0;
146  Use["CutsPCA"] = 0;
147  Use["CutsGA"] = 0;
148  Use["CutsSA"] = 0;
149  //
150  // --- 1-dimensional likelihood ("naive Bayes estimator")
151  Use["Likelihood"] = 0;
152  Use["LikelihoodD"] = 0; // the "D" extension indicates decorrelated input variables (see option strings)
153  Use["LikelihoodPCA"] = 0; // the "PCA" extension indicates PCA-transformed input variables (see option strings)
154  Use["LikelihoodKDE"] = 0;
155  Use["LikelihoodMIX"] = 0;
156  //
157  // --- Mutidimensional likelihood and Nearest-Neighbour methods
158  Use["PDERS"] = 0;
159  Use["PDERSD"] = 0;
160  Use["PDERSPCA"] = 0;
161  Use["PDEFoam"] = 0;
162  Use["PDEFoamBoost"] = 0; // uses generalised MVA method boosting
163  // if (myMethodList == "KNN") {
164  // Use["KNN"] = 1; // k-nearest neighbour method
165  // } else {
166  Use["KNN"] = 0; // k-nearest neighbour method
167  // }
168  //
169  // --- Linear Discriminant Analysis
170  Use["LD"] = 0; // Linear Discriminant identical to Fisher
171  Use["Fisher"] = 0;
172  Use["FisherG"] = 0;
173  Use["BoostedFisher"] = 0; // uses generalised MVA method boosting
174  Use["HMatrix"] = 0;
175  //
176  // --- Function Discriminant analysis
177  Use["FDA_GA"] = 0; // minimisation of user-defined function using Genetics Algorithm
178  Use["FDA_SA"] = 0;
179  Use["FDA_MC"] = 0;
180  Use["FDA_MT"] = 0;
181  Use["FDA_GAMT"] = 0;
182  Use["FDA_MCMT"] = 0;
183  //
184  // --- Neural Networks (all are feed-forward Multilayer Perceptrons)
185  // if (myMethodList == "MLP") {
186  Use["MLP"] = 0; // Recommended ANN
187  // } else {
188  // Use["MLP"] = 0; // Recommended ANN
189  // }
190  Use["MLPBFGS"] = 0; // Recommended ANN with optional training method
191  Use["MLPBNN"] = 0; // Recommended ANN with BFGS training method and bayesian regulator
192  Use["CFMlpANN"] = 0; // Depreciated ANN from ALEPH
193  // if (myMethodList == "TMlpANN") {
194  // Use["TMlpANN"] = 1; // ROOT's own ANN
195  // } else {
196  Use["TMlpANN"] = 1; // ROOT's own ANN
197  // }
198  //
199  // --- Support Vector Machine
200  Use["SVM"] = 0;
201  //
202  // --- Boosted Decision Trees
203  // if (myMethodList == "BDT") {
204  // Use["BDT"] = 1; // uses Adaptive Boost
205  // } else {
206  Use["BDT"] = 0; // uses Adaptive Boost
207  // }
208  Use["BDTG"] = 0; // uses Gradient Boost
209  Use["BDTB"] = 0; // uses Bagging
210  Use["BDTD"] = 0; // decorrelation + Adaptive Boost
211  Use["BDTF"] = 0; // allow usage of fisher discriminant for node splitting
212  //
213  // --- Friedman's RuleFit method, ie, an optimised series of cuts ("rules")
214  Use["RuleFit"] = 0;
215 
216 }
217 
219  std::cout << "IN internalTMVASetup" << std::endl;
220  std::cout << std::endl;
221  std::cout << "==> Start TMVAClassification" << std::endl;
222 
223  // Select methods (don't look at this code - not of interest)
224  if (myMethodList != "") {
225  for (std::map<std::string, int>::iterator it = Use.begin(); it != Use.end(); it++) it->second = 0;
226 
227  std::vector<TString> mlist = TMVA::gTools().SplitString(myMethodList, ',');
228  if (mlist.size() > 1) {
229  std::cout << "Only one MVA method may be specified at a time";
230  // exit;
231  }
232 
233  for (UInt_t i = 0; i < mlist.size(); i++) {
234  std::string regMethod(mlist[i]);
235 
236  if (Use.find(regMethod) == Use.end()) {
237  std::cout << "Method \"" << regMethod << "\" not known in TMVA under this name. Choose among the following:" << std::endl;
238  for (std::map<std::string, int>::iterator it = Use.begin(); it != Use.end(); it++) std::cout << it->first << " ";
239  std::cout << std::endl;
240  return;
241  }
242  Use[regMethod] = 1;
243  }
244  }
245 
246 }
247 
249  std::cout << "IN preTrainingSetup" << std::endl;
250  // --- Here the preparation phase begins
251 
252  // Create a ROOT output file where TMVA will store ntuples, histograms, etc.
253  // TString outfileName( "TMVA.root" );
254  // TFile*
255  outputFile = TFile::Open(TMVAEIDTraining::outfileName, "RECREATE");
256 
257  // Create the factory object. Later you can choose the methods
258  // whose performance you'd like to investigate. The factory is
259  // the only TMVA object you have to interact with
260  //
261  // The first argument is the base of the name of all the
262  // weightfiles in the directory weight/
263  //
264  // The second argument is the output file for the training results
265  // All TMVA output can be suppressed by removing the "!" (not) in
266  // front of the "Silent" argument in the option string
267  // TMVA::Factory *#if ROOT_VERSION_CODE >= ROOT_VERSION(6,8,0)
268  factory = new TMVA::Factory("TMVAClassification", outputFile,
269  "!V:!Silent:Color:DrawProgressBar:Transformations=I;D;P;G,D:AnalysisType=Classification");
270 
271 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,8,0)
272  dataloader = new TMVA::DataLoader("eid_dataset");
273 #else
275 #endif
276 
277  // If you wish to modify default settings
278  // (please check "src/Config.h" to see all available global options)
279  // (TMVA::gConfig().GetVariablePlotting()).fTimesRMS = 8.0;
280  // (TMVA::gConfig().GetIONames()).fWeightFileDir = "myWeightDirectory";
281 
282  // Define the input variables that shall be used for the MVA training
283  // note that you may also use variable expressions, such as: "3*var1/var2*abs(var3)"
284  // [all types of expressions that can also be parsed by TTree::Draw( "expression" )]
285 
286  dataloader->AddVariable("egLLL", 'F'); //"spec1 := var1*2", "Spectator 1", "units", 'F'
287  dataloader->AddVariable("egLLT", 'F');
288  dataloader->AddVariable("emuLLL", 'F');
289  dataloader->AddVariable("emuLLT", 'F');
290  dataloader->AddVariable("epi0LLL", 'F');
291  dataloader->AddVariable("epi0LLT", 'F');
292  dataloader->AddVariable("epLLL", 'F');
293  dataloader->AddVariable("epLLT", 'F');
294  dataloader->AddVariable("enLLL", 'F');
295  dataloader->AddVariable("enLLT", 'F');
296  dataloader->AddVariable("epiLLL", 'F');
297  dataloader->AddVariable("epiLLT", 'F');
298  dataloader->AddVariable("gap", 'F');
299  dataloader->AddVariable("pi0mass", 'F');
300  dataloader->AddVariable("vtxgev", 'F');
301  dataloader->AddVariable("shE", 'F');
302  std::cout << "### useENu = " << useENu << std::endl;
303  if (useENu) {
304  dataloader->AddVariable("nueRecEnergy", 'F');
305  }
306  // dataloader->AddVariable("costheta", 'F');
307  // dataloader->AddVariable("length", 'F');
308 
309 
310 
311  // You can add so-called "Spectator variables", which are not used in the MVA training,
312  // but will appear in the final "TestTree" produced by TMVA. This TestTree will contain the
313  // input variables, the response values of all trained MVAs, and the spectator variables
314  dataloader->AddSpectator("sbtype", 'I');
315  dataloader->AddSpectator("specType", 'I');
316  // dataloader->AddSpectator("spec2 := var1*3", "Spectator 2", "units", 'F');
317 
318 }
319 
321  std::cout << "IN setupTrainingAndTestTrees" << std::endl;
322  // global event weights per tree (see below for setting event-wise weights)
323 // Double_t signalWeight = 1.0;
324 // Double_t backgroundWeight = 1.0;
325 
326  // You can add an arbitrary number of signal or background trees
327  TCut signalCut = "sbtype==1"; // how to identify signal events
328  TCut backgrCut = "sbtype==0"; // how to identify background events
329  dataloader->SetInputTrees(trainingSample, signalCut, backgrCut);
330 
331 
332  // To give different trees for training and testing, do as follows:
333  // dataloader->AddSignalTree( signalTrainingTree, signalTrainWeight, "Training" );
334  // dataloader->AddSignalTree( signalTestTree, signalTestWeight, "Test" );
335 
336  // Use the following code instead of the above two or four lines to add signal and background
337  // training and test events "by hand"
338  // NOTE that in this case one should not give expressions (such as "var1+var2") in the input
339  // variable definition, but simply compute the expression before adding the event
340  //
341  // // --- begin ----------------------------------------------------------
342  // std::vector<Double_t> vars( 4 ); // vector has size of number of input variables
343  // Float_t treevars[4], weight;
344  //
345  // // Signal
346  // for (UInt_t ivar=0; ivar<4; ivar++) signal->SetBranchAddress( Form( "var%i", ivar+1 ), &(treevars[ivar]) );
347  // for (UInt_t i=0; i<signal->GetEntries(); i++) {
348  // signal->GetEntry(i);
349  // for (UInt_t ivar=0; ivar<4; ivar++) vars[ivar] = treevars[ivar];
350  // // add training and test events; here: first half is training, second is testing
351  // // note that the weight can also be event-wise
352  // if (i < signal->GetEntries()/2.0) dataloader->AddSignalTrainingEvent( vars, signalWeight );
353  // else dataloader->AddSignalTestEvent ( vars, signalWeight );
354  // }
355  //
356  // // Background (has event weights)
357  // background->SetBranchAddress( "weight", &weight );
358  // for (UInt_t ivar=0; ivar<4; ivar++) background->SetBranchAddress( Form( "var%i", ivar+1 ), &(treevars[ivar]) );
359  // for (UInt_t i=0; i<background->GetEntries(); i++) {
360  // background->GetEntry(i);
361  // for (UInt_t ivar=0; ivar<4; ivar++) vars[ivar] = treevars[ivar];
362  // // add training and test events; here: first half is training, second is testing
363  // // note that the weight can also be event-wise
364  // if (i < background->GetEntries()/2) dataloader->AddBackgroundTrainingEvent( vars, backgroundWeight*weight );
365  // else dataloader->AddBackgroundTestEvent ( vars, backgroundWeight*weight );
366  // }
367  // --- end ------------------------------------------------------------
368  //
369  // --- end of tree registration
370 
371  // Set individual event weights (the variables must exist in the original TTree)
372  // for signal : dataloader->SetSignalWeightExpression ("weight1*weight2");
373  // for background: dataloader->SetBackgroundWeightExpression("weight1*weight2");
374  // dataloader->SetBackgroundWeightExpression("weight");
375 
376  // Apply additional cuts on the signal and background samples (can be different)
377  TCut mycuts = ""; // for example: TCut mycuts = "abs(var1)<0.5 && abs(var2-0.5)<1";
378  TCut mycutb = ""; // for example: TCut mycutb = "abs(var1)<0.5";
379 
380 
381 
382 
383  // Tell the factory how to use the training and testing events
384  //
385  // If no numbers of events are given, half of the events in the tree are used
386  // for training, and the other half for testing:
387  // dataloader->PrepareTrainingAndTestTree( mycut, "SplitMode=random:!V" );
388  // To also specify the number of testing events, use:
389  // dataloader->PrepareTrainingAndTestTree( mycut,
390  // "NSigTrain=3000:NBkgTrain=3000:NSigTest=3000:NBkgTest=3000:SplitMode=Random:!V" );
391  // dataloader->PrepareTrainingAndTestTree(mycuts, mycutb,
392  // "nTrain_Signal=0:nTrain_Background=0:SplitMode=Random:NormMode=NumEvents:!V");
393  TString optionsString = TString::Format(
394  "nTrain_Signal=%d:nTrain_Background=%d:SplitMode=Alternate:NormMode=NumEvents:!V", 0, 0);
395 
396  dataloader->PrepareTrainingAndTestTree(mycuts, mycutb,
397  // "nTrain_Signal=0:nTrain_Background=0:nTest_Signal=0:nTest_Background=0:SplitMode=Random:NormMode=NumEvents:!V");
398  // "nTrain_Signal=0:nTrain_Background=0:SplitMode=Alternate:NormMode=NumEvents:!V");
399  optionsString.Data());
400  // "nTrain_Signal=0:nTrain_Background=0:SplitMode=Alternate:NormMode=NumEvents:!V");
401 }
402 
404  std::cout << "IN bookTMVAMethods" << std::endl;
405  // ---- Book MVA methods
406  //
407  // Please lookup the various method configuration options in the corresponding cxx files, eg:
408  // src/MethoCuts.cxx, etc, or here: http://tmva.sourceforge.net/optionRef.html
409  // it is possible to preset ranges in the option string in which the cut optimisation should be done:
410  // "...:CutRangeMin[2]=-1:CutRangeMax[2]=1"...", where [2] is the third input variable
411 
412  // Cut optimisation
413  if (Use["Cuts"])
414 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,8,0)
415  factory->BookMethod(dataloader, TMVA::Types::kCuts, "Cuts",
416  "!H:!V:FitMethod=MC:EffSel:SampleSize=200000:VarProp=FSmart");
417 #else
418  factory->BookMethod(TMVA::Types::kCuts, "Cuts",
419  "!H:!V:FitMethod=MC:EffSel:SampleSize=200000:VarProp=FSmart");
420 #endif
421  if (Use["CutsD"])
422 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,8,0)
423  factory->BookMethod(dataloader, TMVA::Types::kCuts, "CutsD",
424  "!H:!V:FitMethod=MC:EffSel:SampleSize=200000:VarProp=FSmart:VarTransform=Decorrelate");
425 #else
426  factory->BookMethod(TMVA::Types::kCuts, "CutsD",
427  "!H:!V:FitMethod=MC:EffSel:SampleSize=200000:VarProp=FSmart:VarTransform=Decorrelate");
428 #endif
429 
430  if (Use["CutsPCA"])
431 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,8,0)
432  factory->BookMethod(dataloader, TMVA::Types::kCuts, "CutsPCA",
433  "!H:!V:FitMethod=MC:EffSel:SampleSize=200000:VarProp=FSmart:VarTransform=PCA");
434 #else
435  factory->BookMethod(TMVA::Types::kCuts, "CutsPCA",
436  "!H:!V:FitMethod=MC:EffSel:SampleSize=200000:VarProp=FSmart:VarTransform=PCA");
437 #endif
438 
439  if (Use["CutsGA"])
440 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,8,0)
441  factory->BookMethod(dataloader, TMVA::Types::kCuts, "CutsGA",
442  "H:!V:FitMethod=GA:CutRangeMin[0]=-10:CutRangeMax[0]=10:VarProp[1]=FMax:EffSel:Steps=30:Cycles=3:PopSize=400:SC_steps=10:SC_rate=5:SC_factor=0.95");
443 #else
444  factory->BookMethod(TMVA::Types::kCuts, "CutsGA",
445  "H:!V:FitMethod=GA:CutRangeMin[0]=-10:CutRangeMax[0]=10:VarProp[1]=FMax:EffSel:Steps=30:Cycles=3:PopSize=400:SC_steps=10:SC_rate=5:SC_factor=0.95");
446 #endif
447 
448  if (Use["CutsSA"])
449 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,8,0)
450  factory->BookMethod(dataloader, TMVA::Types::kCuts, "CutsSA",
451  "!H:!V:FitMethod=SA:EffSel:MaxCalls=150000:KernelTemp=IncAdaptive:InitialTemp=1e+6:MinTemp=1e-6:Eps=1e-10:UseDefaultScale");
452 #else
453  factory->BookMethod(TMVA::Types::kCuts, "CutsSA",
454  "!H:!V:FitMethod=SA:EffSel:MaxCalls=150000:KernelTemp=IncAdaptive:InitialTemp=1e+6:MinTemp=1e-6:Eps=1e-10:UseDefaultScale");
455 #endif
456 
457  // Likelihood ("naive Bayes estimator")
458  if (Use["Likelihood"])
459 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,8,0)
460  factory->BookMethod(dataloader, TMVA::Types::kLikelihood, "Likelihood",
461  "H:!V:TransformOutput:PDFInterpol=Spline2:NSmoothSig[0]=20:NSmoothBkg[0]=20:NSmoothBkg[1]=10:NSmooth=1:NAvEvtPerBin=50");
462 #else
463  factory->BookMethod(TMVA::Types::kLikelihood, "Likelihood",
464  "H:!V:TransformOutput:PDFInterpol=Spline2:NSmoothSig[0]=20:NSmoothBkg[0]=20:NSmoothBkg[1]=10:NSmooth=1:NAvEvtPerBin=50");
465 #endif
466 
467  // Decorrelated likelihood
468  if (Use["LikelihoodD"])
469 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,8,0)
470  factory->BookMethod(dataloader, TMVA::Types::kLikelihood, "LikelihoodD",
471  "!H:!V:TransformOutput:PDFInterpol=Spline2:NSmoothSig[0]=20:NSmoothBkg[0]=20:NSmooth=5:NAvEvtPerBin=50:VarTransform=Decorrelate");
472 #else
473  factory->BookMethod(TMVA::Types::kLikelihood, "LikelihoodD",
474  "!H:!V:TransformOutput:PDFInterpol=Spline2:NSmoothSig[0]=20:NSmoothBkg[0]=20:NSmooth=5:NAvEvtPerBin=50:VarTransform=Decorrelate");
475 #endif
476 
477  // PCA-transformed likelihood
478  if (Use["LikelihoodPCA"])
479 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,8,0)
480  factory->BookMethod(dataloader, TMVA::Types::kLikelihood, "LikelihoodPCA",
481  "!H:!V:!TransformOutput:PDFInterpol=Spline2:NSmoothSig[0]=20:NSmoothBkg[0]=20:NSmooth=5:NAvEvtPerBin=50:VarTransform=PCA");
482 #else
483  factory->BookMethod(TMVA::Types::kLikelihood, "LikelihoodPCA",
484  "!H:!V:!TransformOutput:PDFInterpol=Spline2:NSmoothSig[0]=20:NSmoothBkg[0]=20:NSmooth=5:NAvEvtPerBin=50:VarTransform=PCA");
485 #endif
486 
487  // Use a kernel density estimator to approximate the PDFs
488  if (Use["LikelihoodKDE"])
489 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,8,0)
490  factory->BookMethod(dataloader, TMVA::Types::kLikelihood, "LikelihoodKDE",
491  "!H:!V:!TransformOutput:PDFInterpol=KDE:KDEtype=Gauss:KDEiter=Adaptive:KDEFineFactor=0.3:KDEborder=None:NAvEvtPerBin=50");
492 #else
493  factory->BookMethod(TMVA::Types::kLikelihood, "LikelihoodKDE",
494  "!H:!V:!TransformOutput:PDFInterpol=KDE:KDEtype=Gauss:KDEiter=Adaptive:KDEFineFactor=0.3:KDEborder=None:NAvEvtPerBin=50");
495 #endif
496  // Use a variable-dependent mix of splines and kernel density estimator
497  if (Use["LikelihoodMIX"])
498 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,8,0)
499  factory->BookMethod(dataloader, TMVA::Types::kLikelihood, "LikelihoodMIX",
500  "!H:!V:!TransformOutput:PDFInterpolSig[0]=KDE:PDFInterpolBkg[0]=KDE:PDFInterpolSig[1]=KDE:PDFInterpolBkg[1]=KDE:PDFInterpolSig[2]=Spline2:PDFInterpolBkg[2]=Spline2:PDFInterpolSig[3]=Spline2:PDFInterpolBkg[3]=Spline2:KDEtype=Gauss:KDEiter=Nonadaptive:KDEborder=None:NAvEvtPerBin=50");
501 #else
502  factory->BookMethod(TMVA::Types::kLikelihood, "LikelihoodMIX",
503  "!H:!V:!TransformOutput:PDFInterpolSig[0]=KDE:PDFInterpolBkg[0]=KDE:PDFInterpolSig[1]=KDE:PDFInterpolBkg[1]=KDE:PDFInterpolSig[2]=Spline2:PDFInterpolBkg[2]=Spline2:PDFInterpolSig[3]=Spline2:PDFInterpolBkg[3]=Spline2:KDEtype=Gauss:KDEiter=Nonadaptive:KDEborder=None:NAvEvtPerBin=50");
504 #endif
505 
506  // Test the multi-dimensional probability density estimator
507  // here are the options strings for the MinMax and RMS methods, respectively:
508  // "!H:!V:VolumeRangeMode=MinMax:DeltaFrac=0.2:KernelEstimator=Gauss:GaussSigma=0.3" );
509  // "!H:!V:VolumeRangeMode=RMS:DeltaFrac=3:KernelEstimator=Gauss:GaussSigma=0.3" );
510  if (Use["PDERS"]){
511 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,8,0)
512  factory->BookMethod(dataloader, TMVA::Types::kPDERS, "PDERS",
513  "!H:!V:NormTree=T:VolumeRangeMode=Adaptive:KernelEstimator=Gauss:GaussSigma=0.3:NEventsMin=400:NEventsMax=600");
514 #else
515  factory->BookMethod(TMVA::Types::kPDERS, "PDERS",
516  "!H:!V:NormTree=T:VolumeRangeMode=Adaptive:KernelEstimator=Gauss:GaussSigma=0.3:NEventsMin=400:NEventsMax=600");
517 #endif
518  }
519  if (Use["PDERSD"]){
520 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,8,0)
521  factory->BookMethod(dataloader, TMVA::Types::kPDERS, "PDERSD",
522  "!H:!V:VolumeRangeMode=Adaptive:KernelEstimator=Gauss:GaussSigma=0.3:NEventsMin=400:NEventsMax=600:VarTransform=Decorrelate");
523 #else
524  factory->BookMethod(TMVA::Types::kPDERS, "PDERSD",
525  "!H:!V:VolumeRangeMode=Adaptive:KernelEstimator=Gauss:GaussSigma=0.3:NEventsMin=400:NEventsMax=600:VarTransform=Decorrelate");
526 #endif
527  }
528  if (Use["PDERSPCA"]){
529 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,8,0)
530  factory->BookMethod(dataloader, TMVA::Types::kPDERS, "PDERSPCA",
531  "!H:!V:VolumeRangeMode=Adaptive:KernelEstimator=Gauss:GaussSigma=0.3:NEventsMin=400:NEventsMax=600:VarTransform=PCA");
532 #else
533  factory->BookMethod(TMVA::Types::kPDERS, "PDERSPCA",
534  "!H:!V:VolumeRangeMode=Adaptive:KernelEstimator=Gauss:GaussSigma=0.3:NEventsMin=400:NEventsMax=600:VarTransform=PCA");
535 #endif
536  }
537 
538  // Multi-dimensional likelihood estimator using self-adapting phase-space binning
539  if (Use["PDEFoam"]){
540 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,8,0)
541  factory->BookMethod(dataloader, TMVA::Types::kPDEFoam, "PDEFoam",
542  "!H:!V:SigBgSeparate=F:TailCut=0.001:VolFrac=0.0666:nActiveCells=500:nSampl=2000:nBin=5:Nmin=100:Kernel=None:Compress=T");
543 #else
544  factory->BookMethod(TMVA::Types::kPDEFoam, "PDEFoam",
545  "!H:!V:SigBgSeparate=F:TailCut=0.001:VolFrac=0.0666:nActiveCells=500:nSampl=2000:nBin=5:Nmin=100:Kernel=None:Compress=T");
546 #endif
547  }
548  if (Use["PDEFoamBoost"]){
549 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,8,0)
550  factory->BookMethod(dataloader, TMVA::Types::kPDEFoam, "PDEFoamBoost",
551  "!H:!V:Boost_Num=30:Boost_Transform=linear:SigBgSeparate=F:MaxDepth=4:UseYesNoCell=T:DTLogic=MisClassificationError:FillFoamWithOrigWeights=F:TailCut=0:nActiveCells=500:nBin=20:Nmin=400:Kernel=None:Compress=T");
552 #else
553  factory->BookMethod(TMVA::Types::kPDEFoam, "PDEFoamBoost",
554  "!H:!V:Boost_Num=30:Boost_Transform=linear:SigBgSeparate=F:MaxDepth=4:UseYesNoCell=T:DTLogic=MisClassificationError:FillFoamWithOrigWeights=F:TailCut=0:nActiveCells=500:nBin=20:Nmin=400:Kernel=None:Compress=T");
555 #endif
556  }
557 
558  // K-Nearest Neighbour classifier (KNN)
559  if (Use["KNN"])
560 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,8,0)
561  factory->BookMethod(dataloader, TMVA::Types::kKNN, "KNN",
562  "H:CreateMVAPdfs:nkNN=20:ScaleFrac=0.8:SigmaFact=1.0:Kernel=Gaus:UseKernel=F:UseWeight=T:!Trim");
563 #else
564  factory->BookMethod(TMVA::Types::kKNN, "KNN",
565  "H:CreateMVAPdfs:nkNN=20:ScaleFrac=0.8:SigmaFact=1.0:Kernel=Gaus:UseKernel=F:UseWeight=T:!Trim");
566 #endif
567 
568  // H-Matrix (chi2-squared) method
569  if (Use["HMatrix"])
570 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,8,0)
571  factory->BookMethod(dataloader, TMVA::Types::kHMatrix, "HMatrix", "!H:!V:VarTransform=None");
572 #else
573  factory->BookMethod(TMVA::Types::kHMatrix, "HMatrix", "!H:!V:VarTransform=None");
574 #endif
575 
576  // Linear discriminant (same as Fisher discriminant)
577  if (Use["LD"])
578 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,8,0)
579  factory->BookMethod(dataloader, TMVA::Types::kLD, "LD", "H:!V:VarTransform=None:CreateMVAPdfs:PDFInterpolMVAPdf=Spline2:NbinsMVAPdf=50:NsmoothMVAPdf=10");
580 #else
581  factory->BookMethod(TMVA::Types::kLD, "LD", "H:!V:VarTransform=None:CreateMVAPdfs:PDFInterpolMVAPdf=Spline2:NbinsMVAPdf=50:NsmoothMVAPdf=10");
582 #endif
583 
584  // Fisher discriminant (same as LD)
585  if (Use["Fisher"])
586 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,8,0)
587  factory->BookMethod(dataloader, TMVA::Types::kFisher, "Fisher", "H:!V:Fisher:VarTransform=None:CreateMVAPdfs:PDFInterpolMVAPdf=Spline2:NbinsMVAPdf=50:NsmoothMVAPdf=10");
588 #else
589  factory->BookMethod(TMVA::Types::kFisher, "Fisher", "H:!V:Fisher:VarTransform=None:CreateMVAPdfs:PDFInterpolMVAPdf=Spline2:NbinsMVAPdf=50:NsmoothMVAPdf=10");
590 #endif
591 
592  // Fisher with Gauss-transformed input variables
593  if (Use["FisherG"])
594 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,8,0)
595  factory->BookMethod(dataloader, TMVA::Types::kFisher, "FisherG", "H:!V:VarTransform=Gauss");
596 #else
597  factory->BookMethod(TMVA::Types::kFisher, "FisherG", "H:!V:VarTransform=Gauss");
598 #endif
599 
600  // Composite classifier: ensemble (tree) of boosted Fisher classifiers
601  if (Use["BoostedFisher"])
602 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,8,0)
603  factory->BookMethod(dataloader, TMVA::Types::kFisher, "BoostedFisher",
604  "H:!V:Boost_Num=20:Boost_Transform=log:Boost_Type=AdaBoost:Boost_AdaBoostBeta=0.2:!Boost_DetailedMonitoring");
605 #else
606  factory->BookMethod(TMVA::Types::kFisher, "BoostedFisher",
607  "H:!V:Boost_Num=20:Boost_Transform=log:Boost_Type=AdaBoost:Boost_AdaBoostBeta=0.2:!Boost_DetailedMonitoring");
608 #endif
609 
610  // Function discrimination analysis (FDA) -- test of various fitters - the recommended one is Minuit (or GA or SA)
611  if (Use["FDA_MC"]){
612 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,8,0)
613  factory->BookMethod(dataloader, TMVA::Types::kFDA, "FDA_MC",
614  "H:!V:Formula=(0)+(1)*x0+(2)*x1+(3)*x2+(4)*x3:ParRanges=(-1,1);(-10,10);(-10,10);(-10,10);(-10,10):FitMethod=MC:SampleSize=100000:Sigma=0.1");
615 #else
616  factory->BookMethod(TMVA::Types::kFDA, "FDA_MC",
617  "H:!V:Formula=(0)+(1)*x0+(2)*x1+(3)*x2+(4)*x3:ParRanges=(-1,1);(-10,10);(-10,10);(-10,10);(-10,10):FitMethod=MC:SampleSize=100000:Sigma=0.1");
618 #endif
619  }
620  if (Use["FDA_GA"]){
621  // can also use Simulated Annealing (SA) algorithm (see Cuts_SA options])
622 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,8,0)
623  factory->BookMethod(dataloader, TMVA::Types::kFDA, "FDA_GA",
624  "H:!V:Formula=(0)+(1)*x0+(2)*x1+(3)*x2+(4)*x3:ParRanges=(-1,1);(-10,10);(-10,10);(-10,10);(-10,10):FitMethod=GA:PopSize=300:Cycles=3:Steps=20:Trim=True:SaveBestGen=1");
625 #else
626  factory->BookMethod(TMVA::Types::kFDA, "FDA_GA",
627  "H:!V:Formula=(0)+(1)*x0+(2)*x1+(3)*x2+(4)*x3:ParRanges=(-1,1);(-10,10);(-10,10);(-10,10);(-10,10):FitMethod=GA:PopSize=300:Cycles=3:Steps=20:Trim=True:SaveBestGen=1");
628 #endif
629  }
630  if (Use["FDA_SA"]){
631  // can also use Simulated Annealing (SA) algorithm (see Cuts_SA options])
632 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,8,0)
633  factory->BookMethod(dataloader, TMVA::Types::kFDA, "FDA_SA",
634  "H:!V:Formula=(0)+(1)*x0+(2)*x1+(3)*x2+(4)*x3:ParRanges=(-1,1);(-10,10);(-10,10);(-10,10);(-10,10):FitMethod=SA:MaxCalls=15000:KernelTemp=IncAdaptive:InitialTemp=1e+6:MinTemp=1e-6:Eps=1e-10:UseDefaultScale");
635 #else
636  factory->BookMethod(TMVA::Types::kFDA, "FDA_SA",
637  "H:!V:Formula=(0)+(1)*x0+(2)*x1+(3)*x2+(4)*x3:ParRanges=(-1,1);(-10,10);(-10,10);(-10,10);(-10,10):FitMethod=SA:MaxCalls=15000:KernelTemp=IncAdaptive:InitialTemp=1e+6:MinTemp=1e-6:Eps=1e-10:UseDefaultScale");
638 #endif
639  }
640  if (Use["FDA_MT"]){
641 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,8,0)
642  factory->BookMethod(dataloader, TMVA::Types::kFDA, "FDA_MT",
643  "H:!V:Formula=(0)+(1)*x0+(2)*x1+(3)*x2+(4)*x3:ParRanges=(-1,1);(-10,10);(-10,10);(-10,10);(-10,10):FitMethod=MINUIT:ErrorLevel=1:PrintLevel=-1:FitStrategy=2:UseImprove:UseMinos:SetBatch");
644 #else
645  factory->BookMethod(TMVA::Types::kFDA, "FDA_MT",
646  "H:!V:Formula=(0)+(1)*x0+(2)*x1+(3)*x2+(4)*x3:ParRanges=(-1,1);(-10,10);(-10,10);(-10,10);(-10,10):FitMethod=MINUIT:ErrorLevel=1:PrintLevel=-1:FitStrategy=2:UseImprove:UseMinos:SetBatch");
647 #endif
648  }
649  if (Use["FDA_GAMT"]){
650 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,8,0)
651  factory->BookMethod(dataloader, TMVA::Types::kFDA, "FDA_GAMT",
652  "H:!V:Formula=(0)+(1)*x0+(2)*x1+(3)*x2+(4)*x3:ParRanges=(-1,1);(-10,10);(-10,10);(-10,10);(-10,10):FitMethod=GA:Converger=MINUIT:ErrorLevel=1:PrintLevel=-1:FitStrategy=0:!UseImprove:!UseMinos:SetBatch:Cycles=1:PopSize=5:Steps=5:Trim");
653 #else
654  factory->BookMethod(TMVA::Types::kFDA, "FDA_GAMT",
655  "H:!V:Formula=(0)+(1)*x0+(2)*x1+(3)*x2+(4)*x3:ParRanges=(-1,1);(-10,10);(-10,10);(-10,10);(-10,10):FitMethod=GA:Converger=MINUIT:ErrorLevel=1:PrintLevel=-1:FitStrategy=0:!UseImprove:!UseMinos:SetBatch:Cycles=1:PopSize=5:Steps=5:Trim");
656 #endif
657  }
658  if (Use["FDA_MCMT"]){
659 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,8,0)
660  factory->BookMethod(dataloader, TMVA::Types::kFDA, "FDA_MCMT",
661  "H:!V:Formula=(0)+(1)*x0+(2)*x1+(3)*x2+(4)*x3:ParRanges=(-1,1);(-10,10);(-10,10);(-10,10);(-10,10):FitMethod=MC:Converger=MINUIT:ErrorLevel=1:PrintLevel=-1:FitStrategy=0:!UseImprove:!UseMinos:SetBatch:SampleSize=20");
662 #else
663  factory->BookMethod(TMVA::Types::kFDA, "FDA_MCMT",
664  "H:!V:Formula=(0)+(1)*x0+(2)*x1+(3)*x2+(4)*x3:ParRanges=(-1,1);(-10,10);(-10,10);(-10,10);(-10,10):FitMethod=MC:Converger=MINUIT:ErrorLevel=1:PrintLevel=-1:FitStrategy=0:!UseImprove:!UseMinos:SetBatch:SampleSize=20");
665 #endif
666  }
667  // TMVA ANN: MLP (recommended ANN) -- all ANNs in TMVA are Multilayer Perceptrons
668  if (Use["MLP"]){
669 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,8,0)
670  // factory->BookMethod(TMVA::Types::kMLP, "MLP",
671  // "H:!V:CreateMVAPdfs:NeuronType=tanh:VarTransform=N:NCycles=600:HiddenLayers=N+5:TestRate=5:!UseRegulator");
672  factory->BookMethod(dataloader, TMVA::Types::kMLP, "MLP",
673  "H:!V:CreateMVAPdfs:NeuronType=tanh:TrainingMethod=BFGS:VarTransform=N_AllClasses:NCycles=600:HiddenLayers=22,12,6:TestRate=10:!UseRegulator");
674 #else
675  factory->BookMethod(TMVA::Types::kMLP, "MLP",
676  "H:!V:CreateMVAPdfs:NeuronType=tanh:TrainingMethod=BFGS:VarTransform=N_AllClasses:NCycles=600:HiddenLayers=22,12,6:TestRate=10:!UseRegulator");
677 #endif
678  }
679  if (Use["MLPBFGS"]){
680 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,8,0)
681  factory->BookMethod(dataloader, TMVA::Types::kMLP, "MLPBFGS", "H:!V:NeuronType=tanh:VarTransform=N:NCycles=600:HiddenLayers=N+5:TestRate=5:TrainingMethod=BFGS:!UseRegulator");
682 #else
683  factory->BookMethod(TMVA::Types::kMLP, "MLPBFGS", "H:!V:NeuronType=tanh:VarTransform=N:NCycles=600:HiddenLayers=N+5:TestRate=5:TrainingMethod=BFGS:!UseRegulator");
684 #endif
685  }
686  if (Use["MLPBNN"]){
687 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,8,0)
688  factory->BookMethod(dataloader, TMVA::Types::kMLP, "MLPBNN", "H:!V:NeuronType=tanh:VarTransform=N:NCycles=600:HiddenLayers=N+5:TestRate=5:TrainingMethod=BFGS:UseRegulator"); // BFGS training with bayesian regulators
689 #else
690  factory->BookMethod(TMVA::Types::kMLP, "MLPBNN", "H:!V:NeuronType=tanh:VarTransform=N:NCycles=600:HiddenLayers=N+5:TestRate=5:TrainingMethod=BFGS:UseRegulator"); // BFGS training with bayesian regulators
691 #endif
692  }
693 
694  // CF(Clermont-Ferrand)ANN
695  if (Use["CFMlpANN"])
696 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,8,0)
697  factory->BookMethod(dataloader, TMVA::Types::kCFMlpANN, "CFMlpANN", "!H:!V:NCycles=2000:HiddenLayers=N+1,N"); // n_cycles:#nodes:#nodes:...
698 #else
699  factory->BookMethod(TMVA::Types::kCFMlpANN, "CFMlpANN", "!H:!V:NCycles=2000:HiddenLayers=N+1,N"); // n_cycles:#nodes:#nodes:...
700 #endif
701 
702  // Tmlp(Root)ANN
703  if (Use["TMlpANN"])
704 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,8,0)
705  // factory->BookMethod(TMVA::Types::kTMlpANN, "TMlpANN", "H:!V:NCycles=200:HiddenLayers=N+1,N:LearningMethod=BFGS:ValidationFraction=0.3"); // n_cycles:#nodes:#nodes:...
706  // factory->BookMethod(TMVA::Types::kTMlpANN, "TMlpANN",
707  // "H:!V:NCycles=600:HiddenLayers=N+1,N:LearningMethod=BFGS:ValidationFraction=0.5"); // n_cycles:#nodes:#nodes:...
708  // factory->BookMethod(TMVA::Types::kTMlpANN, "TMlpANN",
709  // "H:!V:NCycles=600:HiddenLayers=22,12,6:LearningMethod=BFGS:ValidationFraction=0.5:VarTransform=N_AllClasses"); // n_cycles:#nodes:#nodes:...
710  // factory->BookMethod(TMVA::Types::kTMlpANN, "TMlpANN",
711  // "H:V:NCycles=600:HiddenLayers=22,12,6:LearningMethod=BFGS:ValidationFraction=0.0:VarTransform=N_AllClasses"); // n_cycles:#nodes:#nodes:...
712  factory->BookMethod(dataloader, TMVA::Types::kTMlpANN, "TMlpANN",
713  "H:V:NCycles=600:HiddenLayers=22,12,6:LearningMethod=BFGS:ValidationFraction=0.0:VarTransform=N_AllClasses"); // n_cycles:#nodes:#nodes:...
714 #else
715  factory->BookMethod(TMVA::Types::kTMlpANN, "TMlpANN",
716  "H:V:NCycles=600:HiddenLayers=22,12,6:LearningMethod=BFGS:ValidationFraction=0.0:VarTransform=N_AllClasses"); // n_cycles:#nodes:#nodes:...
717 #endif
718 
719  // Support Vector Machine
720  if (Use["SVM"])
721 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,8,0)
722  factory->BookMethod(dataloader, TMVA::Types::kSVM, "SVM", "Gamma=0.25:Tol=0.001:VarTransform=Norm");
723 #else
724  factory->BookMethod(TMVA::Types::kSVM, "SVM", "Gamma=0.25:Tol=0.001:VarTransform=Norm");
725 #endif
726 
727  // Boosted Decision Trees
728  // Gradient Boost
729  if (Use["BDTG"]){
730 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,8,0)
731  factory->BookMethod(dataloader, TMVA::Types::kBDT, "BDTG",
732  "!H:!V:NTrees=1000:BoostType=Grad:Shrinkage=0.10:UseBaggedGrad:GradBaggingFraction=0.5:nCuts=20:NNodesMax=5");
733 #else
734  factory->BookMethod(TMVA::Types::kBDT, "BDTG",
735  "!H:!V:NTrees=1000:BoostType=Grad:Shrinkage=0.10:UseBaggedGrad:GradBaggingFraction=0.5:nCuts=20:NNodesMax=5");
736 #endif
737 }
738 if (Use["BDT"]){
739  // Adaptive Boost
740 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,8,0)
741  factory->BookMethod(dataloader, TMVA::Types::kBDT, "BDT",
742  "!H:!V:NTrees=850:nEventsMin=150:MaxDepth=3:BoostType=AdaBoost:AdaBoostBeta=0.5:SeparationType=GiniIndex:nCuts=20:PruneMethod=NoPruning");
743 #else
744  factory->BookMethod(TMVA::Types::kBDT, "BDT",
745  "!H:!V:NTrees=850:nEventsMin=150:MaxDepth=3:BoostType=AdaBoost:AdaBoostBeta=0.5:SeparationType=GiniIndex:nCuts=20:PruneMethod=NoPruning");
746 #endif
747 }
748 
749 if (Use["BDTB"]){
750  // Bagging
751 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,8,0)
752  factory->BookMethod(dataloader, TMVA::Types::kBDT, "BDTB",
753  "!H:!V:NTrees=400:BoostType=Bagging:SeparationType=GiniIndex:nCuts=20:PruneMethod=NoPruning");
754 #else
755  factory->BookMethod(TMVA::Types::kBDT, "BDTB",
756  "!H:!V:NTrees=400:BoostType=Bagging:SeparationType=GiniIndex:nCuts=20:PruneMethod=NoPruning");
757 #endif
758 }
759 if (Use["BDTD"]){
760  // Decorrelation + Adaptive Boost
761 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,8,0)
762  factory->BookMethod(dataloader, TMVA::Types::kBDT, "BDTD",
763  "!H:!V:NTrees=400:nEventsMin=400:MaxDepth=3:BoostType=AdaBoost:SeparationType=GiniIndex:nCuts=20:PruneMethod=NoPruning:VarTransform=Decorrelate");
764 #else
765  factory->BookMethod(TMVA::Types::kBDT, "BDTD",
766  "!H:!V:NTrees=400:nEventsMin=400:MaxDepth=3:BoostType=AdaBoost:SeparationType=GiniIndex:nCuts=20:PruneMethod=NoPruning:VarTransform=Decorrelate");
767 #endif
768 }
769 if (Use["BDTF"]){
770  // Allow Using Fisher discriminant in node splitting for (strong) linearly correlated variables
771 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,8,0)
772  factory->BookMethod(dataloader, TMVA::Types::kBDT, "BDTMitFisher",
773  "!H:!V:NTrees=50:nEventsMin=150:UseFisherCuts:MaxDepth=3:BoostType=AdaBoost:AdaBoostBeta=0.5:SeparationType=GiniIndex:nCuts=20:PruneMethod=NoPruning");
774 #else
775  factory->BookMethod(TMVA::Types::kBDT, "BDTMitFisher",
776  "!H:!V:NTrees=50:nEventsMin=150:UseFisherCuts:MaxDepth=3:BoostType=AdaBoost:AdaBoostBeta=0.5:SeparationType=GiniIndex:nCuts=20:PruneMethod=NoPruning");
777 #endif
778 }
779  // RuleFit -- TMVA implementation of Friedman's method
780 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,8,0)
781 factory->BookMethod(dataloader, TMVA::Types::kRuleFit, "RuleFit",
782  "H:!V:RuleFitModule=RFTMVA:Model=ModRuleLinear:MinImp=0.001:RuleMinDist=0.001:NTrees=20:fEventsMin=0.01:fEventsMax=0.5:GDTau=-1.0:GDTauPrec=0.01:GDStep=0.01:GDNSteps=10000:GDErrScale=1.02");
783 #else
784 factory->BookMethod(TMVA::Types::kRuleFit, "RuleFit",
785  "H:!V:RuleFitModule=RFTMVA:Model=ModRuleLinear:MinImp=0.001:RuleMinDist=0.001:NTrees=20:fEventsMin=0.01:fEventsMax=0.5:GDTau=-1.0:GDTauPrec=0.01:GDStep=0.01:GDNSteps=10000:GDErrScale=1.02");
786 #endif
787 methodsBooked = true;
788 }
789 
791  std::cout << "IN trainTestAndEvaluateTMVAMethods" << std::endl;
792  // For an example of the category classifier usage, see: TMVAClassificationCategory
793 
794  // --------------------------------------------------------------------------------------------------
795 
796  // ---- Now you can optimize the setting (configuration) of the MVAs using the set of training events
797 
798  // factory->OptimizeAllMethods("SigEffAt001","Scan");
799  // factory->OptimizeAllMethods("ROCIntegral","GA");
800 
801  // --------------------------------------------------------------------------------------------------
802 
803  // ---- Now you can tell the factory to train, test, and evaluate the MVAs
804 
805  std::cout << "PRE-TRAIN: trainingSample->GetEntries() = " << trainingSample->GetEntries() << std::endl;
806 
807  // Train MVAs using the set of training events
808  factory->TrainAllMethods();
809 
810  // ---- Evaluate all MVAs using the set of test events
811  factory->TestAllMethods();
812 
813  // ----- Evaluate and compare performance of all configured MVAs
814  factory->EvaluateAllMethods();
815 }
816 
818 
819  std::cout << "IN saveOutput" << std::endl;
820 
821  // Save the output
822  outputFile->Close();
823 
824  std::cout << "==> Wrote root file: " << outputFile->GetName() << std::endl;
825  std::cout << "==> TMVAClassification is done!" << std::endl;
826 
827 }
828 
830  std::cout << "IN preTrainHistograms" << std::endl;
831  TCut MCNUECCQE = "evtTrueNuPdg==12&&(evtRwRand[0]>0&&evtRwRand[0]<evtRwProb[0])";
832  TCut MCBKG = "(evtRwRand[2]>0&&evtRwRand[2]<evtRwProb[2])";
833  htruenue = new TH1F("htruenue", "htruenue", 100, 0, 5);
834  htruenupi0 = new TH1F("htruenupi0", "htruenupi0", 100, 0, 5);
835  signalIn->Draw("evtTrueNuEnergy>>htruenue", MCNUECCQE);
836  backgroundIn->Draw("evtTrueNuEnergy>>htruenupi0", MCBKG);
837 
838 // double nsig = htruenue->GetEntries();
839 // double nbkg = htruenupi0->GetEntries();
840 
841  TCanvas *c0 = new TCanvas("c0", "c0", 200, 10, 1200, 600);
842  c0->Divide(2, 1);
843  c0->cd(1);
844  htruenue->Draw();
845  c0->cd(2);
846  htruenupi0->Draw();
847  c0->Update();
848  c0->SaveAs(plotDir + "JMCanvas0.png");
849  c0->SaveAs(plotDir + "JMCanvas0.eps");
850 }
851 
853  std::cout << "--------------------------------------------" << std::endl;
854  std::cout << "egLLL/T" << egLLL << "," << egLLT << std::endl;
855  std::cout << "emuLLL/T" << emuLLL << "," << emuLLT << std::endl;
856  std::cout << "epi0LLL/T" << epi0LLL << "," << epi0LLT << std::endl;
857  std::cout << "epLLL/T" << epLLL << "," << epLLT << std::endl;
858  std::cout << "enLLL/T" << enLLL << "," << enLLT << std::endl;
859  std::cout << "epiLLL/T" << epiLLL << "," << epiLLT << std::endl;
860  std::cout << "gap" << gap << std::endl;
861  std::cout << "pi0mass" << pi0mass << std::endl;
862  std::cout << "vtxgev" << vtxgev << std::endl;
863  std::cout << "shE" << shE << std::endl;
864  if (useENu) {
865  std::cout << "nueRecEnergy" << nueRecEnergy << std::endl;
866  }
867  std::cout << "costheta" << costheta << std::endl;
868  std::cout << "length" << length << std::endl;
869  std::cout << "type" << type << std::endl;
870 }
871 
873 
874  if (step == CSVOperation::INIT) {
875  std::cout << "Dumping training sample to CSV File:" << cSVFileName << std::endl;
876  lineBufferSize = 0;
877  outfile.open((cSVFileName.Data()));
878  // Header
879  outfile << "egLLL\t" << "egLLT\t" << "emuLLL\t" << "emuLLT\t"
880  << "epi0LLL\t" << "epi0LLT\t" << "epLLL\t" << "epLLT\t"
881  << "enLLL\t" << "enLLT\t" << "epiLLL\t" << "epiLLT\t"
882  << "gap\t" << "pi0mass\t" << "vtxgev\t" << "shE\t"
883  << "nueRecEnergy\t" << "costheta\t" << "length\t" << "type\t"
884  << "specType\t" << "reWeight" << "\n";
885  } else if (step == CSVOperation::FILL) {
886 
887  std::string sigback;
888  if (type == 0) {
889  sigback = "Background";
890  } else {
891  sigback = "Signal";
892  }
893  outfile << egLLL << "\t" << egLLT << "\t" << emuLLL << "\t" << emuLLT << "\t"
894  << epi0LLL << "\t" << epi0LLT << "\t" << epLLL << "\t" << epLLT << "\t"
895  << enLLL << "\t" << enLLT << "\t" << epiLLL << "\t" << epiLLT << "\t"
896  << gap << "\t" << pi0mass << "\t" << vtxgev << "\t" << shE << "\t"
897  << nueRecEnergy << "\t" << costheta << "\t" << length << "\t" << sigback << "\t"
898  << specType << "\t" << reWeight << "\n";
899 
900 
901  } else if (step == CSVOperation::CLOSE) {
902  outfile.close();
903  }
904 }
905 
906 /**
907  *
908  * Transform input trees to internal trees used for training.
909  *
910  * @param fillTrainingSample
911  */
912 void TMVAEIDTraining::transformInputTrees(bool fillTrainingSample) {
913  std::cout << "IN transformInputTrees" << std::endl;
914 
915 // double Elow = 1.0;
916 // double Ehigh = 3.0;
917  double rsig = 1.0;
918  double rbg = 0.9;
919 
922 
924 
925  std::cout << "fillTrainingSample = " << fillTrainingSample << std::endl;
926  if (fillTrainingSample) {
927  trainingSample = new TTree("TrainingSample", "Filtered Monte Carlo Events");
928  trainingSample->Branch("egLLL", &egLLL, "egLLL/F");
929  trainingSample->Branch("egLLT", &egLLT, "egLLT/F");
930  trainingSample->Branch("emuLLL", &emuLLL, "emuLLL/F");
931  trainingSample->Branch("emuLLT", &emuLLT, "emuLLT/F");
932  trainingSample->Branch("epi0LLL", &epi0LLL, "epi0LLL/F");
933  trainingSample->Branch("epi0LLT", &epi0LLT, "epi0LLT/F");
934  trainingSample->Branch("epLLL", &epLLL, "epLLL/F");
935  trainingSample->Branch("epLLT", &epLLT, "epLLT/F");
936  trainingSample->Branch("enLLL", &enLLL, "enLLL/F");
937  trainingSample->Branch("enLLT", &enLLT, "enLLT/F");
938  trainingSample->Branch("epiLLL", &epiLLL, "epiLLL/F");
939  trainingSample->Branch("epiLLT", &epiLLT, "epiLLT/F");
940  trainingSample->Branch("gap", &gap, "gap/F");
941  trainingSample->Branch("pi0mass", &pi0mass, "pi0mass/F");
942  trainingSample->Branch("vtxgev", &vtxgev, "vtxgev/F");
943  trainingSample->Branch("shE", &shE, "shE/F");
944  trainingSample->Branch("nueRecEnergy", &nueRecEnergy, "nueRecEnergy/F");
945  trainingSample->Branch("costheta", &costheta, "costheta/F");
946  trainingSample->Branch("length", &length, "length/F");
947  trainingSample->Branch("sbtype", &type, "sbtype/I");
948  trainingSample->Branch("specType", &specType, "specType/C");
949 
950  }
951 
952 
953 
962 
964  normSig = 0, normBkg = 0, totBkg = 0;
965  //
966  // Signal
967  //
968  type = 1;
969  std::cout << "signalIn->GetEntries() = " << signalIn->GetEntries() << std::endl;
970  for (int i = 0; i < signalIn->GetEntries(); i++) {
971 
972  signalIn->GetEntry(i);
973  //
974  // If not QE interaction AND NuE AND not selected by reweight for oscillation then skip event
975  //
976  if (!(evtTrueNuCCNC == kQE && evtTrueNuPdg == kIDNuE &&
977  (evtRwRand[0] > 0 && evtRwRand[0] < evtRwProb[0]))) continue;
978  if (i % 2 == 1) numSignalCorrectProcess++;
979  //
980  // If leading shower is not and electron or positron then skip event
981  //
982  if (fabs(evtSh1TruePdg) != kIDElectron) continue;
983  //
984  // If the angle between the leading shower and the true Nue is greater than 30 degrees,
985  // (ie., shower not matched to true MC NuE in MC) then skip
986  //
987  if (evtSh1TrueDang > 30) continue;
988  //
989  // If leading shower energy not positive, then skip
990  //
991  if (!(evtSh1Energy > 0.0))continue;
992 
993  if (i % 2 == 1) numSignalCorrectShower++;
994  //
995  // Cut out shower candidates that are obvious muons
996  //
997  if (evtSh1DedxLLL[0] - evtSh1DedxLLL[2]<-0.2)continue;
998  if (i % 2 == 1) numSignalCorrectShowerNotObviousMuon++;
999 
1000  //
1001  // Cut events to get desired proportion of signal
1002  //
1003  if ((double) i > (double) signalIn->GetEntries() * rsig)continue;
1004 
1005  if (fabs(evtSh1TruePdg) != 11) {
1006  specType = "BeamNuE\0";
1007  if (evtRwRand[0] < evtRwProb[0]) {
1008  reWeight = 1.0;
1009  } else {
1010  reWeight = 0.0;
1011  }
1012  } else {
1013  specType = "Signal\0";
1014  if (evtRwRand[0] < evtRwProb[0]) {
1015  reWeight = 1.0;
1016  } else {
1017  reWeight = 0.0;
1018  }
1019  }
1020  // std::cout << "SIG specType = " << specType << std::endl;
1021 
1024  if (fillTrainingSample) trainingSample->Fill();
1025  // printTrainingVars();
1026 
1027  normSig++;
1028  }
1029  std::cout << "normSig = " << normSig << std::endl;
1030  //
1031  // Background
1032  //
1033  type = 0;
1034  std::cout << "backgroundIn->GetEntries() = " << backgroundIn->GetEntries() << std::endl;
1035  for (int i = 0; i < backgroundIn->GetEntries(); i++) {
1036 
1037 
1038 
1039 // if (iEvt % 2 == 0)continue;
1040 // backgroundIn->GetEntry(iEvt);
1041 // if (!((fabs(evtTrueNuPdg) == 14 && evtRwRand[2] > 0 && evtRwRand[2] < evtRwProb[2])
1042 // || (fabs(evtTrueNuPdg) == 12 && evtRwRand[1] > 0 && evtRwRand[1] < evtRwProb[1])))continue;
1043 // if (!(evtSh1NPlane>-15 && evtSh1Energy > 0.0))continue;
1044 // if (evtSh1DedxLLL[0] - evtSh1DedxLLL[2]<-0.2)continue;
1045 // if (!((evtSh1Energy + 0.282525 + 1.0766 * (evtEtot - evtSh1Energy)) / 1.015 > Elow
1046 // && (evtSh1Energy + 0.282525 + 1.0766 * (evtEtot - evtSh1Energy)) / 1.015 < Ehigh))continue;
1047 
1048 
1049 
1050  backgroundIn->GetEntry(i);
1051  //
1052  // If not selected for NuMu OR NuE oscillation reweight then skip
1053  //
1054  if (!((fabs(evtTrueNuPdg) == kIDNuMu && evtRwRand[2] > 0 && evtRwRand[2] < evtRwProb[2])
1055  || (fabs(evtTrueNuPdg) == kIDNuE && evtRwRand[1] > 0 && evtRwRand[1] < evtRwProb[1])))continue;
1056  if (i % 2 == 1) numBackgroundNuMuOrNuE0++;
1057  //
1058  // If not NuMu or NuE event
1059  //
1060  if (!((fabs(evtTrueNuPdg) == kIDNuMu) || (fabs(evtTrueNuPdg) == kIDNuE)))continue;
1061  if (i % 2 == 1 && evtTrueNuCCNC == 0 && fabs(evtTrueNuPdg) == kIDNuMu
1062  && evtRwRand[2] > 0 && evtRwRand[2] < evtRwProb[2]) numBackgroundNumuCC0++;
1063  if (i % 2 == 1 && evtTrueNuCCNC == 1 && ((fabs(evtTrueNuPdg) == kIDNuMu)
1064  || (fabs(evtTrueNuPdg) == kIDNuE && evtRwRand[1] > 0 && evtRwRand[1] < evtRwProb[1]))) numBackgroundNC0++;
1065  if (i % 2 == 1 && evtTrueNuCCNC == 0 && fabs(evtTrueNuPdg) == kIDNuE
1066  && evtRwRand[1] > 0 && evtRwRand[1] < evtRwProb[1]) numBackgroundNueCC0++;
1067  //
1068  // Leading shower energy > 0, number of planes > -15
1069  //
1070  if (!(evtSh1NPlane>-15 && evtSh1Energy > 0.0))continue;
1071  if (i % 2 == 1) numBackgroundPassingFirstShCuts++;
1072 
1073  if (i % 2 == 1 && evtTrueNuCCNC == 0 && fabs(evtTrueNuPdg) == kIDNuMu
1075  if (i % 2 == 1 && evtTrueNuCCNC == 1 && ((fabs(evtTrueNuPdg) == kIDNuMu
1076  && evtRwRand[2] > 0 && evtRwRand[2] < evtRwProb[2])
1077  || (fabs(evtTrueNuPdg) == kIDNuE && evtRwRand[1] > 0 && evtRwRand[1] < evtRwProb[1])))
1079  if (i % 2 == 1 && evtTrueNuCCNC == 0 && fabs(evtTrueNuPdg) == kIDNuE && evtRwRand[1] > 0
1081  //
1082  // Cut out shower candidates that are obvious muons
1083  //
1084  if (evtSh1DedxLLL[0] - evtSh1DedxLLL[2]<-0.2) continue;
1085 
1087  if (i % 2 == 1 && evtTrueNuCCNC == 0 && fabs(evtTrueNuPdg) == kIDNuMu && evtRwRand[2] > 0
1089  if (i % 2 == 1 && evtTrueNuCCNC == 1 && ((fabs(evtTrueNuPdg) == kIDNuMu && evtRwRand[2] > 0
1090  && evtRwRand[2] < evtRwProb[2]) || (fabs(evtTrueNuPdg) == kIDNuE && evtRwRand[1] > 0
1092  if (i % 2 == 1 && evtTrueNuCCNC == 0 && fabs(evtTrueNuPdg) == kIDNuE && evtRwRand[1] > 0
1094  // //
1095  // // Skip if quasi-elastic NuE
1096  // //
1097  // if (evtTrueNuCCNC == kQE && fabs(evtTrueNuPdg) == kIDNuE) continue;
1098  // // nb3++;
1099 
1100  //
1101  // Cut events to get desired proportion of background
1102  //
1103  if ((double) i > (double) backgroundIn->GetEntries() * rbg)continue;
1104 
1105 
1106  // For CSV file, write out CCQE NuE (beam Nue) before skipping over
1107  if (evtTrueNuCCNC == kQE && fabs(evtTrueNuPdg) == kIDNuE) {
1108  specType = "CCQEE\0";
1109  if (evtRwRand[1] < evtRwProb[1]) {
1110  reWeight = 1.0;
1111  } else {
1112  reWeight = 0.0;
1113  }
1114  }
1115 
1116  //
1117  // Skip if quasi-elastic NuE
1118  //
1119  if (evtTrueNuCCNC == kQE && fabs(evtTrueNuPdg) == kIDNuE) continue;
1120  // nb3++;
1121 
1122  if (evtTrueNuCCNC == 1 && ((fabs(evtTrueNuPdg) == kIDNuMu)
1123  || (fabs(evtTrueNuPdg) == kIDNuE))) {
1124  specType = "NCMuOrE\0";
1125  if (fabs(evtTrueNuPdg) == kIDNuE) {
1126  if (evtRwRand[1] < evtRwProb[1] ) {
1127  reWeight = 1.0;
1128  } else {
1129  reWeight = 0.0;
1130  }
1131  } else if (fabs(evtTrueNuPdg) == kIDNuMu) {
1132  if (evtRwRand[2] < evtRwProb[2]) {
1133  reWeight = 1.0;
1134  } else {
1135  reWeight = 0.0;
1136  }
1137  }
1138 
1139  } else if (evtTrueNuCCNC == kQE && fabs(evtTrueNuPdg) == kIDNuMu) {
1140  specType = "CCQEMu\0";
1141  if (evtRwRand[2] < evtRwProb[2]) {
1142  reWeight = 1.0;
1143  } else {
1144  reWeight = 0.0;
1145  }
1146 
1147  }
1148  // std::cout << "BG specType = " << specType << std::endl;
1151 
1152  if (fillTrainingSample) trainingSample->Fill();
1153  normBkg++;
1154  }
1155  std::cout << "normBkg = " << normBkg << std::endl;
1157 
1158  std::cout << "Nsig, Nbig " << normSig << " " << normBkg << " " << std::endl;
1159  scalepotsig = 18.0e+20 / ((numSignalCorrectProcess / (2.87128765282905034e-02 * 6938.035))*2.99649003626026520e+21); //2.38399959425089536e+19);//NDOS
1160  scalepotNumuCC = 18.0e+20 / ((numBackgroundNumuCC0 / (6938.035 - 199.21 - 5511.0658))*3.02657807308939830e+21); //2.38399959425089536e+19);//NDOS
1161  scalepotNueCC = 18.0e+20 / ((numBackgroundNueCC0 / 82.02)*3.02657807308939830e+21); //2.38399959425089536e+19);//NDOS
1162  scalepotNC = 18.0e+20 / ((numBackgroundNC0 / 2664.42968)*3.02657807308939830e+21); //2.38399959425089536e+19);//NDOS
1163  //
1164  // Write scale factors out to a CSV file to be used in R
1165  //
1166  std::ofstream potScaleFile;
1167  TString potScaleFileName = "POTScales.csv";
1168  potScaleFile.open((potScaleFileName.Data()));
1169  // Header
1170  potScaleFile << "scalepotsig\t" << "scalepotNumuCC\t" << "scalepotNueCC\t" << "scalepotNC" << std::endl;
1171  potScaleFile << scalepotsig << "\t" << scalepotNumuCC << "\t" << scalepotNueCC << "\t" << scalepotNC << std::endl;
1172  potScaleFile.close();
1173 
1174 
1175 
1176  trainingSample->Write("trainingSample");
1177 }
1178 
1179 void TMVAEIDTraining::setTreeBranchAddresses(TTree* treeIncoming) {
1180 
1181  treeIncoming->SetBranchAddress("evtTrueNuCCNC", &evtTrueNuCCNC);
1182  treeIncoming->SetBranchAddress("evtTrueNuMode", &evtTrueNuMode);
1183  treeIncoming->SetBranchAddress("evtTrueNuPdg", &evtTrueNuPdg);
1184  treeIncoming->SetBranchAddress("evtTrueNuEnergy", &evtTrueNuEnergy);
1185  treeIncoming->SetBranchAddress("evtTrueNuIsFid", &evtTrueNuIsFid);
1186  treeIncoming->SetBranchAddress("evtTrueNuP4[4]", evtTrueNuP4);
1187  treeIncoming->SetBranchAddress("evtTrueNuVtx[3]", evtTrueNuVtx);
1188  treeIncoming->SetBranchAddress("evtOsTag", &evtOsTag);
1189  treeIncoming->SetBranchAddress("evtOsProb[20]", evtOsProb);
1190  treeIncoming->SetBranchAddress("evtOsRand", &evtOsRand);
1191  treeIncoming->SetBranchAddress("evtRwProb[20]", evtRwProb);
1192  treeIncoming->SetBranchAddress("evtRwRand[20]", evtRwRand);
1193  treeIncoming->SetBranchAddress("evtTruePdg[20]", evtTruePdg);
1194  treeIncoming->SetBranchAddress("evtTrueMother[20]", evtTrueMother);
1195  treeIncoming->SetBranchAddress("evtTrueEnergy[20]", evtTrueEnergy);
1196  treeIncoming->SetBranchAddress("evtSh1TruePdg", &evtSh1TruePdg);
1197  treeIncoming->SetBranchAddress("evtSh1TrueDang", &evtSh1TrueDang);
1198  treeIncoming->SetBranchAddress("evtSh1TrueEDang", &evtSh1TrueEDang);
1199  treeIncoming->SetBranchAddress("evtSh1Energy", &evtSh1Energy);
1200  treeIncoming->SetBranchAddress("evtSh1Energy0", &evtSh1Energy0);
1201  treeIncoming->SetBranchAddress("evtTotalRecoGeV", &evtTotalRecoGeV);
1202  treeIncoming->SetBranchAddress("evtSh1SliceGeV", &evtSh1SliceGeV);
1203  treeIncoming->SetBranchAddress("evtSheEnergy", &evtSheEnergy);
1204  treeIncoming->SetBranchAddress("evtSh1Gap", &evtSh1Gap);
1205  treeIncoming->SetBranchAddress("evtSh1IsFid", &evtSh1IsFid);
1206  treeIncoming->SetBranchAddress("evtSh1P4[4]", evtSh1P4);
1207  treeIncoming->SetBranchAddress("evtSh1Start[3]", evtSh1Start);
1208  treeIncoming->SetBranchAddress("evtSh1DedxLLL[20]", evtSh1DedxLLL);
1209  treeIncoming->SetBranchAddress("evtSh1DedxLLT[20]", evtSh1DedxLLT);
1210  treeIncoming->SetBranchAddress("evtSh1NPlane", &evtSh1NPlane);
1211  treeIncoming->SetBranchAddress("evtSh1NMIPPlane", &evtSh1NMIPPlane);
1212  treeIncoming->SetBranchAddress("evtSh1ContStartPlane", &evtSh1ContStartPlane);
1213  treeIncoming->SetBranchAddress("evtEtot", &evtEtot);
1214  treeIncoming->SetBranchAddress("evtSh1SumPt", &evtSh1SumPt);
1215  treeIncoming->SetBranchAddress("evtNTrk", &evtNTrk);
1216  treeIncoming->SetBranchAddress("evtSh1Enue", &evtSh1Enue);
1217  treeIncoming->SetBranchAddress("evtSh2Energy", &evtSh2Energy);
1218  treeIncoming->SetBranchAddress("evtPi0Mgg", &evtPi0Mgg);
1219  treeIncoming->SetBranchAddress("evtSh1Pi0Mgg", &evtSh1Pi0Mgg);
1220  treeIncoming->SetBranchAddress("evtSh1VtxGeV", &evtSh1VtxGeV);
1221  treeIncoming->SetBranchAddress("evtSh1CosTheta", &evtSh1CosTheta);
1222  treeIncoming->SetBranchAddress("evtSh1Length", &evtSh1Length);
1223  treeIncoming->SetBranchAddress("evtPi0StartDist", &evtPi0StartDist);
1224  treeIncoming->SetBranchAddress("evtPi0LineDist", &evtPi0LineDist);
1225 }
1226 
1228 
1229  egLLL = (float) evtSh1DedxLLL[0] - (float) evtSh1DedxLLL[1];
1230  egLLT = (float) evtSh1DedxLLT[0] - (float) evtSh1DedxLLT[1];
1231  emuLLL = (float) evtSh1DedxLLL[0] - (float) evtSh1DedxLLL[2];
1232  emuLLT = (float) evtSh1DedxLLT[0] - (float) evtSh1DedxLLT[2];
1233  epi0LLL = (float) evtSh1DedxLLL[0] - (float) evtSh1DedxLLL[3];
1234  epi0LLT = (float) evtSh1DedxLLT[0] - (float) evtSh1DedxLLT[3];
1235  epLLL = (float) evtSh1DedxLLL[0] - (float) evtSh1DedxLLL[5];
1236  epLLT = (float) evtSh1DedxLLT[0] - (float) evtSh1DedxLLT[5];
1237  enLLL = (float) evtSh1DedxLLL[0] - (float) evtSh1DedxLLL[6];
1238  enLLT = (float) evtSh1DedxLLT[0] - (float) evtSh1DedxLLT[6];
1239  epiLLL = (float) evtSh1DedxLLL[0] - (float) evtSh1DedxLLL[7];
1240  epiLLT = (float) evtSh1DedxLLT[0] - (float) evtSh1DedxLLT[7];
1241  gap = (float) evtSh1Gap;
1242  pi0mass = TMath::Max((float) evtSh1Pi0Mgg, (float) 0.0);
1244  // costheta = evtSh1SumPt/evtSh1SliceGeV;
1246  length = (float) evtSh1Length / (float) evtSh1Energy;
1247  shE = (float) evtSh1Energy / (float) evtSh1SliceGeV;
1248  nueRecEnergy = ((float) evtSh1Energy + 0.282525 + 1.0766 * ((float) evtEtot - (float) evtSh1Energy));
1249 
1250  ntrk = (float) evtNTrk;
1253 }
1254 
1256  // Launch the GUI for the root macros
1257  // if (!gROOT->IsBatch()) TMVAGui(outfileName);
1258 }
1259 
1261 
1262  // currentAlgorithm = "MLP method";
1263  // currentAlgorithm = "TMlpANN method";
1264  // currentAlgorithm = "KNN method";
1265  // currentAlgorithm = "BDT method";
1266 
1268 
1269  // TMVAEIDTraining::defineEvaluationTree();
1270 
1271  double Elow = 1.0;
1272  double Ehigh = 3.0;
1273  double mVAPlotLow = -1.5;
1274  double mVAPlotHigh = 1.5;
1275  int nBinsANN = 50;
1276  std::cout << std::endl;
1277  std::cout << "==> Start TMVAClassificationApplication" << std::endl;
1278 
1279  if (!methodsBooked) {
1281  }
1282 
1283 
1284  // Select methods (don't look at this code - not of interest)
1285  if (myMethodList != "") {
1286  for (std::map<std::string, int>::iterator it = Use.begin(); it != Use.end(); it++) it->second = 0;
1287 
1288  std::vector<TString> mlist = TMVA::gTools().SplitString(myMethodList, ',');
1289  for (UInt_t i = 0; i < mlist.size(); i++) {
1290  std::string regMethod(mlist[i]);
1291 
1292  if (Use.find(regMethod) == Use.end()) {
1293  std::cout << "Method \"" << regMethod
1294  << "\" not known in TMVA under this name. Choose among the following:" << std::endl;
1295  for (std::map<std::string, int>::iterator it = Use.begin(); it != Use.end(); it++) {
1296  std::cout << it->first << " ";
1297  }
1298  std::cout << std::endl;
1299  return;
1300  }
1301  Use[regMethod] = 1;
1302  currentAlgorithm = regMethod + " method";
1303  if (regMethod == "TMlpANN") {
1304  nBinsANN = 50;
1305  mVAPlotLow = -0.5;
1306  mVAPlotHigh = 1.5;
1307  } else if (regMethod == "MLP") {
1308  nBinsANN = 50;
1309  mVAPlotLow = -0.5;
1310  mVAPlotHigh = 1.5;
1311  } else if (regMethod == "BDT") {
1312  nBinsANN = 100;
1313  mVAPlotLow = -1.0;
1314  mVAPlotHigh = 1.0;
1315  } else if (regMethod == "KNN") {
1316  nBinsANN = 100;
1317  mVAPlotLow = -1.5;
1318  mVAPlotHigh = 1.5;
1319  }
1320  }
1321  }
1322 
1323 
1324  // double rsig = 1.0;
1325  // double rbg = 0.9;
1326  // double trueNuE, etot, enue, sh1Energy, ANN, wSig, wBkg;
1327 
1328  // Use the NN to plot the results for each sample
1329  // This will give approx. the same result as DrawNetwork.
1330  // All entries are used, while DrawNetwork focuses on
1331  // the test sample. Also the xaxis range is manually set.
1332  bg = new TH1F("bgh", "Background", nBinsANN, mVAPlotLow, mVAPlotHigh);
1333  bgNumuCC = new TH1F("bgNumuCC", "NuMu CC Background", nBinsANN, mVAPlotLow, mVAPlotHigh);
1334  bgNC = new TH1F("bgNC", "Neutral Current Background", nBinsANN, mVAPlotLow, mVAPlotHigh);
1335  bgNueCC = new TH1F("bgNueCC", "Beam NuE Charged Current Background", nBinsANN, mVAPlotLow, mVAPlotHigh);
1336 
1337  sig = new TH1F("sigh", "Signal Nue CC QE", nBinsANN, mVAPlotLow, mVAPlotHigh);
1338  sigbg = new TH1F("sigbgh", "Mis-reconstructed Signal \"Background\"", nBinsANN, mVAPlotLow, mVAPlotHigh);
1339  TCanvas* mlpa_canvas2 = new TCanvas("mlpa_canvas2", "Network analysis 2", 10, 10, 1200, 600);
1340  TLine *line1, *line2;
1341 
1342  //Int_t type;
1343 
1344 
1345  // --------------------------------------------------------------------------------------------------
1346 
1347  // --- Create the Reader object
1348  std::cout << "Defining MTVA::Reader object..." << std::endl;
1349  TMVA::Reader reader("!Color:!Silent");
1350  std::cout << "Adding variables to Reader..." << std::endl;
1351 
1352 
1353 
1354  reader.AddVariable("egLLL", &egLLL);
1355  reader.AddVariable("egLLT", &egLLT);
1356  reader.AddVariable("emuLLL", &emuLLL);
1357  reader.AddVariable("emuLLT", &emuLLT);
1358  reader.AddVariable("epi0LLL", &epi0LLL);
1359  reader.AddVariable("epi0LLT", &epi0LLT);
1360  reader.AddVariable("epLLL", &epLLL);
1361  reader.AddVariable("epLLT", &epLLT);
1362  reader.AddVariable("enLLL", &enLLL);
1363  reader.AddVariable("enLLT", &enLLT);
1364  reader.AddVariable("epiLLL", &epiLLL);
1365  reader.AddVariable("epiLLT", &epiLLT);
1366  reader.AddVariable("gap", &gap);
1367  reader.AddVariable("pi0mass", &pi0mass);
1368  reader.AddVariable("vtxgev", &vtxgev);
1369  reader.AddVariable("shE", &shE);
1370  if (useENu) {
1371  reader.AddVariable("nueRecEnergy", &nueRecEnergy);
1372  }
1373  reader.AddSpectator("sbtype", &type);
1374  // reader.AddSpectator("specType", &specType);
1375  // reader.AddSpectator("reWeight", &reWeight);
1376 
1377  // --- Book the MVA methods
1378 
1379  // TString dir = "weights/";
1380  TString prefix = "TMVAClassification";
1381 
1382  // Book method(s)
1383  std::cout << "Booking TMVA methods..." << std::endl;
1384  for (std::map<std::string, int>::iterator it = Use.begin(); it != Use.end(); it++) {
1385  if (it->second) {
1386  TString methodName = TString(it->first) + TString(" method");
1387  TString weightfile = weightDir + prefix + TString("_") + TString(it->first) + TString(".weights.xml");
1388  reader.BookMVA(methodName, weightfile);
1389  }
1390  }
1391 
1392 
1393 
1394 
1397 
1398 
1399  std::cout << "Looping over background input tree to fill hists...Algorithm = " << currentAlgorithm << std::endl;
1400  std::cout << "backgroundIn->GetEntries() = " << backgroundIn->GetEntries() << std::endl;
1401  int bgcnt = 0;
1402  for (int iEvt = 0; iEvt < backgroundIn->GetEntries(); iEvt++) {
1403  // if((double)i>(double)backgroundIn->GetEntries()*rbg)continue;
1404  if (iEvt % 2 == 0)continue;
1405  // if(i>36383)continue;
1406  backgroundIn->GetEntry(iEvt);
1407  if (!((fabs(evtTrueNuPdg) == 14 && evtRwRand[2] > 0 && evtRwRand[2] < evtRwProb[2])
1408  || (fabs(evtTrueNuPdg) == 12 && evtRwRand[1] > 0 && evtRwRand[1] < evtRwProb[1])))continue;
1409  if (!(evtSh1NPlane>-15 && evtSh1Energy > 0.0))continue;
1410  if (evtSh1DedxLLL[0] - evtSh1DedxLLL[2]<-0.2)continue;
1411  if (!((evtSh1Energy + 0.282525 + 1.0766 * (evtEtot - evtSh1Energy)) / 1.015 > Elow
1412  && (evtSh1Energy + 0.282525 + 1.0766 * (evtEtot - evtSh1Energy)) / 1.015 < Ehigh))continue;
1413 
1414  //
1415  // Fill training variables
1416  //
1418 
1419  if (evtTrueNuCCNC == 0 && fabs(evtTrueNuPdg) == 14 && evtRwRand[2] > 0 && evtRwRand[2] < evtRwProb[2]) {
1420  bgNumuCC->Fill(reader.EvaluateMVA(currentAlgorithm), scalepotNumuCC);
1421  bg->Fill(reader.EvaluateMVA(currentAlgorithm), scalepotNumuCC);
1422  }
1423  if (evtTrueNuCCNC == 1 && ((fabs(evtTrueNuPdg) == 14 && evtRwRand[2] > 0 && evtRwRand[2] < evtRwProb[2])
1424  || (fabs(evtTrueNuPdg) == 12 && evtRwRand[1] > 0 && evtRwRand[1] < evtRwProb[1]))) {
1425  bgNC->Fill(reader.EvaluateMVA(currentAlgorithm), scalepotNC);
1426  bg->Fill(reader.EvaluateMVA(currentAlgorithm), scalepotNC);
1427  }
1428  if (evtTrueNuCCNC == 0 && fabs(evtTrueNuPdg) == 12 && evtRwRand[1] > 0 && evtRwRand[1] < evtRwProb[1]) {
1429  bgNueCC->Fill(reader.EvaluateMVA(currentAlgorithm), scalepotNueCC);
1430  bg->Fill(reader.EvaluateMVA(currentAlgorithm), scalepotNueCC);
1431  }
1432  bgcnt++;
1433  }
1434  std::cout << "bgcnt = " << bgcnt << std::endl;
1435 
1436  std::cout << "Looping over signal input tree to fill hists..." << std::endl;
1437  std::cout << "signalIn->GetEntries() = " << signalIn->GetEntries() << std::endl;
1438  int sgcnt = 0;
1439  for (int iSigEvt = 0; iSigEvt < signalIn->GetEntries(); iSigEvt++) {
1440  if (iSigEvt % 2 == 0)continue;
1441  signalIn->GetEntry(iSigEvt);
1442  if (!(evtTrueNuCCNC == 0 && evtTrueNuPdg == 12 && (evtRwRand[0] > 0 && evtRwRand[0] < evtRwProb[0])))continue;
1443  // if(evtSh1TrueDang>10||fabs(evtSh1TruePdg)!=11)continue;
1444  if (!(evtSh1NPlane>-15 && evtSh1Energy > 0.0))continue;
1445  if (evtSh1DedxLLL[0] - evtSh1DedxLLL[2]<-0.2)continue;
1446  if (!((evtSh1Energy + 0.282525 + 1.0766 * (evtEtot - evtSh1Energy)) / 1.015 > Elow
1447  && (evtSh1Energy + 0.282525 + 1.0766 * (evtEtot - evtSh1Energy)) / 1.015 < Ehigh))continue;
1448 
1449  // fill training variables from Tree entry
1451  sig->Fill(reader.EvaluateMVA(currentAlgorithm), scalepotsig);
1452  if (fabs(evtSh1TruePdg) != 11)sigbg->Fill(reader.EvaluateMVA(currentAlgorithm), scalepotsig);
1453  sgcnt++;
1454  }
1455  std::cout << "sgcnt = " << sgcnt << std::endl;
1456 
1457 
1458  hfom = new TH1F("hfom", "Figure of Merit", nBinsANN, mVAPlotLow, mVAPlotHigh);
1459  hfomBeamNueBG = new TH1F("hfomBeamNueBG", "Figure of Merit, Beam Nue as BG", nBinsANN, mVAPlotLow, mVAPlotHigh);
1460  heff = new TH1F("heff", "Signal Efficiency", nBinsANN, mVAPlotLow, mVAPlotHigh);
1461  hbgeff = new TH1F("hbgeff", "Mis-reconstructed Signal \"Background\" Efficiency", nBinsANN, mVAPlotLow, mVAPlotHigh);
1462  hbgNumuCCeff = new TH1F("hbgNumuCCeff", "NuMu CC Background Efficiency", nBinsANN, mVAPlotLow, mVAPlotHigh);
1463  hbgNCeff = new TH1F("hbgNCeff", "Neutral Current Background Efficiency", nBinsANN, mVAPlotLow, mVAPlotHigh);
1464  hbgNueCCeff = new TH1F("hbgNueCCeff", "Beam NuE Charged Current Background Efficiency", nBinsANN, mVAPlotLow, mVAPlotHigh);
1465  for (int iBin = 0; iBin < nBinsANN; iBin++) {
1466  double sumnbg = bg->Integral(iBin + 1, nBinsANN);
1467  double sumnbgNueCC = bgNueCC->Integral(iBin + 1, nBinsANN);
1468  double sumnsig = sig->Integral(iBin + 1, nBinsANN);
1469  double eff = double(sig->Integral(iBin + 1, nBinsANN)) / double(sig->Integral(1, nBinsANN));
1470  double bgeff = double(bg->Integral(iBin + 1, nBinsANN)) / double(bg->Integral(1, nBinsANN));
1471  double bgNumuCCeff = double(bgNumuCC->Integral(iBin + 1, nBinsANN)) / double(bgNumuCC->Integral(1, nBinsANN));
1472  double bgNCeff = double(bgNC->Integral(iBin + 1, nBinsANN)) / double(bgNC->Integral(1, nBinsANN));
1473  double bgNueCCeff = double(bgNueCC->Integral(iBin + 1, nBinsANN)) / double(bgNueCC->Integral(1, nBinsANN));
1474  heff->SetBinContent(iBin + 1, eff);
1475  hbgeff->SetBinContent(iBin + 1, bgeff);
1476  hbgNumuCCeff->SetBinContent(iBin + 1, bgNumuCCeff);
1477  hbgNCeff->SetBinContent(iBin + 1, bgNCeff);
1478  hbgNueCCeff->SetBinContent(iBin + 1, bgNueCCeff);
1479  //
1480  // Fill hfom
1481  //
1482  if (sumnsig + sumnbg - sumnbgNueCC < 0.000001) {
1483  hfom->SetBinContent(iBin + 1, 0);
1484  continue;
1485  }
1486  double significance = 0;
1487  if (sumnsig + sumnbg - sumnbgNueCC > 0) significance = double(sumnsig) / sqrt(sumnsig + sumnbg - sumnbgNueCC);
1488  std::cout << "FOM, sig, bkg " << significance << " " << sumnsig << " " << sumnbg << " " << std::endl;
1489  hfom->SetBinContent(iBin + 1, significance);
1490  //
1491  // Fill hfomBeamNueBG
1492  //
1493  if (sumnsig + sumnbg < 0.000001) {
1494  hfomBeamNueBG->SetBinContent(iBin + 1, 0);
1495  continue;
1496  }
1497  double significanceV2 = 0;
1498 
1499  if (sumnsig + sumnbg > 0) significanceV2 = double(sumnsig) / sqrt(sumnsig + sumnbg);
1500  std::cout << "FOM(Beam Nue as BG), sig, bkg " << significanceV2 << " " << sumnsig << " " << sumnbg << " " << std::endl;
1501  hfomBeamNueBG->SetBinContent(iBin + 1, significanceV2);
1502  }
1503 
1504  bg->SetLineColor(kBlue);
1505  sig->SetLineColor(kRed);
1506  sig->SetFillStyle(3003);
1507  sig->SetFillColor(kRed);
1508 
1509  bg->Draw();
1510  sig->Draw("same");
1511  sigbg->Draw("same");
1512  // sig->DrawNormalized("same", sig->GetEntries()*scalesig);
1513  TLegend *legend = new TLegend(.75, .80, .95, .95);
1514  legend->AddEntry(bg, "Background (nhc)");
1515  legend->AddEntry(sig, "Signal (Nue_CCQE)");
1516  legend->Draw();
1517  mlpa_canvas2->cd(0);
1518  mlpa_canvas2->Update();
1519  mlpa_canvas2->SaveAs(plotDir + "JMCanvas2.png");
1520  mlpa_canvas2->SaveAs(plotDir + "JMCanvas2.eps");
1521 
1522  std::cout << "scalepotsig " << scalepotsig << std::endl;
1523  std::cout << "scalepotNumuCC " << scalepotNumuCC << std::endl;
1524  std::cout << "scalepotNC " << scalepotNC << std::endl;
1525  std::cout << "scalepotNueCC " << scalepotNueCC << std::endl;
1526 
1527  TCanvas* mlpa_canvas3 = new TCanvas("mlpa_canvas3", "Network analysis 3", 10, 10, 800, 400);
1528  mlpa_canvas3->Divide(2, 1);
1529  mlpa_canvas3->cd(1);
1530  Double_t maxFOM = hfom->GetBinContent(hfom->GetMaximumBin());
1531  Double_t maxFOMx = hfom->GetBinCenter(hfom->GetMaximumBin());
1532  hfom->GetYaxis()->SetTitle("FOM");
1533  legend = new TLegend(.75, .80, .95, .95);
1534  legend->AddEntry(hfom, TString::Format("Max FOM = %4.4f at %4.2f", maxFOM, maxFOMx));
1535  hfom->GetXaxis()->SetTitle("ANN");
1536  // hfom->GetXaxis()->SetRangeUser(0.4, 1.2);
1537  hfom->Draw("l");
1538  hfomBeamNueBG->SetLineColor(kRed);
1539  hfomBeamNueBG->Draw("lsame");
1540  line1 = new TLine(hfom->GetBinCenter(hfom->GetMaximumBin()), 0, hfom->GetBinCenter(hfom->GetMaximumBin()), hfom->GetMaximum());
1541  line1->SetLineWidth(4);
1542  line1->SetLineColor(kRed);
1543  line1->Draw("same");
1544 
1545 
1546  mlpa_canvas3->Update();
1547 
1548  mlpa_canvas3->cd(2);
1549  heff->GetYaxis()->SetTitle("PID Eff.");
1550  heff->GetXaxis()->SetTitle("ANN");
1551  heff->GetYaxis()->SetRangeUser(0, 1);
1552  // heff->GetXaxis()->SetRangeUser(0.4, 1.2);
1553  heff->Draw("l");
1554  line2 = new TLine(hfom->GetBinCenter(hfom->GetMaximumBin()), 0, hfom->GetBinCenter(hfom->GetMaximumBin()), 1.0);
1555  line2->SetLineWidth(4);
1556  line2->SetLineColor(kRed);
1557  line2->Draw("same");
1558 
1559  mlpa_canvas3->Update();
1560  mlpa_canvas3->SaveAs(plotDir + "JMCanvas3.png");
1561  mlpa_canvas3->SaveAs(plotDir + "JMCanvas3.eps");
1562 
1563  TCanvas* mlpa_canvas4 = new TCanvas("mlpa_canvas4", "Network analysis 4", 10, 10, 600, 400);
1564  mlpa_canvas4->Divide(1, 1);
1565  mlpa_canvas4->cd(1);
1566  bg->GetYaxis()->SetTitle("Events");
1567  bg->GetXaxis()->SetTitle("ANN");
1568  bg->SetLineColor(kBlue);
1569  bgNumuCC->SetLineColor(9);
1570  bgNC->SetLineColor(7);
1571  bgNueCC->SetLineColor(6);
1572  sig->SetLineColor(kRed);
1573  sig->SetFillStyle(3003);
1574  sig->SetFillColor(kRed);
1575 
1576  bg->Draw();
1577  bgNumuCC->Draw("same");
1578  bgNC->Draw("same");
1579  bgNueCC->Draw("same");
1580  sig->Draw("same");
1581  // sig->DrawNormalized("same", sig->GetEntries()*scalesig);
1582  legend = new TLegend(.75, .80, .95, .95);
1583  legend->AddEntry(bg, "Background (nhc)");
1584  legend->AddEntry(bgNumuCC, "Beam #nu_{#mu} CC");
1585  std::cout << "bgNC = " << bgNC << std::endl;
1586  legend->AddEntry(bgNC, "Beam NC");
1587  legend->AddEntry(bgNueCC, "Beam #nu_{e} CC");
1588  legend->AddEntry(sig, "Signal #nu_{e} CC");
1589  legend->Draw();
1590  mlpa_canvas4->Update();
1591  mlpa_canvas4->SaveAs(plotDir + "JMCanvas4.png");
1592  mlpa_canvas4->SaveAs(plotDir + "JMCanvas4.eps");
1593 
1594  // Get elapsed time
1595  sw.Stop();
1596  std::cout << "--- End of event loop: ";
1597  sw.Print();
1598 
1600 
1601  std::cout << "ns0,ns1,ns2,ns3,scalepotsig " << numSignalCorrectProcess << " " << numSignalCorrectShower << " " << numSignalCorrectShowerNotObviousMuon << " " << ns3 << " " << scalepotsig << std::endl;
1602  double Dns0 = double(numSignalCorrectProcess * scalepotsig);
1603 // double Dns1 = double(numSignalCorrectShower * scalepotsig);
1604 // double Dns2 = double(numSignalCorrectShowerNotObviousMuon * scalepotsig);
1605  // double Dns3 = double(ns3 * scalepotsig);
1606 
1607 // double DnbNumuCC0 = double(numBackgroundNumuCC0 * scalepotNumuCC);
1608 // double DnbNumuCC1 = double(numBackgroundPassingFirstShCutsNumuCC1 * scalepotNumuCC);
1609 // double DnbNumuCC2 = double(numBackgroundPassingNotObviousMuonsShCutsNumuCC2 * scalepotNumuCC);
1610  // double DnbNumuCC3 = double(nbNumuCC3 * scalepotNumuCC);
1611 
1612 // double DnbNC0 = double(numBackgroundNC0 * scalepotNC);
1613 // double DnbNC1 = double(numBackgroundPassingFirstShCutsNC1 * scalepotNC);
1614 // double DnbNC2 = double(numBackgroundPassingNotObviousMuonsShCutsNC2 * scalepotNC);
1615  // double DnbNC3 = double(nbNC3 * scalepotNC);
1616 
1617 // double DnbNueCC0 = double(numBackgroundNueCC0 * scalepotNueCC);
1618 // double DnbNueCC1 = double(numBackgroundPassingFirstShCutsNueCC1 * scalepotNueCC);
1619 // double DnbNueCC2 = double(numBackgroundPassingNotObviousMuonsShCutsNueCC2 * scalepotNueCC);
1620  // double DnbNueCC3 = double(nbNueCC3 * scalepotNueCC);
1621 
1622 // double Dnb0 = double(DnbNumuCC0 + DnbNC0 + DnbNueCC0);
1623 // double Dnb1 = double(DnbNumuCC1 + DnbNC1 + DnbNueCC1);
1624 // double Dnb2 = double(DnbNumuCC2 + DnbNC2 + DnbNueCC2);
1625  // double Dnb3 = double(DnbNumuCC3 + DnbNC3 + DnbNueCC3);
1626 
1627  // double scalepot = 18.0e+20/((nbNumuCC0/7175.2)*2.8123e+21);//2.38399959425089536e+19);//NDOS
1628  std::cout << "ANN = " << hfom->GetBinCenter(hfom->GetMaximumBin()) << std::endl;
1629  std::cout << "FOM hist = " << hfom->GetBinContent(hfom->GetMaximumBin()) << std::endl;
1630  // std::cout << "Dns3 = "<< Dns3 << std::endl;
1631  // std::cout << "Dnb3 = "<< Dnb3 << std::endl;
1632  std::cout << "Dns0 = " << Dns0 << std::endl;
1633  // std::cout << "FOM = " << double(Dns3) / sqrt(Dns3 + Dnb3) << std::endl;
1634 
1635  // std::cout << "Eff. = " << double(Dns3) / double(Dns0) << std::endl;
1637  std::cout << "nb2, nb3 " << numBackgroundPassingNotObviousMuonsShCuts2 << " " << nb3 << std::endl;
1638 
1639  // std::cout << "BG Eff . = " << double(Dnb3) / double(Dnb0) << std::endl;
1640  // std::cout << "BG NumuCC Eff. = " << double(DnbNumuCC3) / double(DnbNumuCC0) << std::endl;
1641  // std::cout << "BG NC Eff. = " << double(DnbNC3) / double(DnbNC0) << std::endl;
1642  // std::cout << "BG NueCC Eff. = " << double(DnbNueCC3) / double(DnbNueCC0) << std::endl;
1643 
1644 }
1645 
1647 
1648 
1649  TFile f(histDir + "annDiagnosticHists.root", "recreate");
1650  hbgeff->Write("hbgeff");
1651  hbgNumuCCeff->Write("hbgNumuCCeff");
1652  hbgNCeff->Write("hbgNCeff");
1653  hbgNueCCeff->Write("hbgNueCCeff");
1654  heff->Write("heff");
1655  hfom->Write("hfom");
1656  // hfomBeamNueBG->Write("hfomBeamNueBG");
1657  sigbg->Write("sigbg");
1658  bg->Write("bg");
1659  bgNumuCC->Write("bgNumuCC");
1660  bgNC->Write("bgNC");
1661  bgNueCC->Write("bgNueCC");
1662  sig->Write("sig");
1663  // annMC->Write("annMC");
1664  f.Close();
1665 
1666 }
void dumpTrainingSampleToCSV(CSVOperation operation, TString cSVFileName)
std::map< std::string, int > Use
double evtTrueEnergy[20]
set< int >::iterator it
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
virtual ~TMVAEIDTraining()
int numSignalCorrectShowerNotObviousMuon
T sqrt(T number)
Definition: d0nt_math.hpp:156
TMVA::Reader * reader
void setBackgroundFilesUntransformed(std::vector< TString > inputFileNameIn)
void setupTrainingAndTestTrees()
int numBackgroundPassingFirstShCuts
int numBackgroundPassingNotObviousMuonsShCuts2
TString currentAlgorithm
void runTraining(TString cSVFileNameIn)
double evtTrueNuP4[4]
void setTreeBranchAddresses(TTree *treeIncoming)
int numBackgroundPassingFirstShCutsNumuCC1
TMVA::Factory * factory
int numBackgroundPassingNotObviousMuonsShCutsNC2
int numBackgroundPassingNotObviousMuonsShCutsNueCC2
double evtTrueNuVtx[3]
std::ofstream outfile
double evtSh1DedxLLL[20]
double evtRwRand[20]
void trainTestAndEvaluateTMVAMethods()
std::string specType
double evtRwProb[20]
TMVAEIDTraining(TString jobNameIn, TString outPathIn, TString myMethodListIn, TString useENuIn)
const XML_Char * prefix
Definition: expat.h:380
double evtOsProb[20]
OStream cout
Definition: OStream.cxx:6
double evtSh1DedxLLT[20]
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
double evtSh1Start[3]
int numBackgroundPassingFirstShCutsNueCC1
int numBackgroundPassingNotObviousMuonsShCutsNumuCC2
int numBackgroundPassingFirstShCutsNC1
void setSignalFilesUntransformed(std::vector< TString > inputFileNameIn)
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
void transformInputTrees(bool fillTrainingSample)
Float_t e
Definition: plot.C:35
TMVA::DataLoader * dataloader
double evtSh1ContStartPlane