//CHOOSE  1keV binning for consistency

#include <iostream>
#include "TCanvas.h"
#include "TFitResult.h"
#include "TDirectory.h"
#include "TAxis.h"
#include "TH1D.h"
#include "TAxis.h"
#include "TStyle.h"  //for gStyle decleration
#include <fstream>
#include "TF1.h"
#include "TFormula.h"
#include "TMath.h"
#include "Math/DistFuncMathCore.h"
#include "TROOT.h"
#include "Math/Minimizer.h"
#include "Minuit2/MnUserParameterState.h"
#include "TMatrixD.h"


using namespace std;
//FUNCTION DECLARATION
void FitIntegration(TH1D *h1, int pol_degree, double n, double BinWidth,  double  xl, double  xr);
TH1D * CreateHisto();

// MAIN function
int DigitalIntegrationProcedures(){

delete gROOT->FindObject("hist");
delete gROOT->FindObject("histo");
delete gROOT->FindObject("bcg_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 = 400.;
double UL = 530. ; //462.; // Keep that in mind: "bin width" = 1.
double n = UL-LL+1.; // number of bins
double BinWidth= 1.;

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();
FitIntegration(h1, 1, n,  BinWidth, LL, UL);  

c1->Modified();
c1->Update();

return 0;
}

//FUNCTION DEFINITIONS
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); 
fstream file;
file.open("/mnt/c/root-6.24.06-build/macros/Co60.txt", ios::in); //ON UBUNTU
///file.open("C:\\root-6.24.06-build\\macros\\Co60.txt", ios::in); //ON WINDOWS
double Counts;
int channel=0;
while (file >> Counts) {	
// cout << channel << "   " << Counts << endl;
hist->Fill(channel,Counts);
channel++;
}
hist->Sumw2(kFALSE); 
hist->ResetStats();
file.close();
//c1->Update();
return hist;
}


void FitIntegration(TH1D* h1, int pol_degree, double n, double BinWidth, double  xl, double  xr){

cout<<"***********************************************************************************" <<endl;
cout<<"This is -FitIntegration- function" <<endl;
ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2");
//Int_t n = 140; // number of bins
//Double_t xl = 400., xr = xl + n; // make sure "bin width" = 1.
TF1 * total;
 if (pol_degree==1) {
total = new TF1("total", "gausn(0) + gausn(3) + ((x < [6]) * [8] * (x - [6])) + [7]", xl, xr);
total->SetParNames("Area_1", "Mean_1", "Sigma_1", // gausn(0)
            "Area_2", "Mean_2", "Sigma_2", // gausn(3)
            "Break_Point", "constant", "slope"); // "broken line" background (6)  y= mx+ constant
total->SetParameters(5.e4, 430., 10., // gausn(0)
                     4.e4, 490., 10., // gausn(3)
                     460.,  167. , -2.22); // "broken line" background (6) 
}
else if (pol_degree==2) {
//“polN” A polynomial of degree N, where N is a number between 0 and 9: f(x) = p0 + p1*x + p2*x2 +...
total = new TF1("total", "gausn(0) + gausn(3) + pol2(6)", xl, xr);
total->SetParNames("Area_1", "Mean_1", "Sigma_1", // gausn(0)
            "Area_2", "Mean_2", "Sigma_2", // gausn(3)
            "c", "b", "a"); // pol2(6)
total->SetParameters(5.e4, 430., 10., // gausn(0)
                     4.e4, 490., 10., // gausn(3)
                     400., 0.4, -0.002); // pol2(6) = f(x) = p0 + p1*x + p2*x2
}
else {cout << "check your pol degree" << endl; }

total->SetLineColor(kRed); total->SetNpx(n); // e.g., (n) or (3 * n);
//total->Print();
//h->FillRandom("total", Int_t(total->Integral(xl, xr) + 0.5));
h1= CreateHisto(); 
h1->SetMarkerStyle(21);
h1->SetMarkerSize(0.8);
h1->SetStats(0);
h1->SetFillStyle(0);
h1->SetLineColor(kBlack); 
h1->GetXaxis()->SetRange(xl,xr);
h1->Draw("");  //bar
//h1->Fit(total, "L"); // e.g., "" or "L"

//***************************************************************************************
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_g2    = cov_total.GetSub(3, 5, 3, 5);
TMatrixD cov_bcg   = cov_total.GetSub(6, 8, 6, 8); //pol2(6) or "broken line" (6)

Double_t *p_total  = total->GetParameters(); // p_total is a pointer which HOLDS 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_g2     = p_total + 3; // gausn(3)
Double_t *p_bcg    = p_total + 6; // pol2(6) or "broken line" (6)
total->SetLineColor(kRed); total->Draw("same");

// & array[1]=array+1   //array[1]=*(array+1)
//*************************************************************************
TF1 *g1 = new TF1("g1", "gausn", xl, xr); 
g1->SetNpx(total->GetNpx());
g1->SetParameters(total->GetParameters() + 0); // gausn(0)  
g1->SetLineColor(kGreen);  g1->Draw("same");

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

TF1 *g2 = new TF1("g2", "gausn", xl, xr); 
g2->SetNpx(total->GetNpx());
g2->SetParameters(total->GetParameters() + 3); // gausn(3).
g2->SetLineColor(kBlue); g2->Draw("same");
//*************************************************************************

TF1 *bcg;
if (pol_degree==1) {bcg = new TF1("bcg", "((x < [0]) * [2] * (x - [0])) + [1]", xl, xr);}
else if (pol_degree==2){ bcg = new TF1("bcg", "pol2", xl, xr);}
else {cout << "check your pol degree" << endl; }
bcg->SetNpx(total->GetNpx());
bcg->SetParameters(total->GetParameters() + 6); // any background (6)
bcg->SetLineColor(kBlack); bcg->Draw("same");
//ROI =Range of interest in a peak
double x1_g1= 400. ;
double x2_g1= 460. ;
double x1_g2= 460. ;
double x2_g2= 524. ;
//////Calculating the correct bin limits  /////////////////////////////////////////
TAxis *Xaxis = h1->GetXaxis();
cout <<"*****************" << h1->GetBinContent( (int)x1_g1) << " " << h1->GetBinContent((int)x2_g1) 
<< " " << h1->GetBinContent((int)x2_g2) << " " << endl;

double  g1_X_LBE = Xaxis->GetBinLowEdge(Xaxis->FindBin(x1_g1) ); //lower bin edge value on x1
double  g1_X_UBE = Xaxis->GetBinUpEdge(Xaxis->FindBin(x2_g1) ); //upper bin edge value on x1


double  g2_X_LBE = Xaxis->GetBinLowEdge(Xaxis->FindBin(x1_g2) ); //lower bin edge value on x1
double  g2_X_UBE = Xaxis->GetBinUpEdge(Xaxis->FindBin(x2_g2) ); //upper bin edge value on x1
//*****************************************************************************************
printf("\nIntegral range on X axis:");
cout<< "\nxl_g1: " << x1_g1 << " xr_g1: " <<  x2_g1 << "\ng1_X_LowerBinEdge: " <<  g1_X_LBE<< " g1_X_UpperBinEdge: " <<  g1_X_UBE << endl;
Double_t err_g1; cout << "raw histogram integral: " << h1->IntegralAndError((int)x1_g1, (int)x2_g1 , err_g1) << " +- " << err_g1 << endl;
///////////////////////////////////////////////////////////////////////
cout<< "\nxl_g2: " << x1_g2 << " xr_g2: " <<  x2_g2 << "\ng2_X_LowerBinEdge: " <<  g2_X_LBE<< " g2_X_UpperBinEdge: " <<  g2_X_UBE << endl;
Double_t err_g2; cout << "raw histogram integral: " << h1->IntegralAndError((int)x1_g2, (int)x2_g2 , err_g2) << " +- " << err_g2 << endl;


cout<<"********************************PEAK 1***************************************************"<< endl;
////Gross AREA from fitted function ///////////////////////////////////
double Gross_Area_g1    = total->Integral( g1_X_LBE, g1_X_UBE) / BinWidth ;// ConfidenceLimitCorrection;  //virtual Double_t Integral (Double_t a, Double_t b, Double_t epsrel=1.e-12)	
double Error_gross_g1   = total->IntegralError(g1_X_LBE, g1_X_UBE , p_total, cov_total.GetMatrixArray(),1e-7) / BinWidth;
cout <<"Gross_area_g1: "<< Gross_Area_g1 <<  " ± " <<  Error_gross_g1  << endl;
////Net AREA from fitted function ///////////////////////////////////
double Net_Area_g1 = g1->Integral( g1_X_LBE, g1_X_UBE) / BinWidth ;// ConfidenceLimitCorrection;  //virtual Double_t Integral (Double_t a, Double_t b, Double_t epsrel=1.e-12)	
double Error_g1    = g1->IntegralError(g1_X_LBE, g1_X_UBE , p_g1, cov_g1.GetMatrixArray(),1e-7) / BinWidth;
cout <<"Net_area_g1: "<< Net_Area_g1 << " ± " << Error_g1 << endl;
//Background AREA from fitted function //////////////////////////////////////////////////////////////////////////////////////////
//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       = bcg->Integral(g1_X_LBE, g1_X_UBE) / BinWidth ;	//bcg->Integral(g1_X_LBE, g1_X_UBE, p_bcg) / BinWidth ;	
double bcg_Error_g1 = bcg->IntegralError( g1_X_LBE, g1_X_UBE , p_bcg , cov_bcg.GetMatrixArray(),1e-7) / BinWidth;
cout <<"Bcg_area_g1: "<< bcg_g1 << " ± " << bcg_Error_g1 << endl; 

cout<<"********************************PEAK 2***************************************************"<< endl;

////Gross AREA from fitted function ///////////////////////////////////
double Gross_Area_g2    = total->Integral( g2_X_LBE, g2_X_UBE) / BinWidth ;// ConfidenceLimitCorrection;  //virtual Double_t Integral (Double_t a, Double_t b, Double_t epsrel=1.e-12)	
double Error_gross_g2   = total->IntegralError(g2_X_LBE, g2_X_UBE , p_total, cov_total.GetMatrixArray(),1e-7) / BinWidth;
cout <<"Gross_area_g2: "<< Gross_Area_g2 <<  " ± " <<  Error_gross_g2  << endl;
////Net AREA from fitted function ///////////////////////////////////
double Net_Area_g2 = g2->Integral( g2_X_LBE, g2_X_UBE) / BinWidth ;// ConfidenceLimitCorrection;  //virtual Double_t Integral (Double_t a, Double_t b, Double_t epsrel=1.e-12)	
double Error_g2    = g2->IntegralError(g2_X_LBE, g2_X_UBE , p_g2, cov_g2.GetMatrixArray(),1e-7) / BinWidth;
cout <<"Net_area_g2: "<< Net_Area_g2 << " ± " << Error_g2 << endl;
//Background AREA from fitted function //////////////////////////////////////////////////////////////////////////////////////////
//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_g2       = bcg->Integral(g2_X_LBE, g2_X_UBE) / BinWidth ;	//bcg->Integral(g1_X_LBE, g1_X_UBE, p_bcg) / BinWidth ;	
double bcg_Error_g2 = bcg->IntegralError( g2_X_LBE, g2_X_UBE , p_bcg , cov_bcg.GetMatrixArray(),1e-7) / BinWidth;
cout <<"Bcg_area_g2: "<< bcg_g2 << " ± " << bcg_Error_g2 << endl; 

cout<<"******************************** THE END ***************************************************"<< endl;

}