// @(#)root/tmva $Id$
/**********************************************************************************
 * Project   : TMVA - a Root-integrated toolkit for multivariate data analysis    *
 * Package   : TMVA                                                               *
 * Root Macro: TMVARegression                                                     *
 *                                                                                *
 * This macro provides examples for the training and testing of the               *
 * TMVA classifiers.                                                              *
 *                                                                                *
 * As input data is used a toy-MC sample consisting of four Gaussian-distributed  *
 * and linearly correlated input variables.                                       *
 *                                                                                *
 * The methods to be used can be switched on and off by means of booleans, or     *
 * via the prompt command, for example:                                           *
 *                                                                                *
 *    root -l TMVARegression.C\(\"LD,MLP\"\)                                      *
 *                                                                                *
 * (note that the backslashes are mandatory)                                      *
 * If no method given, a default set is used.                                     *
 *                                                                                *
 * The output file "TMVAReg.root" can be analysed with the use of dedicated       *
 * macros (simply say: root -l <macro.C>), which can be conveniently              *
 * invoked through a GUI that will appear at the end of the run of this macro.    *
 **********************************************************************************/

#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/Tools.h"
#include "TMVA/Factory.h"
#include "TMVA/TMVARegGui.h"


using namespace TMVA;
   
void TMVARegression( TString myMethodList = "" ) 
{
   // The explicit loading of the shared libTMVA is done in TMVAlogon.C, defined in .rootrc
   // if you use your private .rootrc, or run from a different directory, please copy the 
   // corresponding lines from .rootrc

   // methods to be processed can be given as an argument; use format:
   //
   // mylinux~> root -l TMVARegression.C\(\"myMethod1,myMethod2,myMethod3\"\)
   //

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



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

   // --- Mutidimensional likelihood and Nearest-Neighbour methods
   Use["PDERS"]           = 0;
   Use["PDEFoam"]         = 0;//1; 
   Use["KNN"]             = 0;//1;
   // 
   // --- Linear Discriminant Analysis
   Use["LD"]		        = 0;//1;
   // 
   // --- Function Discriminant analysis
   Use["FDA_GA"]          = 0;//1;
   Use["FDA_MC"]          = 0;
   Use["FDA_MT"]          = 0;
   Use["FDA_GAMT"]        = 0;
   // 
   // --- Neural Network
   Use["MLP"]             = 0;//1; 
   // 
   // --- Support Vector Machine 
   Use["SVM"]             = 0;
   // 
   // --- Boosted Decision Trees
   Use["BDT"]             = 0;//1;//0;
   Use["BDTG"]            = 1;//1;
   // ---------------------------------------------------------------

   std::cout << std::endl;
   std::cout << "==> Start TMVARegression" << 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;
      }
   }

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

   // --- Here the preparation phase begins

   // Create a new root output file
   TString outfileName( "TMVAReg.root" );
   TFile* outputFile = TFile::Open( outfileName, "RECREATE" );

   // Create the factory object. Later you can choose the methods
   // whose performance you'd like to investigate. The factory will
   // then run the performance analysis for you.
   //
   // 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( "TMVARegression", outputFile, 
                                               "!V:!Silent:Color:DrawProgressBar" );

   // If you wish to modify default settings 
   // (please check "src/Config.h" to see all available global options)
   //    (TMVA::gConfig().GetVariablePlotting()).fTimesRMS = 8.0;
   //    (TMVA::gConfig().GetIONames()).fWeightFileDir = "myWeightDirectory";

   // 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" )]
   //   factory->AddVariable( "var1", "Variable 1", "units", 'F' );
   ///   factory->AddVariable( "var2", "Variable 2", "units", 'F' );
   cout<<"prima add var "<<endl;
   //      factory->AddVariable( "MET_px", "Variable 1", "units", 'F' );
   //   factory->AddVariable( "(BL_LepPhi[0])", "lep 1", "units", 'F' );
   //   factory->AddVariable( "(BL_LepPhi[1])", "lep 1", "units", 'F' );
   //   factory->AddVariable( "(BL_LepPhi[2])", "lep 1", "units", 'F' );
   //   factory->AddVariable( "MET_SoftTrk_px", "PV Trk Soft Term", "units", 'F' );
   //   factory->AddVariable( "lep0x : = BL_LepPt[0]*cos(BL_LepPhi[0])", "lep 1", "units", 'F' );

   
   factory->AddVariable( "lep0x : = BL_LepPt[0]*cos(BL_LepPhi[0])", "lep0x", "units", 'F' );
   factory->AddVariable( "lep1x := BL_LepPt[1]*cos(BL_LepPhi[1])", "lep1x", "units", 'F' );
   factory->AddVariable( "lep2x := BL_LepPt[2]*cos(BL_LepPhi[2])", "lep2x", "units", 'F' );
   
   factory->AddVariable( "metx_rest := MET_px - MET_SoftTrk_px + BL_LepPt[0]*cos(BL_LepPhi[0])+BL_LepPt[1]*cos(BL_LepPhi[1])+BL_LepPt[2]*cos(BL_LepPhi[2])", "Rest METx", "units", 'F' );
   factory->AddVariable( "MET_SoftTrk_px", "PV Trk Soft Term", "units", 'F' );
   
   //   factory->AddVariable( "MET_px", "MET Term", "units", 'F' );
   factory->AddVariable("ptall := BL_LepPt[0]+BL_LepPt[1]+BL_LepPt[2]");
   //   factory->AddVariable( "met_softtrk_phi:=atan2(MET_SoftTrk_py,MET_SoftTrk_px)", "PV Trk Soft Term", "units", 'F' );

   //   factory->AddVariable("met := sqrt(MET_px+MET_px + MET_py*MET_py)");
   //   factory->AddVariable("year");
   //   factory->AddVariable("nBaseLep");
   //   factory->AddVariable( "MET_px", "MET Term", "units", 'F' );
   //   factory->AddVariable( "MET_TrkMET_px", "MET Term", "units", 'F' );


   factory->AddVariable( "lep0y := BL_LepPt[0]*sin(BL_LepPhi[0])", "lep0y", "units", 'F' );
   factory->AddVariable( "lep1y := BL_LepPt[1]*sin(BL_LepPhi[1])", "lep1y", "units", 'F' );
   factory->AddVariable( "lep2y := BL_LepPt[2]*sin(BL_LepPhi[2])", "lep2y", "units", 'F' );
   //   factory->AddVariable( "MET_py - MET_SoftTrk_py + BL_LepPt[0]*sin(BL_LepPhi[0])", "PV Trk Soft Term", "units", 'F' );
   factory->AddVariable( "mety_rest := MET_py - MET_SoftTrk_py + BL_LepPt[0]*sin(BL_LepPhi[0])+BL_LepPt[1]*sin(BL_LepPhi[1])+BL_LepPt[2]*sin(BL_LepPhi[2])", "Rest METy", "units", 'F' );
   factory->AddVariable( "MET_SoftTrk_py", "PV Trk Soft Term", "units", 'F' );

   //   factory->AddVariable( "MET_py", "MET Term", "units", 'F' );
   //   factory->AddVariable( "met := sqrt(MET_py*MET_py+MET_px*MET_px)", "MET Term" "units", 'F' );
   //   factory->AddVariable( "met_phi := atan2(MET_py,MET_px)", "MET Term", "units", 'F' );
   //   factory->AddVariable( "BL_LepPhi[0]", "lep pt10", "units", 'F' );
   //   factory->AddVariable( "BL_LepPhi[1]", "lep pt10", "units", 'F' );
   //   factory->AddVariable( "BL_LepPhi[2]", "lep pt10", "units", 'F' );
   //   factory->AddVariable( "cos_met_phi := cos(atan2(MET_py,MET_px))", "MET Term", "units", 'F' );
   /*
   factory->AddVariable( "BL_LepPt[0]", "lep pt10", "units", 'F' );
   factory->AddVariable( "BL_LepEta[0]", "lep eta0", "units", 'F' );
   factory->AddVariable( "BL_LepPhi[0]", "lep phi", "units", 'F' );

   factory->AddVariable( "BL_LepPt[1]", "lep pt1", "units", 'F' );
   factory->AddVariable( "BL_LepEta[1]", "lep eta1", "units", 'F' );
   factory->AddVariable( "BL_LepPhi[1]", "lep phi1", "units", 'F' );
   factory->AddVariable( "BL_LepPt[2]", "lep pt2", "units", 'F' );
   factory->AddVariable( "BL_LepEta[2]", "lep eta2", "units", 'F' );
   factory->AddVariable( "BL_LepPhi[2]", "lep phi2", "units", 'F' );
   */
   //   factory->AddVariable( "met_softtrk_phi := atan2(MET_SoftTrk_py,MET_SoftTrk_px)", "MET Term", "units", 'F' );
   //   factory->AddVariable( "met_rest_phi := atan2(MET_py - MET_SoftTrk_py + BL_LepPt[0]*sin(BL_LepPhi[0])+BL_LepPt[1]*sin(BL_LepPhi[1])+BL_LepPt[2]*sin(BL_LepPhi[2]),MET_px - MET_SoftTrk_px + BL_LepPt[0]*cos(BL_LepPhi[0])+BL_LepPt[1]*cos(BL_LepPhi[1])+BL_LepPt[2]*cos(BL_LepPhi[2]))", "PV Trk Soft Term", "units", 'F' );

   //   factory->AddVariable( "met_softtrk := sqrt(MET_SoftTrk_px*MET_SoftTrk_px + MET_SoftTrk_py*MET_SoftTrk_py)", "MET Term", "units", 'F' );
   //   factory->AddVariable( "MET_TrkMET_py", "MET Term", "units", 'F' );

   //   factory->AddVariable( "BL_LepPt[0]*sinh(BL_LepEta[0])", "lepz 1", "units", 'F' );
   //   factory->AddVariable( "BL_LepPt[1]*sinh(BL_LepEta[1])", "lepz 2", "units", 'F' );
   //   factory->AddVariable( "BL_LepPt[2]*sinh(BL_LepEta[2])", "lepz 3", "units", 'F' );

   //   factory->AddVariable( "MET_py", "MET Term", "units", 'F' );
   //factory->AddVariable( "averageIntPerXing","mu","units",'F');
   //factory->AddVariable( "nVtx","npv","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
   //   factory->AddSpectator( "MET_phi:=atan2(MET_py,MET_px)",  "Spectator 1", "units", 'F' );
   //   factory->AddSpectator( "cos_metlep0phi_reco := (BL_LepPt[0]*cos(BL_LepPhi[0])*MET_px + BL_LepPt[0]*sin(BL_LepPhi[0])*MET_py)/(sqrt(MET_px*MET_px + MET_py*MET_py)*sqrt(BL_LepPt[0]*cos(BL_LepPhi[0])*BL_LepPt[0]*cos(BL_LepPhi[0]) + BL_LepPt[0]*sin(BL_LepPhi[0])*BL_LepPt[0]*sin(BL_LepPhi[0])))","cos","units",'F');
   factory->AddSpectator( "MET_px",  "Spectator 1", "units", 'F' );
   factory->AddSpectator( "MET_py",  "Spectator 1", "units", 'F' );
   
   factory->AddSpectator( "MET_px_truth","sp x","units",'F'); 
   factory->AddSpectator( "MET_py_truth","sp y","units",'F'); 
   factory->AddSpectator( "MET_pz_truth","sp z","units",'F'); 

   /*
   factory->AddSpectator( "lep0_z :=  BL_LepPt[0]*sinh(BL_LepEta[0]) ",  "Spectator 1", "units", 'F' );
   factory->AddSpectator( "lep1_z :=  BL_LepPt[1]*sinh(BL_LepEta[1]) ",  "Spectator 1", "units", 'F' );
   factory->AddSpectator( "lep2_z :=  BL_LepPt[2]*sinh(BL_LepEta[2]) ",  "Spectator 1", "units", 'F' );

   factory->AddSpectator( "lep0_p :=  BL_LepPt[0]*cosh(BL_LepEta[0]) ",  "Spectator 1", "units", 'F' );
   factory->AddSpectator( "lep1_p :=  BL_LepPt[1]*cosh(BL_LepEta[1]) ",  "Spectator 1", "units", 'F' );
   factory->AddSpectator( "lep2_p :=  BL_LepPt[2]*cosh(BL_LepEta[2]) ",  "Spectator 1", "units", 'F' );
   */


   //   factory->AddSpectator( "invmass_truthmet := fabs(sqrt(  ( 0.*BL_LepPt[0]*cosh(BL_LepEta[0]) +  BL_LepPt[1]*cosh(BL_LepEta[1]) +  BL_LepPt[2]*cosh(BL_LepEta[2]) + sqrt(MET_px_truth*MET_px_truth + MET_py_truth*MET_py_truth + MET_pz_truth*MET_pz_truth) )* ( 0.*BL_LepPt[0]*cosh(BL_LepEta[0]) +  BL_LepPt[1]*cosh(BL_LepEta[1]) +  BL_LepPt[2]*cosh(BL_LepEta[2]) + sqrt(MET_px_truth*MET_px_truth + MET_py_truth*MET_py_truth + MET_pz_truth*MET_pz_truth)) -  ( 0.*BL_LepPt[0]*cos(BL_LepPhi[0]) +  BL_LepPt[1]*cos(BL_LepPhi[1]) +   BL_LepPt[2]*cos(BL_LepPhi[2]) + MET_px_truth )*( 0.*BL_LepPt[0]*cos(BL_LepPhi[0]) +  BL_LepPt[1]*cos(BL_LepPhi[1]) +   BL_LepPt[2]*cos(BL_LepPhi[2]) + MET_px_truth ) -( 0.*BL_LepPt[0]*sin(BL_LepPhi[0]) +  BL_LepPt[1]*sin(BL_LepPhi[1]) +   BL_LepPt[2]*sin(BL_LepPhi[2]) + MET_py_truth )*( 0.*BL_LepPt[0]*sin(BL_LepPhi[0]) +  BL_LepPt[1]*sin(BL_LepPhi[1]) +   BL_LepPt[2]*sin(BL_LepPhi[2]) + MET_py_truth ) - (0.*BL_LepPt[0]*sinh(BL_LepEta[0]) + BL_LepPt[1]*sinh(BL_LepEta[1]) + BL_LepPt[2]*sinh(BL_LepEta[2]) + MET_pz_truth)*(0.*BL_LepPt[0]*sinh(BL_LepEta[0]) + BL_LepPt[1]*sinh(BL_LepEta[1]) + BL_LepPt[2]*sinh(BL_LepEta[2]) + MET_pz_truth)  )) ",  "Spectator 1", "units", 'F' );

   factory->AddSpectator( "year",  "year", "units", 'I' );

   factory->AddSpectator("EventNumber","eventnumber","units",'F');
   //   factory->AddSpectator( "spec2:=var1*3",  "Spectator 2", "units", 'F' );

   // Add the variable carrying the regression target
   //   factory->AddTarget( "fvalue" ); 
   //   factory->AddTarget( "year"); 
   //   factory->AddTarget( "cos_metphi_truth:=cos(atan2(MET_py_truth,MET_px_truth))"); 
   //factory->AddTarget( "cos_metlep0phi_truth := (BL_LepPt[0]*cos(BL_LepPhi[0])*MET_px_truth + BL_LepPt[0]*sin(BL_LepPhi[0])*MET_py_truth)/(sqrt(MET_px_truth*MET_px_truth + MET_py_truth*MET_py_truth)*sqrt(BL_LepPt[0]*cos(BL_LepPhi[0])*BL_LepPt[0]*cos(BL_LepPhi[0]) + BL_LepPt[0]*sin(BL_LepPhi[0])*BL_LepPt[0]*sin(BL_LepPhi[0])))");
   //   factory->AddTarget( "metphi_truth:=(atan2(MET_py_truth,MET_px_truth))"); 
   //   factory->AddTarget("up := (lep1x*MET_px_truth+ lep1y*MET_py_truth)/(sqrt(MET_px_truth*MET_px_truth+MET_py_truth*MET_pt_truth)*lep1_p")
   //   factory->AddTarget( "MET_px_truth"); 
   //   factory->AddTarget( "MET_py_truth"); 
   //      factory->AddTarget( "MET_pz_truth" ); 

   
   //   factory->AddTarget( "pxnu_HNL"); 
   //   factory->AddTarget( "pynu_HNL"); 
   //   factory->AddTarget( "pznu_HNL"); 

   //         factory->AddTarget( "fabs(MET_py_truth + MET_px_truth)<10" ); 
   //factory->AddTarget( "sqrt(  ( BL_LepPt[0]*cosh(BL_LepEta[0]) +  BL_LepPt[1]*cosh(BL_LepEta[1]) +  BL_LepPt[2]*cosh(BL_LepEta[2]))* ( BL_LepPt[0]*cosh(BL_LepEta[0]) +  BL_LepPt[1]*cosh(BL_LepEta[1]) +  BL_LepPt[2]*cosh(BL_LepEta[2]))  )  <2" ); 

   
   //   factory->AddTarget( "cond := fabs(sqrt(  ( BL_LepPt[0]*cosh(BL_LepEta[0]) +  BL_LepPt[1]*cosh(BL_LepEta[1]) +  BL_LepPt[2]*cosh(BL_LepEta[2]) + sqrt(MET_px_truth*MET_px_truth + MET_py_truth*MET_py_truth + MET_pz_truth*MET_pz_truth) )* ( BL_LepPt[0]*cosh(BL_LepEta[0]) +  BL_LepPt[1]*cosh(BL_LepEta[1]) +  BL_LepPt[2]*cosh(BL_LepEta[2]) + sqrt(MET_px_truth*MET_px_truth + MET_py_truth*MET_py_truth + MET_pz_truth*MET_pz_truth)) -  ( BL_LepPt[0]*cos(BL_LepPhi[0]) +  BL_LepPt[1]*cos(BL_LepPhi[1]) +   BL_LepPt[2]*cos(BL_LepPhi[2]) + MET_px_truth )*( BL_LepPt[0]*cos(BL_LepPhi[0]) +  BL_LepPt[1]*cos(BL_LepPhi[1]) +   BL_LepPt[2]*cos(BL_LepPhi[2]) + MET_px_truth ) -( BL_LepPt[0]*sin(BL_LepPhi[0]) +  BL_LepPt[1]*sin(BL_LepPhi[1]) +   BL_LepPt[2]*sin(BL_LepPhi[2]) + MET_py_truth )*( BL_LepPt[0]*sin(BL_LepPhi[0]) +  BL_LepPt[1]*sin(BL_LepPhi[1]) +   BL_LepPt[2]*sin(BL_LepPhi[2]) + MET_py_truth ) - (BL_LepPt[0]*sinh(BL_LepEta[0]) + BL_LepPt[1]*sinh(BL_LepEta[1]) + BL_LepPt[2]*sinh(BL_LepEta[2]) + MET_pz_truth)*(BL_LepPt[0]*sinh(BL_LepEta[0]) + BL_LepPt[1]*sinh(BL_LepEta[1]) + BL_LepPt[2]*sinh(BL_LepEta[2]) + MET_pz_truth)  ))");

   //   factory->AddTarget( "cond := fabs(sqrt(  ( BL_LepPt[0]*cosh(BL_LepEta[0]) +  BL_LepPt[1]*cosh(BL_LepEta[1]) +  BL_LepPt[2]*cosh(BL_LepEta[2]) + sqrt(pxnu_HNL*pxnu_HNL + pynu_HNL*pynu_HNL + pznu_HNL*pznu_HNL) )* ( BL_LepPt[0]*cosh(BL_LepEta[0]) +  BL_LepPt[1]*cosh(BL_LepEta[1]) +  BL_LepPt[2]*cosh(BL_LepEta[2]) + sqrt(pxnu_HNL*pxnu_HNL + pynu_HNL*pynu_HNL + pznu_HNL*pznu_HNL)) -  ( BL_LepPt[0]*cos(BL_LepPhi[0]) +  BL_LepPt[1]*cos(BL_LepPhi[1]) +   BL_LepPt[2]*cos(BL_LepPhi[2]) + pxnu_HNL )*( BL_LepPt[0]*cos(BL_LepPhi[0]) +  BL_LepPt[1]*cos(BL_LepPhi[1]) +   BL_LepPt[2]*cos(BL_LepPhi[2]) + pxnu_HNL ) -( BL_LepPt[0]*sin(BL_LepPhi[0]) +  BL_LepPt[1]*sin(BL_LepPhi[1]) +   BL_LepPt[2]*sin(BL_LepPhi[2]) + pynu_HNL )*( BL_LepPt[0]*sin(BL_LepPhi[0]) +  BL_LepPt[1]*sin(BL_LepPhi[1]) +   BL_LepPt[2]*sin(BL_LepPhi[2]) + pynu_HNL ) - (BL_LepPt[0]*sinh(BL_LepEta[0]) + BL_LepPt[1]*sinh(BL_LepEta[1]) + BL_LepPt[2]*sinh(BL_LepEta[2]) + pznu_HNL)*(BL_LepPt[0]*sinh(BL_LepEta[0]) + BL_LepPt[1]*sinh(BL_LepEta[1]) + BL_LepPt[2]*sinh(BL_LepEta[2]) + pznu_HNL)  ))");//,"","",30.,130.);

      //      factory->AddTarget( "invmass := fabs(sqrt(  ( 0.*BL_LepPt[0]*cosh(BL_LepEta[0]) +  BL_LepPt[1]*cosh(BL_LepEta[1]) +  BL_LepPt[2]*cosh(BL_LepEta[2]) + sqrt(pxnu_HNL*pxnu_HNL + pynu_HNL*pynu_HNL + pznu_HNL*pznu_HNL) )* ( 0.*BL_LepPt[0]*cosh(BL_LepEta[0]) +  BL_LepPt[1]*cosh(BL_LepEta[1]) +  BL_LepPt[2]*cosh(BL_LepEta[2]) + sqrt(pxnu_HNL*pxnu_HNL + pynu_HNL*pynu_HNL + pznu_HNL*pznu_HNL)) -  ( 0.*BL_LepPt[0]*cos(BL_LepPhi[0]) +  BL_LepPt[1]*cos(BL_LepPhi[1]) +   BL_LepPt[2]*cos(BL_LepPhi[2]) + pxnu_HNL )*( 0.*BL_LepPt[0]*cos(BL_LepPhi[0]) +  BL_LepPt[1]*cos(BL_LepPhi[1]) +   BL_LepPt[2]*cos(BL_LepPhi[2]) + pxnu_HNL ) -( 0.*BL_LepPt[0]*sin(BL_LepPhi[0]) +  BL_LepPt[1]*sin(BL_LepPhi[1]) +   BL_LepPt[2]*sin(BL_LepPhi[2]) + pynu_HNL )*( 0.*BL_LepPt[0]*sin(BL_LepPhi[0]) +  BL_LepPt[1]*sin(BL_LepPhi[1]) +   BL_LepPt[2]*sin(BL_LepPhi[2]) + pynu_HNL ) - (0.*BL_LepPt[0]*sinh(BL_LepEta[0]) + BL_LepPt[1]*sinh(BL_LepEta[1]) + BL_LepPt[2]*sinh(BL_LepEta[2]) + pznu_HNL)*(0.*BL_LepPt[0]*sinh(BL_LepEta[0]) + BL_LepPt[1]*sinh(BL_LepEta[1]) + BL_LepPt[2]*sinh(BL_LepEta[2]) + pznu_HNL)  ))");//,"",'B' ,0.5,1.5);  // prima volta era 50

   factory->AddTarget( "invmass := fabs(sqrt(  ( 0.*BL_LepPt[0]*cosh(BL_LepEta[0]) +  BL_LepPt[1]*cosh(BL_LepEta[1]) +  BL_LepPt[2]*cosh(BL_LepEta[2]) + sqrt(MET_px_truth*MET_px_truth + MET_py_truth*MET_py_truth + MET_pz_truth*MET_pz_truth) )* ( 0.*BL_LepPt[0]*cosh(BL_LepEta[0]) +  BL_LepPt[1]*cosh(BL_LepEta[1]) +  BL_LepPt[2]*cosh(BL_LepEta[2]) + sqrt(MET_px_truth*MET_px_truth + MET_py_truth*MET_py_truth + MET_pz_truth*MET_pz_truth)) -  ( 0.*BL_LepPt[0]*cos(BL_LepPhi[0]) +  BL_LepPt[1]*cos(BL_LepPhi[1]) +   BL_LepPt[2]*cos(BL_LepPhi[2]) + MET_px_truth )*( 0.*BL_LepPt[0]*cos(BL_LepPhi[0]) +  BL_LepPt[1]*cos(BL_LepPhi[1]) +   BL_LepPt[2]*cos(BL_LepPhi[2]) + MET_px_truth ) -( 0.*BL_LepPt[0]*sin(BL_LepPhi[0]) +  BL_LepPt[1]*sin(BL_LepPhi[1]) +   BL_LepPt[2]*sin(BL_LepPhi[2]) + MET_py_truth )*( 0.*BL_LepPt[0]*sin(BL_LepPhi[0]) +  BL_LepPt[1]*sin(BL_LepPhi[1]) +   BL_LepPt[2]*sin(BL_LepPhi[2]) + MET_py_truth ) - (0.*BL_LepPt[0]*sinh(BL_LepEta[0]) + BL_LepPt[1]*sinh(BL_LepEta[1]) + BL_LepPt[2]*sinh(BL_LepEta[2]) + MET_pz_truth)*(0.*BL_LepPt[0]*sinh(BL_LepEta[0]) + BL_LepPt[1]*sinh(BL_LepEta[1]) + BL_LepPt[2]*sinh(BL_LepEta[2]) + MET_pz_truth)  ))");//,"",'B' ,0.5,1.5);  // prima volta era 50

   //   factory->AddTarget( "invmass := fabs(sqrt(  ( 0.*BL_LepPt[0]*cosh(BL_LepEta[0]) +  BL_LepPt[1]*cosh(BL_LepEta[1]) +  BL_LepPt[2]*cosh(BL_LepEta[2]) + sqrt(MET_px*MET_px + MET_py*MET_py + MET_pz_truth*MET_pz_truth) )* ( 0.*BL_LepPt[0]*cosh(BL_LepEta[0]) +  BL_LepPt[1]*cosh(BL_LepEta[1]) +  BL_LepPt[2]*cosh(BL_LepEta[2]) + sqrt(MET_px*MET_px + MET_py*MET_py + MET_pz_truth*MET_pz_truth)) -  ( 0.*BL_LepPt[0]*cos(BL_LepPhi[0]) +  BL_LepPt[1]*cos(BL_LepPhi[1]) +   BL_LepPt[2]*cos(BL_LepPhi[2]) + MET_px )*( 0.*BL_LepPt[0]*cos(BL_LepPhi[0]) +  BL_LepPt[1]*cos(BL_LepPhi[1]) +   BL_LepPt[2]*cos(BL_LepPhi[2]) + MET_px ) -( 0.*BL_LepPt[0]*sin(BL_LepPhi[0]) +  BL_LepPt[1]*sin(BL_LepPhi[1]) +   BL_LepPt[2]*sin(BL_LepPhi[2]) + MET_py )*( 0.*BL_LepPt[0]*sin(BL_LepPhi[0]) +  BL_LepPt[1]*sin(BL_LepPhi[1]) +   BL_LepPt[2]*sin(BL_LepPhi[2]) + MET_py ) - (0.*BL_LepPt[0]*sinh(BL_LepEta[0]) + BL_LepPt[1]*sinh(BL_LepEta[1]) + BL_LepPt[2]*sinh(BL_LepEta[2]) + MET_pz_truth)*(0.*BL_LepPt[0]*sinh(BL_LepEta[0]) + BL_LepPt[1]*sinh(BL_LepEta[1]) + BL_LepPt[2]*sinh(BL_LepEta[2]) + MET_pz_truth)  ))");//,"",'B' ,0.5,1.5);  // prima volta era 50


   // Elep = 
   // It is also possible to declare additional targets for multi-dimensional regression, ie:
   // -- factory->AddTarget( "fvalue2" );
   // BUT: this is currently ONLY implemented for MLP

   // Read training and test data (see TMVAClassification for reading ASCII files)
   // load the signal and background event samples from ROOT trees
   TFile *input(0);
   TFile *input2(0);
   TFile *input3(0);
   TFile *input4(0);
   //   TString fname = "../ntuples/WmuHNL50_5G_ntp.root";//"./tmva_reg_example.root";
   //   TString fname = "../ntuples/Signal_Wmu/WmuHNL50_10G_ntp_v1.root";
   //   TString fname = "../ntuples/Signal_Wmu/WmuHNL50_10G_ntp_metnu.root";
   //   TString fname2 = "../ntuples/Signal_Wmu/WmuHNL50_30G_ntp_v1.root";
   TString fname4 = "../ntuples/Signal_Wmu/WmuHNL50_50G_ntp.root";
   TString fname3 = "../ntuples/Signal_Wmu/WmuHNL50_30G_ntp_v1.root";
   //   TString fname = "../ntuples/Signal_Wmu/WmuHNL50_50G_ntp.root";
   //   TString fname = "../ntuples/Signal_Wmu/WmuHNL50_10G_ntp.root";
   TString fname = "../ntuples/Signal_Wmu_withxsec/WmuHNL50_10G_ntp.root";
   //   TString fname = "../ntuples/Signal_Wmu_withXsec/WmuHNL50_10G_ntp.root";
   //   TString fname = "small.root";
   //      TString fname = "../ntuples/Wjets_withxsec/Wjets_merged.root";
   TString fname2 = "../ntuples/Wjets_withxsec/Wjets_merged.root";
   if (!gSystem->AccessPathName( fname )) 
      input = TFile::Open( fname ); // check if file in local directory exists
   else 
      input = TFile::Open( "http://root.cern.ch/files/tmva_reg_example.root" ); // if not: download from ROOT server
   
   input2 = TFile::Open( fname2 ); // check if file in local directory exists
   input3 = TFile::Open( fname3 ); // check if file in local directory exists
   input4 = TFile::Open( fname4 ); // check if file in local directory exists
   if (!input) {
      std::cout << "ERROR: could not open data file" << std::endl;
      exit(1);
   }
   std::cout << "--- TMVARegression           : Using input file: " << input->GetName() << std::endl;

   // --- Register the regression tree

   //   TTree *regTree = (TTree*)input->Get("TreeR");
   TTree *regTree = (TTree*)input->Get("HFntupleNONE");
   TTree *regTree2 = (TTree*)input2->Get("HFntupleNONE");
   TTree *regTree3 = (TTree*)input3->Get("HFntupleNONE");
   TTree *regTree4 = (TTree*)input4->Get("HFntupleNONE");

   // global event weights per tree (see below for setting event-wise weights)
   Double_t regWeight  = 1.0;   
   Double_t regWeight2  = 1.0;   
   Double_t regWeight3  = 1.0;   
   Double_t regWeight4  = 1.0;   
  
   // You can add an arbitrary number of regression trees
   factory->AddRegressionTree( regTree,  regWeight );
   //   factory->AddRegressionTree( regTree2, regWeight2 );
   //   factory->AddRegressionTree( regTree3, regWeight3 );
   //   factory->AddRegressionTree( regTree4, regWeight4 );

   // This would set individual event weights (the variables defined in the 
   // expression need to exist in the original TTree)
   //   factory->SetWeightExpression( "var1", "Regression" );
   //   factory->SetWeightExpression( "XSecWeight","Regression");
   // Apply additional cuts on the signal and background samples (can be different)
   //   TCut mycut = "BL_pass3L==1"; // for example: TCut mycut = "abs(var1)<0.5 && abs(var2-0.5)<1";
   //   TCut mycut = "BL_pass3L==1&&nJets==0"; // for example: TCut mycut = "abs(var1)<0.5 && abs(var2-0.5)<1";
   //TCut mycut = "BL_pass3L ==1 && fabs(BL_LepCharge[1]+BL_LepCharge[2])==0";
   //     TCut mycut = "BL_pass3L ==1 && fabs(BL_LepCharge[1]+BL_LepCharge[2])==0 & pxnu_HNL!=0.";
   //     TCut mycut = "BL_pass3L ==1 && fabs(BL_LepCharge[1]+BL_LepCharge[2])==0 & pxnu_HNL!=0. && fabs(BL_LepPhi[1]-BL_LepPhi[2])<1.0";
   //      TCut mycut = "BL_pass3L ==1 && fabs(BL_LepCharge[1]+BL_LepCharge[2])==0  && fabs(BL_LepPhi[1]-BL_LepPhi[2])<1.0&&nJets==0";
   //   TCut mycut = "BL_pass3L ==1 &&nJets==0";
   TCut mycut = "BL_pass3L ==1 && fabs(BL_LepCharge[1]+BL_LepCharge[2])==0 & nJets==0";
   //   TCut mycut = "BL_pass3L ==1 && fabs(BL_LepCharge[1]+BL_LepCharge[2])==0 && sqrt(  ( 0.*BL_LepPt[0]*cosh(BL_LepEta[0]) +  BL_LepPt[1]*cosh(BL_LepEta[1]) +  BL_LepPt[2]*cosh(BL_LepEta[2]) + sqrt(pxnu_HNL*pxnu_HNL + pynu_HNL*pynu_HNL + pznu_HNL*pznu_HNL) )* ( 0.*BL_LepPt[0]*cosh(BL_LepEta[0]) +  BL_LepPt[1]*cosh(BL_LepEta[1]) +  BL_LepPt[2]*cosh(BL_LepEta[2]) + sqrt(pxnu_HNL*pxnu_HNL + pynu_HNL*pynu_HNL + pznu_HNL*pznu_HNL)) -  ( 0.*BL_LepPt[0]*cos(BL_LepPhi[0]) +  BL_LepPt[1]*cos(BL_LepPhi[1]) +   BL_LepPt[2]*cos(BL_LepPhi[2]) + pxnu_HNL )*( 0.*BL_LepPt[0]*cos(BL_LepPhi[0]) +  BL_LepPt[1]*cos(BL_LepPhi[1]) +   BL_LepPt[2]*cos(BL_LepPhi[2]) + pxnu_HNL ) -( 0.*BL_LepPt[0]*sin(BL_LepPhi[0]) +  BL_LepPt[1]*sin(BL_LepPhi[1]) +   BL_LepPt[2]*sin(BL_LepPhi[2]) + pynu_HNL )*( 0.*BL_LepPt[0]*sin(BL_LepPhi[0]) +  BL_LepPt[1]*sin(BL_LepPhi[1]) +   BL_LepPt[2]*sin(BL_LepPhi[2]) + pynu_HNL ) - (0.*BL_LepPt[0]*sinh(BL_LepEta[0]) + BL_LepPt[1]*sinh(BL_LepEta[1]) + BL_LepPt[2]*sinh(BL_LepEta[2]) + pznu_HNL)*(0.*BL_LepPt[0]*sinh(BL_LepEta[0]) + BL_LepPt[1]*sinh(BL_LepEta[1]) + BL_LepPt[2]*sinh(BL_LepEta[2]) + pznu_HNL)  )<40";

   // tell the factory to use all remaining events in the trees after training for testing:
   factory->PrepareTrainingAndTestTree( mycut, 
					//					"nTrain_Regression=1000:nTest_Regression=0:SplitMode=Random:NormMode=NumEvents:!V" ); //def
					"nTrain_Regression=1000:nTest_Regression=0:SplitMode=Random:NormMode=NumEvents:!V" ); //def
					//					"nTrain_Regression=3000:nTest_Regression=0:SplitMode=Random:NormMode=NumEvents:!V" );
   // factory->PrepareTrainingAndTestTree( mycut, 
   //                                      "nTrain_Regression=0:nTest_Regression=0:SplitMode=Random:NormMode=NumEvents:!V" );

   // If no numbers of events are given, half of the events in the tree are used 
   // for training, and the other half for testing:
   //   factory->PrepareTrainingAndTestTree( mycut, "SplitMode=random:!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

   // PDE - RS method
   if (Use["PDERS"])
      factory->BookMethod( TMVA::Types::kPDERS, "PDERS", 
                           "!H:!V:NormTree=T:VolumeRangeMode=Adaptive:KernelEstimator=Gauss:GaussSigma=0.3:NEventsMin=40:NEventsMax=60:VarTransform=None" );
   // And the options strings for the MinMax and RMS methods, respectively:
   //      "!H:!V:VolumeRangeMode=MinMax:DeltaFrac=0.2:KernelEstimator=Gauss:GaussSigma=0.3" );   
   //      "!H:!V:VolumeRangeMode=RMS:DeltaFrac=3:KernelEstimator=Gauss:GaussSigma=0.3" );   

   if (Use["PDEFoam"])
       factory->BookMethod( TMVA::Types::kPDEFoam, "PDEFoam", 
			    "!H:!V:MultiTargetRegression=F:TargetSelection=Mpv:TailCut=0.001:VolFrac=0.0666:nActiveCells=500:nSampl=2000:nBin=5:Compress=T:Kernel=None:Nmin=10:VarTransform=None" );

   // K-Nearest Neighbour classifier (KNN)
   if (Use["KNN"])
      factory->BookMethod( TMVA::Types::kKNN, "KNN", 
                           "nkNN=20:ScaleFrac=0.8:SigmaFact=1.0:Kernel=Gaus:UseKernel=F:UseWeight=T:!Trim" );

   // Linear discriminant
   if (Use["LD"])
      factory->BookMethod( TMVA::Types::kLD, "LD", 
                           "!H:!V:VarTransform=None" );

	// Function discrimination analysis (FDA) -- test of various fitters - the recommended one is Minuit (or GA or SA)
   if (Use["FDA_MC"]) 
      factory->BookMethod( TMVA::Types::kFDA, "FDA_MC",
                          "!H:!V:Formula=(0)+(1)*x0+(2)*x1:ParRanges=(-100,100);(-100,100);(-100,100):FitMethod=MC:SampleSize=100000:Sigma=0.1:VarTransform=D" );
   
   if (Use["FDA_GA"]) // can also use Simulated Annealing (SA) algorithm (see Cuts_SA options) .. the formula of this example is good for parabolas
      factory->BookMethod( TMVA::Types::kFDA, "FDA_GA",
                           "!H:!V:Formula=(0)+(1)*x0+(2)*x1:ParRanges=(-100,100);(-100,100);(-100,100):FitMethod=GA:PopSize=100:Cycles=3:Steps=30:Trim=True:SaveBestGen=1:VarTransform=Norm" );

   if (Use["FDA_MT"]) 
      factory->BookMethod( TMVA::Types::kFDA, "FDA_MT",
                           "!H:!V:Formula=(0)+(1)*x0+(2)*x1:ParRanges=(-100,100);(-100,100);(-100,100);(-10,10):FitMethod=MINUIT:ErrorLevel=1:PrintLevel=-1:FitStrategy=2:UseImprove:UseMinos:SetBatch" );

   if (Use["FDA_GAMT"]) 
      factory->BookMethod( TMVA::Types::kFDA, "FDA_GAMT",
                           "!H:!V:Formula=(0)+(1)*x0+(2)*x1:ParRanges=(-100,100);(-100,100);(-100,100):FitMethod=GA:Converger=MINUIT:ErrorLevel=1:PrintLevel=-1:FitStrategy=0:!UseImprove:!UseMinos:SetBatch:Cycles=1:PopSize=5:Steps=5:Trim" );

   // Neural network (MLP)
   if (Use["MLP"])
     factory->BookMethod( TMVA::Types::kMLP, "MLP", "!H:!V:VarTransform=Norm:NeuronType=tanh:NCycles=20000:HiddenLayers=N+20:TestRate=6:TrainingMethod=BFGS:Sampling=0.3:SamplingEpoch=0.8:ConvergenceImprove=1e-6:ConvergenceTests=15:!UseRegulator" ); //def
     //          factory->BookMethod( TMVA::Types::kMLP, "MLP", "!H:!V:VarTransform=Norm:NeuronType=tanh:NCycles=20000:HiddenLayers=N+20:TestRate=6:TrainingMethod=BFGS:Sampling=0.3:SamplingEpoch=0.8:ConvergenceImprove=1e-10:ConvergenceTests=200:!UseRegulator" ); //def
     //     factory->BookMethod( TMVA::Types::kMLP, "MLP", "!H:!V:VarTransform=Norm:NeuronType=tanh:NCycles=20000:HiddenLayers=N+10,N:TestRate=6:TrainingMethod=BFGS:Sampling=0.3:SamplingEpoch=0.8:ConvergenceImprove=1e-3:ConvergenceTests=30:!UseRegulator" ); //def
     //      factory->BookMethod( TMVA::Types::kMLP, "MLP", "!H:!V:VarTransform=Norm:NeuronType=tanh:NCycles=20000:HiddenLayers=N+10:TestRate=6:TrainingMethod=BFGS:Sampling=0.3:SamplingEpoch=0.8:ConvergenceImprove=1e-8:ConvergenceTests=15:!UseRegulator" ); //def
     //     factory->BookMethod( TMVA::Types::kMLP, "MLP", "H:!V:NeuronType=tanh:VarTransform=N:NCycles=600:HiddenLayers=N+5:TestRate=5:!UseRegulator" ); //from tmvaclassification

   // Support Vector Machine
   if (Use["SVM"])
      factory->BookMethod( TMVA::Types::kSVM, "SVM", "Gamma=0.25:Tol=0.001:VarTransform=Norm" );

   // Boosted Decision Trees
   if (Use["BDT"])
     
     factory->BookMethod( TMVA::Types::kBDT, "BDT",
			  //			  "!H:!V:NTrees=100:MinNodeSize=1.0%:BoostType=AdaBoostR2:SeparationType=RegressionVariance:nCuts=20:PruneMethod=CostComplexity:PruneStrength=30" ); //def
    "!H:!V:NTrees=400:MinNodeSize=1.0%:BoostType=AdaBoostR2:SeparationType=RegressionVariance:nCuts=20:PruneMethod=CostComplexity:PruneStrength=30" ); //def
   
   /*     factory->BookMethod( 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["BDTG"])
     factory->BookMethod( TMVA::Types::kBDT, "BDTG",
			  //			  "!H:!V:NTrees=2000::BoostType=Grad:Shrinkage=0.1:UseBaggedBoost:BaggedSampleFraction=0.5:nCuts=20:MaxDepth=3:MaxDepth=4" );
			  "!H:!V:NTrees=1000::BoostType=Grad:Shrinkage=0.1:UseBaggedBoost:BaggedSampleFraction=0.5:nCuts=20:MaxDepth=3:MaxDepth=4" ); //from keras example
   // --------------------------------------------------------------------------------------------------

   // ---- 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 << "==> Wrote root file: " << outputFile->GetName() << std::endl;
   std::cout << "==> TMVARegression is done!" << std::endl;      

   delete factory;

   // Launch the GUI for the root macros
   if (!gROOT->IsBatch()) TMVA::TMVARegGui( outfileName );
}

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;
   }
   TMVARegression(methodList);
   return 0;
}