//#include "TH1D.h"
//#include "TCanvas.h"
//gROOTsadasfaf->SetOptFit(1111);

gSystem->AddIncludePath("/home/perthi/aliroot-current");

#include "TSpectrum.h"

//#include  "/home/perthi/aliroot-current/HLT/CALO/AliHLTCaloConstants.h"
#include  "AliHLTCaloConstants.h"

using namespace EMCAL;

#define NALGORITHMS 3

TH2D *map[3];


#define TARGET_EDGE 375  //6GeV / 16 MeV

void calib()
{
  map[0] = new TH2D(  "map chisquare",  "map chisquare",       NZROWSMOD, 0, NZROWSMOD -1, NXCOLUMNSMOD, 0, NXCOLUMNSMOD -1 );
  map[1] = new TH2D("map Peak-Finder", "map Peak-Finder;x;y",  NZROWSMOD, 0, NZROWSMOD -1, NXCOLUMNSMOD, 0, NXCOLUMNSMOD -1 );
  map[2] = new TH2D("map Crude", "map Crude;x;y",              NZROWSMOD, 0, NZROWSMOD -1, NXCOLUMNSMOD, 0, NXCOLUMNSMOD -1 );
  gStyle->SetOptFit(1111);
  char rawanalyzers[5][256];
  
  for(int i=0; i< 3; i++)
    {
      for( int j=0; j < 256; j++ )
	{
	  rawanalyzers[i][j] = 0;
	}
    }

  sprintf(&rawanalyzers[0][0], "%s", "Chi Square Fit");
  sprintf(&rawanalyzers[1][0], "%s", "Peak-Finder");
  sprintf(&rawanalyzers[2][0], "%s", "Crude");
  
  for(int i=0; i< 3; i++ )
    {
      cout << __FILE__ << ":" << __LINE__ << "analyzer[" << i << "] = " << rawanalyzers[i] << endl;
    }
  
  TFile *f=new TFile("pscalib-6GeV.root", "read");
  TFile *fout=new TFile("6GeV-fitted.root", "recreate");
  // char hname[4096] = {0};
  // TSpectrum *s= new TSpectrum(5);
  
  TCanvas *c = new TCanvas("test", "test"); 
  c->SetLogy();

  for(int i=0; i < 3; i++)
    {
      for( int row=0; row < 8; row ++ )
	{
	  for( int col=0; col < 8; col ++ )
	   {
	     char hname[4096] = {0}; 
	     sprintf(hname, "%s, row %d, col %d", &rawanalyzers[i][0], ( row ), ( col + 40 ) );
	     
	     TH1D *h = (TH1D*)f->Get(hname);
	    
	     //	     map->SetBinContent( col +40 +1, row +1 );
 
	     if(h != 0)
	       {
		 //	 h->GetYaxis()->SetLog();
		
		 //	 h->GetCanvas();
		 //	 h->SetOptFit(1);
		 cout << __FILE__ << ":" << __LINE__ << "INFO, FOUND !!! histogram with name:  " << hname  << endl;
		 h->Rebin(8);
		 TSpectrum *s= new TSpectrum(100);
		 int n= s->Search(h);
		 Float_t *xpeaks = s->GetPositionX();
		 Float_t tmpmax = 0;

		 for( int k=0; k < n; k++ )
		   {
		     if(  xpeaks[k] > tmpmax )
		       {
			 tmpmax = xpeaks[k];
		       }
		   }
				 
		 delete s;
		 h->Fit("gaus","","",tmpmax -30, tmpmax + 60  );
		 
		 TF1 *formual=h->GetFunction("gaus");
		 
		 cout << "!!!!!!!! P0 ="<< formual->GetParameter(0) 
		      <<  ";  P1 ="  << formual->GetParameter(1) 
		      << "; P2="<<  formual->GetParameter(2)  << endl;

		 map[i]->SetBinContent( col +40 +1, row +1, formual->GetParameter(1) );
		 
		 //SetBinContent(Int_t binx, Int_t biny, Double_t content)
		 fout->cd();
		 h->Write();
	       }
	     else
	       {
		 cout << __FILE__ << ":" << __LINE__ << "ERROR, COULD NOT FIND !!! histogram with name:  " << hname  << endl;
	       }
	   }

	}
      map[i]->Write();
    }
  
  
  TH1D *calibfactors = new TH1D("test", "test", 99, 0, 100 );
  calibfactors->Draw();
  f->Close();
  fout->Close();
  //  delete s;
  delete f;
  delete fout; 
}