#ifndef __GEOMETRY__
#define __GEOMETRY__

#include <iostream>
#include <vector>
#include <map>
#include <stdexcept>

#include "nr/prob.hh"
#include "trigger/util.hh"


/**
 * \file
 * Geometry package.
 */

namespace GEOMETRY {

  using namespace UTIL;


  /**
   * 3D space point with time information.
   */
  class point : public position {
  public:
    /**
     * Default constructor
     */
    point() :
      position(),
      t_(0)
    {}

    /**
     * Constructor
     * \param x      starting x position [m]
     * \param y      starting y position [m]
     * \param z      starting z position [m]
     * \param t      time at position [ns]
     */
    point(const t_pos&  x,
	  const t_pos&  y,
	  const t_pos&  z,
	  const t_time& t) :
      position(x,y,z),
      t_(t)
    {}

    /**
     * Constructor
     * \param pos    starting position [m]
     * \param t0     time at starting position [ns]
     */
    point(const position& pos, const t_time& t) :
      position(pos),
      t_(t)
    {}

    /**
     * Constructor
     * \param in     input stream
     */
    point(std::istream& in) :
      position(),
      t_(0)
    {
      read(in);
    }

    /**
     * Rotate point
     * \param  R  rotation matrix
     * \return    rotated point
     */
    point& rotate(const rotation& R)
    {
      position::rotate(R);
      return *this;
    }

    /**
     * Rotate back point
     * \param  R  rotation matrix
     * \return    rotated point
     */
    point& rotate_back(const rotation& R)
    {
      position::rotate_back(R);
      return *this;
    }

    /** assignment operator */
    point& operator=(const double& c) 
    { 
      vector3d::operator=(c);
      t_ = c;
      return *this;
    }

    /** add operator */
    point& operator+=(const position& o) 
    { 
      vector3d::operator+=(o);
      return *this;
    }

    /** subtract operator */
    point& operator-=(const position& o) 
    { 
      vector3d::operator-=(o);
      return *this;
    }

    /** add operator */
    point& operator+=(const point& o) 
    { 
      vector3d::operator+=(o);
      t_ += o.t();
      return *this;
    }

    /** subtract operator */
    point& operator-=(const point& o) 
    { 
      vector3d::operator-=(o);
      t_ -= o.t();
      return *this;
    }

    /** scale operator */
    point& operator*=(const double& c) 
    { 
      vector3d::operator*=(c);
      t_ *= c;
      return *this;
    }

    /**
     * Get time at starting position
     * \return     time at starting position [ns]
     */
    const t_time& t() const { return t_; }

    /**
     * read from input stream
     * \param  in  input stream
     * \return     input stream
     */
    std::istream& read(std::istream& in)
    {
      position::read(in);
      in >> t_;
      return in;
    }

    /**
     * write to output stream
     * \param  out output stream
     * \return     output stream
     */
    std::ostream& write(std::ostream& out) const
    {
      position::write(out);
      out << ' ';
      out << t_;
      return out;
    }
    
  protected:
    t_time t_;
  };


  /**
   * Hit
   */
  class hit : public point {
  public:
    /**
     * Default constructor
     */
    hit() :
      point(),
      npe_(0),
      sigma_(0)
    {}

    /**
     * Constructor
     * \param  x      x position [m]
     * \param  y      y position [m]
     * \param  z      z position [m]
     * \param  t      time of hit [ns]
     * \param  npe    charge of hit [pe]
     * \param  sigma  time resolution [ns]
     */
    hit(const double& x,
	const double& y,
	const double& z,
	const double& t,
	const double& npe,
	const double& sigma) :
      point(x, y, z, t),
      npe_(npe),
      sigma_(sigma)
    {}

    /**
     * Constructor
     * \param pos    starting position [m]
     * \param t      time at position [ns]
     * \param npe    charge of hit [pe]
     * \param sigma  error on time [ns]
     */
    hit(const position& pos, const t_time& t, const double& npe, const t_time& sigma) :
      point(pos,t),
      npe_(npe),
      sigma_(sigma)
    {}

    /**
     * Constructor
     * \param in     input stream
     */
    hit(std::istream& in) :
      point(),
      npe_(0),
      sigma_(0)
    {
      read(in);
    }

    /**
     * Get charge of hit
     * \return     charge of hit [pe]
     */
    const double& npe() const { return npe_; }

    /**
     * Get error on time
     * \return     error on time [ns]
     */
    const t_time& sigma() const { return sigma_; }

    /**
     * read from input stream
     * \param  in  input stream
     * \return     input stream
     */
    std::istream& read(std::istream& in)
    {
      point::read(in);
      in >> npe_ >> sigma_;
      return in;
    }

    /**
     * write to output stream
     * \param  out output stream
     * \return     output stream
     */
    std::ostream& write(std::ostream& out) const
    {
      point::write(out);
      out << ' ';
      out << npe_;
      out << ' ';
      out << sigma_;
      return out;
    }

  protected:
    double npe_;
    t_time sigma_;
  };


  /**
   * 3D axis
   */
  class axis : public direction {
  public:
    /**
     * Default constructor
     */
    axis() :
      direction(),
      a_(0),
      b_(0)
    {}

    /**
     * Constructor
     * \param  a       transverse position theta^
     * \param  b       transverse position phi^
     * \param  theta   polar angle [rad]
     * \param  phi     azimuth angle [rad]
     */
    axis(const t_pos& a,
	 const t_pos& b,
	 const t_dir& theta,
	 const t_dir& phi) :
      direction(theta,phi),
      a_(a),
      b_(b)
    {}

    /**
     * Constructor
     * \param  dir     direction
     */
    axis(const direction& dir) :
      direction(dir)
    {
      a_ = 0;
      b_ = 0;
    }

    /**
     * Constructor
     * \param  pos     3D position
     * \param  dir     direction
     */
    axis(const position&  pos,
	 const direction& dir) :
      direction(dir)
    {
      position o(pos);
      // BvR
      /*
      rotation R(dir);
      
      o.rotate(R);
      */
      // BvR

      a_ = o.x();
      b_ = o.y();
    }

    /**
     * Rotate axis
     * \param  R  rotation matrix
     * \return    rotated axis
     */
    axis& rotate(const rotation& R)
    {
      direction::rotate(R);
      //a_.rotate(R);
      //b_.rotate(R);
      return *this;
    }

    /**
     * Rotate back axis
     * \param  R  rotation matrix
     * \return    rotated axis
     */
    axis& rotate_back(const rotation& R)
    {
      direction::rotate_back(R);
      //a_.rotate_back(R);
      //b_.rotate_back(R);
      return *this;
    }

    /** 
     * assignment operator
     * \param  pos     3D position
     */
    axis& operator=(const position& pos)
    {
      position o(pos);
      rotation R(*this);
      
      o.rotate(R);

      a_ = o.x();
      b_ = o.y();

      return *this;
    }

    /** assignment operator */
    axis& operator=(const double& c) 
    { 
      a_     = c;
      b_     = c;
      theta_ = c;
      phi_   = c;
      return *this;
    }
 
    /** offset operator */
    axis& operator+=(position o)
    {
      rotation R(*this);
      o.rotate(R);
      a_ += o.x();
      b_ += o.y();
      return *this;
    }

    /** offset operator */
    axis& operator-=(position o)
    {
      rotation R(*this);
      o.rotate(R);
      a_ -= o.x();
      b_ -= o.y();
      return *this;
    }

    /** add operator */
    axis& operator+=(const axis& o) 
    { 
      theta_ += o.theta();
      phi_   += o.phi();
      a_     += o.a();
      b_     += o.b();
      return *this;
    }

    /** subtract operator */
    axis& operator-=(const axis& o) 
    { 
      theta_ -= o.theta();
      phi_   -= o.phi();
      a_     -= o.a();
      b_     -= o.b();
      return *this;
    }

    /** scale operator */
    axis& operator*=(const double& c) 
    { 
      theta_ *= c;
      phi_   *= c;
      a_     *= c;
      b_     *= c;
      return *this;
    }

    const t_pos& a() const { return a_; }    //!< transverse position theta^
    const t_pos& b() const { return b_; }    //!< transverse position phi^

    t_pos x() const { return (t_pos) ( a_*cos(theta_)*cos(phi_) - b_*sin(phi_)); }  //!< X-position
    t_pos y() const { return (t_pos) ( a_*cos(theta_)*sin(phi_) + b_*cos(phi_)); }  //!< Y-position
    t_pos z() const { return (t_pos) (-a_*sin(theta_)); }                           //!< Z-position

    /**
     * angle between directions
     * \param  a direction
     * \param  b direction
     * \return   angle [deg]
     */
    double angle(direction a, direction b)
    {
      rotation R(*this);

      a.rotate(R);
      b.rotate(R);

      return (a.theta() - b.theta())*180/PI;
    }

    /**
     * read from input stream
     * \param  in  input stream
     * \return     input stream
     */
    std::istream& read(std::istream& in)
    {
      in >> a_ >> b_ >> theta_ >> phi_;
      return in;
    }

    /**
     * write to output stream
     * \param  out output stream
     * \return     output stream
     */
    std::ostream& write(std::ostream& out) const
    {
      out << a_ << ' ' << b_ << ' ' << theta_ << ' ' << phi_;
      return out;
    }

  protected:
    t_pos a_;    //!< transverse position theta^
    t_pos b_;    //!< transverse position phi^
  };


  /**
   * 3D plane
   */
  class plane : public direction {
  public:
    /**
     * Default constructor
     */
    plane() :
      direction(),
      c_(0)
    {}

    /**
     * Constructor
     * \param  c       longitudinal position
     * \param  theta   polar angle [rad]
     * \param  phi     azimuth angle [rad]
     */
    plane(const t_pos& c,
	  const t_dir& theta,
	  const t_dir& phi) :
      direction(theta,phi),
      c_(c)
    {}

    /**
     * Constructor
     * \param  pos     3D position
     * \param  dir     direction
     */
    plane(const position&  pos,
	  const direction& dir) :
      direction(dir)
    {
      position o(pos);
      rotation R(*this);
      
      o.rotate(R);

      c_ = o.z();
    }

    const t_pos& c() const { return c_; }    //!< longitudinal position

    t_pos x() const { return (t_pos) (c_*dx()); }  //!< X-position
    t_pos y() const { return (t_pos) (c_*dy()); }  //!< Y-position
    t_pos z() const { return (t_pos) (c_*dz()); }  //!< Z-position

    /**
     * angle between directions
     * \param  a direction
     * \param  b direction
     * \return   angle [deg]
     */
    double angle(direction a, direction b)
    {
      rotation R(*this);

      a.rotate(R);
      b.rotate(R);

      double dot_ = std::max(dot(a,b), dot(a,direction(PI - b.theta(), b.phi())));

      return acos(dot_)*180/PI;
    }

    /**
     * Rotate plane
     * \param  R  rotation matrix
     * \return    rotated plane
     */
    plane& rotate(const rotation& R)
    {
      direction::rotate(R);
      return *this;
    }

    /**
     * Rotate back plane
     * \param  R  rotation matrix
     * \return    rotated plane
     */
    plane& rotate_back(const rotation& R)
    {
      direction::rotate_back(R);
      return *this;
    }

    /**
     * read from input stream
     * \param  in  input stream
     * \return     input stream
     */
    std::istream& read(std::istream& in)
    {
      direction::read(in);
      in >> c_;
      return in;
    }

    /**
     * write to output stream
     * \param  out output stream
     * \return     output stream
     */
    std::ostream& write(std::ostream& out) const
    {
      direction::write(out);
      out << ' ' << c();
      return out;
    }

  protected:
    t_pos c_;    //!< longitudinal position
  };


  /**
   * 1D track
   */
  typedef point track1Z;


  /**
   * 3D track
   */
  class track : public axis {
  public:
    /**
     * Default constructor
     */
    track() :
      axis(),
      t0_(0)
    {}

    /**
     * Constructor
     * \param  a       a position [m]
     * \param  b       b position [m]
     * \param  t0      time at position [ns]
     * \param  theta   polar angle [rad]
     * \param  phi     azimuth angle [rad]
     */
    track(const t_pos&  a,
	  const t_pos&  b,
	  const t_time& t0,
	  const t_dir&  theta,
	  const t_dir&  phi) :
      axis(a,b,theta,phi),
      t0_(t0)
    {}

    /**
     * Constructor
     * \param  pos_    position
     * \param  t0      time
     * \param  dir_    direction
     */
    track(const position&  pos_,
	  const t_time&    t0,
	  const direction& dir_) :
      axis(pos_,dir_),
      t0_(t0)
    {
      // correct t0
      // BvR
      /*
      position o(pos_);
      rotation R(dir_);
      
      o.rotate(R);

      t0_ = t0  -  o.z() * C_INVERSE; 
      */
      // BvR
    }

    /**
     * Constructor
     * \param  point_  3D point
     * \param  dir_    direction
     */
    track(const point&     point_,
	  const direction& dir_) :
      axis(point_,dir_),
      t0_(0)
    {
      // correct t0
      // BvR
      /*
      position o(point_);
      rotation R(dir_);
      
      o.rotate(R);

      t0_ = point_.t()  -  o.z() * C_INVERSE; 
      */
      // BvR
    }
   
    /**
     * Rotate track
     * \param  R  rotation matrix
     * \return    rotated track
     */
    track& rotate(const rotation& R)
    {
      axis::rotate(R);
      return *this;
    }

    /**
     * Rotate back track
     * \param  R  rotation matrix
     * \return    rotated track
     */
    track& rotate_back(const rotation& R)
    {
      axis::rotate_back(R);
      return *this;
    }
 
    /** assignment operator */
    track& operator=(const double& c) 
    { 
      axis::operator=(c);
      t0_ = c;
      return *this;
    }

    /** offset operator */
    track& operator+=(position o)
    {
      rotation R(*this);
      o.rotate(R);
      a_  += o.x();
      b_  += o.y();
      t0_ -= o.z() * C_INVERSE;
      return *this;
    }

    /** offset operator */
    track& operator-=(position o)
    {
      rotation R(*this);
      o.rotate(R);
      a_  -= o.x();
      b_  -= o.y();
      t0_ += o.z() * C_INVERSE;
      return *this;
    }

    /** add operator */
    track& operator+=(const track& o) 
    { 
      axis::operator+=(o);
      t0_ += o.t0();
      return *this;
    }

    /** subtract operator */
    track& operator-=(const track& o) 
    { 
      axis::operator-=(o);
      t0_ -= o.t0();
      return *this;
    }

    /** scale operator */
    track& operator*=(const double& c) 
    { 
      axis::operator*=(c);
      t0_ *= c;
      return *this;
    }

    const t_time& t0() const { return t0_; }   //!< time at starting position [ns]

    /**
     * read from input stream
     * \param  in  input stream
     * \return     input stream
     */
    std::istream& read(std::istream& in)
    {
      axis::read(in);
      in >> t0_;
      return in;
    }

    /**
     * write to output stream
     * \param  out output stream
     * \return     output stream
     */
    std::ostream& write(std::ostream& out) const
    {
      axis::write(out);
      out << ' ' << t0();
      return out;
    }

  protected:
    t_time t0_;
  };


  /**
   * 3D track with velocity
   */
  class vtrack : public track {
  public:
    /**
     * Default constructor
     */
    vtrack() :
      track(),
      beta_(0)
    {}

    /**
     * Constructor
     * \param  a       a position [m]
     * \param  b       b position [m]
     * \param  t0      time at position [ns]
     * \param  theta   polar angle [rad]
     * \param  phi     azimuth angle [rad]
     * \param  beta    velocity
     */
    vtrack(const t_pos&  a,
	   const t_pos&  b,
	   const t_time& t0,
	   const t_dir&  theta,
	   const t_dir&  phi,
	   const float&  beta) :
      track(a,b,t0,theta,phi),
      beta_(beta)
    {}

    /**
     * Constructor
     * \param  pos_    position
     * \param  t0      time
     * \param  dir_    direction
     * \param  beta    velocity
     */
    vtrack(const position&  pos_,
	   const t_time&    t0,
	   const direction& dir_,
	   const float&     beta) :
      track(pos_,t0,dir_),
      beta_(beta)
    {
      // correct t0

      position o(pos_);
      /* // BvR
      rotation R(dir_);
      
      o.rotate(R);*/ // BvR

      t0_ = t0  -  o.z() * C_INVERSE / beta;
    }

    /**
     * Constructor
     * \param  point_  3D point
     * \param  dir_    direction
     * \param  beta    velocity
     */
    vtrack(const point&     point_,
	   const direction& dir_,
	   const float&     beta) :
      track(point_,dir_),
      beta_(beta)
    {
      // correct t0

      position o(point_);
      /* // BvR
      rotation R(dir_);
      
      o.rotate(R);*/ // BvR

      t0_ = point_.t()  -  o.z() * C_INVERSE / beta;
    }
 
    /** assignment operator */
    vtrack& operator=(const double& c) 
    { 
      track::operator=(c);
      beta_ = c;
      return *this;
    }

    /** add operator */
    vtrack& operator+=(const vtrack& o) 
    { 
      track::operator+=(o);
      beta_ += o.beta();
      return *this;
    }

    /** subtract operator */
    vtrack& operator-=(const vtrack& o) 
    { 
      track::operator-=(o);
      beta_ -= o.beta();
      return *this;
    }

    /** scale operator */
    vtrack& operator*=(const double& c) 
    { 
      track::operator*=(c);
      beta_ *= c;
      return *this;
    }

    const float& beta() const { return beta_; }   //!< velocity

    /**
     * read from input stream
     * \param  in  input stream
     * \return     input stream
     */
    std::istream& read(std::istream& in)
    {
      track::read(in);
      in >> beta_;
      return in;
    }

    /**
     * write to output stream
     * \param  out output stream
     * \return     output stream
     */
    std::ostream& write(std::ostream& out) const
    {
      track::write(out);
      out << ' ' << beta();
      return out;
    }

  protected:
    float beta_;
  };


  /** subtract operator */
  static inline point operator-(const point& a, const point& b) 
  { 
    return point(a) -= b;
  }

  /** subtract operator */
  static inline axis operator-(const axis& a, const axis& b) 
  { 
    return axis(a) -= b;
  }

  /** subtract operator */
  static inline track operator-(const track& a, const track& b) 
  { 
    return track(a) -= b;
  }


  /** add operator */
  static inline point operator+(const point& a, const point& b) 
  { 
    return point(a) += b;
  }

  /** add operator */
  static inline axis operator+(const axis& a, const axis& b) 
  { 
    return axis(a) += b;
  }

  /** add operator */
  static inline track operator+(const track& a, const track& b) 
  { 
    return track(a) += b;
  }


  /** scale operator */
  static inline point operator*(const point& a, const double& c) 
  { 
    return point(a) *= c;
  }

  /** scale operator */
  static inline axis operator*(const axis& a, const double& c) 
  { 
    return axis(a) *= c;
  }

  /** scale operator */
  static inline track operator*(const track& a, const double& c) 
  { 
    return track(a) *= c;
  }


  /** scale operator */
  static inline point operator*(const double& c, const point& a)
  { 
    return point(a) *= c;
  }

  /** scale operator */
  static inline axis operator*(const double& c, const axis& a) 
  { 
    return axis(a) *= c;
  }

  /** scale operator */
  static inline track operator*(const double& c, const track& a) 
  { 
    return track(a) *= c;
  }


  /**
   * Get symetry axis.
   * \param begin_     begin of hit data
   * \param end_       end of hit data
   * \param sigma      resolution
   * \return           pair<probability,axis>
   */
  template<class hitIterator>
  static std::pair<double,axis> getAxis(hitIterator begin_, 
					hitIterator end_,
					const double& sigma)
  {
    double x,  y,  z;
    double xx, yy, zz;
    double xz, yz, xy;
      
    if (std::distance(begin_,end_) < 3)
      throw std::invalid_argument("getAxis(): not enough data points");

    x  = y  = z  = 0;
    xx = yy = zz = 0;
    xz = yz = xy = 0;
      
    for (hitIterator i = begin_; i != end_; ++i) {	
      x  += i->x();
      y  += i->y();
      z  += i->z();
	
      xz += i->x() * i->z();
      yz += i->y() * i->z();
      xy += i->x() * i->y();
	
      xx += i->x() * i->x();
      yy += i->y() * i->y();
      zz += i->z() * i->z();
    }
      
      
    const double weight = 1.0 / (double) std::distance(begin_,end_);
      

    vector3d v3[] = {
      vector3d(xx - weight*x*x, xy - weight*x*y, xz - weight*x*z),
      vector3d(xy - weight*x*y, yy - weight*y*y, yz - weight*y*z),
      vector3d(xz - weight*x*z, yz - weight*y*z, zz - weight*z*z)
    };
      

    const int ndf = std::distance(begin_,end_) - 4;


    // principal component

    std::multimap<double, axis, std::greater<double> > zmap;

    for (int i = 0; i != sizeof(v3)/sizeof(v3[0]); ++i) {

      if (v3[i].length() != 0) {

	direction dir(v3[i]);

	position  pos(weight*(x  -  dir.dx() * (i == 0 ? 0 : (i == 1 ? y : z) )),
		      weight*(y  -  dir.dy() * (i == 0 ? x : (i == 1 ? 0 : z) )),
		      weight*(z  -  dir.dz() * (i == 0 ? x : (i == 1 ? y : 0) )));

	axis fit(pos,dir);
      
	rotation R(dir);
      
	double da, db;
	double chi2 = 0;
	  
	for (hitIterator i = begin_; i != end_; ++i) {
	    
	  position o(i->x(),
		     i->y(),
		     i->z());
	    
	  o.rotate(R);
	    
	  da = o.x() - fit.a();
	  db = o.y() - fit.b();
	    
	  chi2 += (da*da + db*db) / (sigma*sigma);
	}

	try {
	  zmap.insert(std::make_pair(NR::prob(chi2,ndf),fit));
	}
	catch(std::exception& error) {
	}
      }
    }

    if (zmap.size() != 0)
      return *(zmap.begin());
    else
      return std::make_pair(0,axis());
  }


  /**
   * Get symetry plane.
   * \param begin_     begin of hit data
   * \param end_       end of hit data
   * \param sigma      resolution
   * \return           pair<chi2,plane>
   */
  template<class hitIterator>
  static std::pair<double,plane> getPlane(hitIterator begin_, 
					  hitIterator end_,
					  const double& sigma)
  {
    double x = 0;
    double y = 0;
    double z = 0;
    
    if (std::distance(begin_,end_) < 3)
      throw std::invalid_argument("getPlane(): not enough data points");

    for (hitIterator i = begin_; i != end_; ++i) {
      x += i->x();
      y += i->y();
      z += i->z();
    }
    
    vector3d pos(x,y,z);

    pos /= (double) std::distance(begin_,end_);


    vector3d v;
    
    for (hitIterator i = begin_; i != end_; ++i) {
      for (hitIterator j = i; ++j != end_; ) {

	vector3d v1(i->x() - pos.x(),
		    i->y() - pos.y(),
		    i->z() - pos.z());
	
	vector3d v2(j->x() - pos.x(),
		    j->y() - pos.y(),
		    j->z() - pos.z());
	
	vector3d r = rot(v1,v2);

	if (dot(v,r) > 0)
	  v += r;
	else
	  v -= r;
      }
    }
    
    direction dir(v);

    plane fit(pos,dir);

    rotation R(dir);

    double residual;
    double chi2 = 0;

    for (hitIterator i = begin_; i != end_; ++i) {
   
      vector3d o(i->x(),
		 i->y(),
		 i->z());

      o.rotate(R);

      residual = (o.z() - fit.c()) / sigma;

      chi2 += residual*residual;
    }

    const int ndf = std::distance(begin_,end_) - 3;

    try {
      return std::make_pair(NR::prob(chi2,ndf), fit);
    }
    catch(std::exception&) {
      return std::make_pair(0, plane());
    }
  }


  /**
   * Input operator for point object.
   * \param  in   input stream
   * \param  o    point
   * \return      input stream
   */
  static inline std::istream& operator>>(std::istream& in, point& o)
  {
    return o.read(in);
  }

  /**
   * Input operator for hit object.
   * \param  in   input stream
   * \param  o    hit
   * \return      input stream
   */
  static inline std::istream& operator>>(std::istream& in, hit& o)
  {
    return o.read(in);
  }

  /**
   * Input operator for axis object.
   * \param  in   input stream
   * \param  o    axis
   * \return      input stream
   */
  static inline std::istream& operator>>(std::istream& in, axis& o)
  {
    return o.read(in);
  }

  /**
   * Input operator for plane object.
   * \param  in   input stream
   * \param  o    plane
   * \return      input stream
   */
  static inline std::istream& operator>>(std::istream& in, plane& o)
  {
    return o.read(in);
  }

  /**
   * Input operator for track object.
   * \param  in   input stream
   * \param  o    track
   * \return      input stream
   */
  static inline std::istream& operator>>(std::istream& in, track& o)
  {
    return o.read(in);
  }

  /**
   * Output operator for point object.
   * \param  out  output stream
   * \param  o    point
   * \return      output stream
   */
  static inline std::ostream& operator<<(std::ostream& out, const point& o)
  {
    return o.write(out);
  }

  /**
   * Output operator for hit object.
   * \param  out  output stream
   * \param  o    hit
   * \return      output stream
   */
  static inline std::ostream& operator<<(std::ostream& out, const hit& o)
  {
    return o.write(out);
  }

  /**
   * Output operator for axis object.
   * \param  out  output stream
   * \param  o    axis
   * \return      output stream
   */
  static inline std::ostream& operator<<(std::ostream& out, const axis& o)
  {
    return o.write(out);
  }

  /**
   * Output operator for plane object.
   * \param  out  output stream
   * \param  o    plane
   * \return      output stream
   */
  static inline std::ostream& operator<<(std::ostream& out, const plane& o)
  {
    return o.write(out);
  }

  /**
   * Output operator for track object.
   * \param  out  output stream
   * \param  o    track
   * \return      output stream
   */
  static inline std::ostream& operator<<(std::ostream& out, const track& o)
  {
    return o.write(out);
  }
}

#endif