//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"
#include "Math/MinimizerOptions.h"
#include "TMatrixD.h"

using namespace std;

//FUNCTION DECLARATIONS
//variables declared inside function are called local variables.
Double_t GaussianPeak(Double_t *x, Double_t *p_total);
Double_t background(Double_t *x, Double_t *p_total);
Double_t totalFcn(Double_t *x, Double_t *p_total) ;
TH1D *CreateHisto();
void FitIntegration(TH1D *h1, double n, double BinWidth,  double  xl, double  xr, double &p1,  double &p2);

//*******************************************************************************************************************************************************
//*******************************************************************************************************************************************************
int test1(){
delete gROOT->FindObject("hist");
delete gROOT->FindObject("histo");
delete gROOT->FindObject("h1");
delete gROOT->FindObject("c1");
delete gROOT->FindObject("h2");
delete gROOT->FindObject("c2");
 //TCanvas::DrawFrame();
////////////Input to initialize the methods /////////
double LL = 402.;
double UL = 461.; // Keep that in mind: "bin width" = 1.
double n = UL-LL+1.; // number of bins
double BinWidth= (UL-LL+1)/n;
double mean_value,sigma_value;

cout<<"n: " << n << endl;
/////////////FUNCTION CALLINGS //////////////////////////////////////////////////
TCanvas *c1 = new TCanvas("c1", "c1", 800,400);
//c1->Divide(2,1);
//c1->cd(1);
gStyle->SetPalette(1);
gStyle->SetOptStat("inM");
TH1D *h1= CreateHisto(); // TH1F *CreateHisto();// TH1F (const char *name, const char *title, Int_t nbinsx, Double_t xlow, Double_t xup);
//when I say content(1), it means the content in the very 1st bin even if the numbering starts from 0
//cout<< "\nchecking: "<<LL <<" bin content: " << h1->GetBinContent(LL) <<", " << UL << " .bin content: " << h1->GetBinContent(UL) << endl;
//cout << "total: " << h1->GetBinContent(LL) +h1->GetBinContent(UL)  << endl;
h1->SetMarkerStyle(21);
h1->SetMarkerSize(0.8);
h1->SetStats(0);
h1->SetFillStyle(0);
h1->SetLineColor(kBlack); 
h1->GetXaxis()->SetRange(LL,UL);
h1->Draw(""); //bar
/////////////FUNCTION CALLINGS //////////////////////////////////////////////////
FitIntegration(h1, n, BinWidth, LL, UL, mean_value, sigma_value); 

c1->Update();

return 0;
}

//*******************************************************************************************************************************************************
//*******************************************************************************************************************************************************

//FUNCTION DEFINITIONS
// Gaussian Peak function
Double_t GaussianPeak(Double_t *x, Double_t *p_total) {
  return p_total[0]*exp(-0.5*TMath::Power((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 peak function
Double_t totalFcn(Double_t *x, Double_t *p_total) {
  return   GaussianPeak(x, p_total )+background(x, &p_total[3]) ;//p_total= &p_total[0]
}


TH1D *CreateHisto(){

// 	TH1F (const char *name, const char *title, Int_t nbinsx, Double_t xlow, Double_t xup)
TH1D *hist = new TH1D("histo","Gaussian Peak and Background",1024,0,1024); // you need to draw the histogram only 1 time in the analysis

ifstream file;
file.open("Co60.txt", ios::in); 
double Counts;
int channel=0;
while (file >> Counts) {	
  // cout << channel << "   " << Counts << endl;
  hist->Fill(channel,Counts);
  channel++;
}
hist->Sumw2(kFALSE); 
hist->ResetStats();
  //cout<< "checking: 1. bin content: " << hist->GetBinContent(0) << " last bin content: " << hist->GetBinContent(1023) << endl; //wrong
  //cout<< "checking: 1. bin content: " << hist->GetBinContent(1) << " last bin content: " << hist->GetBinContent(1024) << endl; //right

//cout <<"Total Chanel Number: " << channel  << endl;
/* Convention for numbering bins
For all histogram types: nbins, xlow, xup
bin = 0;       underflow bin
bin = 1;       first bin with low-edge xlow INCLUDED
bin = nbins;   last bin with upper-edge xup EXCLUDED
bin = nbins+1; overflow bin
*/

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

void FitIntegration(TH1D *h1, double n, double BinWidth, double  xl, double  xr, double &p1, double &p2){

cout<<"***********************************************************************************" <<endl;
cout<<"This is -FitIntegration- function" <<endl;
ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2");

//BACKGROUND from TPA method same as y=ax+b fit
TAxis *Xaxis = h1->GetXaxis();
double Bcg_x1_LowE= Xaxis->GetBinLowEdge(xl );
double Bcg_x1_UpE = Xaxis->GetBinUpEdge(xl); // finding the lower bin edge value on X axis = Bcg_x1_LowE +1 ;
double Bcg_x2_LowE= Xaxis->GetBinLowEdge(xr ); 
double Bcg_x2_UpE = Xaxis->GetBinUpEdge(xr ); //upper bin edge value  = //Bcg_x2_LowE +1 
cout << Bcg_x1_LowE << "  "<<Bcg_x1_UpE << " " << Bcg_x2_LowE << " " << Bcg_x2_UpE << endl;
double Bcg_y1= h1->GetBinContent(xl); //bin width =1 , so low edge and upper edge difference is also 1.
double Bcg_y2= h1->GetBinContent(xr); 
double Bcg_y1_error= h1->GetBinError(xl);
double Bcg_y2_error= h1->GetBinError(xr);
double Bcg_slope= (Bcg_y1-Bcg_y2)/(Bcg_x1_LowE-Bcg_x2_UpE);
double Bcg_a= Bcg_y1 - (Bcg_slope* Bcg_x1_LowE) ;  //if y=a+bx, then a= y-bx.
double Bcg_b= Bcg_slope;  //y=a+bx
cout << "pol1 parameters: a= " << Bcg_a << "  b= " << Bcg_b << endl;
//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
total->SetParNames("Height_1", "Mean_1", "Sigma_1", "a", "b"); // pol1(6) = a + bx =[0]+x*[1]
total->SetNpx(5 * n); // e.g., (n) or (3 * n);
// f->Print();
total->SetParameters(2184., 433., 10., Bcg_a, Bcg_b); //a+bx
// 1589.44 , -3.11111 ); // 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, 0., 2000.); //for pol
///total->SetParLimits(4, -10., 0. );  //for pol 
//total->SetParLimits(5, -0.5 , -0.205 ); //for pol

//NOT= If you fix the parameter, you get zero for the error.
//total->FixParameter(0, 2184.); //height
//total->FixParameter(1, 433. ); //mean
//total->FixParameter(2, 10.  ); //sigma
//total->FixParameter(3, Bcg_a); //for pol
//total->FixParameter(4, Bcg_b); //for pol
//*************************************************************************


// Fitting and retrieving the covariance matrix
//For fit options: https://root.cern.ch/root/htmldoc/guides/users-guide/FittingHistograms.html
//When fitting (make “background” parameters “free”), you can try to play with fit options: "EMRS", "PEMRS", "WEMRS", "LEMRS" ("L" only for histograms)
TFitResultPtr r_total = h1->Fit(total, "EMRS"); // 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_bcg   = cov_total.GetSub(3, 4, 3, 4); //pol1(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; // gausn(0)  p_total= p_total[0], address on the left = address on the right
Double_t *p_bcg    = p_total + 3; // pol1(6) or "broken line" (6)
// & array[1]=array+1   //array[1]=*(array+1)
//*************************************************************************
TF1 *g1 = new TF1("g1", GaussianPeak, xl, xr, 3); 
g1->SetNpx(total->GetNpx());
g1->SetParameters(p_total + 0);
//dereferencing a pointer gives the value where it points to.
double p0= *(p_total + 0); //height  _Gauss 1    
       p1= *(p_total + 1); //mean
       p2= *(p_total + 2); //sigma
//************************************************************************* 
TF1 *bcg = new TF1("bcg", background, xl, xr, 2); //pol2(3)
bcg->SetNpx(total->GetNpx());  
bcg->SetParameters(p_total + 3); // pol1(6)  a + bx =[0]+x[1]
//dereferencing a pointer gives the value where it points to.
double p3= *(p_total + 3); // Pol1  a+bx, a
double p4= *(p_total + 4); // Pol1  a+bx, b
//Drawing the results  **************************************************** 
total->SetLineColor(kRed); total->Draw("same");
g1->SetLineColor(kGreen);  g1->Draw("same");
bcg->SetLineColor(kBlue);  bcg->Draw("same");
///////////////////////////////////////////////////////

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;

Double_t err; cout << "raw histogram integral: " << h1->IntegralAndError(xl, xr, err) << " +- " << err << endl;

////Gross AREA from fitted function ///////////////////////////////////
double GrossArea1 = total->Integral( g1_X_LBE, g1_X_UBE) / BinWidth ;
double Error_GA1  = total->IntegralError(g1_X_LBE, g1_X_UBE , p_total, cov_total.GetMatrixArray(), 1e-7) / BinWidth;
cout <<"Gross_area: "<< GrossArea1 << " +- "<< Error_GA1 << endl;

//Background AREA from fitted function //////////////////////////////////////////////////////////////////////////////////////////
// bcg->SetParameters(p_bcg);
double bcg_g1       = bcg->Integral(g1_X_LBE, g1_X_UBE) / BinWidth ;	
double bcg_g1_error = bcg->IntegralError( g1_X_LBE, g1_X_UBE , p_bcg , cov_bcg.GetMatrixArray(), 1e-7) / BinWidth;
cout <<"Bcg_area: "<< bcg_g1 << " +- "<< bcg_g1_error << endl; //WRONG!!!!

////////NET AREA from fitted function/////////////////////////////////////////////////////////////////////////////////////////////////////////
// g1->SetParameters(p_g1);
double NetArea      = g1->Integral(g1_X_LBE, g1_X_UBE) / BinWidth ;	
double NetAreaError = g1->IntegralError( g1_X_LBE, g1_X_UBE , p_g1, cov_g1.GetMatrixArray(), 1e-7) / BinWidth;

cout << "Net_area: "  << NetArea << " +- "  <<  NetAreaError <<  endl;
cout<<"***********************************************************************************" <<endl;

}