#include <cstdlib>
#include <iostream>
#include <map>
#include <string>

#include "TChain.h"
#include "TFile.h"
#include "TTree.h"
#include "TString.h"
#include "TObjString.h"
#include "TSystem.h"
#include "TROOT.h"

#include "TMVA/Factory.h"
#include "TMVA/DataLoader.h"
#include "TMVA/Tools.h"
#include "TMVA/TMVAGui.h"

int TMVAClassification_Bmesonfinal(TString myMethodList = "")
{

   // This loads the library
   TMVA::Tools::Instance();

   // Default MVA methods to be trained + tested
   std::map<std::string,int> Use;

   //
   // Neural Networks (all are feed-forward Multilayer Perceptrons)
   Use["MLP"]             = 1; // Recommended ANN
   Use["MLPBFGS"]         = 0; // Recommended ANN with optional training method
   Use["MLPBNN"]          = 0; // Recommended ANN with BFGS training method and bayesian regulator
//   Use["CFMlpANN"]        = 0; // Depreciated ANN from ALEPH
//   Use["TMlpANN"]         = 0; // ROOT's own ANN
//   Use["DNN_GPU"]         = 0; // CUDA-accelerated DNN training.
//   Use["DNN_CPU"]         = 0; // Multi-core accelerated DNN.

   // Boosted Decision Trees
   Use["BDT"]             = 1; // uses Adaptive Boost
   Use["BDTG"]            = 0; // uses Gradient Boost
//   Use["BDTB"]            = 0; // uses Bagging
//   Use["BDTD"]            = 0; // decorrelation + Adaptive Boost
//   Use["BDTF"]            = 0; // allow usage of fisher discriminant for node splitting

   // ------------------------------------------------------------------------------------------------------------------------------------
   //Loop to use the selected methods

   std::cout << std::endl;
   std::cout << "==> Start TMVAClassification_Bmesonfinal" << std::endl;

   // Select methods (don't look at this code - not of interest)
   if (myMethodList != "") {
      for (std::map<std::string,int>::iterator it = Use.begin(); it != Use.end(); it++) it->second = 0;

      std::vector<TString> mlist = TMVA::gTools().SplitString( myMethodList, ',' );
      for (UInt_t i=0; i<mlist.size(); i++) {
         std::string regMethod(mlist[i]);

         if (Use.find(regMethod) == Use.end()) {
            std::cout << "Method \"" << regMethod << "\" not known in TMVA under this name. Choose among the following:" << std::endl;
            for (std::map<std::string,int>::iterator it = Use.begin(); it != Use.end(); it++) std::cout << it->first << " ";
            std::cout << std::endl;
            return 1;
         }
         Use[regMethod] = 1;
      }
   }

   // ----------------------------------------------------------------------------------------------------------------------------------
   // Register the training and test trees

   //Monte Carlo data input
   TFile* signalfile = new TFile("/home/karen/Desktop/B_JPsi_03082018/Rootuple_BstoJpsiK_2017_MC_MiniAODSIM_1.root");  //open the file
   TTree* signaltree = (TTree*)signalfile->Get("rootuple/ntuple");
//   std::cout << "--- TMVAClassification   : Using as signal file: " << signalfile->GetName() << std::endl;
   

//   TChain *dataMC = new TChain("ntuple", "");
//   dataMC->Add("/home/karen/Desktop/babymeson/Rootuple_BstoJpsiK_2017_MC_MiniAODSIM_1.root/rootuple/ntuple");
//   dataMC->Add("/home/karen/Desktop/babymeson/Rootuple_BstoJpsiK_2017_MC_MiniAODSIM_1.root/rootuple/ntuple");
//   dataMC->Add("/home/karen/Desktop/babymeson/Rootuple_BstoJpsiK_2017_MC_MiniAODSIM_1.root/rootuple/ntuple");

//   TTree *signaltree = (TTree*)dataMC;


   //Real data input
   TFile* datafile = new TFile("/home/karen/Desktop/B_JPsi_03082018/Rootuple_Run2017F-v1_0000_x0.root");  //open the file
   TTree* datatree = (TTree*)datafile->Get("rootuple/ntuple");
//   std::cout << "--- TMVAClassification   : Using as real data file: " <<datafile->GetName() << std::endl;  
   //------------------------------------------------------------------------------------------------------------------------------------

   // Create a ROOT output file where TMVA will store ntuples, histograms, etc.
   TString outfileName( "TMVAClassification_Bmesonfinal.root" );
   TFile* outputFile = TFile::Open( outfileName, "RECREATE" );
//   std::cout << "--- TMVAClassification   : Creating as resulting file: " << outputFile->GetName() << std::endl;

   //------------------------------------------------------------------------------------------------------------------------------------
   // Create the factory object. Later you can choose the methods
   // whose performance you'd like to investigate. The factory is
   // the only TMVA object you have to interact with
   //
   // The first argument is the base of the name of all the
   // weightfiles in the directory weight/
   //
   // The second argument is the output file for the training results
   // All TMVA output can be suppressed by removing the "!" (not) in
   // front of the "Silent" argument in the option string
   TMVA::Factory *factory = new TMVA::Factory( "TMVAClassification_Bmesonfinal", outputFile,
                                               "!V:!Silent:Color:DrawProgressBar:Transformations=G:AnalysisType=Classification" );

   TMVA::DataLoader *dataloader=new TMVA::DataLoader("Dataset_BJpsi");

   //------------------------------------------------------------------------------------------------------------------------------------
   // Define the input variables that shall be used for the MVA training
   // note that you may also use variable expressions, such as: "3*var1/var2*abs(var3)"
   // [all types of expressions that can also be parsed by TTree::Draw( "expression" )]
//   dataloader->AddVariable( "myvar1 := var1+var2", 'F' );
//   dataloader->AddVariable( "myvar2 := var1-var2", "Expression 2", "", 'F' );
//   dataloader->AddVariable( "var3",                "Variable 3", "units", 'F' );
//   dataloader->AddVariable( "var4",                "Variable 4", "units", 'F' );

   dataloader->AddVariable( "nB", "nB", "units", 'i' );
   dataloader->AddVariable( "nMu", "nMu", "units", 'i' );
//   dataloader->AddVariable( "B_mass", "B_mass", "units", 'F' );
   dataloader->AddVariable( "B_px", "B_px", "units", 'F' );
   dataloader->AddVariable( "B_py", "B_py", "units", 'F' );
   dataloader->AddVariable( "B_pz", "B_pz", "units", 'F' );

//   dataloader->AddVariable( "B_kaon_mass", "B_kaon_mass", "units", 'F' );
   dataloader->AddVariable( "B_kaon_px", "B_kaon_px", "units", 'F' );
   dataloader->AddVariable( "B_kaon_py", "B_kaon_py", "units", 'F' );
   dataloader->AddVariable( "B_kaon_pz", "B_kaon_pz", "units", 'F' );
  
//   dataloader->AddVariable( "B_kaon_parentId1", "B_kaon_parentId1", "units", 'F' );

   dataloader->AddVariable( "B_J_mass", "B_J_mass", "units", 'F' );
   dataloader->AddVariable( "B_J_px", "B_J_px", "units", 'F' );
   dataloader->AddVariable( "B_J_py", "B_J_py", "units", 'F' );
   dataloader->AddVariable( "B_J_pz", "B_J_pz", "units", 'F' );

   dataloader->AddVariable( "B_kaon_px_track", "B_kaon_px_track", "units", 'F' );
   dataloader->AddVariable( "B_kaon_py_track", "B_kaon_py_track", "units", 'F' );
   dataloader->AddVariable( "B_kaon_pz_track", "B_kaon_pz_track", "units", 'F' );
//   dataloader->AddVariable( "B_kaon_charge", "B_kaon_charge", "units", 'F' );
//   dataloader->AddVariable( "B_kaon_pId1", "B_kaon_pId1", "units", 'F' ); //constant

   dataloader->AddVariable( "B_J_px1", "B_J_px1", "units", 'F' );
   dataloader->AddVariable( "B_J_py1", "B_J_py1", "units", 'F' );
   dataloader->AddVariable( "B_J_pz1", "B_J_pz1", "units", 'F' );
//   dataloader->AddVariable( "B_J_charge1", "B_J_charge1", "units", 'F' ); //constant


   dataloader->AddVariable( "B_J_px2", "B_J_px2", "units", 'F' );
   dataloader->AddVariable( "B_J_py2", "B_J_py2", "units", 'F' );
   dataloader->AddVariable( "B_J_pz2", "B_J_pz2", "units", 'F' );
//   dataloader->AddVariable( "B_J_charge2", "B_J_charge2", "units", 'F' ); //constant

   dataloader->AddVariable( "B_chi2", "B_chi2", "units", 'F' );
   dataloader->AddVariable( "B_J_chi2", "B_J_chi2", "units", 'F' );
   dataloader->AddVariable( "B_Prob", "B_Prob", "units", 'F' );
//   dataloader->AddVariable( "B_J_prob", "B_J_Prob", "units", 'F' ); //not valid formula

   dataloader->AddVariable( "B_DecayVtxX", "B_DecayVtxX", "units", 'F' );
   dataloader->AddVariable( "B_DecayVtxY", "B_DecayVtxY", "units", 'F' );
   dataloader->AddVariable( "B_DecayVtxZ", "B_DecayVtxZ", "units", 'F' );
   dataloader->AddVariable( "B_DecayVtxXE", "B_DecayVtxXE", "units", 'F' );
   dataloader->AddVariable( "B_DecayVtxYE", "B_DecayVtxYE", "units", 'F' );
   dataloader->AddVariable( "B_DecayVtxZE", "B_DecayVtxZE", "units", 'F' );
   dataloader->AddVariable( "B_DecayVtxXYE", "B_DecayVtxXYE", "units", 'F' );
   dataloader->AddVariable( "B_DecayVtxXZE", "B_DecayVtxXZE", "units", 'F' );
   dataloader->AddVariable( "B_DecayVtxYZE", "B_DecayVtxYZE", "units", 'F' );

//Works but ANN too complicated
//   dataloader->AddVariable( "B_J_DecayVtxX", "B_J_DecayVtxX", "units", 'F' );
//   dataloader->AddVariable( "B_J_DecayVtxY", "B_J_DecayVtxY", "units", 'F' );
//   dataloader->AddVariable( "B_J_DecayVtxZ", "B_J_DecayYVtxZ", "units", 'F' );
//   dataloader->AddVariable( "B_J_DecayVtxXE", "B_J_DecayVtxXE", "units", 'F' );
//   dataloader->AddVariable( "B_J_DecayVtxYE", "B_J_DecayVtxYE", "units", 'F' );
//   dataloader->AddVariable( "B_J_DecayVtxZE", "B_J_DecayVtxZE", "units", 'F' );
//   dataloader->AddVariable( "B_J_DecayVtxXYE", "B_J_DecayVtxXYE", "units", 'F' );
//   dataloader->AddVariable( "B_J_DecayVtxXZE", "B_J_DecayVtxXZE", "units", 'F' );
//   dataloader->AddVariable( "B_J_DecayVtxYZE", "B_J_DecayVtxYZE", "units", 'F' );

//Works but ANN too complicated
//   dataloader->AddVariable( "priVtxX", "priVtxX", "units", 'F' );
//   dataloader->AddVariable( "priVtxY", "priVtxY", "units", 'F' );
//   dataloader->AddVariable( "priVtxZ", "priVtxZ", "units", 'F' );
//   dataloader->AddVariable( "priVtxXE", "priVtxXE", "units", 'F' );
//   dataloader->AddVariable( "priVtxYE", "priVtxYE", "units", 'F' );
//   dataloader->AddVariable( "priVtxZE", "priVtxZE", "units", 'F' );
//   dataloader->AddVariable( "priVtxXYE", "priVtxXYE", "units", 'F' );
//   dataloader->AddVariable( "priVtxXZE", "priVtxXZE", "units", 'F' );
//   dataloader->AddVariable( "priVtxYZE", "priVtxYZE", "units", 'F' );
//   dataloader->AddVariable( "priVtxCL", "priVtxCL", "units", 'F' );


//   dataloader->AddVariable( "priVtxIPX", "priVtxIPX", "units", 'F' );  //not valid formula
//   dataloader->AddVariable( "priVtxIPY", "priVtxIPY", "units", 'F' );  //not valid formula
//   dataloader->AddVariable( "priVtxIPZ", "priVtxIPZ", "units", 'F' );  //not valid formula
//   dataloader->AddVariable( "priVtxIPXE", "priVtxIPXE", "units", 'F' );  //not valid formula
//   dataloader->AddVariable( "priVtxIPYE", "priVtxIPYE", "units", 'F' );  //not valid formula
//   dataloader->AddVariable( "priVtxIPZE", "priVtxIPZE", "units", 'F' );  //not valid formula
//   dataloader->AddVariable( "priVtxIPXYE", "priVtxIPXYE", "units", 'F' );  //not valid formula
//   dataloader->AddVariable( "priVtxIPXZE", "priVtxIPXZE", "units", 'F' );  //not valid formula
//   dataloader->AddVariable( "priVtxIPYZE", "priVtxIPYZE", "units", 'F' );  //not valid formula
//   dataloader->AddVariable( "priVtxIPCL", "priVtxIPCL", "units", 'F' );  //not valid formula


//   dataloader->AddVariable( "PVXBS", "PVXBS", "units", 'F' );  //problemas de dimensionalidad
//   dataloader->AddVariable( "PVYBS", "PVYBS", "units", 'F' );
//   dataloader->AddVariable( "PVZBS", "PVZBS", "units", 'F' );
//   dataloader->AddVariable( "PVXBSE", "PVXBSE", "units", 'F' );
//   dataloader->AddVariable( "PVYBSE", "PVYBSE", "units", 'F' );
//   dataloader->AddVariable( "PVZBSE", "PVZBSE", "units", 'F' );
//   dataloader->AddVariable( "PVXYBSE", "PVXYBSE", "units", 'F' );
//   dataloader->AddVariable( "PVXZBSE", "PVXZBSE", "units", 'F' );  //constant
//   dataloader->AddVariable( "PVYZBSE", "PVYZBSE", "units", 'F' );  //constant


   dataloader->AddVariable( "vertex_id", "vertex_id", "units", 'F' );
//   dataloader->AddVariable( "nVtx", "nVtx", "units", 'F' );      
//   dataloader->AddVariable( "run", "run", "units", 'F' );  
//   dataloader->AddVariable( "event", "event", "units", 'F' );    
//   dataloader->AddVariable( "lumiblock", "lumiblock", "units", 'F' );   
//   dataloader->AddVariable( "nTgrL", "nTgrL", "units", 'F' );  
//   dataloader->AddVariable( "triggersL", "triggersL", "units", 'F' );    

//-----------Works but ANN too complicated---------------------------
//   dataloader->AddVariable( "kaondxy", "kaondxy", "units", 'F' );
//   dataloader->AddVariable( "kaondz", "kaondz", "units", 'F' );
//   dataloader->AddVariable( "kaondxy_e", "kaondxy", "units", 'F' );
//   dataloader->AddVariable( "kaondz_e", "kaondz_e", "units", 'F' ); 
//-------------------------------------------------------------------------

//-----------------Signal and background get -nan----------------------
//   dataloader->AddVariable( "mumC2", "mumC2", "units", 'F' );
//   dataloader->AddVariable( "mumCat", "mumCat", "units", 'F' );  
//   dataloader->AddVariable( "mumAngT", "mumAngT", "units", 'F' );
   dataloader->AddVariable( "mumNHits", "mumNHits", "units", 'F' );  //Works with only Gaussian transformation

//   dataloader->AddVariable( "mupC2", "mupC2", "units", 'F' );
//   dataloader->AddVariable( "mupCat", "mupCat", "units", 'F' );  
//   dataloader->AddVariable( "mupAngT", "mupAngT", "units", 'F' );
   dataloader->AddVariable( "mupNHits", "mupNHits", "units", 'F' );  //Works with only Gaussian transformation
//--------------------------------------------------------------------


   dataloader->AddVariable( "mumdxy", "mumdxy", "units", 'F' );
   dataloader->AddVariable( "mupdxy", "mupdxy", "units", 'F' );  
   dataloader->AddVariable( "mumdz", "mumdz", "units", 'F' );
   dataloader->AddVariable( "mupdz", "mupdz", "units", 'F' );
   dataloader->AddVariable( "muon_dca", "muon_dca", "units", 'F' );

// <GetSeparation> signal and background histograms have different or invalid dimensions
//   dataloader->AddVariable( "tri_Dim25", "tri_Dim25", "units", 'F' );
//   dataloader->AddVariable( "tri_Dim20", "tri_Dim20", "units", 'F' );
//   dataloader->AddVariable( "tri_JpsiTk", "tri_JpsiTk", "units", 'F' );
//   dataloader->AddVariable( "tri_DoubleMu43Jpsi", "tri_DoubleMu43Jpsi", "units", 'F' );


//   dataloader->AddVariable( "mu1soft", "mu1soft", "units", 'F' );
//   dataloader->AddVariable( "mu2soft", "mu2soft", "units", 'F' );
//   dataloader->AddVariable( "mu1tight", "mu1tight", "units", 'F' );
//   dataloader->AddVariable( "mu2tight", "mu2tight", "units", 'F' );
//   dataloader->AddVariable( "mu1PF", "mu1PF", "units", 'F' );
//   dataloader->AddVariable( "mu2PF", "mu2PF", "units", 'F' );
//   dataloader->AddVariable( "mu1loose", "mu1loose", "units", 'F' );
//   dataloader->AddVariable( "mu2loose", "mu2loose", "units", 'F' );

   // You can add so-called "Spectator variables", which are not used in the MVA training,
   // but will appear in the final "TestTree" produced by TMVA. This TestTree will contain the
   // input variables, the response values of all trained MVAs, and the spectator variables

//   dataloader->AddSpectator( "spec1 := var1*2",  "Spectator 1", "units", 'F' );
//   dataloader->AddSpectator( "spec2 := var1*3",  "Spectator 2", "units", 'F' );
//   dataloader->AddSpectator( "B_kaon_pT:= (B_kaon_px)^2",  "B_kaon_pT", "units", 'F' );


   //------------------------------------------------------------------------------------------------------------------------------------
   // global event weights per tree (see below for setting event-wise weights)
   Double_t signalWeight     = 1.0;
   Double_t backgroundWeight = 1.0;
  
   // You can add an arbitrary number of signal or background trees  //essential part of classification
   dataloader->AddSignalTree    ( signaltree,     signalWeight );
   dataloader->AddBackgroundTree( datatree,   backgroundWeight );
   //------------------------------------------------------------------------------------------------------------------------------------

   // Apply additional cuts on the signal and background samples (can be different)
   TCut mycuts = "B_J_mass==3.097 || B_J_chi2>0.1 || mumNHits>5 || mupNHits>5"; // for example: TCut mycuts = "abs(var1)<0.5 && abs(var2-0.5)<1";
   TCut mycutb = ""; // for example: TCut mycutb = "abs(var1)<0.5";

   // Tell the dataloader how to use the training and testing events
   //
   // If no numbers of events are given, half of the events in the tree are used
   // for training, and the other half for testing:
   //
   //    dataloader->PrepareTrainingAndTestTree( mycut, "SplitMode=random:!V" );
   //
   // To also specify the number of testing events, use:
   //
   //    dataloader->PrepareTrainingAndTestTree( mycut,
   //         "NSigTrain=3000:NBkgTrain=3000:NSigTest=3000:NBkgTest=3000:SplitMode=Random:!V" );
   dataloader->PrepareTrainingAndTestTree( mycuts, mycutb,
                                        "nTrain_Signal=1000:nTrain_Background=1000:SplitMode=Random:NormMode=NumEvents:!V" );

   //------------------------------------------------------------------------------------------------------------------------------------
   // ### Book MVA methods
   //
   // Please lookup the various method configuration options in the corresponding cxx files, eg:
   // src/MethoCuts.cxx, etc, or here: http://tmva.sourceforge.net/optionRef.html
   // it is possible to preset ranges in the option string in which the cut optimisation should be done:
   // "...:CutRangeMin[2]=-1:CutRangeMax[2]=1"...", where [2] is the third input variable

   // TMVA ANN: MLP (recommended ANN) -- all ANNs in TMVA are Multilayer Perceptrons
   if (Use["MLP"])
      factory->BookMethod( dataloader, TMVA::Types::kMLP, "MLP", "H:!V:NeuronType=tanh:VarTransform=N:NCycles=600:HiddenLayers=N+5:TestRate=5:!UseRegulator" );

   if (Use["MLPBFGS"])
      factory->BookMethod( dataloader, TMVA::Types::kMLP, "MLPBFGS", "H:!V:NeuronType=tanh:VarTransform=N:NCycles=600:HiddenLayers=N+5:TestRate=5:TrainingMethod=BFGS:!UseRegulator" );

   if (Use["MLPBNN"])
      factory->BookMethod( dataloader, TMVA::Types::kMLP, "MLPBNN", "H:!V:NeuronType=tanh:VarTransform=N:NCycles=60:HiddenLayers=N+5:TestRate=5:TrainingMethod=BFGS:UseRegulator" ); // BFGS training with bayesian regulators


   // Multi-architecture DNN implementation.
   if (Use["DNN_CPU"] or Use["DNN_GPU"]) {
      // General layout.
      TString layoutString ("Layout=TANH|128,TANH|128,TANH|128,LINEAR");

      // Training strategies.
      TString training0("LearningRate=1e-1,Momentum=0.9,Repetitions=1,"
                        "ConvergenceSteps=20,BatchSize=256,TestRepetitions=10,"
                        "WeightDecay=1e-4,Regularization=L2,"
                        "DropConfig=0.0+0.5+0.5+0.5, Multithreading=True");
      TString training1("LearningRate=1e-2,Momentum=0.9,Repetitions=1,"
                        "ConvergenceSteps=20,BatchSize=256,TestRepetitions=10,"
                        "WeightDecay=1e-4,Regularization=L2,"
                        "DropConfig=0.0+0.0+0.0+0.0, Multithreading=True");
      TString training2("LearningRate=1e-3,Momentum=0.0,Repetitions=1,"
                        "ConvergenceSteps=20,BatchSize=256,TestRepetitions=10,"
                        "WeightDecay=1e-4,Regularization=L2,"
                        "DropConfig=0.0+0.0+0.0+0.0, Multithreading=True");
      TString trainingStrategyString ("TrainingStrategy=");
      trainingStrategyString += training0 + "|" + training1 + "|" + training2;

      // General Options.
      TString dnnOptions ("!H:V:ErrorStrategy=CROSSENTROPY:VarTransform=N:"
                          "WeightInitialization=XAVIERUNIFORM");
      dnnOptions.Append (":"); dnnOptions.Append (layoutString);
      dnnOptions.Append (":"); dnnOptions.Append (trainingStrategyString);

      // Cuda implementation.
      if (Use["DNN_GPU"]) {
         TString gpuOptions = dnnOptions + ":Architecture=GPU";
         factory->BookMethod(dataloader, TMVA::Types::kDNN, "DNN_GPU", gpuOptions);
      }
      // Multi-core CPU implementation.
      if (Use["DNN_CPU"]) {
         TString cpuOptions = dnnOptions + ":Architecture=CPU";
         factory->BookMethod(dataloader, TMVA::Types::kDNN, "DNN_CPU", cpuOptions);
      }
   }

   // CF(Clermont-Ferrand)ANN
   if (Use["CFMlpANN"])
      factory->BookMethod( dataloader, TMVA::Types::kCFMlpANN, "CFMlpANN", "!H:!V:NCycles=200:HiddenLayers=N+1,N"  ); // n_cycles:#nodes:#nodes:...

   // Tmlp(Root)ANN
   if (Use["TMlpANN"])
      factory->BookMethod( dataloader, TMVA::Types::kTMlpANN, "TMlpANN", "!H:!V:NCycles=200:HiddenLayers=N+1,N:LearningMethod=BFGS:ValidationFraction=0.3"  ); // n_cycles:#nodes:#nodes:...



   // Boosted Decision Trees
   if (Use["BDTG"]) // Gradient Boost
      factory->BookMethod( dataloader, TMVA::Types::kBDT, "BDTG",
                           "!H:!V:NTrees=1000:MinNodeSize=2.5%:BoostType=Grad:Shrinkage=0.10:UseBaggedBoost:BaggedSampleFraction=0.5:nCuts=20:MaxDepth=2" );

   if (Use["BDT"])  // Adaptive Boost
      factory->BookMethod( dataloader, TMVA::Types::kBDT, "BDT",
                           "!H:!V:NTrees=850:MinNodeSize=2.5%:MaxDepth=3:BoostType=AdaBoost:AdaBoostBeta=0.5:UseBaggedBoost:BaggedSampleFraction=0.5:SeparationType=GiniIndex:nCuts=20" );

   if (Use["BDTB"]) // Bagging
      factory->BookMethod( dataloader, TMVA::Types::kBDT, "BDTB",
                           "!H:!V:NTrees=400:BoostType=Bagging:SeparationType=GiniIndex:nCuts=20" );

   if (Use["BDTD"]) // Decorrelation + Adaptive Boost
      factory->BookMethod( dataloader, TMVA::Types::kBDT, "BDTD",
                           "!H:!V:NTrees=400:MinNodeSize=5%:MaxDepth=3:BoostType=AdaBoost:SeparationType=GiniIndex:nCuts=20:VarTransform=Decorrelate" );

   if (Use["BDTF"])  // Allow Using Fisher discriminant in node splitting for (strong) linearly correlated variables
      factory->BookMethod( dataloader, TMVA::Types::kBDT, "BDTF",
                           "!H:!V:NTrees=50:MinNodeSize=2.5%:UseFisherCuts:MaxDepth=3:BoostType=AdaBoost:AdaBoostBeta=0.5:SeparationType=GiniIndex:nCuts=20" );
//-----------------------------------------------------------------------------------------------------------------------------------------

   // Now you can tell the factory to train, test, and evaluate the MVAs
   //
   // Train MVAs using the set of training events
   factory->TrainAllMethods();

   // Evaluate all MVAs using the set of test events
   factory->TestAllMethods();

   // Evaluate and compare performance of all configured MVAs
   factory->EvaluateAllMethods();

   // --------------------------------------------------------------

   // Save the output
   outputFile->Close();

   std::cout << "==> TMVAClassification   : Creating as resulting root file: " << outputFile->GetName() << std::endl;
//   std::cout << "==> Wrote root file: " << outputFile->GetName() << std::endl;
   std::cout << "==> TMVAClassification_Bmesonfinal is done!" << std::endl;

   delete factory;
   delete dataloader;
   // Launch the GUI for the root macros
   if (!gROOT->IsBatch()) TMVA::TMVAGui( outfileName );

   return 0;
}

int main( int argc, char** argv )
{
   // Select methods (don't look at this code - not of interest)
   TString methodList;
   for (int i=1; i<argc; i++) {
      TString regMethod(argv[i]);
      if(regMethod=="-b" || regMethod=="--batch") continue;
      if (!methodList.IsNull()) methodList += TString(",");
      methodList += regMethod;
   }
   return TMVAClassification_Bmesonfinal(methodList);
}