#include "TCanvas.h"
#include "TH1.h"
#include "TF1.h"
#include "TRandom.h"
#include "TSpectrum.h"
#include "TVirtualFitter.h"
   
Int_t npeaks = 30;

Double_t fpeaks(Double_t *x, Double_t *par) {
//   Double_t result = par[0] + par[1]*x[0];
  Double_t result =0.0;
  
   for (Int_t p=0;p<npeaks;p++) {
      Double_t norm  = par[3*p];
      Double_t mean  = par[3*p+1];
      Double_t sigma = par[3*p+2];
      result += norm*TMath::Gaus(x[0],mean,sigma);
   }
   return result;
}


void d1_fit(Int_t np=10) {

   npeaks = np;
   
   Double_t par[3000];
   
   TFile *hfile = gROOT->FindObject("test.root"); if (hfile) hfile->Close();
   hfile = new TFile("test.root","READ","Test Calibration");
   
   TCanvas *c1 = new TCanvas("c1","c1",10,10,700,700);
   c1->Divide(1,2);
   c1->cd(1);
  
   h1->GetXaxis()->SetRange(2,500);
   h1->Draw(); // TH1F *h = new TH1F("h","test",1000,0,2000);
   TH1F *h2 = (TH1F*)h1->Clone("h2");
 
   TSpectrum *s = new TSpectrum(npeaks);
   Int_t nfound = s->Search(h2,1,"");
   
   printf("Found %d candidate peaks to fit\n",nfound);
   npeaks = nfound;
   
   c1->Update();
   c1->cd(2); 
   
   Float_t *xpeaks = s->GetPositionX();
   Float_t *ypeaks = s->GetPositionY();
   
   for (int p = 0; p < nfound; p++) {
      par[3*p]   =  xpeaks[p];
      par[3*p+1] =  ypeaks[p];
      par[3*p+2] =  1.0;
   //   cout <<  par[3*p] << " " << par[3*p+1] <<" "<< par[3*p+2] <<endl;
   }
   
   printf("Now fitting: Be patient\n");
   
   TF1 *fit = new TF1("fit",fpeaks,2,500,3*npeaks);
   TVirtualFitter::Fitter(h2,3*npeaks); 
   fit->SetParameters(par);
   fit->SetNpx(100);
   h2->Fit("fit");             
}
// t->Draw()
// TH1F *h1 = new TH1F("h1","h1 title",2000,0,2000);("epsd >> h1"); 
// h1->GetXaxis()->SetRange(50,500) ; h1->Draw();