JMARKOV Namespace Reference

Classes
class	JGenerator
	Abstract interface for the generation of points in 3D space. More...
class	JUniformGenerator
	Implementation of the JGenerator interface. More...
class	JBallGenerator
	Implementation of the JGenerator interface. More...
class	JExpRsqInvGenerator
	Implementation of the JGenerator interface. More...
class	JMagicalDistribution
	The 'magical distributions' are a class of distributions. More...
class	JSingularityGenerator
	Implementation of the JGenerator interface. More...
class	JExponentialGenerator
	Implementation of the JGenerator interface. More...
class	JGaussianGenerator
	Implementation of the JGenerator interface. More...
class	JCombinedGenerator
	Implementation of the JGenerator interface. More...
class	JTripleGenerator
	Implementation of the JGenerator interface. More...
class	JShiftedGenerator
	Implementation of the JGenerator interface. More...
class	JHistGenerator
	Implementation of the JGenerator interface. More...
class	JSphereGenerator
	Implementation of the JGenerator interface. More...
class	J3DhistGenerator
	Implementation of the JGenerator interface. More...
class	JMarkovIntegrator
	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. More...
class	JMarkovUniformIntegrator
	In this implementation of the JMarkovIntegrator interface, the sample distribution is a flat distribution. More...
class	JMarkovEnsembleIntegrator
	Abstract base class for implementations of the JMarkovIntegrator interface, where the sample distribution is based on histograms filled from an ensemble of representative paths. More...
class	JMarkovEnsembleIntegrator1D
	Implementation of JMarkovEnsembleIntegrator interface with 1D histograms. More...
class	JMarkovEnsembleIntegrator3D
	This implementation of the JMarkovIntegrator interface generates 'natural' looking paths that might sample the phase space well in some cases. More...
class	JSourceTargetIntegrator
	In this implementation of the JMarkovIntegrator interface, the sample distribution is built up out of three components: More...
class	JExperimentalIntegrator
	In this implementation of the JMarkovIntegrator interface, the sample distribution is built up out of correlated path vertices. More...
class	JMarkovPathGenerator
	The JMarkovPathGenerator generates ensembles of photon paths using a Markov Chain Monte Carlo (MCMC). More...
class	JMarkovPhotonTracker
	The JMarkovPhotonTracker generates ensembles of photon paths by tracking photons from a source to a target. More...
class	JPhotonPath
	A photon path. More...
class	JPhotonPathReader
class	JPhotonPathWriter
class	JScatteringModel
	Virtual base class for a scattering model. More...
class	JSourceModel
	Virtual base class for a light source. More...
class	JIsotropicSource
	Implementation of the JSourceModel class that represents an isotropic source. More...
class	JDirectedSource
	Implementation of the JSourceModel class that represents a directed source with a flat intensity distribution. More...
class	JTargetModel
	Virtual base class for a light detector ("photon target"). More...
class	JPMTTarget
	Implementation of the JTargetModel class that represents a single PMT on a DOM. More...
class	JIsotropicTarget
	Implementation of the JTargetModel class that represents a spherically symmetric target. More...
class	JCosineTarget
	Implementation of the JTargetModel class that represents a directed target with a cos(theta) acceptance. More...
class	JHenyeyGreensteinScattering
	Implementation of the JScatteringModel interface with scattering according to the Henyey-Greenstein function. More...
class	JIsotropicScattering
	Implementation of the JScatteringModel interface with isotropic scattering. More...
class	JRayleighScattering
	Implementation of the JScatteringModel interface with Rayleigh scattering. More...
class	JCombinedScattering
	Implementation of the JScatteringModel interface with that combines two scattering models into one effective model. More...

Functions
double	getPhotonPathProbabilityDensity (JPhotonPath &p, JSourceModel *src, JScatteringModel *sm, JTargetModel *trg, double lambda_abs)
	Return the probability density for a photon path with the given ingredients. More...

Detailed Description

Author

mjongen

Function Documentation

double JMARKOV::getPhotonPathProbabilityDensity	(	JPhotonPath &	p,
		JSourceModel *	src,
		JScatteringModel *	sm,
		JTargetModel *	trg,
		double	lambda_abs
	)

Return the probability density for a photon path with the given ingredients.

These are

• a model of the source direction distribution
• a scattering model
• a target model
• the absorption length (note that this is allowed to be +INFTY)

Multiply by dV1, dV2, etc. and the target cross-section sigma to get a probability.

The first (last) vertex of the photon path is used as the source (target) position.

Definition at line 582 of file JScatteringModel.hh.

                                                                                                                                          {
    //const bool verbose = true ;

    // calculate inverse scattering and absorption lengths
    const double lambda_scat_inv = 1.0/sm->getScatteringLength() ;
    double lambda_abs_inv = 0 ;
    // if the absorption length is not infinite
    if( lambda_abs == lambda_abs ) lambda_abs_inv = 1.0/lambda_abs ;

    double rho = 1 ;
    const int nvert = p.n   ; // number of vertices in the path
    const int nscat = p.n-2 ; // number of scattering vertices in the path

    //if( verbose ) cout << "$ rho = " << rho << " (before calculating path segment lengths)" << endl ;

    // calculate path segment lengths
    // index i corresponds to the path segment connecting 
    // vertex i and vertex i+1
    double Ltotal = 0 ;
    vector<double> lengths(nvert-1) ; 
    for( int i=0 ; i<nvert-1 ; ++i ) {
      lengths[i] = (p[i+1]-p[i]).getLength() ;
      Ltotal += lengths[i] ;
    }

    //if( verbose ) cout << "$ Total length = " << Ltotal << " m" << endl ;

    // cosine of scattering angles (angle between previous segment and new segment)
    // index i is the cosine of the angle of segment i+1 w.r.t segment i
    vector<double> cts(nscat) ;

    // scattering angles
    for( int i=0 ; i<nscat ; ++i ) {
      JVersor3D seg1(p[i+1]-p[i]) ;
      JVersor3D seg2(p[i+2]-p[i+1]) ;
      cts[i] = seg1.getDot(seg2) ;
    }

    // term for absorption length
    // (=probability to not be absorbed within Ltotal)
    rho *= exp(-Ltotal*lambda_abs_inv ) ;

    //if( verbose ) cout << "$ rho = " << rho << " (after absorption length)" << endl ;

    // term for emission profile
    rho *= src->getEmissionProbability( JVersor3D(p[1]-p[0]) ) ;

    //if( verbose ) cout << "$ rho = " << rho << " (after emission probability)" << endl ;

    // term for target efficiency
    if( trg->getRadius()>0 ) {
      // finite target (spherical), the efficiency is a measure of how good the
      // target is at a given spot
      rho *= trg->getEfficiency( JVersor3D(p[nscat+1]-trg->getPosition()) ) ;
      // this factor accounts for the angle of incidence
      double ct = JVersor3D(p[nscat]-p[nscat+1]).getDot( JVersor3D(p[nscat+1]-trg->getPosition()) ) ;
      rho *= max(0.0,ct) ;

      //if( verbose ) cout << "$ rho = " << rho << " (after target efficiency and angle of incidence)" << endl ;

      // check whether the photon path hits the sphere
      // we do not need to consider the last vertex, because if that intersects the sphere
      // the angle of incidence factor above would be 0.
      if( p.n>2 ) {
        JPhotonPath pshort(p) ;
        pshort.resize( pshort.size()-1 ) ;
        --pshort.n ;
        if( pshort.hitsSphere(trg->getPosition(),trg->getRadius()) ) {
          /*cout << "Path hits sphere at " 
               << JPosition3D(pshort.getSphereHitPosition(trg->getPosition(),trg->getRadius())-trg->getPosition()) 
               << endl ;*/
          rho *= 0 ;
        }
      }
    } else {
      rho *= trg->getEfficiency( p[nscat+1]-p[nscat] ) ;
    }

    //if( verbose ) cout << "$ rho = " << rho << " (after target efficiency)" << endl ;

    // terms for scattering directions
    for( int i=0 ; i<nscat ; ++i ) {
      rho *= sm->getScatteringProbability(cts[i]) ;
    }

    //if( verbose ) cout << "$ rho = " << rho << " (after scattering directions)" << endl ;

    // phase space terms for scattering to each of the scattering 
    // vertices
    // this is the probability/dV to scatter into a volume element 
    // dV at a given distance
    // (although it is missing the direction-dependent term, which is
    // the emission probability or scattering probability)
    for( int i=0; i<nscat; ++i ) {
      // reminder: lengths[i] is the length from vertex i to vertex i+1
      rho *= lambda_scat_inv / (lengths[i]*lengths[i]) * exp(-lengths[i]*lambda_scat_inv) ;
    }

    //if( verbose ) cout << "$ rho = " << rho << " (after scattering phase space term)" << endl ;

    // phase space term for scattering to the target
    // this is the probability/dsigma to hit a target at a given
    // distance, without scattering along the way
    rho *= exp(-lengths[nscat]*lambda_scat_inv) / ( lengths[nscat]*lengths[nscat] ) ;

    //if( verbose ) cout << "$ rho = " << rho << " (after target phase space term)" << endl ;

    // check for NaN
    if( rho!=rho ) {
      cerr << "NaN for probability density" << endl ;
      cerr << "rho = " << rho << endl ;
      exit(1) ;
    }

    //if( verbose ) cout << "$ rho = " << rho << endl ;

    return rho ;
  }