#include <iostream>
#include <iomanip>
#include <string>
#include <math>

#include "TH2F.h"
#include "TH1F.h"
#include "TGraph.h"
#include "TF2.h"
#include "TF1.h"

class TwoDHisto{
 private:
  TH2F** h;
  int n;
  float * emin;
  float * emax;
  
 public:
  TwoDHisto(int n);
  ~TwoDHisto();
  float GetEmin(int i);
  float GetEmax(int i);
  float * Getemin();
  float * Getemax();
  void SetEmin(int i,float v);
  void SetEmin(int i,float v);
  TH2F * getHisto(int i);
  void fill(int x,float energy,float t);
  void draw(int x);
  void setbins(int x,int nbins,float xmin,float xmax);
  void fit(int x,char * funct,float * quad,float * lin,float * pos);
}; 

/*******************************************************************************/

TwoDHisto::TwoDHisto(int n){
  if(h)delete [] h;
  h = new TH2F*[n];
  emin = new float [n];
  if(emin==NULL)cout << "Didn't manage to allocate memory for emin!" <<endl;
  emax = new float [n];
  if(emax==NULL)cout << "Didn't manage to allocate memory for emax!" <<endl;
  char name[64];
  for(int i=0;i < n;i++){
    
    sprintf(name,"%d", i);
    //    h[i] = new TH2F(name,"2Dhisto",int(emax[i]-emin[i]),emin[i],emax[i],500,0.0,10000.0);
    h[i] = new TH2F(name,"2Dhisto",200,900,1000,500,0.0,10000.0);
  }    


}


TwoDHisto::~TwoDHisto(){
  
  delete[] h;
  delete[] emin;
  delete[] emax;
}

float TwoDHisto::GetEmin(int i){
  return emin[i];
}

float TwoDHisto::GetEmax(int i){
  return emax[i];
}

float * TwoDHisto::Getemin(){
  return emin;
}

float * TwoDHisto::Getemax(){
  return emax;
}

void TwoDHisto::SetEmin(int i,float v){
  emin[i] = v;
}

void TwoDHisto::SetEmax(int i,float v){
  emax[i] = v;
}

TH2F * TwoDHisto::getHisto(int i){
  return h[i];
}

void TwoDHisto::fill(int x,float energy,float t){
  h[x]->Fill(energy,t);
}

void TwoDHisto::draw(int x){
  h[x]->Draw();
}

void TwoDHisto::setbins(int x,int nbinsx,float xmin,float xmax,int nbinsy,float ymin,float ymax){
  h[x]->SetBins(nbinsx,xmin,xmax,nbinsy,ymin,ymax);
}  

void TwoDHisto::fit(int x,char * funct){
  h[x]->Fit(funct,"RNV");
  //"R" := Use the range specified in the function range
  //"N" := Do not store the graphics function, do not draw
  //"0" := Do not plot the result of the fit
  //"Q" := Quiet mode
  //"V" := Verbose mode
/*   quad[x]=f->GetParameter(1); */
/*   lin[x]=f->GetParameter(2); */
/*   pos[x]=f->GetParameter(3); */
}

/******************************************************************************************/

void twodfit(char *filename){
  
  float t = 0.0;
  float t_abs = 0.0;
  float energy = 0.0;
  
  float chisq;
  float p1,p2,p3,p4,p5;
  float fpar1,fpar2,fpar3,fpar4,fpar5,fpar6;
  float fqp1,fqp2,fqp3,fqp4,fqp5,fqp6;
  int temp;
  float tempmin,tempmax;
  
  string proc;
  ofstream output;
  output.open("parmeters.txt");
  ifstream input;
  input.open(filename);

  cout << "---------------------------" << endl;
  cout << "HOW MANY PEAKS?" << endl;
  cout << "# of peaks: ";
  cin >> temp;
  cout << "RANGES" << endl;
  
  cout << "n=" << temp << endl;

  const int n = temp;
  
  //  float *emin;
  float *pos;
  float *quad1;
  float *lin1;
  float *quad2;
  float *lin2;
  float *position;

 
  quad1 = new float [n];
  if(quad1==NULL)cout << "Didn't manage to allocate memory for quad1!" <<endl;
  lin1 = new float [n];
  if(lin1==NULL)cout << "Didn't manage to allocate memory for lin1!" <<endl;
  quad2 = new float [n];
  if(quad2==NULL)cout << "Didn't manage to allocate memory for quad2!" <<endl;
  lin2 = new float [n];
  if(lin2==NULL)cout << "Didn't manage to allocate memory for lin2!" <<endl;
  position = new float [n];
  if(position==NULL)cout << "Didn't manage to allocate memory for position!" <<endl;
  pos = new float [n];
  if(pos==NULL)cout << "Didn't manage to allocate memory for pos!" <<endl;
  
  TwoDHisto * h = new TwoDHisto(n);

  cout << "Memomry allocated!" << endl;
  for(int i = 0;i<n;i++){
    cout << "Range #" << i+1 << "(min): ";
    cin >> tempmin;
    h->SetEmin(i,tempmin);
    cout << h->GetEmin(i) << endl;
    cout << "Range #" << i+1 << "(max): ";
    cin >> tempmax;
    h->SetEmax(i,tempmax);
    cout << h->GetEmax(i) << endl;
    cout << "Position #" << i+1 << ": ";
    cin >> position[i];
    h->setbins(i,int(h->GetEmax(i)-(h->GetEmin(i))),h->GetEmin(i),h->GetEmax(i),1000,0.0,10000);
  }

  cout << "Making histos..." << endl;
  cout << "Proceed (y/n)? : ";
  cin >> proc; 
  if(proc == "y"){
 
    cout << "Proceeding..." << endl;
    TF2 * f2 = (TF2*)gROOT->FindObject("f2");
    if(f2)delete f2;
    f2 = new TF2("f2","([5]/[4])*exp(-(x-([1]*y*y+[2]*y+[3]))*(x-([1]*y*y+[2]*y+[3]))*0.5/([4]*[4]))",0,10000,0,10000);
    
    
    for(int j = 0;j<n;j++){
      
      cout << "j = " << j << endl;

      //f2->SetRange(h->GetEmin(j),h->GetEmax(j),0,10000);
      p1 = 1e-08;
      p2 = 1e-05;
      p3 = position[j];
      p4 = 10;
      p5 = 100;
      
      f2->SetParameter(1,p1);
      f2->SetParameter(2,p2);
      f2->SetParameter(3,p3);
      f2->SetParameter(4,p4);
      f2->SetParameter(5,p5);
      
      f2->SetParLimits(1,-1e-07, 1e-07);
      f2->SetParLimits(2,-1e-04,1e-04);
      f2->SetParLimits(3,position[j]-10,position[j]+10);
      f2->SetParLimits(4,0,200);
      f2->SetParLimits(5,0,1000);
      
      while(1) {
	input >> t >> t_abs >> energy;
	if (!input.good()) break;
	if((h->GetEmin(j) <= energy) && (h->GetEmax(j) >= energy)){
	  h->fill(j,energy,t);
	}
      }
      cout << "Histos filled!" <<  endl;


      cout << "Fitting #" << j+1 << "..." << endl; 
      
      (h->getHisto(j))->Fit("f2","NL"); 
      
      cout << "Fitting #" << j+1 << "ready!" << endl;
      
      cout << "Chisquared = " << (f2->GetChisquare()/f2->GetNDF()) << endl;




      output << setw(5) << f2->GetParameter(1) << setw(15) << f2->GetParameter(2) << setw(15) << f2->GetParameter(3) << endl;

      //(f2->GetHistogram())->Draw("surf2");
      

    }
  }

  delete quad1;
  delete lin1;
  delete quad2;
  delete lin2;
  delete position;
  delete pos;

  output.close();
  input.close();
}