//
//  main.cpp
//  uclatest
//
//  Created by Haochen Yan on 16/11/18.
//  Copyright © 2016年 Haochen Yan. All rights reserved.
//


/*
#include "TCanvas.h"
#include "TMath.h"
#include "TH1.h"
#include "TF1.h"
#include "TRandom.h"
#include "TSpectrum.h"
#include "TVirtualFitter.h"


    
Int_t npeaks = 30;
Double_t fpeaks(Double_t *x, Double_t *par) {
    Double_t result = par[0] + par[1]*x[0];
    for (Int_t p=0;p<npeaks;p++) {
        Double_t norm  = par[3*p+2];
        Double_t mean  = par[3*p+3];
        Double_t sigma = par[3*p+4];
        result += norm*TMath::Gaus(x[0],mean,sigma);
    }
    return result;
}
void newtry3(Int_t np=10)
{
    npeaks = TMath::Abs(np);
    TH1F *h = new TH1F("h","test",500,0,1000);
    //generate n peaks at random
    Double_t par[3000];
    par[0] = 0.8;
    par[1] = -0.6/1000;
    Int_t p;
    for (p=0;p<npeaks;p++) {
        par[3*p+2] = 1;
        par[3*p+3] = 10+gRandom->Rndm()*980;
        par[3*p+4] = 3+2*gRandom->Rndm();
    }
    TF1 *f = new TF1("f",fpeaks,0,1000,2+3*npeaks);
    f->SetNpx(1000);
    f->SetParameters(par);
    TCanvas *c1 = new TCanvas("c1","c1",10,10,1000,900);
    c1->Divide(1,2);
    c1->cd(1);
    h->FillRandom("f",200000);
    h->Draw();
    TH1F *h2 = (TH1F*)h->Clone("h2");
    //Use TSpectrum to find the peak candidates
    TSpectrum *s = new TSpectrum(2*npeaks);
    Int_t nfound = s->Search(h,2,"",0.10);
    printf("Found %d candidate peaks to fit\n",nfound);
    //Estimate background using TSpectrum::Background
    TH1 *hb = s->Background(h,20,"same");
    if (hb) c1->Update();
    if (np <0) return;
    
    //estimate linear background using a fitting method
    c1->cd(2);
    TF1 *fline = new TF1("fline","pol1",0,1000);
    h->Fit("fline","qn");
    //Loop on all found peaks. Eliminate peaks at the background level
    par[0] = fline->GetParameter(0);
    par[1] = fline->GetParameter(1);
    npeaks = 0;
    Double_t *xpeaks = s->GetPositionX();
    for (p=0;p<nfound;p++) {
        Double_t xp = xpeaks[p];
        Int_t bin = h->GetXaxis()->FindBin(xp);
        Double_t yp = h->GetBinContent(bin);
        if (yp-TMath::Sqrt(yp) < fline->Eval(xp)) continue;
        par[3*npeaks+2] = yp;
        par[3*npeaks+3] = xp;
        par[3*npeaks+4] = 3;
        npeaks++;
    }
    printf("Found %d useful peaks to fit\n",npeaks);
    printf("Now fitting: Be patient\n");
    TF1 *fit = new TF1("fit",fpeaks,0,1000,2+3*npeaks);
    //we may have more than the default 25 parameters
    TVirtualFitter::Fitter(h2,10+3*npeaks);
    fit->SetParameters(par);
    fit->SetNpx(1000);
    h2->Fit("fit");
}
*/
#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;
    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;
     Double_t** source = new Double_t*[nbinsx];
    for (int i=0;i<nbinsx;i++)
    {
        source[i]=new Double_t[nbinsy];
    }
    Double_t** dest = new Double_t*[nbinsx];
    for (int i=0;i<nbinsx;i++)
    {
        dest[i]=new Double_t[nbinsy];
    }
    
  
    /*
    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);
    
    */
    
    
    
    
    
    const int n=396,m=279;
    float content;
    //TH2F * myhist = new TH2F("myhist", "myhist", n, 0, n, m, 0,m);
    TH2F * myhist2 = new TH2F("myhist2", "myhist2", n, 0, n, m, 0,m);
    
    ifstream fin("image65_2.txt");
    for(int i=0;i<n;i++)
        
    {
        for(int j=0;j<m;j++)
            
        {
            fin>>content;
            
            
            
            
            
            if(content>=0.6||content<=0.45)
            {
                
                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, 0.1, 0.4, kTRUE, 1, kFALSE, 1);
  
    //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();
    myhist2->Draw("CONT2");
    
    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();
}