/*
 * This macro is working perfectly as expected!
 * Ajay Y. Deo
 * Department of Physics
 * IIT Roorkee, India
 *
 * July 5, 2018
 *
 */

#include <TFile.h>
#include <TH1.h>
#include <TF1.h>
#include <TH2.h>
#include <TCanvas.h>
#include "TSpectrum.h"
#include "TVirtualFitter.h"
#include "TLatex.h"

#define CALIBRATION 1   // For 0.5 keV/ch
#define __PEAKS_C_FIT_AREAS__ 1 /* fit peaks' areas */

/*=====================================================*/
char ParNam1[10], ParNam2[10], ParNam3[10];

TH1 *hb;
Int_t npeaks = 30;

Double_t fpeaks(Double_t *x, Double_t *par) {
   Double_t xx = x[0];
   Int_t bin = hb->GetXaxis()->FindBin(xx);
   Double_t result = hb->GetBinContent(bin);
   for (Int_t p=0;p<npeaks;p++) {
      Double_t norm  = par[3*p+0]; // "height" or "area"
      Double_t mean  = par[3*p+1];
      Double_t sigma = par[3*p+2];
#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
pfit_bg(Char_t *hist_name, Int_t blo, Int_t bhi)
{
Int_t xlo, xhi;
Double_t par[300] = {0.0};
Double_t xpos[50] = {0.0}, ypos[50] = {0.0};
Int_t p;
TH1D *hist;
Double_t pPosition[30];
Double_t pHeight[30];

#if CALIBRATION == 1
  xlo = 2*blo; xhi = 2*bhi; // For 0.5 keV/ch dispersion
#else
  xlo = blo; xhi = bhi; // For 1 keV/ch dispersion
#endif

  hist = (TH1D *) gROOT->FindObject(hist_name);
  hist->GetXaxis()->SetRange(xlo, xhi);

  TCanvas *c1 = new TCanvas("c1","c1",10,10,1000,900);
  c1->Divide(1,2);
  c1->cd(1);
  hist->Draw();
  gPad->SetLogy();

   TH1F *h2 = (TH1F*)hist->Clone("h2");
   TSpectrum *s = new TSpectrum(2*npeaks);
   Int_t nfound = s->Search(hist,3,"",0.5);
//   Int_t nfound = s->Search(hist,3,"nodraw",0.05);
   printf("Found %d candidate peaks to fit\n",nfound);
//   hb = s->Background(hist,20,"BackIncreasingWindow,BackOrder2,BackSmoothing3,Compton,same");
//   hb = s->Background(hist,20,"BackOrder2,BackSmoothing3,Compton");
   hb = s->Background(hist,20,"Compton,same");

#if ROOT_VERSION_CODE >= ROOT_VERSION(6,00,00)
   Double_t *xfound; // ROOT 6
#else
   Float_t *xfound; // ROOT 5
#endif

  xfound = s->GetPositionX();
   for (p=0;p<nfound;p++) {
      pPosition[p] = xfound[p];
      Int_t xbin = hist->GetXaxis()->FindBin(pPosition[p]);
      pHeight[p] = hist->GetBinContent(xbin);

        TText *t2 = new TText(pPosition[p]-1,pHeight[p]*0.65,Form("%g",pPosition[p])); // For 1keV/ch
//        TText *t2 = new TText(pPosition[p],pHeight[p]*0.65,Form("%g",pPosition[p]*0.5));  // For 0.5 keV/ch
        t2->SetTextAngle(90);
        t2->SetTextColor(kBlue);
        t2->SetTextSize(0.03);
        t2->Draw();

      xpos[p] = pPosition[p];
      ypos[p] = pHeight[p];
   }

   if (hb) c1->Update();
   if (nfound <0) return;

   c1->cd(2);

   npeaks = 0;
#if ROOT_VERSION_CODE >= ROOT_VERSION(6,00,00)
   Double_t *xpeaks; // ROOT 6
#else
   Float_t *xpeaks; // ROOT 5
#endif

   xpeaks = s->GetPositionX();
   for (p=0;p<nfound;p++) {
      Double_t xp = xpeaks[p];
      Int_t bin = hist->GetXaxis()->FindBin(xp);
      Double_t yp = hist->GetBinContent(bin);
      par[3*npeaks+0] = yp; // "height"
      par[3*npeaks+1] = xp; // "mean"
      par[3*npeaks+2] = 3; // "sigma"
#if defined(__PEAKS_C_FIT_AREAS__)
      par[3*npeaks+0] *= par[3*npeaks+2] * (TMath::Sqrt(TMath::TwoPi())); // "area"
#endif /* defined(__PEAKS_C_FIT_AREAS__) */

//      xpos[p] = xp;
//      ypos[p] = yp;
      npeaks++;
   }

   printf("Found %d useful peaks to fit\n",npeaks);
   printf("Now fitting: Be patient\n");
//   TF1 *fit = new TF1("fit",fpeaks,0,1000,2+3*npeaks);
   TF1 *fit = new TF1("fit",fpeaks,xlo,xhi,3*npeaks);

   TVirtualFitter::Fitter(h2,10+3*npeaks);
   fit->SetParameters(par);
   fit->SetNpx(1000);
	for (p=0;p<nfound;p++) {
	   sprintf(ParNam1,"Area[%d]",p+1);
	   sprintf(ParNam2,"Mean[%d]",p+1);
	   sprintf(ParNam3,"Sigma[%d]",p+1);

	   fit->SetParName(3*p+0, ParNam1);
	   fit->SetParName(3*p+1, ParNam2);
	   fit->SetParName(3*p+2, ParNam3);

	   fit->SetParLimits(3*p+2, 1.0, 10.0);
	}
   h2->Fit("fit");

//   for (Int_t i = 0; i < fit->GetNpar(); i++) printf("%s = %g +- %g\n", fit->GetParName(i), fit->GetParameter(i), fit->GetParError(i));

   gPad->SetLogy();

	for(Int_t ip = 0; ip <nfound;ip++){
//	TText *t1 = new TText(xpos[ip]-4,ypos[ip]*0.65,Form("%g",xpos[ip]*0.5));
	TText *t1 = new TText(xpos[ip]-1,ypos[ip]*0.65,Form("%g",xpos[ip]));
	t1->SetTextAngle(90);
	t1->SetTextColor(kBlue);
	t1->SetTextSize(0.03);
	t1->Draw();
	}

}

/*---------------------------------------------------------------*/