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

Int_t npeaks = 8;//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];
    result += norm*TMath::Gaus(x[0],mean,sigma);
  }
  //fpeaks = gaus_1+gaus_2+...+gaus_npeaks + linear_func.
  return result;
}

void peaks(){
  
  //1.---------Reading data file. 
  ifstream in;
  in.open(Form("./data.txt"));
  Float_t x,y;
  TH1F *h = new TH1F("h1","energy spectrum",4096,0,4096);
  for (Long_t i=0; i<4096;i++){
      in >> x >> y;
      h1->Fill(x,y);
  }
  in.close();
  
  //2.--------- Finding peaks to fit in the spectrum 

  TSpectrum *s = new TSpectrum(npeaks);// TSpectrum class finds candidates peaks
  Int_t nfound = s->Search(h ,0 ,"new");
  printf("Found %d candidate peaks to fit\n",nfound);
  
  //3.--------- Background Estimation
  
  TF1 *fline = new TF1("fline","pol1",0,650);// linear function (pol1) defined between [0,650]
   h->Fit("fline","qn");
   
   //4.-------- 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;
   Float_t *xpeaks = s->GetPositionX();//array with X-positions of the centroids found by TSpectrum 
   for (p=0;p<nfound;p++) {
     Float_t xp = xpeaks[p];
     Int_t bin = h->GetXaxis()->FindBin(xp);
     Float_t yp = h->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] = 50; //sigma -I've choosed it as same for all peaks-
     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->SetNpx(1000);
   f->SetParameters(par);


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