#include "TCanvas.h"
#include "TH1.h"
#include "TF1.h"
#include "TRandom.h"
#include "TSpectrum.h"
#include "TVirtualFitter.h"
#include "Riostream.h"

Int_t npeaks = 15;//number of peaks to fit


//0.------Definition of the function to be fitted to the histogram
Double_t fpeaks(Double_t *x, Double_t *par) {
  Double_t result = par[0] + par[1]*x[0];
  for (Int_t p=0;p<npeaks;p++) {
    Double_t norm  = par[3*p+2];
    Double_t mean  = par[3*p+3];
    Double_t sigma = par[3*p+4];
#if defined(__PEAKS_C_FIT_AREAS__)
      norm /= sigma * (TMath::Sqrt(TMath::TwoPi())); // "area"
#endif /* defined(__PEAKS_C_FIT_AREAS__) */
    result += norm*TMath::Gaus(x[0],mean,sigma);
  }
  return result;
}

void peaks(Int_t np=15){
  
 //1. ------- Path do Data file for the histogram
  
  // Set path to root file
	std::string myPath = "../Analysis/";

	// Set name of root file
	std::string myFileName = "Out_run_1.root";

	// Open the file containing the tree
	TFile *myFile = TFile::Open(std::string(myPath+myFileName).c_str());

	// Create a TTreeReader for the tree by passing the TTree's name and the TFile it is in
	TTreeReader myReader("Secondaries", myFile);
	
	//The branch "ParticleID" contains a vector<int>; access them as vID
   TTreeReaderValue<vector<int>> vPDGcode (myReader, "ParticleID"); //New ADD
   
   // A histogram for Particle ID ranging from 1 to 100  with 200 bins
	TH1F* hPDGcode = new TH1F("hPDGcode", ";Number;#entries", 2500, -30, 30); //New ADD
	

	// The branch "ParticleEnergy" contains a vector<double>; access them as vEnergies
	TTreeReaderValue<vector<double>> vEnergies(myReader, "ParticleEnergy");
	
	// Create histograms for the values we read

	// A histogram for the energies ranging from 0 MeV to 20 MeV with 200 bins
	TH1F* hEnergyGamma = new TH1F("hEnergyGamma", ";E in MeV;#entries", 600, 0,10);
	
	TH1F* hEnergyAll = new TH1F("hEnergyAll", ";E in MeV;#entries", 600, 0, 10);
	
	 while (myReader.Next()) {
    int energySize = vEnergies->size();
		
		for(int i = 0; i<energySize; ++i){
			
			double myEnergy= vEnergies->at(i);
			int myPDG= vPDGcode->at(i);
		
			hEnergyAll->Fill(myEnergy);
			if (myPDG==22){
				hEnergyGamma->Fill(myEnergy);
				}
				
				
   }
}   
  //2.--------- Finding peaks to fit in the spectrum 

  TSpectrum *s = new TSpectrum(npeaks);// TSpectrum class finds candidates peaks
  Int_t nfound = s->Search(hEnergyGamma,2,"new",0.10);
  printf("Found %d candidate peaks to fit\n",nfound);
  
  //3.------- Estimate background using TSpectrum::Background
   //TH1 *hb = s->Background(hEnergyGamma,20,"same");
   //if (hb) c1->Update();
  // if (np <0) return;
   
  //4.--------- Background Estimation
  
   TF1 *fline = new TF1("fline","pol1",0,10);// linear function (pol1) defined between [0,1000]
   hEnergyGamma->Fit("fline","qn");
   
   //5.-------- Setting initial parameters

   Double_t par[100];//Array of fitting parameters 
   // Loop on the found peaks for  setting initial parameters. Peaks at the background level will be disregarded from the fitting 
   par[0] = fline->GetParameter(0);
   par[1] = fline->GetParameter(1);
   npeaks = 0;
   Double_t *xpeaks;
   xpeaks = s->GetPositionX();//array with X-positions of the centroids found by TSpectrum 
   for (npeaks=0;npeaks<nfound;npeaks++) {
     Double_t xp = xpeaks[npeaks];
     Int_t bin = hEnergyGamma->GetXaxis()->FindBin(xp);
     Double_t yp = hEnergyGamma->GetBinContent(bin);
     //Condition for valid peaks: height-sqrt(height) has to be greater than the bckgnd:  
     if (yp-TMath::Sqrt(yp) < fline->Eval(xp)) continue;
     par[3*npeaks+2] = yp; //height
     par[3*npeaks+3] = xp; //centroid
     par[3*npeaks+4] = 3; //sigma -I've choosed it as same for all peaks 50
     
#if defined(__PEAKS_C_FIT_AREAS__)
      par[3*npeaks+2] *= par[3*npeaks+4] * (TMath::Sqrt(TMath::TwoPi())); // "area"
#endif /* defined(__PEAKS_C_FIT_AREAS__) */
     npeaks++;
   }
   
   printf("Found %d useful peaks to fit\n",npeaks);
   printf("Parameter No 1,2 ---> linear background Fit\n");
   printf("Parameters No 3+3*n ---> Gauss Constants\n");
   printf("Parameters No 4+3*n ---> Gauss Means\n");
   printf("Parameters No 5+3*n ---> Gauss Sigmas\n");
   printf("Now fitting: Be patient\n");
   
   //5.----------Fitting peaks
   
   TF1 *f = new TF1("f",fpeaks,0,1000,2+3*npeaks);
   f->SetParameters(par);
   f->SetNpx(1000);
   


   TF1 *fit = new TF1("fit",fpeaks,0,10,3*nfound+2);
   
   TVirtualFitter::Fitter(hEnergyGamma,npeaks); 
   fit->SetParameters(par);
   hEnergyGamma->Fit("fit");             
   
   TCanvas * c1 = new TCanvas("c1");
   hEnergyGamma->Draw();    
   fit->Draw("same");
   
}