#include <iostream>
#include <fstream>

#include "TGraph.h"

//Gaussiana 2D
Double_t Gauss2D(Double_t *x, Double_t *par)
{
  Double_t xx =x[0];
  Double_t yy = x[1];
  Double_t f = par[2]*((exp(-0.5*((xx-par[0])/par[1])*((xx-par[0])/par[1])))*exp(-0.5*((yy-par[3])/par[1])*((yy-par[3])/par[1])));
  return f;
}

//Gaussiana 2D -> TF2
TF2 *f_Gauss2D = new TF2("f_Gauss2D",Gauss2D,-15,15,-15,15,4);


//integrale della gaussiana, apertura = variabile, 4 parametetri (stessi della gaussiana)
Double_t I_Gauss2D(Double_t *x, Double_t *par)
{
  f_Gauss2D->SetParameter(0,par[0]); 
  f_Gauss2D->SetParameter(1,par[1]); 
  f_Gauss2D->SetParameter(2,par[2]); 
  f_Gauss2D->SetParameter(3,par[3]); 
  Double_t integral = f_Gauss2D->Integral(-x[0]/2,x[0]/2,-x[0]/2,x[0]/2,1e-4);
  return integral;
}  


//inizio macro
void Fit_slit(){


  //Disegno una copia di Gauss2D a parametri fissati e il suo integrale in funzione dell'apertura della slitta
  TF2 *f1 = new TF2("f1",Gauss2D,-15,15,-15,15,4);
  TF1 *I_f1 = new TF1("I_f1", I_Gauss2D, 0, 15, 4);

  f1->FixParameter(0,0);
  f1->FixParameter(1,2);
  f1->FixParameter(2,3300);
  f1->FixParameter(3,0);

  I_f1->FixParameter(0,0);
  I_f1->FixParameter(1,2);
  I_f1->FixParameter(2,3300);
  I_f1->FixParameter(3,0);


  TCanvas *can = new TCanvas();
  can->Divide(3,2);
  can->cd(1);
  f1->Draw("SURF");
  
  can->cd(2);
  I_f1->Draw();

  Double_t x_8[13] = {0,1,2,3,4,5,6,7,8,9,10,12,14};
  Double_t I_8[13] = {0,11785,39305,67087,85581,94821,98427,99625,99923,99990,100000,100000,100000};  
  Double_t ex_8[13];
  Double_t eI_8[13];
  
  eI_8[0]=10;
  ex_8[0]=0.1;
  
  for(Int_t i=1; i<13; i++)
    {
      ex_8[i]=0.1;
      eI_8[i]=sqrt(I_8[i]);
    }
  
  Double_t x_40[13] = {0,1,2,3,4,5,6,7,8,9,10,12,14};
  Double_t I_40[13] = {0,2629,9908,21048,34079,47669,60426,71513,80359,86948,91725,97001,99075};
  Double_t ex_40[13];
  Double_t eI_40[13];
  
  eI_40[0]=10;
  ex_40[0]=0.1;
  
  for(Int_t i=1; i<13; i++)
    {
      ex_40[i]=0.1;
      eI_40[i]=sqrt(I_40[i]);
    }
  
  TGraphErrors *g_8= new TGraphErrors(13,x_8,I_8,ex_8,eI_8);
  TGraphErrors *g_40= new TGraphErrors(13,x_40,I_40,ex_40,eI_40);

  g_8->GetXaxis()->SetRangeUser(0, 15);
  g_40->GetXaxis()->SetRangeUser(0, 15);
  
  TF1 *I_f1_fit = new TF1("I_f1_fit", I_Gauss2D, 0, 15, 4);
  TF1 *I_f2_fit = new TF1("I_f2_fit", I_Gauss2D, 0, 15, 4);
  TF1 *I_f3_fit = new TF1("I_f3_fit", I_Gauss2D, 0, 15, 4);

  I_f1_fit->SetNpx(10000);
  I_f1_fit->SetParameter(0,0);
  I_f1_fit->SetParLimits(0,-0.1,0.1);
  I_f1_fit->SetParameter(1,1);
  I_f1_fit->SetParLimits(1,0.1,2);
  I_f1_fit->SetParameter(2,15000);
  I_f1_fit->SetParLimits(2,10000,18000);
  I_f1_fit->SetParameter(3,0);
  I_f1_fit->SetParLimits(3,-0.1,0.1);
  //oppure fisso par[0] = par[3] = 0
  I_f1_fit->FixParameter(0,0);
  I_f1_fit->FixParameter(3,0);
  
  I_f2_fit->SetNpx(10000);
  I_f2_fit->SetParameter(0,0);
  I_f2_fit->SetParLimits(0,-0.5,0.5);
  I_f2_fit->SetParameter(1,2.2);
  I_f2_fit->SetParLimits(1,1.5,2.5);
  I_f2_fit->SetParameter(2,3100);
  I_f2_fit->SetParLimits(2,2000,3500);
  I_f2_fit->SetParameter(3,0);
  I_f2_fit->SetParLimits(3,-0.5,0.5);
  //oppure fisso par[0] = par[3] = 0
  I_f2_fit->FixParameter(0,0);
  I_f2_fit->FixParameter(3,0);

  I_f3_fit->SetNpx(10000);
  I_f3_fit->SetParameter(0,0);
  I_f3_fit->SetParLimits(0,-0.5,0.5);
  I_f3_fit->SetParameter(1,1);
  I_f3_fit->SetParLimits(1,0.1,2);
  I_f3_fit->SetParameter(2,6);
  I_f3_fit->SetParLimits(2,2,40);
  I_f3_fit->SetParameter(3,0);
  I_f3_fit->SetParLimits(3,-0.5,0.5);
  //oppure fisso par[0] = par[3] = 0
  I_f3_fit->FixParameter(0,0);
  I_f3_fit->FixParameter(3,0);
  
  can->cd(3);
  g_8->Draw("AP");
  g_8->Fit("I_f1_fit");
  // g_8->Fit("f_real");
  gStyle->SetOptFit(1);
  
  can->cd(4);
  g_40->Draw("AP");
  g_40->Fit("I_f2_fit");

  const int lines = 8;
    Double_t aperture[lines];
    Double_t current[lines];
    Double_t e_aperture[lines];
    Double_t e_current[lines];

     ifstream infile("slit_current_1.75_7.txt");

     for(int i = 0; i < lines; i++) {
        infile >> aperture[i] >> current[i] >> e_aperture[i] >> e_current[i];
        current[i] = current[i]*1E9;
        e_current[i] = e_current[i]*1E9;

     }


  TGraphErrors *first = new TGraphErrors(lines, aperture, current, e_aperture, e_current);

  can->cd(5);
  first->Draw("AP");
  first->Fit("I_f3_fit");

  
}