/// \author Rene Brun
 
#include "TCanvas.h"
#include "TMath.h"
#include "TH1.h"
#include "TF1.h"
#include "TRandom.h"
#include "TSpectrum.h"
#include "TVirtualFitter.h"
 
Int_t npeaks = 30;//.......arbitrary number?
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]; // "height" or "area"
      Double_t mean  = par[3*p+3];
      Double_t sigma = par[3*p+4];

      result += norm*TMath::Gaus(x[0],mean,sigma);
   }
   return result;
}

void gaussian_fit(Int_t np, double sig, double limit, TH1F* h, Double_t *mean, Double_t *sigma, Double_t *height, Int_t *numbpeaks) {
   npeaks = TMath::Abs(np);

   // Generate n peaks at random  Double_t *sigma, Double_t *height,
   Double_t par[3000];
   par[0] = 0;
   par[1] = 0;

   Int_t p;

   for (p=0;p<npeaks;p++) {
      par[3*p+2] = 0; // "height"
      par[3*p+3] = 10+gRandom->Rndm()*980; // "mean" how to know why 10 and 980
      par[3*p+4] = 3+2*gRandom->Rndm(); // "sigma"
   }

   TF1 *f = new TF1("f",fpeaks,0,1000,2+3*npeaks); // 2+3*npeaks ->number of free parameters used by the function
   f->SetNpx(1000);
   f->SetParameters(par);

   TCanvas *c1 = new TCanvas("c1","c1",10,10,1000,900);
   c1->Divide(1,2);
   c1->cd(1);
   h->Draw();
   TH1F *h2 = (TH1F*)h->Clone("h2");

   // Use TSpectrum to find the peak candidates
   TSpectrum *s = new TSpectrum(2*npeaks);
   Int_t nfound = s->Search(h,sig,"",limit);
   
   printf("Found %d candidate peaks to fit\n",nfound);

   // Estimate background using TSpectrum::Background
   TH1 *hb = s->Background(h,10,"same");
   hb->SetLineColor(kBlue);
   if (hb) c1->Update();
   if (np <0) return;
 
   //estimate linear background using a fitting method
   c1->cd(2);
   TF1 *fline = new TF1("fline","pol1",0,1000);
   h->Fit("fline","qn");

   // Loop on all found peaks. Eliminate peaks at the background level
   par[0] = fline->GetParameter(0);
   par[1] = fline->GetParameter(1);
   npeaks = 0;
   Double_t *xpeaks;
   xpeaks = s->GetPositionX();

   for (p=0;p<nfound;p++) 
   {
      Double_t xp = xpeaks[p];
      Int_t bin = h->GetXaxis()->FindBin(xp);
      Double_t yp = h->GetBinContent(bin);
      if (yp-TMath::Sqrt(yp) < fline->Eval(xp)) continue;
      par[3*npeaks+2] = yp; // "height"
      par[3*npeaks+3] = xp; // "mean"
      par[3*npeaks+4] = 3; // "sigma"

      npeaks++;
   }

   printf("Found %d useful peaks to fit\n",npeaks);
   TF1 *fit = new TF1("fit",fpeaks,0,1000,2+3*npeaks);
   TVirtualFitter::Fitter(h2,10+3*npeaks);
   fit->SetParameters(par);
   fit->SetNpx(1000);
   fit->SetLineColor(kBlue);

   h2->Fit("fit");
   //Double_t sigma[npeaks], height[npeaks];
   Int_t i=3, j=0;
   Double_t aux[npeaks], auxsigma[npeaks], auxheight[npeaks];

   while(i<3*npeaks+1)
   {
     mean[j]=fit->GetParameter(i);
     aux[j]=mean[j];
     sigma[j]=fit->GetParameter(i+1);
     auxsigma[j]=sigma[j];
     height[j]=fit->GetParameter(i-1);
     auxheight[j]=height[j];
     i=i+3;
     j=j+1; // j=(i-3)/3
   }

   sort(mean, mean + npeaks);
   int k=0;

   for (i = 0; i < npeaks; ++i)
   {
      for (int j = 0; j < npeaks; ++j)
      {
         if (mean[i]==aux[j])
              {
                 sigma[k]=auxsigma[j];
                 height[k]=auxheight[j];
                 k++;
              } 
      }
      // cout<<"mean "<<mean[i]<<endl;
      // cout<<"sigma "<<sigma[i]<<endl;
      // cout<<"height "<<height[i]<<endl;
   }

   //for (i = 0; i < npeaks; ++i) cout<<"Resolution peak "<< i << ": " << ((sigma[i]*2.355)/mean[i])*100<<endl; 
   // for (i = 0; i < npeaks; ++i) cout<<"mean peak "<< i << ": " << mean[i]<<endl; 
   // for (i = 0; i < npeaks; ++i) cout<<"sigma peak "<< i << ": " << sigma[i]<<endl; 
   // for (i = 0; i < npeaks; ++i) cout<<"height peak "<< i << ": " << height[i]<<endl; 

   numbpeaks[0]=npeaks;
}