#include "TSpectrum2.h"
#include "Riostream.h"
#include "TH2F.h"
#include "TCanvas.h"
#include "TSpectrum2Fit.h"
#include "TFile.h"
void tutor1() {

   Int_t i, j, nfound;

   Int_t xbin, ybin;

   Int_t nbinsx = 500;

   Int_t nbinsy = 500;  

   Float_t xmin  = 1000.;

   Float_t xmax  = 5000.;

   Float_t ymin  = 0.;

   Float_t ymax  = 7000.;

   Float_t ** source = new float *[nbinsx];  

   Float_t ** dest = new float *[nbinsx];     

   for (i=0;i<nbinsx;i++)

                                                source[i]=new float[nbinsy];

   for (i=0;i<nbinsx;i++)

                                                dest[i]=new float[nbinsy];

   TH2F *search = new TH2F("search","High resolution peak searching",nbinsx,xmin,xmax,nbinsy,ymin,ymax);

   TFile *f = new TFile("D:/documents/GANIL/histo.root");

   search=(TH2F*) f->Get("ht;1");
   search->Print("base");

   TCanvas *Searching = new TCanvas("Searching","Two-dimensional fitting using Stiefel-Hestens method",10,10,1000,700);

   TSpectrum2 *s = new TSpectrum2();


   for (i = 0; i < nbinsx; i++){

     for (j = 0; j < nbinsy; j++){

                    source[i][j] = (Float_t)search->GetBinContent(i + 1,j + 1);

                 }

	}  

   nfound = s->SearchHighRes(source, dest, nbinsx, nbinsy, 3, 20, kTRUE, 3, kFALSE, 3);  

   printf("Found %d candidate peaks\n",nfound);

   Bool_t *FixPosX = new Bool_t[nfound];

   Bool_t *FixPosY = new Bool_t[nfound];   

   Bool_t *FixAmp = new Bool_t[nfound];     

   Float_t *PosX = new Float_t[nfound];        

   Float_t *PosY = new Float_t[nfound];

   Float_t *Amp = new Float_t[nfound];     

   Float_t *AmpXY = new Float_t[nfound];        

   PosX = s->GetPositionX();

   PosY = s->GetPositionY();     
   s->Print();
   
   cout<<"positions acquired!"<<endl;
   for(i = 0; i< nfound ; i++){
     cout<<"i1= "<<i<<endl; 

      FixPosX[i] = kFALSE;

      FixPosY[i] = kFALSE;     

      FixAmp[i] = kFALSE;   

      Amp[i] = source[(int)(PosX[i]+0.5)][(int)(PosY[i]+0.5)];      //initial values of peaks amplitudes, input parameters         
      cout<<Amp[i]<<endl;
      AmpXY[i] = 0;
      cout<<"i2= "<<i<<endl;
	  xbin=1+(int)(PosX[i]+0.5);
		ybin=1+(int)(PosY[i]+0.5);
	  PosX[i]=search->GetXaxis()->GetBinCenter(xbin);
	  PosY[i]=search->GetYaxis()->GetBinCenter(ybin);
   }
  
   //filling in the initial estimates of the input parameters

   TSpectrum2Fit *pfit=new TSpectrum2Fit(nfound);

   pfit->SetFitParameters(xmin, xmax-1, ymin, ymax-1, 1000, 0.1, pfit->kFitOptimChiCounts, pfit->kFitAlphaHalving, pfit->kFitPower2, pfit->kFitTaylorOrderFirst);  

   pfit->SetPeakParameters(2, kFALSE, 2, kFALSE, 0, kTRUE, PosX, (Bool_t *) FixPosX, PosY, (Bool_t *) FixPosY, PosX, (Bool_t *) FixPosX, PosY, (Bool_t *) FixPosY, Amp, (Bool_t *) FixAmp, AmpXY, (Bool_t *) FixAmp, AmpXY, (Bool_t *) FixAmp);     

   pfit->SetBackgroundParameters(0, kFALSE, 0, kTRUE, 0, kTRUE);  
	cout <<"fit started!"<<endl;
   pfit->FitStiefel(source);
	cout <<"Fit finished"<<endl;
    for (i = 0; i < nbinsx; i++){

     for (j = 0; j < nbinsy; j++){

                  search->SetBinContent(i + 1, j + 1,source[i][j]);

                 }

   }  

   search->Draw("SURF");
   
}