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

Int_t npeaks = 30;
// fit function def
   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 sigma  = par[3*p+4];
      Double_t mean  = par[3*p+3];
	  result += norm*TMath::Gaus(x[0],mean,sigma);
   }
   return result;
   }



void test1sig ( const string & filename, int n_iter , int max_pic)
{
	
	
//opening test file
	int xa[1024];
	int y[1024];

  FILE *fp = fopen(filename.c_str(),"r");
  if (fp != NULL){
    bool  error = false;

// reading file
    for (int i = 0; i < 1000; i++){
      int n = fscanf (fp, "%d %d", &xa[i], &y[i]);
      if (n != 2){
        printf( "BOF..." );
        i = 1024;
        error = true;
      }
    }

// making a new ROOT histogram
    if (! error){
      TH1I *histo = new TH1I("histo", "My histo", 300, 100., 300.);
      
      for (int i = 0; i < 300; i++){
        histo->SetBinContent(i,y[i]);
      }

// fitting procedure

//making a new canvas
   TCanvas *c1 = new TCanvas("c1","c1",10,10,500,600);
   pad2 = new TPad("pad1","My pad1",0.02,0.02,0.98,0.5,33);
   pad1 = new TPad("pad2","My pad2",0.02,0.5,0.98,0.98,33);
   pad1->Draw();
   pad2->Draw();
   pad1->SetLogy();
   pad2->SetLogy();
   pad1->cd();
   
   TH1F *h2 = (TH1F*)histo->Clone("h2");
   
//resettint and starting counter
	gBenchmark->Reset(); 
	gBenchmark->Start("hsimple"); 

//Use TSpectrum to find the peak candidates
   TSpectrum *s = new TSpectrum(max_pic,2);
   Double_t mymax;
   mymax = histo->GetMaximum();
   mymax=3/mymax;
   Int_t nfound = s->Search(histo,2,"new",mymax);
   s->Background(histo,n_iter,"same");
	for (int i = 0; i < 300; i++){
        h2->SetBinContent(i,h2->GetBinContent(i)- histo->GetBinContent(i));
      }
	
	pad2->cd();
	printf("Found %d candidate peaks to fit...",nfound);
   

//estimate linear background
   Double_t par[3000];
   par[0] = 0.50;
   par[1] = -1/10000;
   
//TF1(const char* name, const char* formula, Double_t xmin = 0, Double_t xmax = 1)
   TF1 *fline = new TF1("fline","pol1",1.,1000);
   //h2->Fit("fline","qn");
//TF1(const char* name, void* fcn, Double_t xmin, Double_t xmax, Int_t npar)
   TF1 *fit = new TF1("fit",fpeaks,0,1000,2+3*npeaks);

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

   Float_t *xpeaks = s->GetPositionX();
   int_t p=0;
   for (p=0;p<nfound;p++) {
      Float_t xp = xpeaks[p];
      Int_t bin = histo->GetXaxis()->FindBin(xp);
      Float_t yp = histo->GetBinContent(bin);
      
	  if (yp-TMath::Sqrt(yp) < fline->Eval(xp)) continue;
	  par[3*npeaks+4]=2;
	  fit->SetParLimits(3*npeaks+4,1.8,3);
	  par[3*npeaks+2] = yp;
      par[3*npeaks+3] = xp;
	  npeaks++;
   }
   
   printf("Found %d useful peaks to fit...",npeaks);
   printf("Now fitting: Be patient");
      
//fit the found peaks with fpeaks=1 gauss function
   
   //TVirtualFitter::Fitter(h2,10+3*npeaks); //we may have more than the default 25 parameters
   fit->SetParameters(par);
   fit->SetNpx(500);
   h2->Fit("fit","B"); 
   gBenchmark->Show("hsimple");
   c1->Update();
    }

//estimation of background
	
	

	
/*	printf("Integral of fit function = %g\n",fit->Integral(par[3]-2*par[4],par[3]+2*par[4]));
	printf("Integral of function = %g\n",histo->Integral(par[3]-2*par[4],par[3]+2*par[4]));
	printf("Integral of back = %g\n",fline->Integral(par[3]-2*par[4],par[3]+2*par[4]));
	*/
  }
  else
    printf( "ERREUR OUVERTURE FICHIER" );
}