//-----------------------------------------------------------
//
// Description:
//      Implementation of class TpcDevmapCyl
//      see TpcDevmapCyl.h for details
//
// Environment:
//      Software developed for the PANDA Detector at FAIR.
//
// Author List:
//      Cristoforo Simonetto    TUM         (original author)
//      Felix Boehmer           TUM
//
//-----------------------------------------------------------

#include <iostream>
#include "TMath.h"
// This Class' Header ------------------
#include "TpcDevmapCyl.h"


// Class Member definitions -----------
//ClassImp(TpcDevmapCylWrapper)

TpcDevmapCyl::TpcDevmapCyl(const char* const fileName, double vdrift)
  :TpcFieldCylGrid<TVector3>(TVector3(0.,0.,0.),
			     0.,0.,0.,0.), floaded(false)
{

  //TODO: redesign, get vDrift from param management
  //TODO: Remove useless parameter
  fvDrift = vdrift;
  loader = new TpcDevmapCylLoader(this, fileName);
  if(loader->load()==0)floaded=true;
  if (!floaded)
    Error("TpcDevmapCyl::TpcDevmapCyl()",
	  "Failed to load field data - check deviation-map file!");
  evalMaxPoint();
  _lengthR=loader->getLengthR();
  _lengthZ=loader->getLengthZ();
  loader->getValues(_rMax,_rMin,_zMax,_zMin,_rWidth,_zWidth);
}

TpcDevmapCyl::TpcDevmapCyl(const char* const fileName)
  :TpcFieldCylGrid<TVector3>(TVector3(0.,0.,0.),
			     0.,0.,0.,0.), floaded(false)
{

  loader = new TpcDevmapCylLoader(this, fileName);
  fvDrift=0;
  if(loader->load()==0)floaded=true;
  if (!floaded)
    Error("TpcDevmapCyl::TpcDevmapCyl()",
	  "Failed to load field data - check deviation-map file!");
  evalMaxPoint();
  _lengthR=loader->getLengthR();
  _lengthZ=loader->getLengthZ();
  loader->getValues(_rMax,_rMin,_zMax,_zMin,_rWidth,_zWidth);

}

TpcDevmapCyl::TpcDevmapCyl(const char* const fileName, double vdrift, bool inv)
  :TpcFieldCylGrid<TVector3>(TVector3(0.,0.,0.),
			     0.,0.,0.,0.),floaded(false)
{
  //TODO: redesign, get vDrift from param management
  //TODO: Remove useless parameter
  fvDrift = vdrift;
  loader = new TpcDevmapCylLoader(this,fileName);
  if (inv)
    {
      if (loader->loadinv()==0) floaded=true;
      if (!floaded)
	Error("TpcDevmapCyl()","Failed to load inverse field data - check deviation-map file!");
      evalMaxPoint();
    }
  else
    {
      //      loader = new TpcDevmapCylLoader(this, fileName);
      if(loader->load()==0)floaded=true;
      if (!floaded)
	Error("TpcDevmapCyl::TpcDevmapCyl()",
	      "Failed to load field data - check deviation-map file!");
      evalMaxPoint(); 
    }
  _lengthR=loader->getLengthR();
  _lengthZ=loader->getLengthZ();
  loader->getValues(_rMax,_rMin,_zMax,_zMin,_rWidth,_zWidth);

}

TpcDevmapCyl::TpcDevmapCyl(const char* const fileName, bool inv)
  :TpcFieldCylGrid<TVector3>(TVector3(0.,0.,0.),
			     0.,0.,0.,0.),floaded(false)
{
  loader = new TpcDevmapCylLoader(this,fileName);
  fvDrift=0;
  if (inv)
    {
      if (loader->loadinv()==0) floaded=true;
      if (!floaded)
	Error("TpcDevmapCyl()","Failed to load inverse field data - check deviation-map file!");
      evalMaxPoint();
    }
  else
    {
      //      loader = new TpcDevmapCylLoader(this, fileName);
      if(loader->load()==0)floaded=true;
      if (!floaded)
	Error("TpcDevmapCyl::TpcDevmapCyl()",
	      "Failed to load field data - check deviation-map file!");
      evalMaxPoint(); 
    }
  _lengthR=loader->getLengthR();
  _lengthZ=loader->getLengthZ();
  loader->getValues(_rMax,_rMin,_zMax,_zMin,_rWidth,_zWidth,fvDrift);

}

TpcDevmapCyl::TpcDevmapCyl() 
  :TpcFieldCylGrid<TVector3>(TVector3(0.,0.,0.),
			     0.,0.,0.,0.), floaded(false)
{;}

void TpcDevmapCyl::load(const char* const fileName)
{
  loader = new TpcDevmapCylLoader(this, fileName);
  if(loader->load()==0)floaded=true;
  if (!floaded)
    Error("TpcDevmapCyl::TpcDevmapCyl()",
	  "Failed to load field data - check deviation-map file!");
  evalMaxPoint();
  _lengthR=loader->getLengthR();
  _lengthZ=loader->getLengthZ();
  loader->getValues(_rMax,_rMin,_zMax,_zMin,_rWidth,_zWidth);

}

void TpcDevmapCyl::loadinv(const char* const fileName)
{
  loader = new TpcDevmapCylLoader(this,fileName);
  if (loader->loadinv()==0) floaded=true;
  if (!floaded)
    Error("TpcDevmapCyl()","Failed to load inverse field data - check deviation-map file!");
  evalMaxPoint();
  _lengthR=loader->getLengthR();
  _lengthZ=loader->getLengthZ();
  loader->getValues(_rMax,_rMin,_zMax,_zMin,_rWidth,_zWidth);
}

TpcDevmapCyl::~TpcDevmapCyl()
{delete loader;}

void
TpcDevmapCyl::print(std::ostream& s) const
{
  s << "TpcDevmapCyl\n";
  TpcFieldCylGrid<TVector3>::print(s); 
}

void TpcDevmapCyl::printc()
{
  TpcFieldCylGrid<TVector3>::print(std::cout); 
}

TVector3
TpcDevmapCyl::value(const TVector3& point) const
{
  if(!floaded)return TVector3(0,0,0);
  //std::cout<<"TpcDevmapCyl::value("<<point.X()<<","<<point.Y()<<","
  //	   <<point.Z()<<")"<<std::endl;
  if (pointOk(point) == true)
  {
    return TpcFieldCylGrid<TVector3>::value(point);    
  }
  else //return value for boundary point next to point 
  {
    //Error("TpcDevmapCyl::value()", "Point outside of tpc volume!");
    TVector3 boundPoint = point - frelPosition;
    double r2 = boundPoint.X()*boundPoint.X()+boundPoint.Y()*boundPoint.Y();
    if (r2 >= maxR()*maxR())
    {
      boundPoint.SetX(boundPoint.X()*maxR()/TMath::Sqrt(r2+1e-10));//1e-10 guarantees, that the point is inside volume
      boundPoint.SetY(boundPoint.Y()*maxR()/TMath::Sqrt(r2+1e-10));
    }
    else if(r2 < minR()*minR())
    {
      boundPoint.SetX(boundPoint.X()*minR()/sqrt(r2-1e-10));
      boundPoint.SetY(boundPoint.Y()*minR()/sqrt(r2-1e-10));
    }
    if (boundPoint.Z() >= maxZ())
      boundPoint.SetZ(maxZ()-1e-10);
    else if (boundPoint.Z() < minZ())
      boundPoint.SetZ(minZ());
    //recursive: guarantees, that the point is inside volume
    return value(boundPoint + frelPosition);
  }
}

TVector3 TpcDevmapCyl::value2(TVector3 point)
{
  return value(point);
}

 TpcDevmapCylWrapper::TpcDevmapCylWrapper(const char* file) 
   :fLoaded(false),
    fDevmap(NULL)
 {
   fDevmap = new TpcDevmapCyl(file);
   if(fDevmap->loaded()) 
     fLoaded = true;
 }

TpcDevmapCylWrapper::TpcDevmapCylWrapper(const char* file, double vdrift, bool inv)
  :fLoaded(false),
   fDevmap(NULL)
{
  fDevmap = new TpcDevmapCyl(file,vdrift,inv);
  if (fDevmap->loaded())
    fLoaded = true;
}

TpcDevmapCylWrapper::TpcDevmapCylWrapper(const char* file, double vdrift)
  :fLoaded(false),
   fDevmap(NULL)
{
  fDevmap = new TpcDevmapCyl(file,vdrift);
  if (fDevmap->loaded())
    fLoaded = true;
}

TpcDevmapCylWrapper::TpcDevmapCylWrapper(const char* file, bool inv)
  :fLoaded(false),
   fDevmap(NULL)
{
  fDevmap = new TpcDevmapCyl(file,inv);
  if (fDevmap->loaded())
    fLoaded = true;
}


 TpcDevmapCylWrapper::~TpcDevmapCylWrapper() 
 {
   if(fDevmap)
     delete fDevmap;
 }

 TVector3
 TpcDevmapCylWrapper::value(const TVector3& vec) const {
   if(fLoaded)
     return fDevmap->value(vec);
   else
     return TVector3(0.,0.,0.);
  
 }