//SetANGLE  "cThetaMax<110"  for 90degree Clover
//SetANGLE  mystring  for 135degree Clover
//CHOOSE  1keV binning for consistency

#include <iostream>
#include "TAttMarker.h"
#include "TCanvas.h"
#include "TLatex.h"
#include "TText.h"
#include "TDirectory.h"
#include "TAxis.h"
#include <string>
#include "TH1F.h"
#include "TF1.h"
#include "TAxis.h"
#include "TH2F.h"
#include "TTree.h"
#include "TFile.h"
#include "TChain.h"
#include "TStyle.h"  //for gStyle decleration
#include "TNtuple.h"
#include <fstream>
#include "TLine.h"
#include "TLegend.h"
#include "TH2F.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 "TMatrixDsymfwd.h" //https://root.cern.ch/doc/master/TMatrixDSymfwd_8h.html#a614c344e23528bf63193e7e369df7b77

TH1F *CreateHisto(int n, double xl, double xr){

TCanvas *c1 = new TCanvas("c1", "c1", 800,400);
c1->cd();
gStyle->SetPalette(1);
gStyle->SetOptStat("inM");

TH1F *hist = new TH1F("hist", "h;Channel;Counts", n, xl, xr); // you need to draw the histogram only 1 time in the analysis
fstream file;
file.open("C:\\Users\\Spyhunter\\Desktop\\Co60.txt", ios::in); 
//file.open("/home/ic00025/Desktop/Co60.txt", ios::in); 
double Counts;
int channel=0;
while(!file.eof()){	
	file >> Counts ;
	hist->Fill(channel,Counts);
	channel++;
	if(file.eof()) break;
	//cout << channel << "   " << Counts << endl;
}
cout <<"Total Chanel Number: " << channel  << endl;
file.close();
c1->Update();
return hist;
}

double Err_f(double x){
//What they call “erf(x)” in the “Table 1.1” and equations “(1.6)”, “(1.7)” and “(1.8)” is:
double y= ROOT::Math::erf(x / TMath::Sqrt(2.)) ;
return y;
}

void pm(){
   cout << "Integral:" << 100<< "±" << 5<< endl;
}

void FitGussianFunctions(TH1F *h1, int n, double  xl, double  xr){
	
ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2");
//double extention =20. ; //For background  // Gilmore book page 112, figure 5.7
double HistoBinning= (xr-xl)/n;

TF1 *total = new TF1("total","gausn(0) + gausn(3) + pol2(6)", xl, xr);
/////////////////////////////////////////////////////////////////////////////////////////////////////////
//TH1F *h = (TH1F*)h1->Clone("h");
  total->SetParNames("Height_1", "Mean_1", "Sigma_1", // gausn(0)
                     "Height_2", "Mean_2", "Sigma_2", // gausn(3)
                      "a", "b", "c"); // pol2(6)
  total->SetNpx(n); // e.g., (n) or (3 * n);
  // f->Print();
  total->SetParameters(21000., 430., 10., // gausn(0)
                       18000., 490., 10., // gausn(3)
                       400., 0.4, -0.002); // pol2(6)

  h1->GetXaxis()->SetRange(380,540); // for user to focus
  
//https://root.cern.ch/doc/master/classTF1.html#a4f798c9626dc2d2268198ee01dc51437     
///////////////////////////////////////////////////////////////////////////////
  TF1 *g1 = new TF1("g1", "gausn", xl, xr); g1->SetNpx(total->GetNpx());
  g1->SetParameters(total->GetParameters() + 0); // gausn(0)

  TF1 *g2 = new TF1("g2", "gausn", xl, xr); g2->SetNpx(total->GetNpx());
  g2->SetParameters(total->GetParameters() + 3); // gausn(3)
  
  TF1 *bg = new TF1("bg", "pol2", xl, xr); bg->SetNpx(total->GetNpx());
  bg->SetParameters(total->GetParameters() + 6); // pol2(6)

////Getting the parameters individually ///////////////////////

double p0= total->GetParameter(0); //Gaus 1
double p1= total->GetParameter(1); 
double p2= total->GetParameter(2); 
double p3= total->GetParameter(3); //Gaus 2
double p4= total->GetParameter(4); 
double p5= total->GetParameter(5); 
double p6= total->GetParameter(6); //Pol2
double p7= total->GetParameter(7);  
double p8= total->GetParameter(8); 
//////////////////////////////////////
double* g1Parameters = new double[3] ;
g1Parameters[0] = p0;
g1Parameters[1] = p1;
g1Parameters[2] = p2;
double* g2Parameters = new double[3] ;
g2Parameters[0] = p3;
g2Parameters[1] = p4;
g2Parameters[2] = p5;
double* bgParameters = new double[3] ;
bgParameters[0] = p6;
bgParameters[1] = p7;
bgParameters[2] = p8;
////////////////////////////////////////////////////////
// TFitResultPtr contains the TFitResult.
TFitResultPtr r0 = hist->Fit(total,"LS", "same", p1-3*p2, p4+3*p5);
TFitResultPtr r1 = hist->Fit(g1,"LS", "same", p1-2*p2, p1+2*p2);
TFitResultPtr r2 = hist->Fit(g2,"LS", "same", p4-2*p5, p4+2*p5);
//TFitResultPtr r3 = hist->Fit(bg,"LS", "same", p1-4*p2, p4+4*p5);
// Access the covariance matrix.
TMatrixDSym covMatrix1 = r1->GetCovarianceMatrix();
TMatrixDSym covMatrix2 = r2->GetCovarianceMatrix();
//TMatrixDSym covMatrix3 = r3->GetCovarianceMatrix();
//Drawing the results ////////////////////////////////////
total->SetLineColor(kRed); total->Draw("same");
g1->SetLineColor(kGreen); g1->Draw("same");
g2->SetLineColor(kBlue);  g2->Draw("same");
//bg->SetLineColor(kBlack); bg->Draw("same");

////////////////////////////////////////////////////////////
//Fİnding the count
double PeakCoverageFactor = 2.326 ;//in terms of sigma. Page 106, table 5.11Practical Gamma Ray Spectroscopy, or k factor on page 115, why they are different ????????
double ConfidenceLimitCorrection= Err_f(PeakCoverageFactor); // Associated degree of 2 tailed confidence 
printf("For probabilty interval alpha: %1.2f\n",    (1.- ConfidenceLimitCorrection)/2.  );  //Table 5.3 Gilmore page 129 in pdf.
printf("For k_alfa factor %1.3f in statistics.\n",   PeakCoverageFactor);
printf("Err_f()=ConfidenceLimitCorrection= %1.4f in 2 tailed confidence limit.",  ConfidenceLimitCorrection);


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

TAxis *Xaxis = h1->GetXaxis();
double xVmin= p1 - (PeakCoverageFactor* p2);
double xVmax= p1 + (PeakCoverageFactor* p2);
double LL= Xaxis->FindBin(xVmin ); // min x value 
double UL= Xaxis->FindBin(xVmax); // max x value

//This integral initially does not cover the 100% of the area, only the one in the selected range
// ConfidenceLimitCorrection;  //virtual Double_t Integral (Double_t a, Double_t b, Double_t epsrel=1.e-12)
double GrossArea1 = g1->Integral( xVmin, xVmax) / HistoBinning ;	
double Error_GA1  = g1->IntegralError(xVmin, xVmax , g1Parameters, covMatrix1.GetMatrixArray()) / HistoBinning;
cout <<"1st peak's " << "Integral: "<< GrossArea1 << " +- "<< Error_GA1 << endl;

double GrossArea2= g2->Integral( xVmin, xVmax) / HistoBinning ;
double Error_GA2  = g2->IntegralError(xVmin, xVmax , g2Parameters, covMatrix2.GetMatrixArray()) / HistoBinning;
cout <<"2nd peak's " << "Integral: "<< GrossArea2 << " +- "<< Error_GA2 << endl;

//double GrossArea3 = bg->Integral( xVmin, xVmax) / HistoBinning ;	
//double Error_GA3 = bg->IntegralError(xVmin, xVmax , bgParameters, covMatrix3.GetMatrixArray()) / HistoBinning;
//cout <<"Background's " << "Integral: "<< GrossArea3 << " +- "<< Error_GA3 << endl;
}




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


int n = 1024; // number of bins
double LL = 0.;
double UL = 1024.; // make sure "bin width" = 1.
////////////////////////////////////////////////
TH1F *h1= CreateHisto( n,  LL,  UL);
////INPUT FILE for initialization //////////////
double slide=22. ; 
////////////////////////////////////////////////
FitGussianFunctions(h1, n, LL, UL);



return 0;
}