#define ROOT6

#include <TSpectrum2.h>
#include "TCanvas.h"
#include "TH2.h"
#include "TMath.h"
#include "TCanvas.h"
#include "TMath.h"
#include "TH1.h"
#include "TFile.h"
#include "TPolyMarker.h"

#include <fstream>

Int_t newtry3(Double_t sig = 2., Double_t thr = 0.01)
{
	Double_t sigma = TMath::Max(sig, 1.);
	Double_t threshold = TMath::Max(thr, 0.000001);
	
	Int_t npeaks = 100;
   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;
#ifdef ROOT6  
   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];
#else
   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];
#endif
   const int n=nbinsx, m=nbinsy;
   float content;
   TFile *f = new TFile("newtry3.root", "RECREATE");
   TH2F * myhist2 = new TH2F("myhist2", "Input spectrum", n, 0, xmax, m, 0,ymax);
   TH2F * result = new TH2F("result", "Convoluted spectrum", n, 0, xmax, m, 0,ymax);
	myhist2->SetStats(0);
	result->SetStats(0);
    ifstream fin("image65_2.txt");
    for (int i=0;i<n;i++) {
      for (int j=0;j<m;j++) {
         fin>>content;
         // cut away noise around x=20
         if (i < 50 ) content = 0.;
         // cut away the plain around z = 0.5
//         if ( content<=0.6 && content > 0.45 ) content = 0.;
         // or one might cut plain and craters, keep the vulcanos only
         if ( content<=0.6 ) content = 0.;
         // make sure all cells are filled
         myhist2->SetBinContent(i+1, j+1, content);
         source[i][j] = content;
         dest[i][j] = 0.;
      }
    }
	fin.close();

   TSpectrum2 *s = new TSpectrum2(npeaks);
   //now the real stuff: Finding the peaks
   Int_t nfound = s->SearchHighRes(source, dest, n,m, sigma, threshold, kTRUE, 1, kFALSE, 1);
//   Int_t nfound = s->Search(myhist2,sigma,"col",threshold);

   TPolyMarker * pm = (TPolyMarker*)myhist2->GetListOfFunctions()->FindObject("TPolyMarker");
   if (pm) {
      myhist2->GetListOfFunctions()->Remove(pm);
      delete pm;
   }
   pm = new TPolyMarker(nfound, s->GetPositionX(), s->GetPositionY());
   myhist2->GetListOfFunctions()->Add(pm);
   pm->SetMarkerStyle(23);
   pm->SetMarkerColor(kRed);
   pm->SetMarkerSize(1.3);
	// fill hist with resulting convoluted spectrum
   for (int i=0;i<n;i++) {
      for (int j=0;j<m;j++) {
         result->SetBinContent(i+1,j+1,dest[i][j]);
      }
   }
   TCanvas *cc = new TCanvas("newtry3", "newtry3", 600, 1200);
	cc->Divide(1,2);
   cc->cd(1);
   myhist2->Draw("colz");
   myhist2->GetXaxis()->SetRangeUser(150., 270.);
   myhist2->GetYaxis()->SetRangeUser(60., 180.);
   cc->cd(2);
   result->Draw("colz");
   result->GetXaxis()->SetRangeUser(150., 270.);
   result->GetYaxis()->SetRangeUser(60., 180.);
	s->Print();
   f->Write();
   return nfound;
}