/**********************************************************************************
 * Project   : TMVA - a Root-integrated toolkit for multivariate data analysis    *
 * Package   : TMVA                                                               *
 * Exectuable: TMVAClassificationApplication                                      *
 *                                                                                *
 * This macro provides a simple example on how to use the trained classifiers     *
 * within an analysis module                                                      *
 **********************************************************************************/

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

#include "TFile.h"
#include "TTree.h"
#include "TString.h"
#include "TSystem.h"
#include "TROOT.h"
#include "TStopwatch.h"

#include "TMVAGui.C"

#if not defined(__CINT__) || defined(__MAKECINT__)
#include "TMVA/Tools.h"
#include "TMVA/Reader.h"
#include "TMVA/MethodCuts.h"
#endif

using namespace TMVA;

void myTMVAClassificationApplication( TString myMethodList = "" ) 
{   
#ifdef __CINT__
   gROOT->ProcessLine( ".O0" ); // turn off optimization in CINT
#endif

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

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

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

 // --- Cut optimisation
 
   // 
   // --- Boosted Decision Trees
   Use["BDT"]             = 1; // uses Adaptive Boost
   
   std::cout << std::endl;
   std::cout << "==> Start TMVAClassificationApplication" << 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 = 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;
         }
         Use[regMethod] = 1;
      }
   }

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

   // --- Create the Reader object

   TMVA::Reader *reader = new TMVA::Reader( "!Color:!Silent" );    

   // Create a set of variables and declare them to the reader
   // - the variable names MUST corresponds in name and type to those given in the weight file(s) used
   Float_t axrec, ayrec, aprec;
   Float_t azrec, athrec, aphrec;
   reader->AddVariable( "axrec", &axrec);
   reader->AddVariable( "ayrec", &ayrec);
   reader->AddVariable( "azrec", &azrec);
   reader->AddVariable( "aprec", &aprec);
   reader->AddVariable( "athrec", &athrec);
   reader->AddVariable( "aphrec", &aphrec);


   // // Spectator variables declared in the training have to be added to the reader, too
   // Float_t spec1,spec2;
   //   Float_t threc;
   //   reader->AddSpectator( "spec1 := athrec",  &threc);
   // reader->AddSpectator( "spec1 := var1*2",   &spec1 );
   // reader->AddSpectator( "spec2 := var1*3",   &spec2 );

   //Float_t Category_cat1, Category_cat2, Category_cat3;
   // if (Use["Category"]){
   //    // Add artificial spectators for distinguishing categories
   //   //   reader->AddSpectator( "Category_cat1 := aprec>1.6", &Category_cat1);
   //   //   reader->AddSpectator( "Category_cat2 := aprec<1.6", &Category_cat2);
   //    // reader->AddSpectator( "Category_cat2 := (var3>0)&&(var4<0)",  &Category_cat2 );
   //    // reader->AddSpectator( "Category_cat3 := (var3>0)&&(var4>=0)", &Category_cat3 );
   // }

   // --- Book the MVA methods

   TString dir    = "weights/";
   TString prefix = "TMVAClassification";

   // Book method(s)
   for (std::map<std::string,int>::iterator it = Use.begin(); it != Use.end(); it++) {
      if (it->second) {
         TString methodName = TString(it->first) + TString(" method");
         TString weightfile = dir + prefix + TString("_") + TString(it->first) + TString(".weights.xml");
         reader->BookMVA( methodName, weightfile ); 
      }
   }
   
   // Book output histograms
   UInt_t nbin = 1000;
   TH1F    *histBdt(0);
   
   if (Use["BDT"])           histBdt     = new TH1F( "MVA_BDT",           "MVA_BDT",           nbin, -1.1, 1.1);
   

   // Prepare input tree (this must be replaced by your data source)
   
   // TString fname = "/panda/myResults/DPMoutputDec2013/FTF/XYCut2/mom_1_5/GoodBadSig_3_1sigma_3_1sigma.root";
 TString fname = "/panda/myResults/DPMoutputMarch2014_BeamProf_newXYCut/mom_1_5/GoodBadSig_3sigma_3sigma_DPMel.root";
   TFile *input = TFile::Open(fname);

   std::cout << "--- TMVAClassificationApp    : Using input file: " << input->GetName() << std::endl;
   
   // --- Event loop

   // Prepare the event tree
   // - here the variable names have to corresponds to your tree
   // - you can use the same variables as above which is slightly faster,
   //   but of course you can use different ones and copy the values inside the event loop
   //
   std::cout << "--- Select signal sample" << std::endl;
   TTree* theTree = (TTree*)input->Get("sigGood");//2591101 events
   //TTree* theTree = (TTree*)input->Get("sigBad");
   //  double f_axrec, f_ayrec, f_azrec, f_aprec, f_athrec, f_aphrec;
  float f_axrec, f_ayrec, f_azrec, f_aprec, f_athrec, f_aphrec;
   theTree->SetBranchAddress( "axrec", &f_axrec );
   theTree->SetBranchAddress( "ayrec", &f_ayrec );
   theTree->SetBranchAddress( "azrec", &f_azrec );
   theTree->SetBranchAddress( "aprec", &f_aprec );
   theTree->SetBranchAddress( "athrec", &f_athrec );
   theTree->SetBranchAddress( "aphrec", &f_aphrec );
  
   std::cout << "--- Processing: " << theTree->GetEntries() << " events" << std::endl;
   TStopwatch sw;
   sw.Start();
   for (Long64_t ievt=0; ievt<theTree->GetEntries();ievt++) {
   //   for (Long64_t ievt=0; ievt<10000;ievt++) {
     axrec = f_axrec;    ayrec = f_ayrec;    azrec = f_azrec;
     aprec = f_aprec;    athrec = f_athrec;    aphrec = f_aphrec;
     // cout<<" "<<f_axrec<<" "<<f_ayrec<<" "<<f_aprec<<endl;
     // cout<<" "<<axrec<<" "<<ayrec<<" "<<aprec<<endl;
      if (ievt%1000 == 0) std::cout << "--- ... Processing event: " << ievt << std::endl;

      theTree->GetEntry(ievt);

      // var1 = userVar1 + userVar2;
      // var2 = userVar1 - userVar2;

      // --- Return the MVA outputs and fill into histograms

      if (Use["BDT"          ])   histBdt    ->Fill( reader->EvaluateMVA( "BDT method"     ) );
    
   }

   // Get elapsed time
   sw.Stop();
   std::cout << "--- End of event loop: "; sw.Print();

   

   // --- Write histograms

   TFile *target  = new TFile( "TMVApp.root","RECREATE" );
   
   if (Use["BDT"          ])   histBdt    ->Write();
  
   std::cout << "--- Created root file: \"TMVApp.root\" containing the MVA output histograms" << std::endl;
  
   delete reader;
    
   std::cout << "==> TMVAClassificationApplication is done!" << endl << std::endl;
}