#include "TSpectrum2.h"
#include "Riostream.h"
#include "TH2F.h"
#include "TCanvas.h"
#include "TSpectrum2Fit.h"
#include "TFile.h"
void fit_peak(Int_t nbinsx,Int_t xmin,Int_t xmax, Int_t nbinsy, Int_t ymin, Int_t ymax){
  TSpectrum2* s= new TSpectrum2(20,1);
  //searching
  TFile* f=new TFile("histo.root");
  TH2F* histo;
  histo=(TH2F*)f->Get("ht;1");
  Int_t val=s->Search(histo,3.,"",0.2);
  Int_t n=s->GetNPeaks();
  cout<<n<<" peaks found!"<<endl;
  ofstream out;
  out.open("peaks.dat");
  /*
  for (Int_t i=0;i<n;i++){
    out<<i+1<<" "<<y[i]<<" "<<x[i]<<" "<<endl;
  }
  s->Print();
  out.close();
  */
  Float_t ** source = new float *[nbinsx];  

  Float_t ** dest = new float *[nbinsx];     
   for (Int_t i=0;i<nbinsx;i++)

                                                source[i]=new float[nbinsy];

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

                                                dest[i]=new float[nbinsy];

   TCanvas* srch=new TCanvas("Search","peak fit",10,10,1000,700);
   Int_t count=0;
   for(Int_t i=0;i<nbinsx;i++){
     for(Int_t j=0;j<nbinsy;j++){
       source[i][j]=histo->GetBinContent(i+1,j+1);
       if(source[i][j]>10){count=count+1;
	 //out<<" ("<<i<<", "<<j<<")=  "<<source[i][j];
       }
     }
     out<<" "<<endl;
   }
   cout<<"source filled!"<<endl;
   cout<<"count= "<<count<<endl;

   Bool_t* FixPosX=new Bool_t[n];
   Bool_t* FixPosY=new Bool_t[n];
   Bool_t* FixAmp= new Bool_t[n];

   Float_t* PosX=new Float_t[n];
   Float_t* PosY=new Float_t[n];
   Float_t* Amp= new Float_t[n];
   Float_t* AmpXY= new Float_t[n];

   Float_t binx, biny;
   cout<<"variables set!"<<endl;

   binx=(xmax-xmin)/nbinsx;
   biny=(ymax- ymin)/nbinsy;
   
   PosX=s->GetPositionX();
   PosY=s->GetPositionY();
   cout<<"positions acquired"<<endl;
   //initial values
   for(Int_t i=0;i<n;i++){
     FixPosX[i]=kFALSE;
     FixPosY[i]=kFALSE;
     FixAmp[i]=kFALSE;
     Amp[i]=source[(int)((PosX[i]-xmin)/binx+0.5)][(int)((PosY[i]-ymin)/biny+0.5)];
     cout<<"Amp"<<i<<" = "<<Amp[i]<<endl;
     cout<<"posx"<<i<<" = "<<PosX[i]<<endl;
     cout<<"binx"<<i<<" = "<<(int)((PosX[i]-xmin)/binx+0.5)<<endl;
     cout<<"posy"<<i<<" = "<<PosY[i]<<endl;
     cout<<"biny"<<i<<" = "<<(int)((PosY[i]-ymin)/biny+0.5)<<endl;
     AmpXY[i]=0;
   }
       cout<<"initial values set!"<<endl;
         
     //fitting
     TSpectrum2Fit *pfit=new TSpectrum2Fit(n);
     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 parameters set!"<<endl;
     pfit->FitAwmi(source);
     cout<<"fit performed!"<<endl;
     for(Int_t i=0; i<nbinsx;i++){
       for(Int_t j=0; j<nbinsx;j++){
	 histo->SetBinContent(i+1,j+1,source[i][j]);
       }
     }
     histo->Draw("SURF");
      
   out.close();
}