#ifndef __JTRACK3D__
#define __JTRACK3D__

#include <istream>
#include <ostream>

#include "JPhysics/JConstants.hh"
#include "JGeometry3D/JAxis3D.hh"
#include "JGeometry3D/JTime.hh"
#include "JGeometry3D/JVertex3D.hh"
#include "JIO/JSerialisable.hh"


/**
 * \author mdejong
 */

namespace JGEOMETRY3D {}
namespace JPP { using namespace JGEOMETRY3D; }

namespace JGEOMETRY3D {

  using JIO::JReader;
  using JIO::JWriter;


  /**
   * 3D track.
   */
  class JTrack3D :
    public JAxis3D,
    public JTime
  {
  public:

    using JTime::getT;
    using JTime::add;
    using JTime::sub;
    using JAxis3D::add;
    using JAxis3D::sub;
    using JAxis3D::getDirection;


    /**
     * Default constructor.
     */
    JTrack3D() :
      JAxis3D(),
      JTime  ()
    {}


    /**
     * Constructor.
     *
     * \param  axis          track axis
     * \param  t0            time at vertex
     */
    JTrack3D(const JAxis3D& axis,
	     const JTime&   t0) :
      JAxis3D(axis),
      JTime  (t0)
    {}


    /**
     * Constructor.
     *
     * \param  pos           position
     * \param  dir           direction
     * \param  t0            time at vertex
     */
    JTrack3D(const JVector3D& pos,
	     const JVersor3D& dir,
	     const JTime&     t0) :
      JAxis3D(pos, dir),
      JTime  (t0)
    {}


    /**
     * Constructor.
     *
     * \param  vertex        vertex
     * \param  dir           direction
     */
    JTrack3D(const JVertex3D& vertex,
	     const JVersor3D& dir) :
      JAxis3D(vertex.getPosition(), dir),
      JTime  (vertex.getT())
    {}


    /**
     * Get vertex of this track.
     *
     * \return               vertex
     */
    JVertex3D getVertex() const
    {
      return JVertex3D(getPosition(), getT());
    }

    

    /**
     * Move vertex along this track with given velocity.
     *
     * \param  step          step
     * \param  velocity      velocity
     */
    void move(const double step, const double velocity)
    {
      JAxis3D::move(step);

      __t += step / velocity;
    }


    /**
     * 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 JPHYSICS;

      JPosition3D D(pos);

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

      const double dz = D.getDot(this->getDirection());
      const double rt = sqrt(D.getLengthSquared() - dz*dz);
      
      return this->getT() + (dz + rt * 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 JPHYSICS;

      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());
    }


    /**
     * Get cosine angle of emission of Cherenkov light on PMT.
     *
     * \param  pos              PMT position
     * \return                  cosine angle of impact
     */
    inline double getDot(const JVector3D& pos) const
    {
      return JDirection3D(pos - getPosition()).getDot(getDirection());
    }


    /**
     * Read track from input.
     *
     * \param  in            input stream
     * \param  track         track
     * \return               input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JTrack3D& track)
    {
      in >> static_cast<JAxis3D&>(track);
      in >> static_cast<JTime&>  (track);

      return in;
    }


    /**
     * Write track to output.
     *
     * \param  out           output stream
     * \param  track         track
     * \return               output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const JTrack3D& track)
    {
      out << static_cast<const JAxis3D&>(track);
      out << ' ';
      out << static_cast<const JTime&>  (track);

      return out;
    }


    /**
     * Read track from input.
     *
     * \param  in            reader
     * \param  track         track
     * \return               reader
     */
    friend inline JReader& operator>>(JReader& in, JTrack3D& track)
    {
      in >> static_cast<JAxis3D&>(track);
      in >> static_cast<JTime&>  (track);

      return in;
    }


    /**
     * Write track to output.
     *
     * \param  out           writer
     * \param  track         track
     * \return               writer
     */
    friend inline JWriter& operator<<(JWriter& out, const JTrack3D& track)
    {
      out << static_cast<const JAxis3D&>(track);
      out << static_cast<const JTime&>  (track);

      return out;
    }
  };
}

#endif