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
      JVersor3D 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( 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