//CHOOSE  1keV binning for consistency

#include <iostream>
#include "TCanvas.h"
//#include "TLatex.h"
//#include "TText.h"
#include "TDirectory.h"
#include "TAxis.h"
//#include <string>
#include "TH1D.h"
#include "TAxis.h"
//#include "TFile.h"
#include "TStyle.h"  //for gStyle decleration
#include <fstream>
//#include "TLine.h"
//#include "TLegend.h"
//#include "AnalyticalIntegrals.h"
#include "TROOT.h"
#include "TF1.h"
#include "TFormula.h"
#include "TMath.h"
#include "Math/DistFuncMathCore.h"
//#include <stdio.h>
#include "TFitResult.h"


TH1D *CreateHisto(int BinWidth,int xlow,int xup){

TCanvas *c1 = new TCanvas("c1", "c1", 800,400);
c1->cd();
gStyle->SetPalette(1);
gStyle->SetOptStat("inM");
// 	TH1F (const char *name, const char *title, Int_t nbinsx, Double_t xlow, Double_t xup)
TH1D *hist = new TH1D("histo","Gaussian Peak on Quadratic Background",(xup-xlow)/BinWidth,xlow,xup); 
// you need to draw the histogram only 1 time in the analysis
fstream file;
file.open("C:\\Users\\ilker\\OneDrive\\OneDrive - harran.edu.tr\\ROOTPRO\\Co60.txt", ios::in); 
double Counts;
int channel=0;
while(!file.eof()){	
	file >> Counts ;
  	//cout << channel << "   " << Counts << endl;
	hist->Fill(channel,Counts);
	channel++;
	if(file.eof()) break;
}
//cout <<"Total Chanel Number: " << channel  << endl;
 

file.close();
c1->Update();
return hist;
}


void pm(){
   cout << "Integral:" << 100<< "±" << 5<< endl;
}
// Gaussian Peak function
Double_t GaussianPeak(Double_t *x, Double_t *p_total) {
  return p_total[0]*exp(-0.5*((x[0]-p_total[1])/p_total[2])**2) ;
}

// Quadratic background function
Double_t background(Double_t *x, Double_t *p_total) {
   return p_total[0] + p_total[1]*x[0] + p_total[2]*x[0]*x[0];
   //return p_total[0] + p_total[1]*x[0] ;
}


// Sum of background and Gauss function
Double_t totalFcn(Double_t *x, Double_t *p_total) {
  return   GaussianPeak(x, &p_total[0] )+background(x, &p_total[3]) ;//p_total= &p_total[0]
}


void FitGussianFunctions(TH1D *h1, double  xl, double  xr){
ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2");
double  n = xr-xl+1.; // number of bins
double BinWidth=  (xr-xl+1.)/n ;
//First we fit the combination of all functions corresponding to signal+background, then we retrieve them back one by one later from the fit parameters we find here.
//TF1 *total = new TF1("total", totalFcn , xl,xr, 5); //for pol1
TF1 *total = new TF1("total", totalFcn , xl,xr, 6); //for pol1

/////////////////////////////////////////////////////////////////////////////////////////////////////////
//TH1F *h = (TH1F*)h1->Clone("h");
  total->SetParNames("Height_1", "Mean_1", "Sigma_1", // gaus(0) // if gausn(0), the names would be "Area_1", "Mean_1", "Sigma_1".
                      "a", "b", "c"); // pol1(6) = a + bx =[0]+x*[1]
  total->SetNpx(n); // e.g., (n) or (3 * n);
  // f->Print();
  total->SetParameters(2184., 433., 10., // gaus(0)
                       //935.134,  -1.573); // pol1(6)  a + bx =[0]+x[1]
                      -36342. , 173.63 , -0.205); //calculated by hand for pol2 , bump up 
                       //2252.5 ,  1.58 , -0.004);  //for pol2  // bump down
  //total->SetParLimits(3,-40000.,-30000.); //for pol
  //total->SetParLimits(4,100.,200.);  //for pol 
  //total->SetParLimits(5, -1., 0. ); //for pol

  total->FixParameter(3, -36342.); //for pol
  total->FixParameter(4, 173.63);  //for pol
  total->FixParameter(5, -0.205);  //for pol

h1->GetXaxis()->SetRange(xl,xr); // for fitting to focus on the limits
// Fitting and retrieving the covariance matrix
TFitResultPtr r_total = h1->Fit(total, "LRS"); // fit should be done using the fit option S.
TMatrixD cov_total = r_total->GetCovarianceMatrix(); //Cov matrix is always square and symmetric matrix
TMatrixD cov_g1 = cov_total.GetSub(0, 2, 0, 2); // gaus(0)
//TMatrixD cov_bg = cov_total.GetSub(3, 4, 3, 4); //pol1(3)
TMatrixD cov_bg = cov_total.GetSub(3, 5, 3, 5); //pol2(3)

////////////// ////////////////////////////////////////////////////////////
Double_t *p_total = total->GetParameters(); // p_total is a pointer which refers to the address of a double variable
Double_t *p_g1 = p_total + 0; // p_total= p_total[0], address on the left = address on the right
Double_t *p_bg = p_total + 3; 
// & array[1]=array+1
//array[1]=*(array+1)

TF1 *g1 = new TF1(Form("g1"),GaussianPeak, xl, xr,3 ); 
g1->SetNpx(total->GetNpx());
g1->SetParameters(p_total + 0); //virtual void TF1::SetParameters	(	const Double_t * 	params	)	
TF1 *bg = new TF1(Form("bg"),background , xl, xr, 3); //pol2(3)
bg->SetNpx(total->GetNpx());  
bg->SetParameters(p_total + 3); // pol1(6)  a + bx =[0]+x[1]
/*
//dereferencing a pointer gives the value where it points to.
double p0= *(p_total + 0); //Gauss 1
double p1= *(p_total + 1); 
double p2= *(p_total + 2); 
double p3= *(p_total + 3); // Pol1
double p4= *(p_total + 4);
double p5= *(p_total + 5); // for pol2
*/
//Drawing the results ////////////////////////////////////
total->SetLineColor(kRed); total->Draw("same");
g1->SetLineColor(kGreen); g1->Draw("same");
bg->SetLineColor(kBlack); bg->Draw("same");
///////////////////////////////////////////////////////

TAxis *Xaxis = h1->GetXaxis();
double  g1_X_LBE = Xaxis->GetBinLowEdge(Xaxis->FindBin(xl) ); // finding the lower bin edge value on X axis
double  g1_X_UBE = Xaxis->GetBinUpEdge(Xaxis->FindBin(xr) ); //upper bin edge value
printf("\nIntegral range on X axis:");
cout<< "\nxl: " << xl << " xr: " <<  xr << "\ng1_X_LowerBinEdge: " <<  g1_X_LBE<< " g1_X_UpperBinEdge: " <<  g1_X_UBE << endl;


//Background AREA from fittedfit //////////////////////////////////////////////////////////////////////////////////////////
//cout<<"check from function: "<< bg->Integral(g1_X_LBE, g1_X_LBE+1, p_bg) / BinWidth << "  and  " <<  bg->Integral(g1_X_UBE, g1_X_UBE+1, p_bg) / BinWidth << endl;
double bcg_g1 = bg->Integral(g1_X_LBE, g1_X_UBE, p_bg) / BinWidth ;	
double bcg_g1_error = bg->IntegralError( g1_X_LBE, g1_X_UBE , p_bg, cov_bg.GetMatrixArray()) / BinWidth;
cout <<"Bcg area from bcg pol. function: "<< bcg_g1 << " +- "<< bcg_g1_error << endl;


/////////////Background AREA from histogram////////////////////////////////////////////////////////////////////////////
int NBins_g1 = Xaxis->FindBin(xr)- Xaxis->FindBin(xl) +1;
//cout<<"check from histogram: "<<  h1->GetBinContent(Xaxis->FindBin(xl)-1) <<"  and  " << h1->GetBinContent(Xaxis->FindBin(xr)+1) << endl;
double BCG_g1 =( h1->GetBinContent(Xaxis->FindBin(xl)-1) + h1->GetBinContent(Xaxis->FindBin(xr)+1) )*0.5*NBins_g1;
double BCG_g1_Error =( 0.5*NBins_g1 )* sqrt ( h1->GetBinContent(Xaxis->FindBin(xl)-1) + h1->GetBinContent(Xaxis->FindBin(xr)+1) ) ;
//cout<<"checking:  "<< double( h1->GetBinContent(410)) << "  " << h1->GetBinError(410) << " " << sqrt(double( h1->GetBinContent(410)))  << endl;
cout <<"BCG Area from histo: " <<BCG_g1  <<" +- " << BCG_g1_Error <<  endl;


}
int FindFittingValues(){
delete gROOT->FindObject("hist");
delete gROOT->FindObject("histo");
delete gROOT->FindObject("h1");
delete gROOT->FindObject("c1");

////////////Input to initialize the methods /////////
double LL = 400.;
double UL = 462.; // Keep that in mind: "bin width" = 1.
/////////////////////////////////////////////////////////////
TH1D *h1= CreateHisto(1, 0, 1024); // TH1F *CreateHisto();// TH1F (const char *name, const char *title, Int_t nbinsx, Double_t xlow, Double_t xup);
h1->SetMarkerStyle(21);
h1->SetMarkerSize(0.8);
h1->SetStats(0);
h1->SetFillStyle(0);
h1->Draw("bar"); //px
FitGussianFunctions(h1, LL, UL);


return 0;
}