#include "Fit/Fitter.h"
#include "Fit/BinData.h"
#include "Fit/Chi2FCN.h"
#include "TList.h"
#include "Math/WrappedMultiTF1.h"
#include "HFitInterface.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TGraphMultiErrors.h"
#include "TFitResult.h"
#include <ROOT/TThreadExecutor.hxx>

using namespace std;

ROOT::Fit::FitResult DoROOTFit(TGraphErrors* g, TF1* ffit);

void test_FitParallel()  {

   ROOT::TThreadExecutor pool; // or ROOT::TProcessExecutor

   // Create the list of graphs:
   vector<TGraphErrors*> m_TGraphErrors;
   ROOT::Fit::FitResult fit_results_tmp[2];
   // graph 1
   const Int_t n = 10;
   Double_t x[n]  = {-0.22, 0.05, 0.25, 0.35, 0.5, 0.61,0.7,0.85,0.89,0.95};
   Double_t y[n]  = {1,2.9,5.6,7.4,9,9.6,8.7,6.3,4.5,1};
   Double_t ex[n] = {0};
   Double_t ey[n] = {.8,.7,.6,.5,.4,.4,.5,.6,.7,.8};

   // graph 2
   Double_t x2[n];
   Double_t y2[n];
   for(Int_t i = 0; i < n; i++ ){
      y2[i] = y[i]*2;
      x2[i] = x[i]+0.2;
   }
   auto gr1 = new TGraphErrors(n,x,y,ex,ey);
   auto gr2 = new TGraphErrors(n,x2,y2,ex,ey);

   m_TGraphErrors.push_back(gr1);
   m_TGraphErrors.push_back(gr2);


   TF1 * f1 = new TF1("f1","pol3", 0, 360);

   cout<<"--------------------> Fit with TThreadExecutor <--------------------------------"<<endl;
   auto fitResults = pool.Map([=](TGraphErrors* a) { return DoROOTFit(a, f1); }, m_TGraphErrors);

   cout<<"--------------------> Graph 1 Fit with ROOT::Fit <--------------------------------"<<endl;
   fit_results_tmp[0] = DoROOTFit(m_TGraphErrors[0], f1);

   cout<<"--------------------> Graph 2 Fit with ROOT::Fit <--------------------------------"<<endl;
   fit_results_tmp[1] = DoROOTFit(m_TGraphErrors[1], f1);

   TMultiGraph* mgr = new TMultiGraph();
   mgr->SetTitle("total cross section + fit; #phi [deg];d^{4}#sigma [nb/GeV^{4}]");
   mgr->Add(m_TGraphErrors[0]);
   mgr->Add(m_TGraphErrors[1]);

   auto c1 = new TCanvas("c1","A Simple Graph with error bars",200,10,700,500);
   mgr->Draw("AP");

}
//___________________________________________________________________________________________________
ROOT::Fit::FitResult DoROOTFit(TGraphErrors* g, TF1* ffit) { //Fit with ROOT::Fit method

	ROOT::Fit::DataOptions opt;
	ROOT::Fit::DataRange range;

	ROOT::Fit::BinData data(opt,range);
	ROOT::Fit::FillData(data, g);

	ROOT::Math::WrappedMultiTF1 fitfunc(*ffit,1);

	ROOT::Fit::Fitter fitter;
	fitter.SetFunction(fitfunc, false);

	ROOT::Math::MinimizerOptions(mopt);

   mopt.SetMinimizerType("Minuit2");
   mopt.SetStrategy(0);
   mopt.SetTolerance(10);

   mopt.SetPrintLevel(1);
	// print the default minimizer option values
   fitter.Config().SetMinimizerOptions(mopt);
   fitter.Fit(data); // chi2 fit'

   ROOT::Fit::FitResult result = fitter.Result();

   result.Print(std::cout);

   ffit ->SetFitResult(result);
   g->GetListOfFunctions()->Add(ffit);

   return result;
}