//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>
/*
double Gauss_PDF(double x) {

	 TF1*f1=new TF1("myf", "[0]*exp(-0.5*((x-[1])/[2])**2)", -4, 4);
	     f1->SetParameter(0, (1./sqrt(2.*TMath::Pi())));
	     f1->SetParameter(1, 0);
         f1->SetParameter(2, 1);
         f1->Draw();
		 double y = f1->Eval(x);

	return y;
}
*/
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.)) ;
//double y= ROOT::Math::chisquared_cdf((x * 1.), 1.) ;

 //TF1*f2 =new TF1("myf2", "ROOT::Math::erf(x / TMath::Sqrt(2.))", -4, 4); //f2->SetParameter(1, 0);
 //f2->Draw();
////////////////////////////////////////////
/////////////////////
/*Try:
Wikipedia → Chi-squared distribution 1
Wolfram MathWorld → Chi-Squared Distribution
*/
return y;
}


void FittingOneGussianFunctions(TH1F *h1, double  LL, double  UL){

//cout <<"old " <<LL << "  " << UL << endl;
/////////////////////////////////////////////////////////////////////////////////
//For pol. background functions only
double extention =20. ; //For background  // Gilmore book page 112, figure 5.7

//cout<<"For 1st Gaussian Limits, LL_BCG= " <<LL_BCG<< " UL_BCG= " <<UL_BCG << endl;
////////////////////Defining Functions //////////////////////////////////////////////////////////////
//The limits are put to get exact function shapes, later they should be change to heir real values by reading the fit parameters.
TF1 *g1    = new TF1("Myg1","gaus", LL, UL);	
TF1 *pol1  = new TF1("Mypol1","pol1",   LL- extention, UL+extention) ;	//POLYNOMIAL BACKGROUND
TF1 *total = new TF1("total","gaus(0)+pol1(3)",  LL- extention, UL+extention) ;//TOTAL FUNCTION
/////////////////////////////////////////////////////////////////////////////////////////////////////////
//TH1F *h = (TH1F*)h1->Clone("h");
h1->SetLineColor(kBlack);
g1->SetLineColor(kBlue);
pol1->SetLineColor(kBlack);
total->SetLineColor(kRed);
h1->Fit(g1,"R");
h1->Fit(pol1,"R+"); 
///////////////////////////////////////////////////////////
// To store the parameter values in to an ARRAY called "par[5]"
double par[5];
g1->GetParameters(&par[0]);
pol1->GetParameters(&par[3]);
total->SetParameters(par);  //first get the parameter from g1 and pol1 , then set them for total function as parameters.
h1->Fit(total, "RNQ+");   ///RBE  //Now, I am able to fit and draw this function since I assigned its values which were taken from the first gaussian and pol1 fits
//double chi2 = fitting->GetChisquare();
// value of the first parameter
double p0 =  par[0];  //height of gaussian
double p1 =  par[1]; //mean value of gaussian
double p2 =  par[2]; //if this is sigma, width of gaussian=8*sigma
cout<<"mean value for this peak= "<< p1 <<endl;
double p3 = par[3];
double p4 = par[4];
 //https://root.cern.ch/doc/master/classTF1.html#a4f798c9626dc2d2268198ee01dc51437     
// error of the first parameter
double e0 = g1->GetParError(0);
double e1 = g1->GetParError(1);
double e2 = g1->GetParError(2);
double e3 = pol1->GetParError(0);
double e4 = pol1->GetParError(1);
//////////////////////////////////
/*
//g1->SetParNames("height","Mean","Sigma");
//g1->SetParameters(100, 15, 2);
g1->SetParameter(0, p0);
g1->SetParameter(1, p1);
g1->SetParameter(2, p2);
pol1->SetParameter(0, p3);
pol1->SetParameter(1, p4);
*/
//For interested gaussian peak ONLY
double* gaussParameters = new double[3] ;
gaussParameters[0] = p0;
gaussParameters[1] = p1;
gaussParameters[2] = p2;
double* MyConvarienceMatrix =new double[9];
MyConvarienceMatrix[0]=e0 ;
MyConvarienceMatrix[1]=0  ;
MyConvarienceMatrix[2]=0  ;
MyConvarienceMatrix[3]=0  ;
MyConvarienceMatrix[4]=e1 ;
MyConvarienceMatrix[5]=0  ;
MyConvarienceMatrix[6]=0  ;
MyConvarienceMatrix[7]=0  ;
MyConvarienceMatrix[8]=e2 ;	

//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);
cout<<"Err_f()=ConfidenceLimitCorrection= " << ConfidenceLimitCorrection <<endl;
/*
     if (PeakCoverageFactor ==3.){ConfidenceLimitCorrection= 0.9973;   }
else if (PeakCoverageFactor ==2.326){ConfidenceLimitCorrection= 0.9799;   }
else if (PeakCoverageFactor ==2.){ConfidenceLimitCorrection= 0.9545;   }
else if (PeakCoverageFactor ==1.645){ConfidenceLimitCorrection= 0.90;  }
else if (PeakCoverageFactor ==1.){ConfidenceLimitCorrection= 0.6828;   }
*/
//TF1 *myf = new TF1("myf", "ROOT::Math::erf(x)", -5., 5.);
//gPad->SetGrid(1, 1);
//myf->Draw();




//updated real values for peaks ,all calculations have to be done with integer bin numbers.
//UPDATED after fitting
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
//cout <<"For x1:" << p1 - (PeakCoverageFactor* p2) <<" to x2: " << p1 + (PeakCoverageFactor* p2) <<   " , updated bins of LL and UL: " <<LL << "  " << UL << endl;

//From the information below, you can implement COVELL, TPA, ... methods in to the code 
//cout <<"!!!!!!!!!!!!!! " << h1->GetXaxis()->FindBin(p1 - (PeakCoverageFactor* p2)) << " to " <<  h1->GetXaxis()->FindBin(p1 + (PeakCoverageFactor* p2)) << endl;
//cout <<"For integral correction from a histogram by using directly the bin content purest range from " << newLL<< " to " << newUL << " for " << PeakCoverageFactor << " sigma" <<endl;
//SET gaussian ranges for drawing
    /*
virtual void 	SetRange (Int_t first=0, Int_t last=0)
 	Set the viewing range for the axis using bin numbers. More...
 
virtual void 	SetRangeUser (Double_t ufirst, Double_t ulast)
 	Set the viewing range for the axis from ufirst to ulast (in user coordinates, that is, the "natural" axis coordinates).
	 */
//hist->GetXaxis()->SetRange(ceil(p1 - (4.* p2)) ,ceil(p1 + (4.* p2))  ); 
g1->SetRange(xVmin , xVmax)  ; g1->Draw("same");  g1->Print("Q");
//pol1->SetRange(xVmin-extention , xVmax+extention)  ;  pol1->Draw("same");  pol1->Print("Q"); 
//total->SetRange(xVmin-extention , xVmax+extention)  ;   //total->Draw("same"); // total->Print("Q"); 
//Q: Quiet mode (minimum printing) &  V: Verbose (detailed) mode 
/*
std::cout <<"!!!!!!!!TESTING2!!!!!:  @channel " <<g1->GetParameter(1) <<", the count is:  "<< g1->Eval(g1->GetParameter(1)) << 
" which equals to " << g1->GetParameter(0) <<" or " << g1->Eval(432.431) << endl;
cout << h1->FindBin(432.431) << endl;

double yunkown= -1000, xunkown= -1000;
for (int i=315; i < 360 ; i++) {
double vy1= h1->GetBinContent(i);
double vy2= h1->GetBinContent(i+1);
if (vy2 > vy1){ cout <<"!!!!!!!!!!!!!     "<<  i << "   "<< vy1 << "   " << vy2  << endl; 
 yunkown= vy1;  xunkown= i;  cout << "xunkown:  " << xunkown << " yunknown " << yunkown << endl;  break;}
}
*/
double yunkown= -1000, xunkown= -1000;
for (int i= 420; i < 465; i++) {
double vy1= g1->Eval(i);
double vy2= g1->Eval(i+1);
if (vy2 > vy1){ cout <<"!!!!!!!!!!!!!     "<<  i << "   "<< vy1 << "   " << vy2  << endl; 
 yunkown= vy1;  xunkown= i;  
 cout << "xunkown:  " << xunkown << " yunknown " << yunkown << endl;  break;}
}

int NofBins = UL-LL+1;
TAxis *Xaxis = h1->GetXaxis();
int binwidth= Xaxis->GetBinWidth(LL);
double TotalCounts;
double TotalError, TotalErrorUpdated;
cout <<"  Bin width:   "<<  binwidth  << endl;
cout <<"NbinsX: " << h1->GetNbinsX() << endl;
	for (int i = LL; i < UL+1 ; i++ ){
    
	//cout<< i << "  " << h1->GetBinContent(i) << endl;
    TotalCounts += h1->GetBinContent(i); ////This integral initially does not cover the 100% of the area, only the one in the selected range
	//TotalError +=h1->GetBinError(i);
	TotalError += std::pow(h1->GetBinError(i), 2.);  // GetBinError(i) = sigma value of the count on that bin.
	}
	TotalErrorUpdated= sqrt(TotalError);

cout <<"From " <<LL << " to  " << UL <<"   Total GetBinContent: " << TotalCounts <<" +- " << TotalErrorUpdated <<  endl;

double myerror; 
//This integral initially does not cover the 100% of the area, only the one in the selected range
double myintegral = h1->IntegralAndError(LL, UL, myerror, ""); // "" ... or ... "width"
//TH1F * h1= (TH1F*)gDirectory->FindObjectAny("h1")
cout << "IntegralAndError(...) = " << myintegral << " +- " << myerror << endl;

//for ( int bin=0 ; bin< h1->GetNbinsX() ; bin++ ){
//printf("!!!!!????!!!!%d\t%.1f\t%.3e\t%.3e\n"  ,bin, h1->GetBinCenter(bin) ,h1->GetBinContent(bin), h1->GetBinError(bin) );   }

/////////////////////////////////////////////////////////////ıntegral from fit parameters ///////////////////////
/////Integrall correction for Fitted function is meaningless, so just use the widest rane to cover 100%
   double HistoBinning=1.;
   int binx1= Xaxis->FindBin(LL);
   double LowEdge1=Xaxis->GetBinLowEdge(binx1);
   double Center1 =Xaxis->GetBinCenter(binx1) ;
   double UpEdge1=Xaxis->GetBinUpEdge(binx1) ;
   //////////////////////////////
   int binx2= Xaxis->FindBin(UL);
   double LowEdge2=Xaxis->GetBinLowEdge(binx2);
   double Center2 =Xaxis->GetBinCenter(binx2) ;
   double UpEdge2=Xaxis->GetBinUpEdge(binx2) ;

//testing function's integral in ROOT fitting
//This integral initially does not cover the 100% of the area, only the one in the selected range
double GrossArea1 = g1->Integral( xVmin, xVmax) ;// ConfidenceLimitCorrection;  //virtual Double_t Integral (Double_t a, Double_t b, Double_t epsrel=1.e-12)	
double Error_GA1  = g1-> IntegralError(xVmin, xVmax , gaussParameters, MyConvarienceMatrix) ; //   / ConfidenceLimitCorrection ; // HistoBinning 
cout <<"1st peak's " << "Integral: "<< GrossArea1 << " +- "<< Error_GA1 << endl;
cout << "check bin numbers: "  <<  LowEdge1 << "  " <<Center1 << " " << UpEdge1 << "  " <<LowEdge2 << "  " <<Center2 << " " << UpEdge2 << endl;
cout <<"For x1:" << xVmin <<" to x2: " << xVmax <<   " , updated bins of LL and UL: " <<LL << "  " << UL << endl;
///////////////////////////////////////////////////////////////////////////////////////
//TLine(Double_t x1,Double_t y1,Double_t x2,Double_t y2) //x1, y1, x2, y2 are the coordinates of the first and the second point.
TLine* threshold_min = new TLine (p1-PeakCoverageFactor*p2,0,  p1-PeakCoverageFactor*p2,h1->GetMaximum());
TLine* threshold_max = new TLine (p1+PeakCoverageFactor*p2,0,  p1+PeakCoverageFactor*p2,h1->GetMaximum());
threshold_min->SetLineColor(kRed);  
threshold_max->SetLineColor(kRed);
threshold_min->Draw();              
threshold_max->Draw();
////////////////////////////////
   int binx1_BCG= Xaxis->FindBin(p1-PeakCoverageFactor*p2);
   double LowEdge1_BCG=Xaxis->GetBinLowEdge(binx1_BCG);
   double Center1_BCG =Xaxis->GetBinCenter(binx1_BCG) ;
   double UpEdge1_BCG=Xaxis->GetBinUpEdge(binx1_BCG) ;
   //////////////////////////////
   int binx2_BCG= Xaxis->FindBin(p1+PeakCoverageFactor*p2);
   double LowEdge2_BCG=Xaxis->GetBinLowEdge(binx2_BCG);
   double Center2_BCG =Xaxis->GetBinCenter(binx2_BCG) ;
   double UpEdge2_BCG=Xaxis->GetBinUpEdge(binx2_BCG) ;
TLine* BCG = new TLine (Center1_BCG, h1->GetBinContent(binx1_BCG), Center2_BCG, h1->GetBinContent(binx2_BCG));
BCG->SetLineColor(kBlack);
BCG->Draw();  
///////////////////////////////////
TText* text_min= new TText(p1-PeakCoverageFactor*p2, g1->Eval(p1)*2.2, Form("%.1f",p1-PeakCoverageFactor*p2 )) ;
TText* text_max= new TText(p1+PeakCoverageFactor*p2, g1->Eval(p1)*2.2 , Form("%.1f", p1+PeakCoverageFactor*p2 ));
text_min->SetTextSize(0.04) ; 
text_max->SetTextSize(0.04) ; 
text_min->Draw();  
text_max->Draw();

//Total background area from the fit has to recover from the ConfidenceLimitCorrection
double BCG_Area=  ( h1->GetBinContent(binx1_BCG)  +  h1->GetBinContent(binx2_BCG) )* 0.5 * (Center2_BCG - Center1_BCG) / ConfidenceLimitCorrection;
double BCG_Area_Error=  sqrt( pow(h1->GetBinError(binx1_BCG), 2.)  +  pow(h1->GetBinError(binx2_BCG), 2.) )* 0.5 * ((Center2_BCG - Center1_BCG) / ConfidenceLimitCorrection ) ;

 cout << "BCG_Area: "  << BCG_Area << " +- "  <<  BCG_Area_Error <<  endl;

 double NetArea= GrossArea1 - BCG_Area ;
 double NetAreaError= sqrt( pow(Error_GA1,2.)+ pow(BCG_Area_Error,2.) );
 cout << "Net_Area: "  << NetArea << " +- "  <<  NetAreaError <<  endl;

//double GrossArea = (  g1->Integral(newLL, newUL) ) / ConfidenceLimitCorrection; 
/*Double_t TH1::Integral(Int_t binx1,Int_t binx2,Option_t *option = "") const;
By default the integral is computed as the sum of bin contents in the range. if option "width" is specified, 
the integral is the sum of the bin contents multiplied by the bin width in x.
Double_t TH1::IntegralAndError	(	Int_t 	binx1,
Int_t 	binx2,
Double_t & 	error,
Option_t * 	option = "" 
)	
Return integral of bin contents in range [binx1,binx2] and its error.

By default the integral is computed as the sum of bin contents in the range. if option "width" is specified, 
the integral is the sum of the bin contents multiplied by the bin width in x.
 the error is computed using error propagation from the bin errors assuming that all the bins are uncorrelated
*/	
 //TF1 , Integral (Double_t a, Double_t b, Double_t epsrel=1.e-12)
/*
virtual Double_t 	Integral (Double_t a, Double_t b, Double_t epsrel=1.e-12)
 	IntegralOneDim or analytical integral. More...
 
virtual Double_t 	IntegralError (Double_t a, Double_t b, const Double_t *params=0, const Double_t *covmat=0, Double_t epsilon=1.E-2)
 	Return Error on Integral of a parametric function between a and b due to the parameter uncertainties and their covariance matrix from the fit.
*/
//https://root.cern.ch/doc/master/classTString.html#aeb44d6183d8b1f1b7090dbd3c93f5e39
//example: Form("%3.2f^{#circ}",GettingProtonAngleInLabFrame->Eval(0))
//////////////Scanning the histogram bin by bin until finding where on X axis I have 5% ratio interms of integral  /////////////////////////////////////////////////////////////////
/*
double xt= -1000; // the threshold
double step = h1->GetBinWidth(2); // bin width of the bin number 2

for(unsigned int i = 0 ; i< 60 ; i++ ){
	double x = UL1 + i* step;
//cout <<"hey:   " <<g1->Eval(300) << "    " <<g1->Eval(330) <<endl;
		if(g1->Eval(x+1) > g1->Eval(x)){
			xt = x;
			break;
		}	
}
cout << "xt: " << xt << endl;
*/
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////


}

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


int FindFittingValues(){

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

TH1F *hist = new TH1F("hist", "histogram",1024.,0.,1024.); // 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();
hist->GetXaxis()->SetTitle("Channel");
hist->GetYaxis()->SetTitle("Counts");
hist->GetXaxis()->SetRange(380,540); // for user to focus
hist->Draw("hist");

/////////////////////////////INPUT FILE for initialization /////////////////////////////////////////////////
double slide=22. ; 
//for peak 1
double Centroid1 =432.; // 313.; //peak 1
double LL1= Centroid1 - slide;
double UL1= Centroid1 + slide; 
//for peak 2
double Centroid2 = 489.; // 354. ; //peak 8
double LL2= Centroid2 - slide;
double UL2= Centroid2 + slide; 
////////////////////////////////////////////////////////////////////////////////////
FittingOneGussianFunctions(hist, LL1, UL1);
FittingOneGussianFunctions(hist, LL2, UL2);

c1->Update();
/*
TCanvas *c2 = new TCanvas("c2", "c2", 800,400);
c2->cd();
gStyle->SetPalette(1);
gStyle->SetOptStat("inM");
double myY1= 2.*Gauss_PDF(1.);
cout<<"Gauss_PDF!!!!: " << myY1 <<endl;
c2->Update();


TCanvas *c3 = new TCanvas("c3", "c3", 800,400);
c3->cd();
gStyle->SetPalette(1);
gStyle->SetOptStat("inM");
double myY2= Err_f(1.);
cout<<"Err_f!!!!: " << myY2 <<endl;

c3->Update();
*/



return 0;
}