#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 <TSystem.h>
#include "TPaletteAxis.h"
#include "TPaveStats.h"
#include "TROOT.h"
#include "TF2.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]*exp(-(rx*rx+ry*ry-rz)/(2.*den));
  return fitval;
   }
else 
  return 0;
}




void run() {

//======================================  STYLE  ===========================================//
   gStyle->SetOptDate(0);
   gROOT->SetStyle("Plain"); // It forces white background everywhere
   gStyle->SetPalette(1, 0);
   gStyle->SetOptStat(2201);
   TH1::AddDirectory(kFALSE);  // <- This avoids the messages "memory leak" when plotting more than one histogram
   TH2::AddDirectory(kFALSE);  // <- This avoids the messages "memory leak" when plotting more than one histogram


// =================================== RETRIVE DATA ==========================================//

       TFile *f = new TFile("Beam_image.root","read");
       TH2D*h2D=(TH2D*)f->Get("image2D_beam_001_01");
       h2D->GetXaxis()->SetTitle("Horizontal axis");
       h2D->GetYaxis()->SetTitle("Vertical axis");
       h2D->GetZaxis()->SetTitle("Counts");
       TCanvas *c1=new TCanvas("c1","Histo");
       c1->cd();
       h2D->Draw("lego2");

       
   
//====================Fitting the histogram =================================================//


        TF2 *g2D=new TF2("g2d",Gaus2D,-4,4,-4,4,6);
        g2D->SetParNames("Const","X_{0}","#sigma_{x}","Y_{0}","#sigma_{y}","#rho");
        g2D->SetParameters(100,-0.1,0.02,-0.1,0.02,0.);
        h2D->Fit("g2d","R0+");
	g2D->Draw("surf same");
       
       cout << "Chisquare given by root " << g2D->GetChisquare()<<endl;
       cout << "NDF given by root " << g2D->GetNDF()<<endl;  
       cout << "Number of free parameters " << g2D->GetNumberFreeParameters()<<endl;
       cout << "Number of x bins " << h2D->GetXaxis()->GetNbins() << endl;
       cout << "Number of y bins " << h2D->GetYaxis()->GetNbins() << endl;
       
       double chi=0;
       int num_points=0;
       double expected=0;
       double observed=0;
       
       for(int m=1;m<h2D->GetXaxis()->GetNbins()+1;m++ )
       {
	 for(int n=1;n<h2D->GetYaxis()->GetNbins()+1;n++ )
	 //if(h_image->GetBinContent(m,n)>0) 
	 {
	 if(g2D->Eval(h2D->GetXaxis()->GetBinCenter(m),h2D->GetYaxis()->GetBinCenter(n))>0)
	 {
	   num_points++;
//	   expected=g2D->Eval(h2D->GetXaxis()->GetBinLowEdge(m),h2D->GetYaxis()->GetBinLowEdge(n));
	   expected=g2D->Eval(h2D->GetXaxis()->GetBinCenter(m),h2D->GetYaxis()->GetBinCenter(n));
           observed=h2D->GetBinContent(m,n);
           chi=chi+( ( (observed-expected)*(observed-expected) )/expected);
	 }
	 if(g2D->Eval(h2D->GetXaxis()->GetBinCenter(m),h2D->GetYaxis()->GetBinCenter(n))==0) cout << m << " "  << n << endl;
	 }
        }
       cout << "Calculated chisquare " << chi << endl;
       cout << "Number of points used to calculate the chisquare " << num_points << endl;
       cout << "Calculated number of degrees of freedom " << num_points-g2D->GetNumberFreeParameters() -1 << endl;

  
}