//-----------------------------------------------------------
// File and Version Information:
// $Id$
//
// Description:
//      Implementation of class TpcPRLookupTable
//      see TpcPRLookupTable.h for details
//
// Environment:
//      Software developed for the PANDA Detector at FAIR.
//
// Author List:
//      Sebastian Neubert    TUM            (original author)
//      Cristoforo Simonetto TUM
//
//
//-----------------------------------------------------------

// This Class' Header ------------------
#include "TpcPRLookupTable.h"

// C/C++ Headers ----------------------
#include <fstream>
#include <sstream>

// Collaborating Class Headers --------
#include "TpcAbsPadShape.h"
#include "TpcGem.h"
#include "TError.h"
#include "TMath.h"
// Class Member definitions -----------


TpcPRLookupTable::TpcPRLookupTable()
  : isbuilt(false),range(0),step(0),nx(0),ny(0),table(0)
{;}


TpcPRLookupTable::~TpcPRLookupTable()
{
  if(table!=0)
  {
    for(int i=0;i<nx;++i)
      delete [] table[i];
    delete [] table;
  }
}

double
TpcPRLookupTable::GetValue(const double dx, const double dy) const{
  if(!isbuilt)return 0;

  int offsetX=0;
  int offsetY=0;

  if( int(nx/2.)!=nx/2. )
	  offsetX=(nx-1)/2;
  else
	  offsetX=nx/2;

  if (int(ny/2.)!=ny/2. )
	  offsetY=(ny-1)/2;
  else
	  offsetY=ny/2;

  int ilx=(int)TMath::Floor(dx/step);
  int ily=(int)TMath::Floor(dy/step);
  int iux=ilx+1;
  int iuy=ily+1;

  ilx+=offsetX;
  iux+=offsetX;
  ily+=offsetY;
  iuy+=offsetY;

  //std::cout<<"offsetX:"<<offsetX<<" OffsetY:"<<offsetY<<" ilx:"<<ilx<<" ily:"<<ily<<" iux:"<<iux<<" iuy:"<<iuy<<" nx:"<<nx<<" ny:"<<ny<<" dx:"<<dx<<" dy:"<<dy<<" step:"<<step<<"\n";

  if(ilx<0 || ilx>=nx){
    //ErrMsg(warning)<<"ilx out of range"<<endmsg;
    return 0;
  }
  if(ily<0 || ily>=ny){
    //ErrMsg(warning)<<"ily out of range"<<endmsg;
    return 0;
  }
  if(iux<0 || iux>=nx){
    //ErrMsg(warning)<<"iux out of range"<<endmsg;
    return 0;
  }
  if(iuy<0 || iuy>=ny){
    //ErrMsg(warning)<<"iuy out of range"<<endmsg;
    return 0;
  }


  /*
  ilx=std::max(ilx,0);
  ilx=std::min(ilx,nx-1);
  ily=std::max(ily,0);
  ily=std::min(ily,ny-1);
  iux=std::max(iux,0);
  iux=std::min(iux,nx-1);
  iuy=std::max(iuy,0);
  iuy=std::min(iuy,ny-1);
   */

  double tab_11=0;
  double tab_12=0;
  double tab_21=0;
  double tab_22=0;


  /*
  if(ilx<0 and iux<0)
	  iux=0;
  if(ilx>=nx and iux>=nx)
	  ilx=nx-1;

  if(ily<0 and iuy<0)
	  iuy=0;
  if(ily>=ny and iuy>=ny)
	  ily=ny-1;
	*/

  if (ilx>=0 and ilx<nx)
  {
	  if(ily>=0 and ily<ny)
		  tab_11=table[ilx][ily];
	  if(iuy>=0 and iuy<ny)
		  tab_12=table[ilx][iuy];
  }
  if (iux>=0 and iux<nx)
  {
	  if(ily>=0 and ily<ny)
		  tab_21=table[iux][ily];
	  if(iuy>=0 and iuy<ny)
		  tab_22=table[iux][iuy];
  }

  double lx=(ilx-offsetX)*step;
  double ly=(ily-offsetY)*step;
  double ux=(iux-offsetX)*step;
  double uy=(iuy-offsetY)*step;

  //std::cout<<"lx="<<lx<<" ly="<<ly<<" ux="<<ux<<" uy="<<uy<<"\n";
  //std::cout<<"ilx:"<<ilx<<" ily:"<<ily<<" iux:"<<iux<<" iuy:"<<iuy<<" dx:"<<dx<<" dy:"<<dy<<" step:"<<step<<"\n";
  //std::cout<<"tab_11="<<tab_11<<" tab_12="<<tab_12<<" tab21="<<tab_21<<" tab_22="<<tab_22<<"\n";

  //bilinear interpolation:
  double t=(dx-lx)/(ux-lx);
  double u=(dy-ly)/(uy-ly);
  double result=0;
  /*
  result+=(1-t)*(1-u)*table[ilx][ily];
  result+=t*(1-u)*table[iux][ily];
  result+=t*u*table[iux][iuy];
  result+=(1-t)*u*table[ilx][iuy];
*/
  result+=(1-t)*(1-u)*tab_11;
  result+=t*(1-u)*tab_21;
  result+=t*u*tab_22;
  result+=(1-t)*u*tab_12;

  //std::cout<<"t="<<t<<" u="<<u<<"\n";
  //std::cout<<"result="<<result<<"\n";
  return result;
}

void 
TpcPRLookupTable::BuildTable(const TpcGem& gem, 
			     const TpcAbsPadShape& pad, 
			     const double frange,      
			     const double fstep,
			     const double fintStep){
  if(table!=0){
    for(int i=0;i<nx;++i){
      delete [] table[i];
    }
    delete [] table;
  }
  range=frange; step=fstep;
  double padwidth;
  double padheight;
  pad.EvalBoundingRect(padwidth, padheight, 0);
  nx=2*(int)ceil( (range+0.5*std::max(padwidth,padheight) ) /step);//# of steps in x direction
  ny=2*(int)ceil( (range+0.5*std::max(padwidth,padheight) ) /step);//# of steps in y
  //build table
  table=new double*[nx];
  for(int i=0;i<nx;++i){
    table[i]=new double[ny];
  }
  std::ostream& s = std::cerr;
  //loop over table
  s<<"Building TPC LookupTable ";
  for(int ix=0;ix<nx;++ix){
    if(ix%10 == 0) s<<".";
    for(int iy=0;iy<ny;++iy){
      double thisx=(ix-nx/2)*step;
      double thisy=(iy-ny/2)*step;
      table[ix][iy]=Integrate(gem,pad,fintStep,thisx,thisy);
    }
  }
  s<<std::endl;
  isbuilt=true;
}

void
TpcPRLookupTable::ReadTable(TString filename)
{
	std::cout<<"TpcPRLookupTable: Reading File: "<<filename<<"\n";
	ifstream input(filename.Data());
	std::ostream& s = std::cerr;
	if(!input)
	{
		std::cout<<"TpcPRLookupTable: could not open par file!\n";
		return;
	}
	std::string line;

	int ix, iy;
	double response;
	int lc=-1;
	s<<"Reading TPC LookupTable ";
	while (getline(input,line))
	{
		lc++;
		if(lc%10 == 0) s<<".";
		std::istringstream lineStream(line);
		if (lc==0)
		{
			if(!(lineStream>>nx>>ny>>step))
			{
				std::string msg("TpcPRLookupTable: error while reading file (first line)!");
				std::cout<<msg<<"\n";
			}
			table=new double*[nx];
			for(int i=0;i<nx;i++)
				table[i]=new double[ny];
			continue;
		}

		if (!(lineStream>>ix>>iy>>response))
		{
			std::string msg("TpcPRLookupTable: error while reading file!");
			std::cout<<msg<<"\n";
		}
		table[ix][iy]=response;
	}
	s<<std::endl;
	isbuilt=true;
}

double 
TpcPRLookupTable::Integrate(const TpcGem& gem, 
			    const TpcAbsPadShape& p,
			    const double intStep,
			    const double dx,      
			    const double dy) const 
{
  double padwidth;
  double padheight;
  p.EvalBoundingRect(padwidth, padheight, 0);
  int nintx=((int)floor(std::max(padwidth,padheight)/intStep));
  int ninty=((int)floor(std::max(padwidth,padheight)/intStep));
  double sum=0;
  double offsetx=0.5*(padwidth-intStep);
  double offsety=0.5*(padheight-intStep);
  //initialize Avalanche
  for(int ix=0;ix<nintx;++ix)
    for(int iy=0;iy<ninty;++iy)
      {
	if (!p.Contains(ix*intStep-offsetx,iy*intStep-offsety)) 
	  continue;
      double thisx=ix*intStep-offsetx-dx;
      double thisy=iy*intStep-offsety-dy;
      double r=sqrt(thisx*thisx+thisy*thisy);
      sum+=gem.cloudShape()->Eval(r)*1/gem.spread()*1/sqrt(2*3.1416); // TODO: cache these values!!!
    }
  
  //normalize to area
  double s2=intStep*intStep;
  sum*=s2;
  return sum;
}

double TpcPRLookupTable::GetLTValue(const int x, const int y) 
{
  if (x<nx and y<ny and x>=0 and y>=0)
    return table[x][y];
  else
    return -1;
} 

double TpcPRLookupTable::GetWeight(double amp, double dirX, double dirY)
{
  double val=1.;
  double X=0;
  double Y=0;

  //std::cout<<"start calculating weight: "<<val<<" "<<amp<<"\n";
  /*
  while (val>amp)
    {
      val=this->GetValue(X,Y);
      std::cout<<"xpos: "<<X<<" Ypos: "<<Y<<" val: "<<val<<"\n";
      X+=0.01*dirX;
      Y+=0.01*dirY;

    
    }
  */

  double xdir=dirX;
  double ydir=dirY;
  double mag=sqrt(dirX*dirX+dirY*dirY);
  /*
  if (mag<=0.15/2.)
  {
	  //std::cout<<"magnitude was smaller 0.15\n";
	  return 1;
  }
  */
  xdir=(xdir*2)/mag;
  ydir=(ydir*2)/mag;

  double scale=1;
  bool isbigger=true;
  int searchCounter=-1;
  //std::cout<<"starting fine search:"<<fabs(val-amp)<<"\n";

  if (amp==0)
	  return 0;

  if (amp==1)
	  amp=0.999;

  val=this->GetValue(X,Y);
  while (fabs((val-amp)/amp)>0.001)
  {
	  searchCounter+=1;
	  if (searchCounter>100)
	  {
		  X=0;
		  Y=0;
		  //std::cout<<"TpcPRLookupTable: weighting: Can't find point "<<amp<<";"<<searchCounter<<"\n";
		  //break;
		  return 0;
	  }
	  if (val>amp)
	  {
		  if (not isbigger)
			  scale/=2;

		  X+=scale*xdir;
		  Y+=scale*ydir;
		  //std::cout<<"bigger\n";
		  isbigger=true;
	  }
	  if (val<amp)
	  {
		  if(isbigger)
			  scale/=2.;
		  X-=scale*xdir;
		  Y-=scale*ydir;
		  //std::cout<<"smaller\n";
		  isbigger=false;
	  }
	  val=this->GetValue(X,Y);


	  /*
	  if (sqrt(X*X+Y*Y)/mag > 2)
	  {
		  std::cout<<"xpos: "<<X<<" Ypos: "<<Y<<" val: "<<val<<" scale: "<<scale<<" amp="<<amp<<" mag="<<mag<<"\n";
		  std::cout<<"retval would be:"<<sqrt(X*X+Y*Y)/mag<<"\n";
	  }*/
  }
  double retval=sqrt(X*X+Y*Y)/mag;


  	  /*
  if (fabs(retval)>1)
  {
	  //std::cout<<"distance was bigger than 1 "<<retval<<"\n";
	  //return 1.+(retval-1.)/10.;
	  return retval;

  }
	*/

  //std::cout<<"Lookup mag="<<mag<<" xdir="<<xdir<<" ydir="<<ydir<<"\n";
  //std::cout<<"Retval:"<<sqrt(X*X+Y*Y)<<"\n";
  //if (retval!=0)
  //    return 1./retval;
  //  else
  //return 1.;
  return retval;// sqrt(X*X+Y*Y);
}