#ifndef H_J_PHOTON_PATH
#define H_J_PHOTON_PATH

/**
 * \author mjongen
 */

// c++ standard library includes
#include<iostream>
#include<sstream>
#include<iomanip>
#include<cmath>
#include<cstdlib>
#include<vector>

// root includes
#include "TObject.h"
#include "TH1F.h"
#include "TRandom3.h"

// JPP includes
#include "JGeometry3D/JPolyline3D.hh"
#include "JGeometry3D/JVersor3D.hh"
#include "JIO/JSerialisable.hh"

namespace JMARKOV {}
namespace JPP { using namespace JMARKOV; }

namespace JMARKOV {
  using namespace std ;
  using namespace JGEOMETRY3D ;
  using namespace JGEOMETRY2D ;

  /**
     A photon path.
     This is basically a JPolyline3D with at least two vertices and a weight.
  **/
  class JPhotonPath : public JPolyline3D, public JIO::JSerialisable {
   public:
    
    /// default constructor
    JPhotonPath() : n(2), weight(1) { resize(n) ; }
    
    /// constructor to actually be used
    JPhotonPath( int _nscat ) : n(_nscat+2), weight(1) { resize(n) ; }
    
    /// get a string with the path vertex positions
    string getString() ;

    /// get the total path length
    double getLength() ;

    void addVertex() {
      JPosition3D pos(0,0,0) ;
      addVertex(pos) ;
    }
    void addVertex( JPosition3D& pos ) ;

    /// Returns whether the photon path intersects a sphere of radius r at position pos
    bool hitsSphere( const JPosition3D& pos, double r ) ;

    /// Returns the position where the photon hits a sphere of radius r at position pos. If the photon does not hit, it returns (0,0,0).
    JPosition3D getSphereHitPosition( const JPosition3D& pos, double r ) ;
  

    /// return all coordinates where the photon path has the given x
    std::vector<JPosition3D> getPointsWithX( double x ) ;

    /// return all coordinates where the photon path has the given y
    std::vector<JPosition3D> getPointsWithY( double y ) ;

    /// return all coordinates where the photon path has the given z
    std::vector<JPosition3D> getPointsWithZ( double z ) ;

    JIO::JReader& read(JIO::JReader& in ) {
      // read number of vertices and adjust size accordingly
      in >> n ;
      if( n>0 ) resize(n) ;
      // read vertices
      for( iterator it=begin() ; it!=end() ; ++it )
	in >> *it ;
      // read weight
      in >> weight ;

      return in ;
    }

    JIO::JWriter& write(JIO::JWriter& out ) const {
      // write size
      out << (int)size() ;
      // write vertices
      for( const_iterator it=begin() ; it!=end() ; ++it )
	out << *it ; //it->getX() << it->getY() << it->getZ() ;
      // write weight
      out << weight ;

      return out ;
    }

    int n ; // the number of vertices in the path
    double weight ;
  } ;

  double JPhotonPath::getLength() {
    iterator v1 = begin() ;
    iterator v2(v1) ;
    ++v2 ;
    double L = 0 ;
    for(  ; v2!=end() ; ++v1,++v2 ) {
      // length of path segment
      L += (*v2-*v1).getLength() ;
    }
    return L ;
  }

  void JPhotonPath::addVertex( JPosition3D& pos ) {
    ++n ;
    push_back(pos) ;
  }

  std::vector<JPosition3D> JPhotonPath::getPointsWithX( double x ) {
    std::vector<JPosition3D> res ;

    // loop over the path segments
    for( unsigned int i=1 ; i<size() ; ++i ) {
      JPosition3D v1( at(i-1) ) ;
      JPosition3D v2( at(i)   ) ;
      if( (v1.getX()>x) != (v2.getX()>x) ) {
	double frac =  (x-v1.getX())/(v2.getX()-v1.getX()) ;
	res.push_back( v1+frac*(v2-v1) ) ;
      }
    }
    return res ;
  }

  string JPhotonPath::getString() {
    ostringstream oss ;
    for( unsigned int i=0 ; i<size() ; ++i ) {
      oss << at(i) ;
      if( i+1!=size() ) oss << ", " ;
    }
    return oss.str() ;
  }

  std::vector<JPosition3D> JPhotonPath::getPointsWithY( double y ) {
    std::vector<JPosition3D> res ;

    // loop over the path segments
    for( unsigned int i=1 ; i<size() ; ++i ) {
      JPosition3D v1( at(i-1) ) ;
      JPosition3D v2( at(i)   ) ;
      if( (v1.getY()>y) != (v2.getY()>y) ) {
	double frac =  (y-v1.getY())/(v2.getY()-v1.getY()) ;
	res.push_back( v1+frac*(v2-v1) ) ;
      }
    }
    return res ;
  }

  std::vector<JPosition3D> JPhotonPath::getPointsWithZ( double z ) {
    std::vector<JPosition3D> res ;

    // loop over the path segments
    for( unsigned int i=1 ; i<size() ; ++i ) {
      JPosition3D v1( at(i-1) ) ;
      JPosition3D v2( at(i)   ) ;
      if( (v1.getZ()>z) != (v2.getZ()>z) ) {
	double frac =  (z-v1.getZ())/(v2.getZ()-v1.getZ()) ;
	res.push_back( v1+frac*(v2-v1) ) ;
      }
    }
    return res ;
  }


  bool JPhotonPath::hitsSphere( const JPosition3D& pos, double r ) {
    // loop over the path segments
    for( unsigned int i=1 ; i<size() ; ++i ) {
      // the current path segment is from x1 to x2

      // vector pointing from x1 to x2
      JPosition3D seg = at(i)-at(i-1) ;
      JVersor3D dir_seg(seg) ; // normalized version of seg
      double Lseg = seg.getLength() ;

      // vector pointing from the center of the sphere to x1
      JPosition3D to_center(at(i-1)-pos) ;
      JVersor3D dir_to_center(to_center) ; // normalized version of to_center
      double r_to_center = to_center.getLength() ;

      // cosine of angle between the two directions
      double ct = dir_seg.getDot(dir_to_center) ;
      
      // if sin^2(theta) > (r_to_center/radius)^2
      if( 1-ct*ct > r*r/(r_to_center*r_to_center) ) {
	// definitely no hit in this segment
	// continue with the next one
	continue ;
      }

      // the segment can be parametrized as x1 + xi * dir_seg
      // where xi is in the range [0,Lseg]
      // there are two possible solutions for xi
      double D = sqrt( r_to_center*r_to_center*(ct*ct-1) + r*r ) ;
      double xi1 = -r_to_center*ct + D ;
      if( xi1>0 && xi1<Lseg ) return true ;
      double xi2 = -r_to_center*ct - D ;
      if( xi2>0 && xi2<Lseg ) return true ;
    }
    return false ;
  }

  JPosition3D JPhotonPath::getSphereHitPosition( const JPosition3D& pos, double r ) {
    // loop over the path segments
    for( unsigned int i=1 ; i<size() ; ++i ) {
      // the current path segment is from x1 to x2

      // vector pointing from x1 to x2
      JPosition3D seg = at(i)-at(i-1) ;
      JVersor3D dir_seg(seg) ; // normalized version of seg
      double Lseg = seg.getLength() ;

      // vector pointing from the center of the sphere to x1
      JPosition3D to_center(at(i-1)-pos) ;
      JVersor3D dir_to_center(to_center) ; // normalized version of to_center
      double r_to_center = to_center.getLength() ;

      // cosine of angle between the two directions
      double ct = dir_seg.getDot(dir_to_center) ;
      
      // if sin^2(theta) > (r_to_center/radius)^2
      if( 1-ct*ct > r*r/(r_to_center*r_to_center) ) {
	// definitely no hit in this segment
	// continue with the next one
	continue ;
      }

      // the segment can be parametrized as x1 + xi * dir_seg
      // where xi is in the range [0,Lseg]
      // there are two possible solutions for xi
      double D = sqrt( r_to_center*r_to_center*(ct*ct-1) + r*r ) ;

      // first xi
      double xi1 = -r_to_center*ct + D ;
      if( xi1<0 || xi1>Lseg ) {
	xi1 = 1.1*Lseg ;
      }

      // second xi
      double xi2 = -r_to_center*ct - D ;
      // if this is a xi within the segment
      if( xi2>0 && xi2<Lseg ) {
	// and it is before xi1, we replace xi1
	if( xi2<xi1 ) xi1 = xi2 ;
      }

      if( xi1>0 && xi1<Lseg ) {
	return at(i-1) + JPosition3D(dir_seg)*xi1 ;
      }
    }
    return JPosition3D(0,0,0) ;
  }


}

#endif