// Example to illustrate the 2-d peak finder (class TSpectrum2).
// This script generates a random number of 2-d gaussian peaks
// The position of the peaks is found via TSpectrum2
// To execute this example, do
//  root > .x peaks2.C  (generate up to 50 peaks by default)
//  root > .x peaks2.C(10) (generate up to 10 peaks)
//  root > .x peaks2.C+(200) (generate up to 200 peaks via ACLIC)
//
// The script will iterate generating a new histogram having
// between 5 and the maximun number of peaks specified.
// Double Click on the bottom right corner of the pad to go to a new spectrum
// To Quit, select the "quit" item in the canvas "File" menu
//
//Author: Rene Brun
   
#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 = 0;
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(0);
   //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;
   Float_t *xpeaks = s->GetPositionX();
   Float_t *ypeaks = s->GetPositionY();
   for (p=0;p<npeaks;p++) {
      for (Int_t 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 (Int_t 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);
   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();
   }
}