#ifndef H_MARKOV_INTEGRATOR
#define H_MARKOV_INTEGRATOR

/**
 * \author mjongen
 */

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

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

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

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

namespace JMARKOV {

  /**
     Abstract base class for calculating the total probability
     (/m^2 target cross-section) for a photon from the source to hit
     the target (with a given, fixed number of scatterings) by Monte 
     Carlo sampling the available nscat*3 dimensional phase space. 

     The sample distribution is implemented in derived classes.

     Random numbers are drawn from gRandom.
  **/
  class JMarkovIntegrator {
    public:
    
    /// standard constructor
    JMarkovIntegrator() {}
    
    /** 
	Integrate with N samples.

	Returns a vector with the contribution to the integral of each sample.
	The mean of those values is the estimate of the result of the integral, 
	while the distribution itself can be used to estimate the stability of 
	the result.
	In this distribution, you want to avoid
	- long tails (because they make the result unstable)
	- small contributions (because it means that parts of the parameter
	  space are being oversampled, so it is less efficient)
	This can be achieved by tuning the sample distribution to the problem 
	at hand.
    **/
    vector<double> integrate( int N, int nscat, JSourceModel* src, JScatteringModel* sm, JTargetModel* trg, double lambda_abs ) ;

    /** 
	Integrate a test function with N samples.

	This can be used as a sanity check for derived classes of 
	JMarkovIntegrator. 

	The integral should yield 1 when the complete relevant part of the
	volume is taken into account.
	If it does not, it may be a sign that the implementation is not
	correct.

	Returns a vector with the contribution to the integral of each sample.
    **/
    vector<double> dummy_integrate( int N, int nscat, JSourceModel* src, JScatteringModel* sm, JTargetModel* trg, double lambda_abs ) ;

   /** 
       Return photon paths generated with the generatePath method.

       This can be used to identify the parts of parameter space that are 
       over- or undersampled in a given problem so that the integrator may 
       be optimized to handle those better.
    **/
    vector<JPhotonPath> get_diagnostic_ensemble( int N, int nscat, JSourceModel* src, JScatteringModel* sm, JTargetModel* trg, double lambda_abs ) ;

    protected:

    /**
       Generate a random photon path with a given number of 
       scatterings

       winv must be set to the inverted probability density
       to generate this particular path.
     **/
    virtual JPhotonPath generatePath( int nscat, double& winv ) {
      // this default implementation assumes that the vertex positions
      // are completely uncorrelated
      JPhotonPath p(nscat) ;
      double _winv = 1 ;
      for( int nv=0 ; nv<nscat+2 ; ++nv ) {
	double part_winv ;
	p[nv] = generatePosition( nscat, nv, part_winv ) ;
	_winv *= part_winv ;
      }
      winv = _winv ;
      return p ;
    }


    /**
       Generate a random position for vertex nv 

       nv = 1 is the first scattering vertex,
       nv = 2 is the second scattering vertex 
       etc.

       So we require  that
       0 < nv <= nscat

       winv has to be set to 1 over the probability density
       for this particular position.
       NOTE: this only works when the positions of each vertex are
       UNCORRELATED to each other (as in the default implementation 
       of generatePath. If they are not, this method can
       be given a dummy implementation and be ignored.
     **/
    virtual JPosition3D generatePosition( int nscat, int nv, double& winv ) = 0 ;

  } ;

  vector<double> JMarkovIntegrator::integrate( int N, int nscat, JSourceModel* src, JScatteringModel* sm, JTargetModel* trg, double lambda_abs ) {
    vector<double> contributions(N,-1) ;

    for( int i=0 ; i<N ; ++i ) {
      double winv ;
      JPhotonPath p = generatePath(nscat,winv) ;

      double rho = getPhotonPathProbabilityDensity(p,src,sm,trg,lambda_abs) ;
      contributions[i] = rho * winv ; 
    }
    return contributions ;
  }

  vector<double> JMarkovIntegrator::dummy_integrate( int N, int nscat, JSourceModel* src, JScatteringModel* sm, JTargetModel* trg, double lambda_abs ) {
    vector<double> contributions(N,-1) ;

    for( int i=0 ; i<N ; ++i ) {
      const double r = 10 ;
      double winv ;
      JPhotonPath p = generatePath(nscat,winv) ;
      double rho = 1.0/(4.0/3.0*M_PI*r*r*r) ;
      if( p[1].getLength()>r ) rho = 0 ;
      contributions[i] = rho * winv ; 
    }
    return contributions ;
  }

  vector<JPhotonPath> JMarkovIntegrator::get_diagnostic_ensemble( int N, int nscat, JSourceModel* src, JScatteringModel* sm, JTargetModel* trg, double lambda_abs ) {
    vector<JPhotonPath> paths ;

    for( int i=0 ; i<N ; ++i ) {
      double winv ;
      JPhotonPath p = generatePath(nscat,winv) ;
      paths.push_back(p) ;
    }
    return paths ;
  }


  /**
     In this implementation of the JMarkovIntegrator interface,
     the sample distribution is a flat distribution.
  **/
  class JMarkovUniformIntegrator : public JMarkovIntegrator {
    public :

    /** 
       constructor. Posmin and posmax contain the minimal and maximal
       values of the coordinate in each direction.
    **/
    JMarkovUniformIntegrator( JPosition3D _posmin, JPosition3D _posmax ) : gen(_posmin,_posmax) {}

    protected:

    virtual JPosition3D generatePosition( int nscat, int nv, double& winv ) {
      JPosition3D pos = gen.getPosition() ;
      winv = 1.0/gen.getWeight(pos) ;
      return pos ;
    }

    JUniformGenerator gen ;
  } ;

  /**
     Abstract base class for implementations of the JMarkovIntegrator interface,
     where the sample distribution is based on histograms filled from an ensemble 
     of representative paths.
  **/
  class JMarkovEnsembleIntegrator : public JMarkovIntegrator {
    public :
    
    // constructor 
    JMarkovEnsembleIntegrator( const vector<JPhotonPath>& _ensemble, JTargetModel* _trg, int _nbinsx, int _nbinsy, int _nbinsz ) : ensemble(_ensemble), nbinsx(_nbinsx), nbinsy(_nbinsy), nbinsz(_nbinsz), trg(_trg) {
      // check that the ensemble size is at least one photon path
      if( ensemble.size()==0 ) {
	cerr << "FATAL ERROR in JMarkov3DensembleIntegrator constructor. Ensemble size = 0." << endl ;
	exit(1) ;
      }

      // check that the position of the first vertex is the same throughout the ensemble
      for( vector<JPhotonPath>::iterator it=ensemble.begin(); it!=ensemble.end(); ++it ) {
	if( (*it)[0].getDistance(ensemble.front()[0]) > 0 ) {
	  cerr << "Fatal error in JMarkovEnsembleIntegrator, first vertex positions are not all the same." << endl ;
	  exit(1) ;
	}
      }
      // set source position based on ensemble
      source_position = ensemble.front()[0] ;
    }

    virtual ~JMarkovEnsembleIntegrator() {} ;

    /// write the histograms filled from the ensemble
    virtual void writeHistograms() = 0 ;

    protected :
    
    /// generate a position for vertex nv for the given number of scatterings, winv will be set to 1/weight
    JPosition3D generatePosition( int nscat, int nv, double& winv ) {
      // source vertex
      if( nv==0 ) {
	winv = 1.0 ;
	return source_position ;
      }
      // target vertex
      if( nv == nscat+1 ) {
	if( nscat+1>(int)target_gens.size() || target_gens[nscat] == NULL ) { 
	  initgenerators(nscat) ; 
	}
	JPosition3D pos = target_gens[nscat]->getPosition() ;
	winv = 1.0/target_gens[nscat]->getWeight(pos) ;
	return pos ;
      } else {
	// scattering vertex
	return generateScatteringVertexPosition(nscat,nv,winv) ;
      }
    }

    virtual JPosition3D generateScatteringVertexPosition(int nscat, int nv, double& winv) = 0 ;

    /// return the internal index for a given number of scatterings and vertex number
    int getIndex( int nv, int nscat ) {
      // indices are distributed as follows
      // [0]: nscat=1, nv=1
      // [1]: nscat=2, nv=1
      // [2]: nscat=2, nv=2
      // [3]: nscat=3, nv=1
      // [4]: nscat=3, nv=2
      // etc.
      return (nscat*(nscat-1))/2 + nv - 1 ;
    }

    /// initialize the position generators for a given number of scatterings (i.e. create histograms and fill them from the ensemble)
    void initgenerators( int nscat ) {
      // target vertex
      init_target_generator(nscat) ;

      // scattering vertices
      init_scattering_vertex_generators(nscat) ;
    }

    /// initialize the generator for a target vertex for a given number of scatterings
    void init_target_generator( int nscat ) {
      if( (int)target_gens.size() <= nscat ) target_gens.resize( nscat+1 ) ;

      if( trg->getRadius()>0 ) {
	// for a finite target, we fill a "sphere generator" from the ensemble
	char hname[200] ;
	sprintf( hname, "htarget_nscat%i", nscat ) ;
	TH2D* ht = new TH2D(hname,";cos(#theta);#phi",100,-1,1,100,-M_PI,M_PI) ;
	for( vector<JPhotonPath>::iterator it=ensemble.begin() ; it!=ensemble.end() ; ++it ) {
	  if( it->n-2 != nscat ) continue ;
	  JDirection3D dir( it->back() - trg->getPosition() ) ;
	  ht->Fill( dir.getDZ(), dir.getPhi() ) ;
	}
	target_gens[nscat] = new JSphereGenerator( trg->getPosition(), trg->getRadius(), ht ) ;
	delete ht ;
      } else {
	// for an infinitesimal target, we create a "generator" that only generates the single point
	target_gens[nscat] = new JSphereGenerator( trg->getPosition() ) ;
      }
    }

    virtual void init_scattering_vertex_generators( int nscat ) = 0 ;

    void free_target_gens() {
      for( vector<JSphereGenerator*>::iterator it=target_gens.begin() ; it!=target_gens.end() ; ++it ) {
	if( *it != NULL ) delete *it ;
      }
    }

    vector<JSphereGenerator*> target_gens ;
    vector<JPhotonPath> ensemble ;
    int nbinsx ;
    int nbinsy ;
    int nbinsz ;
    JTargetModel* trg ;
    JPosition3D source_position ;
  } ;

  /**
     Implementation of JMarkovEnsembleIntegrator interface with 1D histograms.
  **/
  class JMarkovEnsembleIntegrator1D : public JMarkovEnsembleIntegrator {
    public :
    
    // constructor 
    JMarkovEnsembleIntegrator1D( const vector<JPhotonPath>& _ensemble, JTargetModel* _trg, int _nbinsx, int _nbinsy, int _nbinsz ) : JMarkovEnsembleIntegrator(_ensemble,_trg,_nbinsx,_nbinsy,_nbinsz) {}

    ~JMarkovEnsembleIntegrator1D() {
      free_target_gens() ;
      free_scattering_vertex_gens() ;
    }

    void free_scattering_vertex_gens() {
      for( vector<JHistGenerator*>::iterator it=gens.begin() ; it!=gens.end() ; ++it ) {
	if( *it != NULL ) delete *it ;
      }
    }

    void writeHistograms() {
      for( vector<JHistGenerator*>::iterator it=gens.begin() ; it!=gens.end() ; ++it ) {
	if( *it != NULL ) {
	  (*it)->hx->Write() ;
	  (*it)->hy->Write() ;
	  (*it)->hz->Write() ;
	}
      }
      for( vector<JSphereGenerator*>::iterator it=target_gens.begin(); it!=target_gens.end(); ++it ) {
	if( *it != NULL ) {
	  if( (*it)->h != NULL ) (*it)->h->Write() ;
	}
      }
    }

    protected :
    
    void init_scattering_vertex_generators( int nscat ) {
      int ifirst = getIndex(1,nscat) ;
      int ilast  = getIndex(1,nscat+1) ;
      if( (int)gens.size() < ilast ) gens.resize( ilast, NULL ) ;
      for( int nv=1 ; nv<=nscat ; ++nv ) {
	int index = ifirst + nv - 1 ;
	// determine minimal and maximal coordinate values for this vertex
	double xmin =  1.0/0.0 ;
	double xmax = -1.0/0.0 ;
	double ymin =  1.0/0.0 ;
	double ymax = -1.0/0.0 ;
	double zmin =  1.0/0.0 ;
	double zmax = -1.0/0.0 ;
	for( vector<JPhotonPath>::iterator it=ensemble.begin() ; it!=ensemble.end() ; ++it ) {
	  if( it->n-2 != nscat ) continue ;
	  JPosition3D vtx( (*it)[nv] ) ;
	  if( vtx.getX()<xmin ) xmin = vtx.getX() ;
	  if( vtx.getX()>xmax ) xmax = vtx.getX() ;
	  if( vtx.getY()<ymin ) ymin = vtx.getY() ;
	  if( vtx.getY()>ymax ) ymax = vtx.getY() ;
	  if( vtx.getZ()<zmin ) zmin = vtx.getZ() ;
	  if( vtx.getZ()>zmax ) zmax = vtx.getZ() ;
	}
	// allocate histograms
	char hname[200] ;
	sprintf( hname, "hx_nscat%i_vtx%i", nscat, nv ) ;
	TH1F hx(hname,"Vertex position;X",nbinsx,xmin,xmax) ;
	sprintf( hname, "hy_nscat%i_vtx%i", nscat, nv ) ;
	TH1F hy(hname,"Vertex position;Y",nbinsy,ymin,ymax) ;
	sprintf( hname, "hz_nscat%i_vtx%i", nscat, nv ) ;
	TH1F hz(hname,"Vertex position;Z",nbinsz,zmin,zmax) ;
	// fill histograms from ensemble
	for( vector<JPhotonPath>::iterator it=ensemble.begin() ; it!=ensemble.end() ; ++it ) {
	  if( it->n-2 != nscat ) continue ;
	  JPosition3D vtx( (*it)[nv] ) ;
	  hx.Fill( vtx.getX() ) ;
	  hy.Fill( vtx.getY() ) ;
	  hz.Fill( vtx.getZ() ) ;
	}
	// extra: make sure every field is filled at least once
	for( Int_t bin=1 ; bin<=hx.GetNbinsX() ; ++bin ) 
	  hx.AddBinContent(bin,1) ;
	for( Int_t bin=1 ; bin<=hy.GetNbinsX() ; ++bin ) 
	  hy.AddBinContent(bin,1) ;
	for( Int_t bin=1 ; bin<=hz.GetNbinsX() ; ++bin ) 
	  hz.AddBinContent(bin,1) ;
	
	gens[index] = new JHistGenerator( &hx, &hy, &hz ) ;
      }
    }


    JPosition3D generateScatteringVertexPosition(int nscat, int nv, double& winv) {
      int index = getIndex(nv,nscat) ;
      if( index+1>(int)gens.size() || gens[index] == NULL ) initgenerators(nscat) ;

      JPosition3D pos = gens[index]->getPosition() ;
      winv = 1.0/gens[index]->getWeight(pos) ;
      return pos ;
    }

    vector<JHistGenerator*> gens ;
  } ;

  /**
     This implementation of the JMarkovIntegrator interface generates
     'natural' looking paths that might sample the phase space well
     in some cases
  **/
  /*class JMarkovNaturalIntegrator : public JMarkovIntegrator {
    public :

    JMarkovNaturalIntegrator( double _l ) {
      l = _l ; // approximate average path length above the minimum
    }

    ~JMarkovNaturalIntegrator() {} 

    protected :

    /// dummy implementation
    JPosition3D generatePosition( int nscat, int nv, double& winv ) {
      return JPosition3D(0,0,0) ;
    } ;

    /// Generate a 'natural' path
    JPhotonPath generatePath( int nscat, double& winv ) {
      double winvtot = 1 ;
      
      // distance from source to target
      double d = (source_pos-target_pos).getLength() ;

      double g = d / (l+d) ;
      JHenyeyGreensteinScattering sm(100,g) ;

      double lambda = (d+l)/(nscat+1) ;

      JPhotonPath p(nscat) ;
      p.front() = source_pos ;
      p.back()  = target_pos ;

      for( int nv=1 ; nv<=nscat ; ++nv ) {	
	// generate travel length
	double r = gRandom->Exp(lambda) ;
	// generate direction from a HG distribution, relative to THE TARGET.
	JDirection3D dir( p.back()-p[nv-1] ) ;
	JRotation3D rot(dir) ;
	dir = sm.generateDirection() ;
	double ct = dir.getDZ() ;
	dir = dir.rotate_back(rot) ;
	// set vertex
	p[nv] = p[nv-1] + r*JVector3D(dir) ;
	// double 
	winvtot /= ( r * r * sm.getScatteringProbability(ct) * 1.0/lambda*exp(-r/lambda) ) ;
      }

      return p ;
    }

    double l ;
    } ;*/

  /**
     Similar to JMarkovEnsembleIntegrator, but using a 3D histogram for each
     vertex instead of 3 1D histograms.
  **/
  class JMarkovEnsembleIntegrator3D : public JMarkovEnsembleIntegrator {
    public :
    
    // constructor 
    JMarkovEnsembleIntegrator3D( const vector<JPhotonPath>& _ensemble, JTargetModel* _trg, int _nbinsx, int _nbinsy, int _nbinsz, JPosition3D _posmin, JPosition3D _posmax ) : JMarkovEnsembleIntegrator(_ensemble,_trg,_nbinsx,_nbinsy,_nbinsz), posmin(_posmin), posmax(_posmax) {}

    ~JMarkovEnsembleIntegrator3D() {
      free_target_gens() ;
      free_scattering_vertex_gens() ;
    }

    void free_scattering_vertex_gens() {
      // free memory
      for( vector<J3DhistGenerator*>::iterator it=gens.begin() ; it!=gens.end() ; ++it ) {
	if( *it != NULL ) delete *it ;
      }
    }

    void writeHistograms() {
      for( vector<J3DhistGenerator*>::iterator it=gens.begin() ; it!=gens.end() ; ++it ) {
	if( *it != NULL ) {
	  (*it)->h->Write() ;
	  TH2D* htmp ;
	  htmp = (TH2D*)(*it)->h->Project3D("yx") ;
	  htmp->SetOption("colz") ;
	  htmp->SetEntries(10) ;
	  htmp->Write() ;
	  htmp = (TH2D*)(*it)->h->Project3D("yz") ;
	  htmp->SetOption("colz") ;
	  htmp->SetEntries(10) ;
	  htmp->Write() ;
	  htmp = (TH2D*)(*it)->h->Project3D("zx") ;
	  htmp->SetOption("colz") ;
	  htmp->SetEntries(10) ;
	  htmp->Write() ;
	}
      }
    }

    JPosition3D generateScatteringVertexPosition(int nscat, int nv, double& winv) {
      int index = getIndex(nv,nscat) ;
      if( index+1>(int)gens.size() || gens[index] == NULL ) initgenerators(nscat) ;

      JPosition3D pos = gens[index]->getPosition() ;
      winv = 1.0/gens[index]->getWeight(pos) ;
      return pos ;
    }

    protected :
    
    void init_scattering_vertex_generators( int nscat ) {
      int ifirst = getIndex(1,nscat) ;
      int ilast  = getIndex(1,nscat+1) ;
      if( (int)gens.size() < ilast ) gens.resize( ilast, NULL ) ;
      
      for( int nv=1 ; nv<=nscat ; ++nv ) {
	int index = ifirst + nv - 1 ;
	// determine minimal and maximal value
	double xmin =  1.0/0.0 ;
	double xmax = -1.0/0.0 ;
	double ymin =  1.0/0.0 ;
	double ymax = -1.0/0.0 ;
	double zmin =  1.0/0.0 ;
	double zmax = -1.0/0.0 ;
	xmin = posmin.getX() ;
	ymin = posmin.getY() ;
	zmin = posmin.getZ() ;
	xmax = posmax.getX() ;
	ymax = posmax.getY() ;
	zmax = posmax.getZ() ;
	// allocate histogram
	char hname[200] ;
	sprintf( hname, "h3_nscat%i_vtx%i", nscat, nv ) ;
	TH3F h(hname,"Vertex position;X;Y;Z",
	       nbinsx,xmin,xmax,
	       nbinsy,ymin,ymax,
	       nbinsz,zmin,zmax ) ;
	// fill histogram from ensemble
	for( vector<JPhotonPath>::iterator it=ensemble.begin() ; it!=ensemble.end() ; ++it ) {
	  if( it->n-2 != nscat ) continue ;
	  JPosition3D vtx( (*it)[nv] ) ;
	  h.Fill( vtx.getX(), vtx.getY(), vtx.getZ() ) ;
	}
	// extra: make sure every field is filled at least once
	double excuseval = 0.01*h.Integral() / ( h.GetNbinsX() * h.GetNbinsY() * h.GetNbinsZ() ) ;
	for( Int_t xbin=1 ; xbin<=h.GetNbinsX() ; ++xbin ) { 
	  for( Int_t ybin=1 ; ybin<=h.GetNbinsY() ; ++ybin ) { 
	    for( Int_t zbin=1 ; zbin<=h.GetNbinsZ() ; ++zbin ) { 
	      Int_t bin = h.GetBin(xbin,ybin,zbin) ;
	      h.AddBinContent(bin,excuseval) ;
	    }
	  }
	}
	gens[index] = new J3DhistGenerator( &h ) ;
      }
    }

    vector<J3DhistGenerator*> gens ;
    JPosition3D posmin ;
    JPosition3D posmax ;
  } ;

  /**
     In this implementation of the JMarkovIntegrator interface,
     the sample distribution is built up out of three components:
     - a 1/r^2 distribution centered around the source
     - a 1/r^2 distribution centered around the target
     - an exponential distribution around the point directly between the source and the target
  **/
  class JSourceTargetIntegrator : public JMarkovIntegrator {
    public :

    /** 
       constructor.

       Arguments:
       - position of the source
       - position of the target
       - the fraction of points that is drawn from the exponential distribution
       - the minimal radius of the exp. distr.
       - the maximal radius of the exp. distr.
       - the length parameter of the exp. distr.
       - the maximal radius for the 1/r^2 distributions
    **/
    JSourceTargetIntegrator( double Rsing, JPosition3D src_pos, JPosition3D trg_pos, double fraction_exp, double Lexp, double Rexp ) : gen_exp(0,Rexp,Lexp), gen_s1(Rsing,src_pos), gen_s2(Rsing,trg_pos) {
      JPosition3D cnt_pos = 0.5*(src_pos+trg_pos) ;
      gen_exp_shift = new JShiftedGenerator(&gen_exp,cnt_pos) ;
      double fraction_sing = 0.5*(1-fraction_exp) ;
      gen = new JTripleGenerator( fraction_sing, &gen_s1, fraction_sing, &gen_s2, fraction_exp, gen_exp_shift ) ;
    }

    ~JSourceTargetIntegrator() {
      delete gen_exp_shift ;
      delete gen ;
    }

    protected:

    virtual JPosition3D generatePosition( int nscat, int nv, double& winv ) {
      JPosition3D pos = gen->getPosition() ;
      winv = 1.0/gen->getWeight(pos) ;
      return pos ;
    }

    JExponentialGenerator gen_exp ;
    JSingularityGenerator gen_s1 ;
    JSingularityGenerator gen_s2 ;
    JGenerator* gen_exp_shift ;
    JTripleGenerator* gen ;
  } ;

  /**
     In this implementation of the JMarkovIntegrator interface,
     the sample distribution is built up out of correlated path vertices.


     lambda is the effective length scale defined as:

     1/lambda := 1/lambda_scat + 1/lambda_abs
  **/
  class JExperimentalIntegrator : public JMarkovIntegrator {
    public :

    /** 
       constructor.
    **/
    JExperimentalIntegrator( double _lambda ) : gen(_lambda) {}

    ~JExperimentalIntegrator() {}

    protected:

    /// dummy implementation
    JPosition3D generatePosition( int nscat, int nv, double& winv ) {
      return JPosition3D(0,0,0) ;
    } ;

    JPhotonPath generatePath( int nscat, double& winv, JPosition3D source_pos, JPosition3D target_pos ) {
      JPhotonPath p(nscat) ;
      p.front() = source_pos ;
      p.back()  = target_pos ;
      double _winv = 1 ;

      for( int nv=1 ; nv<=nscat ; ++nv ) {
	double __winv ;
	p[nv] = getFirstScatteringVertex(nscat+1-nv,p[nv-1],p.back(),__winv ) ;
	_winv *= __winv ;
      }
      winv = _winv ;
      return p ;
    }

    protected:

    /**
       recursive function to randomly generate the first scattering
       vertex for isotropic scattering from a source to a target
    **/
    JPosition3D getFirstScatteringVertex( int nscat, JPosition3D src_pos, JPosition3D trg_pos, double& winv ) {
      double _winv = 1 ;
      if( nscat == 0 ) {
	winv = 1 ;
	return trg_pos ;
      }
      if( gRandom->Integer(2) == 0 ) {
	src_pos = getFirstScatteringVertex(nscat-1,src_pos,trg_pos,_winv) ;
      }
      JPosition3D pos = gen.getPosition() ;
      winv = 1.0 / gen.getWeight(pos) ;
      return src_pos + pos ;
    }

    JExpRsqInvGenerator gen ;
  } ;

}

#endif