#include <Fit/Fitter.h>
#include <Fit/BinData.h>
#include <Fit/Chi2FCN.h>
#include <TH1.h>
#include <Math/WrappedMultiTF1.h>
#include <HFitInterface.h>
#include <TCanvas.h>
#include <TStyle.h>
#include <iostream>
#include <fstream>

#include "TGraph.h"

using namespace std; 

  
 
// definition of shared parameter
// background function

//we define the functions here (Gauss_2d and Gauss_2d_integral)

Double_t center = 6441.55;

//best parameter for alpha -0.0785
//best center 6441.5

Double_t Gauss2D_1(Double_t *x, Double_t *par)
{
  Double_t xx =x[0] + center;
  Double_t yy = x[1];
  Double_t alpha_left = -0.092;
  Double_t alpha_right = 0.15;
  //Double_t f = 100000/(2*3.1415*par[1]*par[1])*((exp(-0.5*((xx-par[0])/par[1])*((xx-par[0])/par[1])))*exp(-0.5*((yy-par[2])/par[1])*((yy-par[2])/par[1])));
  Double_t f = (par[3]*((exp(-0.5*((xx-par[0])/par[1])*((xx-par[0])/par[1])))*exp(-0.5*((yy-par[2])/par[1])*((yy-par[2])/par[1])))) * (1+alpha_left*(xx-par[0])) * (1 + alpha_right*(par[0] - xx));
  return f;
}

Double_t Gauss2D_5(Double_t *x, Double_t *par)
{
  Double_t xx =x[0] + center;
  Double_t yy = x[1];
  Double_t alpha_left_5 = -0.08;
  Double_t alpha_right_5 = 0.03;
  //Double_t f = 100000/(2*3.1415*par[1]*par[1])*((exp(-0.5*((xx-par[0])/par[1])*((xx-par[0])/par[1])))*exp(-0.5*((yy-par[2])/par[1])*((yy-par[2])/par[1])));
  Double_t f = (par[3]*((exp(-0.5*((xx-par[0])/par[1])*((xx-par[0])/par[1])))*exp(-0.5*((yy-par[2])/par[1])*((yy-par[2])/par[1])))) * (1+alpha_left_5*(xx-par[0])) * (1 + alpha_right_5*(par[0] - xx));
  return f;
}

Double_t Gauss2D_10(Double_t *x, Double_t *par)
{
  Double_t xx =x[0] + center;
  Double_t yy = x[1];
  Double_t alpha_left_10 = -0.086;
  Double_t alpha_right_10 = 0.2;
  //Double_t f = 100000/(2*3.1415*par[1]*par[1])*((exp(-0.5*((xx-par[0])/par[1])*((xx-par[0])/par[1])))*exp(-0.5*((yy-par[2])/par[1])*((yy-par[2])/par[1])));
  Double_t f = (par[3]*((exp(-0.5*((xx-par[0])/par[1])*((xx-par[0])/par[1])))*exp(-0.5*((yy-par[2])/par[1])*((yy-par[2])/par[1]))))* (1+alpha_left_10*(xx-par[0])) * (1 + alpha_right_10*(par[0] - xx));
  return f;
}



TF2 *f_Gauss2D_1 = new TF2("f_Gauss2D",Gauss2D_1,6400,6470,-15,15,4);
TF2 *f_Gauss2D_5 = new TF2("f_Gauss2D",Gauss2D_5,6400,6490,-15,15,4);
TF2 *f_Gauss2D_10 = new TF2("f_Gauss2D",Gauss2D_10,6400,6480,-15,15,4);

Double_t I_Gauss2D_1(Double_t *x, Double_t *par)
{
  f_Gauss2D_1->SetParameter(0,x[0]); //la variabile dell'integrale e' il parametro 0 della gaussiana
  f_Gauss2D_1->SetParameter(1,par[0]); 
  f_Gauss2D_1->SetParameter(2,par[1]);
  f_Gauss2D_1->SetParameter(3,par[2]);
  Double_t integral = f_Gauss2D_1->Integral(-par[3]/2,par[3]/2,-par[3]/2,par[3]/2,1e-4);
  return integral;
}

Double_t I_Gauss2D_5(Double_t *x, Double_t *par)
{
  f_Gauss2D_5->SetParameter(0,x[0]); //la variabile dell'integrale e' il parametro 0 della gaussiana
  f_Gauss2D_5->SetParameter(1,par[0]); 
  f_Gauss2D_5->SetParameter(2,par[1]);
  f_Gauss2D_5->SetParameter(3,par[2]);
  Double_t integral = f_Gauss2D_5->Integral(-par[3]/2,par[3]/2,-par[3]/2,par[3]/2,1e-4);
  return integral;
}

Double_t I_Gauss2D_10(Double_t *x, Double_t *par)
{
  f_Gauss2D_10->SetParameter(0,x[0]); //la variabile dell'integrale e' il parametro 0 della gaussiana
  f_Gauss2D_10->SetParameter(1,par[0]); 
  f_Gauss2D_10->SetParameter(2,par[1]);
  f_Gauss2D_10->SetParameter(3,par[2]);
  Double_t integral = f_Gauss2D_10->Integral(-par[3]/2,par[3]/2,-par[3]/2,par[3]/2,1e-4);
  return integral;
}



//in our case we have different gaussian functions with four shared parameters: sigma, normalization, x0 and y0
//we also define the convolution as the gaussian integrated on a square domain -> the aperture (the limits of the integral) is always different 

//we define the vectors of the parameters -> we can define three vectors: one for the gaussian and two for the integrated gaussians  
int par_1[4] = {
   0, //x0
   1, //sigma
   2, //y0
   3, //normalization
};

// signal + background function
int par_5[4] = {
   0, 
   1, 
   2, 
   4
};

int par_10[4] = {
   0, 
   1, 
   2, 
   5

};




 
// Create the GlobalCHi2 structure


 
struct GlobalChi2 {
   GlobalChi2(ROOT::Math::IMultiGenFunction &f1, ROOT::Math::IMultiGenFunction &f2, ROOT::Math::IMultiGenFunction &f3 ) : fChi2_1(&f1), fChi2_2(&f2), fChi2_3(&f3) {}
 
   // parameter vector is first background (in common 1 and 2)
   // and then is signal (only in 2)

   //I don't know how to do this in our case 

   //It takes a pointer to an array of double values (par) as input. Tha values in the par array are used to set the parameters 
   //of the two chi-square functions. The values of the shared parameters are set using iparB array for the  background function
   //the values of the signal parameters are set using the iparSB array for the signal + bkg function 

   //It returns the sum of the chi-square values calculated using the two functions with the input parameters set to par 
   //The chi-square is defined as the sum of (observed values - expected values)^2 divided for the expected values 
   double operator()(const double *par) const
   {
      double p1[4];
      for (int i = 0; i < 4; ++i)
         p1[i] = par[par_1[i]];
 
      double p2[4];
      for (int i = 0; i < 4; ++i)
         p2[i] = par[par_5[i]];

         double p3[4];
      for (int i = 0; i < 4; ++i)
         p3[i] = par[par_10[i]];
 
 
      return (*fChi2_1)(p1) + (*fChi2_2)(p2) + (*fChi2_3)(p3) ;
   }
 
   const ROOT::Math::IMultiGenFunction *fChi2_1;
   const ROOT::Math::IMultiGenFunction *fChi2_2;
   const ROOT::Math::IMultiGenFunction *fChi2_3;
};
 
void combinedFit_real()
{

  ifstream in_file_1;
  ifstream in_file_5;
  ifstream in_file_10;
  in_file_1.open("2023-03-02_scan_12.txt");
  
  int numLines = 0;
  std::string unused;
  while ( std::getline(in_file_1, unused) )
    ++numLines;
  cout << "numero linee " << numLines << endl;
  in_file_1.close();

  in_file_1.open("2023-03-02_scan_12.txt");
  in_file_5.open("2023-02-28_scan_37.txt");
  in_file_10.open("2023-02-28_scan_48.txt");

  double B[numLines];
  double eB[numLines];
  double i_MiB[numLines];
  double ei_MiB[numLines];
  double B_5[numLines];
  double eB_5[numLines];
  double i_MiB_5[numLines];
  double ei_MiB_5[numLines];
  double B_10[numLines];
  double eB_10[numLines];
  double i_MiB_10[numLines];
  double ei_MiB_10[numLines];
  double dummy; 
 
  ofstream out_file;
  out_file.open("out_slit.txt");

  
  for(int i=0; i<numLines; i++)
    {
      in_file_1 >> B[i] >> eB[i] >> dummy >> dummy >> dummy >> dummy >> i_MiB[i] >> ei_MiB[i];
      in_file_5 >> B_5[i] >> eB_5[i] >> dummy >> dummy >> dummy >> dummy >> i_MiB_5[i] >> ei_MiB_5[i];
      in_file_10 >> B_10[i] >> eB_10[i] >> dummy >> dummy >> dummy >> dummy >> i_MiB_10[i] >> ei_MiB_10[i];

      i_MiB[i] *= 1e9;
      eB[i] += 0.8;
      eB_5[i] += 0.8;
      eB_10[i] += 0.8;
      ei_MiB[i] *= 1e9;
      i_MiB_5[i] *= 1e9;
      ei_MiB_5[i] *= 1e9;
      i_MiB_10[i] *= 1e9;
      ei_MiB_10[i] *= 1e9;

    }
  

  
  TGraphErrors *g1 = new TGraphErrors(numLines,B,i_MiB,eB,ei_MiB);
  TGraphErrors *g5 = new TGraphErrors(numLines,B_5,i_MiB_5,eB_5,ei_MiB_5);
  TGraphErrors *g10 = new TGraphErrors(numLines,B_10,i_MiB_10,eB_10,ei_MiB_10);


  
  TH1F *h_1=new TH1F("h_1","test_1",100,-15,15);
  h_1->Clear();
  h_1->Reset();
  auto nPoints_1 = g1->GetN(); // number of points in your TGraph
  for(int i=0; i < nPoints_1; ++i) {
  double x = g1->GetX()[i];
  double y = g1->GetY()[i];
  h_1->SetBinContent(h_1->FindBin(x),y);
  }

  TH1F *h_5=new TH1F("h_5","test_5",100,-15,15);
  h_5->Clear();
  h_5->Reset();
  auto nPoints_5 = g5->GetN(); // number of points in your TGraph
  for(int i=0; i < nPoints_5; ++i) {
  double x = g5->GetX()[i];
  double y = g5->GetY()[i];
  h_5->SetBinContent(h_5->FindBin(x),y);
  }

  TH1F *h_10=new TH1F("h_10","test_10",100,-15,15);
  h_10->Clear();
  h_10->Reset();
  auto nPoints_10 = g10->GetN(); // number of points in your TGraph
  for(int i=0; i < nPoints_10; ++i) {
  double x = g10->GetX()[i];
  double y = g10->GetY()[i];
  h_10->SetBinContent(h_10->FindBin(x),y);
  }

  

   //we don't need histograms, we have to fit tha data (TGraphErrors)
 
   //TH1D *hB = new TH1D("hB", "histo B", 100, 0, 100);
   //TH1D *hSB = new TH1D("hSB", "histo S+B", 100, 0, 100);
 
   //we create TF1 objects for gaussian and gaussian integrated 

   TF1 *I_t_Gauss2D_1 = new TF1("I_t_Gauss2D",I_Gauss2D_1,6400,6480,4);
   TF1 *I_t_Gauss2D_5 = new TF1("I_t_Gauss2D",I_Gauss2D_5,6400,6490,4);
   TF1 *I_t_Gauss2D_10 = new TF1("I_t_Gauss2D",I_Gauss2D_10,6300,6500,4);
   /*TF1 *fB = new TF1("fB", "expo", 0, 100);
   fB->SetParameters(1, -0.05);
   hB->FillRandom("fB");
 
   TF1 *fS = new TF1("fS", "gaus", 0, 100);
   fS->SetParameters(1, 30, 5);
 
   hSB->FillRandom("fB", 2000);
   hSB->FillRandom("fS", 1000);*/
 
   // perform now global fit
 
   //TF1 *fSB = new TF1("fSB", "expo + gaus(2)", 0, 100);
 
   ROOT::Math::WrappedMultiTF1 w_1(*I_t_Gauss2D_1, 1);
   ROOT::Math::WrappedMultiTF1 w_5(*I_t_Gauss2D_5, 1);
   ROOT::Math::WrappedMultiTF1 w_10(*I_t_Gauss2D_10, 1);
 
   ROOT::Fit::DataOptions opt;
   /*ROOT::Fit::DataRange range_1;
   // set the data range
   range_1.SetRange(-5, 5);
   ROOT::Fit::BinData data_1(opt, range_1);
   ROOT::Fit::FillData(data_1, h_1);
 
   ROOT::Fit::DataRange range_5;
   range_5.SetRange(-5, 5);
   ROOT::Fit::BinData data_5(opt, range_5);
   ROOT::Fit::FillData(data_5, h_5);
 
   ROOT::Fit::DataRange range_10;
   range_10.SetRange(-5, 5);
   ROOT::Fit::BinData data_10(opt, range_10);
   ROOT::Fit::FillData(data_10, h_10);*/

   int n = numLines; 
   ROOT::Fit::BinData data_1(n, B,i_MiB,eB,ei_MiB);
   ROOT::Fit::FillData(data_1, g1, I_t_Gauss2D_1);

   ROOT::Fit::BinData data_5(n, B_5,i_MiB_5,eB_5,ei_MiB_5);
   ROOT::Fit::FillData(data_5, g5, I_t_Gauss2D_5);

   ROOT::Fit::BinData data_10(n, B_10,i_MiB_10,eB_10,ei_MiB_10);
   ROOT::Fit::FillData(data_10, g10, I_t_Gauss2D_10);


 
   ROOT::Fit::Chi2Function chi2_1(data_1, w_1);
   ROOT::Fit::Chi2Function chi2_5(data_5, w_5);
   ROOT::Fit::Chi2Function chi2_10(data_10, w_10);
 
   GlobalChi2 globalChi2(chi2_1, chi2_5, chi2_10);
 
   ROOT::Fit::Fitter fitter;
   
   //total parameters                                  
   const int Npar = 6;
   double par0[Npar] = {3.4, 0, 2, 3, 15.2, 36};
 
   // create before the parameter settings in order to fix or set range on them
   fitter.Config().SetParamsSettings(6, par0);
   // fix 5-th parameter
   /*fitter.Config().ParSettings(3).Fix();
   fitter.Config().ParSettings(2).Fix();
   fitter.Config().ParSettings(4).Fix();  */  
   // set limits on the third and 4-th parameter
   fitter.Config().ParSettings(2).Fix();
   fitter.Config().ParSettings(0).SetLimits(1, 3.5);
   //fitter.Config().ParSettings(2).SetLimits(0.05, 2);
   fitter.Config().ParSettings(3).SetLimits(1,20);
   fitter.Config().ParSettings(4).SetLimits(8,25);
   fitter.Config().ParSettings(5).SetLimits(15,50);

   fitter.Config().MinimizerOptions().SetPrintLevel(0);
   fitter.Config().SetMinimizer("Minuit2", "Migrad");
 
   // fit FCN function directly
   // (specify optionally data size and flag to indicate that is a chi2 fit)
   fitter.FitFCN(6, globalChi2, 0, data_1.Size() + data_5.Size() + data_10.Size(), true);
   ROOT::Fit::FitResult result = fitter.Result();
   result.Print(std::cout);
 
   TCanvas *c1 = new TCanvas();
   c1->Divide(1, 3);
   c1->cd(1);
   gStyle->SetOptFit(1111);
 
   I_t_Gauss2D_1->SetFitResult(result, par_1);
   //fB->SetRange(rangeB().first, rangeB().second);
   I_t_Gauss2D_1->SetLineColor(kBlue);
   g1->GetListOfFunctions()->Add(I_t_Gauss2D_1);
   g1->GetXaxis()->SetTitle("B [Gauss]");
   g1->GetYaxis()->SetTitle("Current [nA]");
   g1->SetTitle("Current vs magnetic field, slit opening 1 mm");
   g1->Draw();
 
   c1->cd(2);
   I_t_Gauss2D_5->SetFitResult(result, par_5);
   //fSB->SetRange(rangeSB().first, rangeSB().second);
   I_t_Gauss2D_5->SetLineColor(kRed);
   g5->GetListOfFunctions()->Add(I_t_Gauss2D_5);
   //g5->GetXaxis()->SetRangeUser(6400,6480);
   g5->GetXaxis()->SetTitle("B [Gauss]");
   g5->GetYaxis()->SetTitle("Current [nA]");
   g5->SetTitle("Current vs magnetic field, slit opening 5 mm");
   g5->Draw();

   c1->cd(3);
   I_t_Gauss2D_10->SetFitResult(result, par_10);
   //fSB->SetRange(rangeSB().first, rangeSB().second);
   I_t_Gauss2D_10->SetLineColor(kGreen);
   g10->GetListOfFunctions()->Add(I_t_Gauss2D_10);
   g10->GetXaxis()->SetTitle("B [Gauss]");
   g10->GetYaxis()->SetTitle("Current [nA]");
   g10->SetTitle("Current vs magnetic field, slit opening 10 mm");
   g10->Draw();



}