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



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



// Class Member definitions -----------

TpcEFieldCyl::TpcEFieldCyl(const char* const fileName)
  : TpcFieldCylGrid<TVector3>(TVector3(0.,0.,0.),
			    0.,0.,0.,0.)
{
  loader = new TpcEFieldCylLoader(this, fileName);
  loader->load();
  if (_pGrid == 0)
    Fatal("TpcEFieldCyl::TpcEFieldCyl()",
	  "Failed to load field data - check E-field file!");
  evalMaxPoint();
}

TpcEFieldCyl::~TpcEFieldCyl()
{delete loader;} //TpcEFieldCyl has ownership of the fieldmap created by
                 //TpcEFieldCylLoader.


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

TVector3
TpcEFieldCyl::value(const TVector3& point) const
{
  if (pointOk(point) == true)
    return TpcFieldCylGrid<TVector3>::value(point);
  else //return value for boundary point next to point 
  {
    Error("TpcEFieldCyl::value()", "Point outside of tpc volume!");
    TVector3 boundPoint = point - _relPosition;
    double r2 = boundPoint.X()*boundPoint.X()+boundPoint.Y()*boundPoint.Y();
    if (r2 >= maxR()*maxR())
    {
      boundPoint.SetX(boundPoint.X()*maxR()/sqrt(r2+1e-10));//1e-10 guarantees, that the point is inside volume
      boundPoint.SetY(boundPoint.Y()*maxR()/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 + _relPosition);
  }
}