//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 "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"




TH1F *Recall_Create_Histogram(string address1){

ifstream in;
in.open(address1);
float Count;
TH1F*h =new TH1F("h", "h" , 2048, 0, 2048);  //TH1F*h =new TH1F(Form("h%d", i), Form("h%d", i) ,1024, 0, 1024);
Int_t nlines= 0;
TFile*f1= new TFile("/home/ic00025/Desktop/AliHoca/EnamicData.root","RECREATE");

while (1) {  //!in.eof
in  >> Count;
if (!in.good()) break;
h->Fill( nlines,Count); 
nlines++;
}

printf("found %dpoints\n",nlines);
in.close();
f1->Write();
//h->Draw("hist"); //"Count:Channel","", ""
return h;
}

void FittingOneGussianFunctions(TH1F *h1, int  LL1, int  UL1, int  extention1_right, int  extention1_left){


/////////////////////////////////////////////////////////////////////////////////
//For pol. background functions only
int LL_BCG=LL1-extention1_left, UL_BCG= UL1+extention1_right; //For background

cout<<"For 1st Gaussian Limits, LL_BCG= " <<LL_BCG<< " UL_BCG= " <<UL_BCG << endl;
////////////////////Defining Functions //////////////////////////////////////////////////////////////
//The limits are put t get exact function shapes, so they might be narrow.
TF1 *g1    = new TF1("g1","gaus", LL1, UL1 );	
TF1 *pol1  = new TF1("pol1","pol1",  LL_BCG, UL_BCG );	//POLYNOMIAL BACKGROUND
TF1 *total = new TF1("total","gaus(0)+pol1(3)",  LL1 , UL1 );//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);

///////////////////////////////////////////////
//SET gaussian ranges for drawing
g1->SetRange(LL_BCG, UL_BCG); 
pol1->SetRange(LL_BCG, UL_BCG); 
total->SetRange(LL_BCG, UL_BCG); 
////////////////////////////////////////
g1->Draw("same");
pol1->Draw("same");  
total->Draw("same");

//Fİnding the count
double PeakCoverageFactor =3.; //in terms of sigma. Page 106, table 5.11Practical Gamma Ray Spectroscopy, or k factor on page 115, why they are different ????????
double PeakCountCorrection;
     if (PeakCoverageFactor ==3.){PeakCountCorrection= 0.999;   }
else if (PeakCoverageFactor ==2.576){PeakCountCorrection= 0.990; } //bckg count yanlış çıkar, UL ve LL limitlerden alınmazsa sayımlar.
else if (PeakCoverageFactor ==2.){PeakCountCorrection= 0.955;   }
else if (PeakCoverageFactor ==1.645){PeakCountCorrection= 0.90;  }
else if (PeakCoverageFactor ==1.){PeakCountCorrection= 0.683;   }

//for peaks
int newLL= ceil(p1 - (PeakCoverageFactor* p2) );
int newUL= ceil(p1 + (PeakCoverageFactor* p2) );
//for background
int newLL_BG= newLL+2;
int newUL_BG= newUL-2;
//////////////////////////////
int NofBins = newUL-newLL+1;
TAxis *axis = h1->GetXaxis();
int binwidth= axis->GetBinWidth(newLL);
cout <<"newLL:  "<< newLL << "	newUL:  "<< newUL <<"  Bin width:   "<<binwidth  << endl;

//fonksiyon fit ettirmeye gerek olmaksızın olan yöntem. Sadece histogram ile işlemler yapılıyor.
double GrossArea= (  h1->Integral(newLL, newUL) ) / PeakCountCorrection;  //Double_t TH1::Integral(Int_t binx1,Int_t binx2,Option_t *option = "") 	
double GrossArea_adjusted = (  h1->Integral(newLL+3, newUL-3) ) ;
///////////////////////////////////////////////////////////
double BCG_LowerSide= h1->Integral(newLL,  newLL_BG ) ;
double BCG_UpperSide= h1->Integral(newUL_BG,newUL ) ;
double BCG_Area = (NofBins)* ((BCG_LowerSide+BCG_UpperSide ) /6. );
double BCG_Area_adjusted= (BCG_Area/ NofBins)* (NofBins-6 );
////////////////////////////////////////////////////////////
double Net_Count = GrossArea-BCG_Area ;
double Net_Count_adjusted = GrossArea_adjusted-BCG_Area_adjusted ;
double AreaControl= total->Integral(newLL_BG, newUL_BG);
cout<< "AreaControl: function  " << AreaControl << "  vs histo "<<Net_Count << "  Net_Count_adjusted :  " << Net_Count_adjusted   <<endl;
cout << "Not:REMEMBER: Bin no 10 is between 9 and 10 on X axis if you compare ROOT integral results and your calculations here."<<endl;
cout<< "Check on X axis "<< newLL <<" ve " << newUL << " for Gross Area:\t"<< GrossArea <<"\tBCG:\t"<<  BCG_Area <<  "\tNet integral:\t" << Net_Count   << endl;

}




int FindFittingValues(){


TCanvas *c1 = new TCanvas("c1", "c1",546,179,2369,1449);
c1->cd();

gStyle->SetPalette(1);
gStyle->SetOptStat("inM");
string address2="EnamicData";
string address1="/home/ic00025/Desktop/AliHoca/"+address2+".txt" ;
Recall_Create_Histogram(address1);
//TH1F* h1 = (TH1F*)gDirectory->FindObjectAny("h1");
//TH1F *h1 = (TH1F*)h2->Clone("h1");
//TFile *file1= TFile::Open("/home/ic00025/Desktop/HarranProgramFile/Co60Source.root");
TFile *file1 = new TFile(Form("/home/ic00025/Desktop/AliHoca/%s.root",address2.c_str()),"NEW");
if ( file1->IsOpen() ) printf("File opened successfully\n");
TH1F*h1= (TH1F*)gDirectory->FindObjectAny("h");
h1->Draw("hist");
file1->Close();
delete file1;
/////////////////////////////INPUT FILE /////////////////////////////////////////////////
double p2,p5; //To recall these values, first define
int Centroid1 = 1502; //peak 7
//int Centroid2 = 1919; //peak 8
int slide= 5 ; //sadece çizim ve fit için. slide indeed should be 3 sigma for %99.7 data.
//Covell method says, take equal number of bins in each side of the centroid of each peak. Then calculate the ratio of this integral over whole, for the compensation
//slide= 3 sigma olmalı
//To find correct height, mean and sigma from Gaussians
int LL1= Centroid1 - slide,  UL1= Centroid1 + slide; //it covers almost %99.7 of the peak range
//Only to draw the gaussian and pol. functions in the whole range. Nothing else.
int  extention1_left = 20  , extention1_right = 20;  //for peak 1
//page 110, practical gamma ray spectroscopy: We can never know which counts within the peak region are due to background and which are the peak counts.
/////////////CONTROLLING /////////////////////////////////////
 if (extention1_right <=1 || extention1_left <=1){cout<<"extention must be greater than 1. !!!!!!!!!!!!!!!!!!!!!!!!!"<<endl;} 
FittingOneGussianFunctions(h1,  LL1, UL1, extention1_right,  extention1_left );

//EXTRA INFORMATION  /////////////////////////////////////////
//FittingOneGussianFunctions(h1,  LL2, UL2, extention2_right,  extention2_left, p2);
//TF1*g1= (TF1*)gDirectory->FindObjectAny("g1");
//TF1*g2= (TF1*)gDirectory->FindObjectAny("g2");






//TH1F *h1 = (TH1F*)h2->Clone("h1");

//how to calculate the integral with functions // same as calculating from histogram
//Important to make this convertion when you got energy on x axis instead of cahnnel number.
//COVELL METHOD - Gilmore sayfa 110
//To see examples, type on a terminal:  grep -r "GetBin(" ${ROOTSYS}/ , star koy
//int bin = h1->GetXaxis()->FindBin(x);
//For investigation
//h1->GetXaxis()->SetRangeUser(LL-1,UL); // Between 75-76 on X axis, bin no: 76, Between 85-86 on X axis, bin no 86
c1->Update();
return 0;
}