//-----------------------------------------------------------
//
// 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>

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


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


TpcDevmapCyl::TpcDevmapCyl(const char* const fileName, double vDrift)
  :TpcFieldCylGrid<TVector3>(TVector3(0.,0.,0.),
			    0.,0.,0.,0.)
{
  _vDrift = vDrift;
  loader = new TpcDevmapCylLoader(this, fileName);
  loader->load();
  if (_pGrid == 0)
    Fatal("TpcDevmapCyl::TpcDevmapCyl()",
	  "Failed to load field data - check deviation-map file!");
  evalMaxPoint();
}

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

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

TVector3
TpcDevmapCyl::value(const TVector3& point) const
{
  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 - _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);
  }
}