#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 Gauss2D(Double_t *x, Double_t *par)
{
  Double_t xx =x[0]+6443.5;
  Double_t yy = x[1];
  //Double_t f = 220/(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])));
  return f;
}


TF2 *f_Gauss2D = new TF2("f_Gauss2D",Gauss2D,6420,6460,-15,15,4);

Double_t I_Gauss2D_1(Double_t *x, Double_t *par)
{
  f_Gauss2D->SetParameter(0,x[0]); //la variabile dell'integrale e' il parametro 0 della gaussiana
  f_Gauss2D->SetParameter(1,par[0]); 
  f_Gauss2D->SetParameter(2,par[1]);
  f_Gauss2D->SetParameter(3,par[2]);
  Double_t integral = f_Gauss2D->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->SetParameter(0,x[0]); //la variabile dell'integrale e' il parametro 0 della gaussiana
  f_Gauss2D->SetParameter(1,par[0]); 
  f_Gauss2D->SetParameter(2,par[1]);
  f_Gauss2D->SetParameter(3,par[2]);
  Double_t integral = f_Gauss2D->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->SetParameter(0,x[0]); //la variabile dell'integrale e' il parametro 0 della gaussiana
  f_Gauss2D->SetParameter(1,par[0]); 
  f_Gauss2D->SetParameter(2,par[1]);
  f_Gauss2D->SetParameter(3,par[2]);
  Double_t integral = f_Gauss2D->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 squared 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, //sigma
   1, //y0
   2, //normalization
   3  // slit's opening
};


int par_5[4] = {
   0, // sigma
   1, // y0
   2, // normalization
   4 // slit's opening 

};

int par_10[4] = {
   0, // sigma
   1, // y0
   2, // normalization
   5 // slit's opening

};




 
// 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){}
 

   //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()
{

  ifstream in_file_1;
  ifstream in_file_5;
  ifstream in_file_10;
  in_file_1.open("2023-03-02_scan_50.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_5.close();
  in_file_10.close();

  in_file_1.open("2023-03-02_scan_50.txt");
  in_file_5.open("2023-03-02_scan_52.txt");
  in_file_10.open("2023-03-02_scan_55.txt");

  double B[numLines];
  //double number[numLines];
   double eB[numLines];    // Error in magnetic field
   double i_DF[numLines];  // Current of our Faraday cup
   double ei_DF[numLines]; // Error in current of our Faraday cup
   double i_MiB[numLines]; // Current of Milan's Faraday cup 
   double ei_MiB[numLines];// Error in current of Milan's Faraday cup
   double dummy;


   double B_5[numLines];
  //double number[numLines];
   double eB_5[numLines];    // Error in magnetic field
   double i_DF_5[numLines];  // Current of our Faraday cup
   double ei_DF_5[numLines]; // Error in current of our Faraday cup
   double i_MiB_5[numLines]; // Current of Milan's Faraday cup 
   double ei_MiB_5[numLines];// Error in current of Milan's Faraday cup

   double B_10[numLines];
  //double number[numLines];
   double eB_10[numLines];    // Error in magnetic field
   double i_DF_10[numLines];  // Current of our Faraday cup
   double ei_DF_10[numLines]; // Error in current of our Faraday cup
   double i_MiB_10[numLines]; // Current of Milan's Faraday cup 
   double ei_MiB_10[numLines];// Error in current of Milan's Faraday cup
  

  


 
  //ofstream out_file;
  //out_file.open("out_slit_1.txt");

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

      i_MiB[i] *= 1e9;
      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",51,6400,6480);
  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",51,6400,6480);
  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",51,6400,6480);
  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,6480,4);
   TF1 *I_t_Gauss2D_10 = new TF1("I_t_Gauss2D",I_Gauss2D_10,6400,6480,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(6420, 6470);
   //y.SetRange(0,200);
   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::DataRange range_5;
   //range_5.SetRange(6420, 6470);
   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::DataRange range_10;
   //range_10.SetRange(6420, 6470);
   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] = {7.44, 0.4306, 0.85, 3.011, 15.24, 37.02};
   double par0[Npar] = {4.706, 0, 1.42, 3.62, 18.1, 43.44};

   //1 mm equivale a circa 3 Gauss 
 
   // 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(0).Fix();
   //fitter.Config().ParSettings(1).Fix();
   //fitter.Config().ParSettings(2).Fix();
   /*fitter.Config().ParSettings(3).Fix();
   fitter.Config().ParSettings(4).Fix();
   fitter.Config().ParSettings(5).Fix();*/

   //fitter.Config().ParSettings(1).Fix();

   // set limits on the third and 4-th parameter
   fitter.Config().ParSettings(0).SetLimits(2, 8);   //sigma -> mi aspetto che sia circa 4.5 (fascio di 2 mm ma misuriamo meno)
   fitter.Config().ParSettings(1).SetLimits(-0.5, 0.5); //y0 -> mi aspetto che sia centrato in 0 
   fitter.Config().ParSettings(2).SetLimits(0.3, 2); //normalizzazione -> mi aspetto che sia circa 1.4
   fitter.Config().ParSettings(3).SetLimits(1,5); //apertura slitta 1 mm -> mi aspetto circa 4 
   fitter.Config().ParSettings(4).SetLimits(10,25); //apertura slitta 5 mm -> mi aspetto circa 18
   fitter.Config().ParSettings(5).SetLimits(20,50); //apertura slitta 12 mm -> mi aspetto circa 44

   /*fitter.Config().ParSettings(0).SetLimits(3, 9);   //sigma 
   fitter.Config().ParSettings(1).SetLimits(-0.5, 0.5); //y0
   fitter.Config().ParSettings(2).SetLimits(0.3, 1.5); //normalizzazione 
   fitter.Config().ParSettings(3).SetLimits(1,10); //apertura slitta 1 mm
   fitter.Config().ParSettings(4).SetLimits(10,30); //apertura slitta 5 mm 
   fitter.Config().ParSettings(5).SetLimits(20,55);*/ //apertura slitta 12 mm 
   
 
   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->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->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->Draw();



}