#include <TSpectrum2.h>
#include "TCanvas.h"
#include "TRandom.h"
#include "TH2.h"
#include "TF2.h"
#include "TMath.h"
#include "TROOT.h"
#include "TCanvas.h"
#include "TMath.h"
#include "TH1.h"
#include "TF1.h"
#include "TRandom.h"

#include "TVirtualFitter.h"

TSpectrum2 *s;

Int_t npeaks = 30;
Double_t fpeaks2(Double_t *x, Double_t *par) {
    Double_t result = 0.1;
    for (Int_t p=0;p<npeaks;p++) {
        Double_t norm   = par[5*p+0];
        Double_t mean1  = par[5*p+1];
        Double_t sigma1 = par[5*p+2];
        Double_t mean2  = par[5*p+3];
        Double_t sigma2 = par[5*p+4];
        result += norm*TMath::Gaus(x[0],mean1,sigma1)*TMath::Gaus(x[1],mean2,sigma2);
    }
    return result;
}
void newtry3()

{
   Int_t nbinsx = 396;
   Int_t nbinsy = 279;
   Double_t xmin   = 0;
   Double_t xmax   = (Double_t)nbinsx;
   Double_t ymin   = 0;
   Double_t ymax   = (Double_t)nbinsy;
   Double_t dx = (xmax-xmin)/nbinsx;
   Double_t dy = (ymax-ymin)/nbinsy;
   //delete h2;

   //generate n peaks at random
   Double_t par[3000];
   Int_t p;
   Float_t** source = new Float_t*[nbinsx];
   for (int i=0;i<nbinsx;i++) source[i]=new Float_t[nbinsy];
   Float_t** dest = new Float_t*[nbinsx];
   for (int i=0;i<nbinsx;i++) dest[i]=new Float_t[nbinsy];

   const int n=396,m=279;
   float content;
   TH2F * myhist2 = new TH2F("myhist2", "myhist2", n, 0, n, m, 0,m);
    TCanvas *c1 = (TCanvas*)gROOT->GetListOfCanvases()->FindObject("c1");
    if (!c1) c1 = new TCanvas("c1","c1",10,10,1000,700);

    ifstream fin("image65_2.txt");
    for (int i=0;i<n;i++) {
      for (int j=0;j<m;j++) {
         fin>>content;
         if (content>=0.61||content<=0.43) {
            myhist2->SetBinContent(i+1,j+1,content);
            source[i][j] = myhist2->GetBinContent(i + 1,j + 1);
         }
      }
    }


   TSpectrum2 *s = new TSpectrum2(200);
   //now the real stuff: Finding the peaks
   Int_t nfound = s->SearchHighRes(source, dest, n,m, 4.2, 3, kTRUE, 1, kFALSE, 1);

   //searching good and ghost peaks (approximation)
   Int_t pf,ngood = 0;
   Float_t *xpeaks = s->GetPositionX();
   Float_t *ypeaks = s->GetPositionY();

   for (p=0;p<npeaks;p++) {
      for (pf=0;pf<nfound;pf++) {
         Double_t diffx = TMath::Abs(xpeaks[pf] - par[5*p+1]);
         Double_t diffy = TMath::Abs(ypeaks[pf] - par[5*p+3]);
         if (diffx < 2*dx && diffy < 2*dy) ngood++;
     }
   }
   if (ngood > nfound) ngood = nfound;

   //Search ghost peaks (approximation)
   Int_t nghost = 0;
   for (pf=0;pf<nfound;pf++) {
      Int_t nf=0;
      for (p=0;p<npeaks;p++) {
         Double_t diffx = TMath::Abs(xpeaks[pf] - par[5*p+1]);
         Double_t diffy = TMath::Abs(ypeaks[pf] - par[5*p+3]);
         if (diffx < 2*dx && diffy < 2*dy) nf++;
      }
      if (nf == 0) nghost++;
   }

   s->Print();
   printf(" Found=%d, Good=%d, Ghost=%d\n",nfound,ngood,nghost);

   TPolyMarker *pm = new TPolyMarker(nfound, xpeaks, ypeaks);
   myhist2->GetListOfFunctions()->Add(pm);
   pm->SetMarkerStyle(23);
   pm->SetMarkerColor(kRed);
   pm->SetMarkerSize(2.);

    myhist2->Draw("COL");
    c1->Update();
}

void peaks2(Int_t maxpeaks=50) {
    s = new TSpectrum2(2*maxpeaks);
    for (int i=0; i<10; ++i) {
        npeaks = (Int_t)gRandom->Uniform(5,maxpeaks);
        findPeak2();
    }
}