#include <cmath>
#include <cstdlib>
#include <sstream>
#include <iomanip>
#include "Cell.h"
#include "TGraph2DErrors.h"

using namespace std;
/*************************************************************************************/
Cell::Cell(double delX, double delY, double rad, double X0, double Y0) : dx(delX),dy(delY),r(rad),x0(X0),y0(Y0){
  InitializeCells();//figure out n_cell and assign
  f = NULL;
}

Cell::Cell(const char* filename){//load from filename
  f = NULL;
  OpenFile(filename,"READ");
  ReadData(filename);
    //set TDirectory(0) to keep read objects in memory after closing file
  CloseFile();
}

Cell::~Cell(){}

void Cell::ResizeVectors(int n){
  xlo.ResizeTo(n);
  xmid.ResizeTo(n);
  xhi.ResizeTo(n);
  ylo.ResizeTo(n);
  ymid.ResizeTo(n);
  yhi.ResizeTo(n);
}
void Cell::InitializeCells(){
  int k;
  double xtemp,ytemp;
  
  //determine number of cells just overtaking circle horizontally and vertically
  nx = 2*ceil((r-dx/2)/dx) + 1;
  ny = 2*ceil((r-dy/2)/dy) + 1;
  ResizeVectors(nx*ny);
  
  //determine resulting bounding square
  X_bound = (dx/2 + (nx-1)/2*dx);
  Y_bound = (dy/2 + (ny-1)/2*dy);

  k = 0;
  ytemp = -Y_bound;
  for(int i=0;i<ny;i++){
    xtemp = - X_bound;
    for(int j=0;j<nx;j++){
      if (InCircle(xtemp,ytemp) && InCircle(xtemp+dx,ytemp) && InCircle(xtemp,ytemp+dy) && InCircle(xtemp+dx, ytemp+dy)){
        //Assign
        xlo[k] = xtemp;
        xmid[k] = xtemp +dx/2;
        xhi[k] = xtemp+dx;
        ylo[k] = ytemp;
        ymid[k] = ytemp +dy/2;
        yhi[k] = ytemp+dy;
        //Iterate
        k++;
      }
      xtemp += dx;
    }
    ytemp += dy;
  }
  n_cell = k;
  ResizeVectors(n_cell); //delete unused space
}
void Cell::SaveToFile(const char* filename){//save to file
  OpenFile(filename,"UPDATE");
  WriteData(filename);
  CloseFile();
}

void Cell::OpenFile(const char* filename, const char* opt){
  if (f!=NULL && f->IsOpen()){
    cout << "A file is already open. Close file this file before opening another." << endl;
  }
  f = new TFile(filename,opt);
  if (f->IsZombie()){
    cout << "ERROR: Could not open file "<<filename<<"with option "<<opt<<"!"<<endl;
    cout << "Check the full path!\n";
	  exit(1);
  }
}
 
void Cell::CloseFile(){
  if (f!=NULL){
    //f->Close();
    delete f;
    f = NULL;
  }
}

void Cell::ReadData(const char* filename){
  if (f->IsZombie()||f==NULL){
    cout << "File"<<filename<<" not open for reading."<<endl;
	  exit(1);
  }
  //Define one element vectors for reading 
  TVectorD dx_t(1),dy_t(1),r_t(1),x0_t(1),y0_t(1),X_bound_t(1),Y_bound_t(1);
  
  dx_t = *(TVectorD*) f->Get("dx");
  dy_t = *(TVectorD*) f->Get("dy");
  r_t = *(TVectorD*) f->Get("r");
  x0_t = *(TVectorD*) f->Get("x0");
  y0_t = *(TVectorD*) f->Get("y0");
  X_bound_t = *(TVectorD*) f->Get("X_bound");
  Y_bound_t = *(TVectorD*) f->Get("Y_bound");
  
  dx = dx_t[0];
  dy = dy_t[0];
  r = r_t[0];
  x0 = x0_t[0];
  y0 = y0_t[0];
  X_bound = X_bound_t[0];
  Y_bound = Y_bound_t[0];
  
  //Resize vectors before assignment
  TVectorD* xlo_ptr =  (TVectorD*) f->Get("xlo");
  n_cell = xlo_ptr->GetNoElements();
  ResizeVectors(n_cell);

    
  xlo = *(TVectorD*) f->Get("xlo");
  xmid = *(TVectorD*) f->Get("xmid");
  xhi = *(TVectorD*) f->Get("xhi");
  ylo = *(TVectorD*) f->Get("ylo");
  ymid = *(TVectorD*) f->Get("ymid");
  yhi = *(TVectorD*) f->Get("yhi");
  
  nx = (int) 2*X_bound/dx;
  ny = (int) 2*Y_bound/dy;
}

void Cell::WriteData(const char* filename){
  if (f->IsZombie()||f==NULL){
    cout << "File"<<filename<<" not open for writing."<<endl;
	  exit(1);
  }
  
  //Define one element vectors for writing
  TVectorD dx_t(1),dy_t(1),r_t(1),x0_t(1),y0_t(1),X_bound_t(1),Y_bound_t(1); 
  dx_t[0] = dx;
  dy_t[0] = dy;
  r_t[0] = r;
  x0_t[0] = x0;
  y0_t[0] = y0;
  X_bound_t[0] = X_bound;
  Y_bound_t[0] = Y_bound;
  
  dx_t.Write("dx",TObject::kOverwrite);
  dy_t.Write("dy",TObject::kOverwrite);
  r_t.Write("r",TObject::kOverwrite);
  x0_t.Write("x0",TObject::kOverwrite);
  y0_t.Write("y0",TObject::kOverwrite);
  X_bound_t.Write("X_bound",TObject::kOverwrite);
  Y_bound_t.Write("Y_bound",TObject::kOverwrite);
   
  xlo.Write("xlo",TObject::kOverwrite);
  xmid.Write("xmid",TObject::kOverwrite);
  xhi.Write("xhi",TObject::kOverwrite);
  ylo.Write("ylo",TObject::kOverwrite);
  ymid.Write("ymid",TObject::kOverwrite);
  yhi.Write("yhi",TObject::kOverwrite);
}
void Cell::PrintClass(){
  xmid.Print();
}

bool Cell::InCircle(double x, double y){
  if ((x-x0)*(x-x0) + (y-y0)*(y-y0) < r*r) return true;
  else return false;
}

//int Cell::GetIndex(double x, double y){}
void Cell::GetRange(int i, double& x1, double& x2, double& y1, double& y2){
  x1 = xlo[i];
  x2 = xhi[i];
  y1 = ylo[i];
  y2 = yhi[i];
}
void Cell::GetCenter(int i, double& x, double& y){
  x = xmid[i];
  y = ymid[i];
}
TCut Cell::GetXYCut(int i,string Xstring, string Ystring){
  TCut cut1, cut2, cut3, cut4, cut;
  double x1,x2,y1,y2;
  GetRange(i,x1,x2,y1,y2);
  ostringstream convert; 
  convert << setprecision(3);
  string dumstr;
 
  convert <<x1;
  dumstr = Xstring + ">" + convert.str();
  cut1 = dumstr.c_str();

  
  convert.str("");
  convert.clear();
  dumstr.clear();  
  convert <<x2;
  dumstr = Xstring + "<" + convert.str();
  cut2 = dumstr.c_str();

  convert.str("");  
  convert.clear();
  dumstr.clear();  
  convert <<y1;
  dumstr = Ystring + ">" + convert.str();
  cut3 = dumstr.c_str();

  convert.str("");  
  convert.clear();
  dumstr.clear();  
  convert <<y2;
  dumstr = Ystring + "<" + convert.str();
  cut4 = dumstr.c_str();
  
  cut = cut1 && cut2 && cut3 && cut4;
  
  return cut;
}
/*************************************************************************************/


//Then define derived class T0_Cell...

/*************************************************************************************/
Q_Cell::Q_Cell(double delX, double delY, double rad, double X0, double Y0): Cell(delX, delY, rad, X0, Y0){
  Q_mean.ResizeTo(n_cell);
  //Q_std.ResizeTo(n_cell);
  Q_err.ResizeTo(n_cell);
}

Q_Cell::Q_Cell(const char* filename): Cell(filename){//load from file
  OpenFile(filename,"READ"); 
//set TDirectory(0) to keep read objects in memory after closing file
  ReadData(filename); //call own method
  CloseFile();
}

Q_Cell::~Q_Cell(){
 // delete Q_graph2D;
}

void Q_Cell::SaveToFile(const char* filename){//save to file
  OpenFile(filename,"UPDATE");
  WriteData(filename);
  CloseFile();
}
void Q_Cell::ReadData(const char* filename){
    Q_mean.ResizeTo(n_cell);
  //Q_std.ResizeTo(n_cell);
  Q_err.ResizeTo(n_cell);
  Q_mean = *(TVectorD*) f->Get("Q_mean");
//  Q_std = *(TVectorD*) f->Get("Q_std");
  Q_err = *(TVectorD*) f->Get("Q_err");

}

void Q_Cell::WriteData(const char* filename){
  Cell::WriteData(filename); //call parent method
  Q_mean.Write("Q_mean",TObject::kOverwrite);
  //Q_std.Write("Q_std",TObject::kOverwrite);
  Q_err.Write("Q_err",TObject::kOverwrite);
  
  TH2D h2D_Q("h2D_Q","MCP Mean Q",nx,-X_bound,X_bound,ny,-Y_bound,Y_bound); 
  TGraph2DErrors graph2D_Q(n_cell);
  graph2D_Q.SetTitle("MCP Mean Q");
  for(int i=0;i<n_cell;i++){
    h2D_Q.Fill(xmid[i],ymid[i],Q_mean[i]);
    graph2D_Q.SetPoint(i,xmid[i],ymid[i],Q_mean[i]);
    graph2D_Q.SetPointError(i,0,0,Q_err[i]);
  }
  h2D_Q.Write("h2D_Q",TObject::kOverwrite);
  graph2D_Q.Write("graph2D_Q",TObject::kOverwrite);
//Q_mean->SetDirectory(0); ....
}

void Q_Cell::SetPoint(int i, double Qmean, double Qerr){
  double x, y;
  GetCenter(i,x,y);
  SetQMean(i,Qmean);
  //SetQStd(i,Qstd);
  SetQErr(i,Qerr);
  //cout<<"n_cell ="<<n_cell<<"\t i= "<<i<<" x = "<<x<<" y = "<<y<<" Q_mean = "<<Qmean<<endl;

} 
/*************************************************************************************/
T0_Cell::T0_Cell(double delX, double delY, double rad, double X0, double Y0): Cell(delX, delY, rad, X0, Y0){
  ResizeVectors(n_cell);
  rel_idx = 0;
}

T0_Cell::T0_Cell(const char* filename): Cell(filename){//load from file
  OpenFile(filename,"READ"); 
//set TDirectory(0) to keep read objects in memory after closing file
  ReadData(filename); //call own method
  CloseFile();
}

T0_Cell::~T0_Cell(){
}

void T0_Cell::ResizeVectors(int n){
  T0_abs.ResizeTo(n);
  T0_abs_err.ResizeTo(n);
  T0_rel.ResizeTo(n);
}
void T0_Cell::SaveToFile(const char* filename){//save to file
  SetRelT0Corr();
  OpenFile(filename,"UPDATE");
  WriteData(filename);
  CloseFile();
}
void T0_Cell::ReadData(const char* filename){
  ResizeVectors(n_cell);
  TVectorD* T0_abs_ptr = NULL;
  T0_abs_ptr = (TVectorD*) f->Get("T0_abs");
  if (T0_abs_ptr ==NULL){
    cout<<"WARNING: File "<<filename<<" does not contain class T0_Cell data."<<endl;
    return;
  }
  T0_abs = *(TVectorD*) f->Get("T0_abs");
  T0_abs_err = *(TVectorD*) f->Get("T0_abs_err");
  T0_rel = *(TVectorD*) f->Get("T0_rel");
  
  //Define one element vectors for reading 
  TVectorD rel_idx_t(1); 
  rel_idx_t = *(TVectorD*) f->Get("rel_idx");  
  rel_idx = rel_idx_t[0];
}

void T0_Cell::WriteData(const char* filename){
  Cell::WriteData(filename); //call parent method
  T0_abs.Write("T0_abs",TObject::kOverwrite);
  T0_abs_err.Write("T0_abs_err",TObject::kOverwrite);
  T0_rel.Write("T0_rel",TObject::kOverwrite);
  
  //Define one element vectors for writing
  TVectorD rel_idx_t(1); 
  rel_idx_t[0] = rel_idx;
  rel_idx_t.Write("rel_idx",TObject::kOverwrite);
  
  
  TH2D h2D_T0rel("h2D_T0rel","Relative T0 Correction",nx,-X_bound,X_bound,ny,-Y_bound,Y_bound); 
  TGraph2DErrors graph2D_T0abs(n_cell);
  graph2D_T0abs.SetTitle("Absolute T0");
  for(int i=0;i<n_cell;i++){
    h2D_T0rel.Fill(xmid[i],ymid[i],T0_rel[i]);
    graph2D_T0abs.SetPoint(i,xmid[i],ymid[i],T0_abs[i]);
    graph2D_T0abs.SetPointError(i,0,0,T0_abs_err[i]);
  }
  h2D_T0rel.Write("h2D_T0rel",TObject::kOverwrite);
  graph2D_T0abs.Write("graph2D_T0abs",TObject::kOverwrite);
//Q_mean->SetDirectory(0); ....
}
void T0_Cell::PrintClass(){
  T0_abs.Print();
  xmid.Print();
  
}
void T0_Cell::SetPoint(int i, double Qmean, double Qerr){
  double x, y;
  GetCenter(i,x,y);
  //SetT0_abs(i,Qmean);
  //SetT0_abs_err(i,Qerr);
  //cout<<"n_cell ="<<n_cell<<"\t i= "<<i<<" x = "<<x<<" y = "<<y<<" Q_mean = "<<Qmean<<endl;
}  
void T0_Cell::SetRelT0Corr(){
  double t0 = GetT0_abs(rel_idx);
  T0_rel = T0_abs;
  T0_rel -= t0;
}