JMarkovPhotonTracker.hh

Go to the documentation of this file.

#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