#ifndef __JLINE3Z__
#define __JLINE3Z__

#include <istream>
#include <ostream>

#include "JMath/JMath.hh"
#include "JPhysics/JConstants.hh"
#include "JGeometry3D/JVector3D.hh"
#include "JGeometry3D/JVersor3Z.hh"
#include "JGeometry3D/JPosition3D.hh"
#include "JGeometry3D/JDirection3D.hh"
#include "JGeometry3D/JTrack3D.hh"
#include "JFit/JLine1Z.hh"


/**
 * \author mdejong
 */

namespace JFIT {}
namespace JPP { using namespace JFIT; }

namespace JFIT {

  using JMATH::JMath;
  using JGEOMETRY3D::JVector3D;
  using JGEOMETRY3D::JVersor3Z;
  using JGEOMETRY3D::JTrack3D;
    

  /**
   * Data structure for fit of straight line in positive z-direction.
   * Note that the position coordinates are defined with respect to the given direction.
   */
  class JLine3Z : 
    public JLine1Z,
    public JVersor3Z,
    public JMath<JLine3Z>
  {
  public:

    using JLine1Z::add;
    using JLine1Z::sub;
    using JLine1Z::getT;
    using JVersor3Z::add;
    using JVersor3Z::sub;
    using JVersor3Z::getDirection;


    /**
     * Default constructor.
     */
    JLine3Z() :
      JLine1Z(),
      JVersor3Z()
    {}


    /**
     * Constructor.
     *
     * \param  line             line
     */
    JLine3Z(const JLine1Z& line) :
      JLine1Z(line),
      JVersor3Z()
    {}


    /**
     * Constructor.
     *
     * \param  line             line
     * \param  dir              direction
     */
    JLine3Z(const JLine1Z&   line,
	    const JVersor3Z& dir) :
      JLine1Z(line),
      JVersor3Z(dir)
    {}


    /**
     * Move vertex along this line with given velocity.
     *
     * \param  step             step
     * \param  velocity         velocity
     */
    void move(const double step, const double velocity)
    {
      __x += step * getDX();
      __y += step * getDY();
      __z += step * getDZ();
      __t += step / velocity;
    }


    /**
     * Set z-position.
     *
     * \param  z                z-position
     * \param  velocity         velocity
     */
    void setZ(const double z, const double velocity)
    {
      move((z - getZ()) / getDZ(), velocity);
    }


    /**
     * Prefix unary minus.
     *
     * \return                  line
     */
    JLine3Z& negate()
    {
      JLine1Z  ::negate();
      JVersor3Z::negate();

      return *this;
    }


    /**
     * Addition operator.
     *
     * \param  value            line
     * \return                  line
     */
    JLine3Z& add(const JLine3Z& value)
    {
      JLine1Z  ::add(value);
      JVersor3Z::add(value);

      return *this;
    }


    /**
     * Subtraction operator.
     *
     * \param  value            line
     * \return                  line
     */
    JLine3Z& sub(const JLine3Z& value)
    {
      JLine1Z  ::sub(value);
      JVersor3Z::sub(value);

      return *this;
    }


    /**
     * Multiplication operator.
     *
     * \param  value            multiplication factor
     * \return                  line
     */
    JLine3Z& mul(const double value)
    {
      JLine1Z  ::mul(value);
      JVersor3Z::mul(value);

      return *this;
    }


    /**
     * Division operator.
     *
     * \param  value            multiplication factor
     * \return                  line
     */
    JLine3Z& div(const double value)
    {
      JLine1Z  ::div(value);
      JVersor3Z::div(value);

      return *this;
    }


    /**
     * Type conversion.
     *
     * \return                  track
     */
    inline operator JTrack3D() const
    {
      return JTrack3D(JAxis3D(this->getPosition(),
			      this->getDirection()),
		      this->getT());
    }


    /**
     * Get distance squared.
     *
     * \param  pos              position
     * \return                  square of distance
     */
    inline double getDistanceSquared(const JVector3D& pos) const
    {
      JPosition3D D(pos);

      D.sub(this->getPosition());

      const double dz = D.getDot(this->getDirection());

      return D.getLengthSquared() - dz*dz;
    }


    /**
     * Get distance.
     *
     * \param  pos              position
     * \return                  distance
     */
    inline double getDistance(const JVector3D& pos) const
    {
      return sqrt(this->getDistanceSquared(pos));
    }


    /**
     * Get arrival time of Cherenkov light at given position.
     *
     * \param  pos              position [m]
     * \return                  time     [ns]
     */
    inline double getT(const JVector3D& pos) const
    {
      using namespace JPP;

      JPosition3D D(pos);

      D.sub(this->getPosition());

      const double dz = D.getDot(this->getDirection());
      const double R  = sqrt(D.getLengthSquared() - dz*dz);
      
      return this->getT() + (dz + R * getKappaC()) * getInverseSpeedOfLight();
    }


    /**
     * Get photon direction of Cherenkov light on PMT.
     *
     * \param  pos              PMT position
     * \return                  direction
     */
    inline JVersor3D getDirection(const JVector3D& pos) const
    {
      using namespace JPP;

      JPosition3D D(pos);

      D.sub(this->getPosition());

      const double dz = D.getDot(this->getDirection());
      const double R  = sqrt(D.getLengthSquared() - dz*dz);

      D.sub(JPosition3D(this->getDirection()) * (dz - R/getTanThetaC()));
      
      return JDirection3D(D);
    }


    /**
     * Get cosine angle of impact of Cherenkov light on PMT.
     *
     * \param  axis             PMT axis
     * \return                  cosine angle of impact
     */
    inline double getDot(const JAxis3D& axis) const
    {
      return getDirection(axis.getPosition()).getDot(axis.getDirection());
    }


    /**
     * Read object from input.
     *
     * \param  in            input stream
     * \param  object        object
     * \return               input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JLine3Z& object)
    {
      in >> static_cast<JLine1Z&>  (object);
      in >> static_cast<JVersor3Z&>(object);

      return in;
    }


    /**
     * Write object to output.
     *
     * \param  out           output stream
     * \param  object        object
     * \return               output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const JLine3Z& object)
    {
      out << static_cast<const JLine1Z&>  (object);
      out << ' ';
      out << static_cast<const JVersor3Z&>(object);

      return out;
    }


    typedef double JLine3Z::*parameter_type;

    static parameter_type pDX() { return &JLine3Z::__dx; }
    static parameter_type pDY() { return &JLine3Z::__dy; }
  };
}

#endif