//CHOOSE  1keV binning for consistency

#include <stdio.h> 
#include<locale.h> 

#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 x1_g1,double x2_g1 , double  xl, double  xr);
TH1D * CreateHisto();
double Err_f(double x);
double Gauss_PDF(double x);

// MAIN function
int TestMatrix1(){

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 x1_g1= 200. ;
double x2_g1= 300. ;

double LL = x1_g1 ;
double UL = x2_g1 ; //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();
int mypol_degree = 1 ;
FitIntegration(h1, mypol_degree , n,  BinWidth,  x1_g1, x2_g1 ,LL, UL);  

c1->Modified();
c1->Update();
return 0;
}



//FUNCTION DEFINITIONS
TH1D *CreateHisto(){
TH1D *hist = new TH1D("histo","Gaussian Peak and Background",1024,0,1024); 
fstream file;
file.open("/mnt/c/root-6.24.06-build/macros/Cs137.txt", ios::in); //ON UBUNTU
double Counts;
int channel=0;
while (file >> Counts) {	
//cout << channel << "   " << Counts << endl;
hist->Fill(channel,Counts);
channel++;
}
hist->Sumw2(kFALSE); 
hist->ResetStats();
cout<< "test: " <<hist->GetBinContent(250) << "  " <<hist->GetBinContent(251) << "  " << endl;
Double_t err_g1; cout << "raw histogram integral: " << hist->IntegralAndError( 250, 251 , err_g1) << " +- " << err_g1 << endl;
file.close();
//c1->Update();
hist->Draw("bar");  //bar

return hist;
}

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

cout<<"***********************************************************************************" <<endl;
cout<<"This is -FitIntegration- function" <<endl;
ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2");
TAxis *Xaxis = h1->GetXaxis();
/*
double x1_g1= 400. ;
double x2_g1= 462. ;
double x1_g2= 462. ;
double x2_g2= 524. ;
*/
double y1_g1= h1->GetBinContent((int)x1_g1);
double y2_g1= h1->GetBinContent((int)x2_g1);
double slope_g1= (y1_g1 - y2_g1) / (x1_g1-x2_g1) ;

cout <<"*****x1_g1= " << x1_g1 << " y1_g1= " << y1_g1 <<  endl;
cout <<"*****x2_g1= " << x2_g1 << " y2_g1= " << y2_g1 <<  endl;
cout <<"*****slope1= "<< slope_g1  << endl;

///////////////////////////////////////////////////////////////////////////////////////

TF1 * total;
 if (pol_degree==1) {
total = new TF1("total", "gausn(0) + pol1(3)", xl , xr);
total->SetParNames( "Area_1", "Mean_1", "Sigma_1", // gausn(0)
                     "b", "a"); // pol1(3)


total->SetParameters(11745., 247., 9.517, // gausn(0)
                    563.116  , -2.292 );  // y=b+ax 

/*
total->FixParameter(6, break_pointX);                 
total->FixParameter(7, break_pointY);                 
total->FixParameter(8,  slope_g1);                 
total->FixParameter(9,  slope_g2);                 
*/
}
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) + pol2(3)", xl, xr);
total->SetParNames("Area_1", "Mean_1", "Sigma_1", // gausn(0)
            "c", "b", "a"); // pol2(6)
total->SetParameters(11745., 247., 9.517, // gausn(0)
                     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); //https://root.cern.ch/doc/master/classTAttFill.html#ATTFILL2
h1->SetLineColor(kBlack); 
h1->GetXaxis()->SetRange(xl,xr);
h1->GetXaxis()->SetTitle("Channel Numbers");
h1->GetYaxis()->SetTitle("Counts");

TFitResultPtr r_total = h1->Fit(total, "EMRS", "L", xl, xr); // 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;
    
 if (pol_degree==1) {
 cov_bcg   = cov_total.GetSub(3, 4, 3, 4); //pol1
}
else if (pol_degree==2) {
 cov_bcg   = cov_total.GetSub(3, 5, 3, 5); //pol2
}
else {cout<<"check if something is wrong!" << endl;}

cov_total.Print();
cov_g1.Print();
cov_bcg.Print();



}