#ifndef H_MARKOV_PHOTON_TRACKER
#define H_MARKOV_PHOTON_TRACKER

/**
 * \author mjongen
 */

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

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

// JPP includes
#include "JMarkov/JScatteringModel.hh"

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

namespace JMARKOV {

  /**
     The JMarkovPhotonTracker generates ensembles of photon paths by
     tracking photons from a source to a target.

     The target is modelled as a sphere with finite radius r.

     Random numbers are drawn from gRandom.
  **/
  class JMarkovPhotonTracker {
    public:
    
    /// standard constructor
    JMarkovPhotonTracker() {}
  
    /** 
	Track n photons. 

	The tracks of all photons that hit the target are returned.
    **/
    std::vector<JPhotonPath> trackPhotons( unsigned long int n, JSourceModel* src, JScatteringModel* sm, JTargetModel* trg, double lambda_abs ) ;

  } ;

  std::vector<JPhotonPath> JMarkovPhotonTracker::trackPhotons( unsigned long int n, JSourceModel* src, JScatteringModel* sm, JTargetModel* trg, double lambda_abs ) {
    const double lambda_scat = sm->getScatteringLength() ;
    std::vector<JPhotonPath> ensemble ;

    for( unsigned long int i=0 ; i<n ; ++i ) {
      // how long until the photon is absorbed
      double Lmax = gRandom->Exp(lambda_abs) ;
      // initial direction
      JDirection3D dir = src->generateDirection() ;

      // current length of the photon path
      double L = 0 ;
      // how often the photon has scattered so far
      int nscat = 0 ;

      // initialize path
      JPhotonPath p(0) ;
      p[0] = src->getPosition() ;
      p[1] = src->getPosition() ; // temporary placeholder

      while( true ) {
	// distance until next scattering (or absorption)
	double stepsize = gRandom->Exp( lambda_scat ) ;
	// whether this is going to be the last step
	bool isAbsorbed = false ;
	if( stepsize > Lmax-L ) {
	  stepsize = Lmax-L ;
	  isAbsorbed = true ;
	}

	// set new vertex position
	p[1+nscat] = p[nscat] + stepsize * JPosition3D(dir) ;

	// check whether the photon hits the target
	JPhotonPath lastvertex(0) ;
	lastvertex[0] = p[p.size()-2] ;
	lastvertex[1] = p[p.size()-1] ;
	if( lastvertex.hitsSphere(trg->getPosition(),trg->getRadius()) ) {
	  const JPosition3D hit_position = lastvertex.getSphereHitPosition( trg->getPosition(), trg->getRadius() ) ;

	  // correct the position of the final vertex
	  p[p.size()-1] = hit_position  ;
	  // weigh the photon path with the target efficiency
	  p.weight = trg->getEfficiency( JDirection3D(p[p.size()-1]-p[p.size()-2]) ) ;
	  // add to the ensemble
	  ensemble.push_back(p) ;
	  break ;
	}

	if( isAbsorbed ) break ;

	L += lastvertex.getLength() ;

	// scatter
	++nscat ;
	JRotation3D rot(dir) ;
	dir = sm->generateDirection() ;
	dir = dir.rotate_back(rot) ;
	p.addVertex() ;
      }
    }

    return ensemble ;
  }

}

#endif