//-----------------------------------------------------------
// 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"
#include "TMath.h"

ClassImp(TpcEFieldCyl)

// Class Member definitions -----------
TpcEFieldCyl::TpcEFieldCyl()
: TpcFieldCylGrid<TVector3>(TVector3(0.,0.,0.),
		0.,0.,0.,0.)
		{
	iter=0;
	verbose=kFALSE;
		}

TpcEFieldCyl::TpcEFieldCyl(const char* const fileName)
: TpcFieldCylGrid<TVector3>(TVector3(0.,0.,0.),
		0.,0.,0.,0.)
		{
	iter=0;
	verbose=kFALSE;

	loader = new TpcEFieldCylLoader(this, fileName);
	loader->load();
	if (fpGrid == 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)
	{
		//      setiter(0);
		return TpcFieldCylGrid<TVector3>::value(point);
	}
	else //return value for boundary point next to point
	{
		TString error="Point outside of tpc volume!";
		error+=" X:";
		error+=point.X();
		error+=" Y:";
		error+=point.Y();
		error+=" Z:";
		error+=point.Z();
		error+=" R:";
		error+=TMath::Sqrt(point.X()*point.X()+point.Y()*point.Y());
		//    error+=" Iteration: ";
		//    error+=iter;
		//    Error("TpcEFieldCyl::value()", "Point outside of tpc volume!");
		if (verbose)
			Error("TpcEFieldCyl::value()", error.Data());
		TVector3 boundPoint = point - frelPosition;
		double r2 = boundPoint.X()*boundPoint.X()+boundPoint.Y()*boundPoint.Y();
		if (verbose)
			std::cout<<"r2: "<<r2<<" boundZ: "<<boundPoint.Z()<<std::endl;
		//    frelPosition.Print();
		//std::cout<<std::endl;
		if (r2 >= maxR()*maxR())
		{
			if (verbose)
				std::cout<<"too big"<<std::endl;
			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())
		{
			if (verbose)
				std::cout<<"too small"<<std::endl;
			boundPoint.SetX(boundPoint.X()*minR()/TMath::Sqrt(r2-1e-10));
			boundPoint.SetY(boundPoint.Y()*minR()/TMath::Sqrt(r2-1e-10));
			if (verbose)
				std::cout<<"new X: "<<boundPoint.X()<<" new Y: "<<boundPoint.Y()<<std::endl;
		}
		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
		//    this->setiter(iter+1);
		return value(boundPoint + frelPosition);
	}
}