#include <cstdlib>
#include <cstdio>
#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <sstream>
#include <cmath>
#include "TColor.h"
#include "TStyle.h"
#include <TFile.h>
#include <TCanvas.h>
#include <TH2F.h>
#include <TH2D.h>
#include <TH3D.h>
#include <TSystem.h>
#include "TPaletteAxis.h"
#include <TRandom3.h>
#include <TF1.h>
#include <TF2.h>
#include <TMinuit.h>
#include "TMatrixD.h"
#include "TMatrixDSym.h"
#include "TFitter.h"
#include "TVirtualFitter.h"
#include "TRegexp.h"
#include "TCanvas.h"
#include "TMath.h"
#include "TROOT.h"
#include "TLegend.h"
using namespace std;




Double_t Gaus2D(Double_t *x, Double_t *par)
{
if(par[2] > 0 && par[4] > 0)
   {
  double rx=0.;
  double ry=0.;
  double rz=0.;
  double den=0;
  rx=(x[0]-par[1])/par[2];
  ry=(x[1]-par[3])/par[4];
  rz=2.*par[5]*rx*ry;
  den= 1.-par[5]*par[5];
  double fitval = par[0]/ (2 * TMath::Pi() * par[2] * par[4] * sqrt(den))*exp(-(rx*rx+ry*ry-rz)/(2.*den));

  return fitval;
   }
else 
  return 0;
}


void fit(){

   TFile *f_image = new TFile("results_from_9072_to_9078.root","read");   
   int bins_per_decade=3;
   int bins= bins_per_decade*8; // 8=number of decades, since you want to cover energies from 10^-2 eV to 10^6 eV
   int intervals=bins+1; 
   Double_t *energy = new Double_t[intervals];       
   for(int i=0;i<intervals;i++) energy[i]= TMath::Power(10,TMath::Log10(1e-02)+(Float_t)i*(TMath::Log10(1e+06)-TMath::Log10(1e-02))/bins);
   double *bifs= new double[bins];
   double *bifsERR= new double[bins];
   double *bifb= new double[bins];
   double *bifbERR= new double[bins];
   
   char histoname[256];  
   for(int j=5;j<6;j++)   
      {  // Loop on j 
   sprintf(histoname,TString::Format("h_det%d",j));  
   TH2D *h_image = (TH2D*)f_image->Get(histoname); 
   
   //============================== FIT WITH A 2D GAUSSIAN =================
   printf("\n");
   printf("\n");
   printf ("\033[34;01m Fitting %s\033[00m\n", histoname);
   printf("\n");
   printf("\n");

   TF2 *g2D=new TF2("g2D",Gaus2D,h_image->GetMean(1)-2.,h_image->GetMean(1)+2.,h_image->GetMean(2)-2.,h_image->GetMean(2)+2.,6);
   g2D->SetParNames("Const","X_{0}","#sigma_{x}","Y_{0}","#sigma_{y}","#rho");
   g2D->SetNpy(1000); 
   g2D->SetParameters(300, 0, 0.8,0,0.8,0);
   double f=h_image->GetBinContent( h_image->GetXaxis()->FindBin(h_image->GetMean(1)),h_image->GetYaxis()->FindBin(h_image->GetMean(2))) ;
   if (f<100)h_image->Fit(g2D,"RLL0");
   if(f>100) h_image->Fit(g2D,"RL0");
   g2D->SetLineColor(2);
   delete h_image;
   g2D->ReleaseParameter(0);
   g2D->ReleaseParameter(1);
   g2D->ReleaseParameter(2);
   g2D->ReleaseParameter(3);
   g2D->ReleaseParameter(4);
   g2D->ReleaseParameter(5);
   delete g2D;
    
      }
      
}