#include "TSpectrum2.h"
#include "TCanvas.h"
#include "TRandom.h"
#include "TH2.h"
#include "TF2.h"
#include "TMath.h"
#include "TROOT.h"

TSpectrum2 *s;
TH2F *h2 = nullptr;
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 findPeak2()
{
   printf("Generating histogram with %d peaks\n", npeaks);
   Int_t nbinsx = 200;
   Int_t nbinsy = 200;
   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;
   h2 = new TH2F("h2", "test", nbinsx, xmin, xmax, nbinsy, ymin, ymax);
   h2->SetStats(false);
   // generate n peaks at random
   Double_t par[3000];
   Int_t p;
   for (p = 0; p < npeaks; p++) {
      par[5 * p + 0] = gRandom->Uniform(0.2, 1);
      par[5 * p + 1] = gRandom->Uniform(xmin, xmax);
      par[5 * p + 2] = gRandom->Uniform(dx, 5 * dx);
      par[5 * p + 3] = gRandom->Uniform(ymin, ymax);
      par[5 * p + 4] = gRandom->Uniform(dy, 5 * dy);
   }
   TF2 *f2 = new TF2("f2", fpeaks2, xmin, xmax, ymin, ymax, 5 * npeaks);
   f2->SetNpx(100);
   f2->SetNpy(100);
   f2->SetParameters(par);
   TCanvas *c1 = (TCanvas *)gROOT->GetListOfCanvases()->FindObject("c1");
   if (!c1)
      c1 = new TCanvas("c1", "c1", 10, 10, 1000, 700);
   h2->FillRandom("f2", 500000);

   // now the real stuff: Finding the peaks
   Int_t nfound = s->Search(h2, 2, "col");

   // searching good and ghost peaks (approximation)
   Int_t pf, ngood = 0;
   Double_t *xpeaks = s->GetPositionX();
   Double_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++;
   }
   c1->Update();

   s->Print();
   printf("Gener=%d, Found=%d, Good=%d, Ghost=%d\n", npeaks, nfound, ngood, nghost);
   if (!gROOT->IsBatch()) {
      printf("\nDouble click in the bottom right corner of the pad to continue\n");
      c1->WaitPrimitive();
   }
}
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();
   }
}