#include <stdio.h>
#include <math.h>

#include <TFile.h>
#include <TMapFile.h>
#include <TH1.h>
#include <TF1.h>
#include <TH2.h>
#include <TCanvas.h>
#include <TSystem.h>
#include <TCutG.h>
#include <TKey.h>
#include <values.h>
#include "TSpectrum.h"
#include "TMath.h"
#include "TVirtualFitter.h"
#include "TAxis.h"
#include "TGaxis.h"
#include "TLatex.h"

/*=====================================================*/
Int_t npeaks = 30;
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];
#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(Char_t *hist_name, Int_t xlo, Int_t xhi)
{

Double_t par[3000];
Int_t p;
TH1D *hist;

  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();

   TH1F *h2 = (TH1F*)hist->Clone("h2");
   TSpectrum *s = new TSpectrum(2*npeaks);
   Int_t nfound = s->Search(hist,2,"",0.10);
   printf("Found %d candidate peaks to fit\n",nfound);
   TH1 *hb = s->Background(hist,20,"same");
   if (hb) c1->Update();
   if (nfound <0) return;

   c1->cd(2);

   TF1 *fline = new TF1("fline","pol1",0,1000);
   hist->Fit("fline","qn");

   par[0] = fline->GetParameter(0);
   par[1] = fline->GetParameter(1);

   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);

	TText *t1 = new TText(xp,yp,Form("%g",xp));
	t1->SetTextAngle(90);
	t1->Draw();
      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"
#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("Now fitting: Be patient\n");
   TF1 *fit = new TF1("fit",fpeaks,0,1000,2+3*npeaks);

   TVirtualFitter::Fitter(h2,10+3*npeaks);
   fit->SetParameters(par);
   fit->SetNpx(1000);
   h2->Fit("fit");


	TText *t = new TText(600,3000,"Ajay");
	t->SetTextAngle(90);
	t->Draw();

}

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