Beamsource - topasmc/dicom-interface GitHub Wiki

Beamsource has a hierarchical structure. It consists of beamlets, each beamlet consists of distribution functions of energy and phase-space variables.

Distributions

All distributions in RTI library is a child class of the pure virtual class

template<typename T, std::size_t M> class pdf_Md

pdf_Md<typename T, size_t M> returns M random variables of type T from its probability distribution function. For example, 1 dimensional distributions in RTI are

template<typename T> class const_1d : public pdf_Md<T,1> //constant distribution
template<typename T> class uni_1d   : public pdf_Md<T,1>   //uniform distribution
template<typename T> class norm_1d  : public pdf_Md<T,1>  //normal distribution

These distributions are used mostly for sampling energy.

For phase-space variables, two 6-dimensional distributions are provided

template<typename T> class php_6d : public pdf_Md<T,6>   //Bivariate gaussian distribution for a spot like beam
template<typename T> class php_6d_fanbeam : public pdf_Md<T,6>   //Bivariate gaussian distribution for a line like beam

For example, a 1D gaussian distribution with mean value 150 and 0.8 can be defined as

rti::distributions::norm_1d<float> pnorm({150},{0.8});
pnorm(); //will return energy value in std::array<T,1>.

Beamlet

Beamlet class consists of PDF of energy and phase-space and coordinate transform.

template <typename T>
class beamlet {
protected:
    rti::distributions::pdf_Md<T,1>* energy  = nullptr;  
    rti::distributions::pdf_Md<T,6>* fluence = nullptr;
    coordinate_transform<T> p_coord;
}

Following code is an example of a phase-space distribution phsp_6d that we samples, $x, xp, y,yp, z, zp$ from. The spot's mean position is (x,y,z) = (0,0,39.0) and its sigma is (0.1, 0.1, 0) meaning that spot-size x and y is 0.1 and 0.1 but no z position uncertainty. The 0,0, and 0 are the mean of direction cosine but actually z' is calculated from $\sqrt{1-xp^2-yp^2}$ so the $zp$ in spot_mean is ignored.

Figure A diagram of beamlet.

std::array<float,6> spot_mean ={0, 0, 39.0,  0, 0, 0};
std::array<float,6> spot_sigm ={0.1, 0.1, 0, 0.01, 0.05, 0};
std::array<float,2> corr ={0.0, 0.0};
rti::distributions::phsp_6d<float> spot_(spot_mean, spot_sigm, corr);

A coordinate system related with gantry, couch, collimator, etc and iso-center shift can be created as follow.

//colli, gantry, couch, iec2dicom -> order is important
std::array<float, 4> angles = {0*deg2rad, 0*deg2rad, 0*deg2rad, 90.0*M_PI/180.0}; 
coordinate_transform<float> coordinate0(angles,{0,0,0});

Finally, a history can be sampled from beamlet as

beamlet<float> beamlet_1(&pnorm, &spot_) ;
beamlet_1.set_coordinate_transform(coordinate0);

beamlet<float> beamlet_1_wo_coordinate(&pnorm, &spot_) ;

auto h1 = beamlet_1();                 //return one history in global coordinate system
auto h2 = boomlet_1_wo_coordinate();   //return one history in local coordinate system.

float energy = std::get<0>(h1); //energy 
rti::vec3<float> pos = std::get<1>(h2); //position vector, rti::vec3
rti::vec3<float> dir = std::get<2>(h2); //direction cosine vector, rti::vec3

Beamsource

Beamsource class is a container for Beamlet instances and interacts with MC engine.

Figure A diagram of beamsource.

template <typename T> 
class beamsource {
public:
    /// Each element is a tuple of beamlet, number of histories, and accumlated histories
    std::vector< std::tuple<beamlet<T>, size_t, size_t> > beamlets; 

    /// A lookup table to map a history to beamlet id.
    std::map<std::size_t, std::uint32_t>   cdf2beamlet; 
}

A characterized boomlet is added to beam object with its number of histories to be simulated as

beamsource <float> beam_0;
beam_0.append_beamlet(beamlet_1, 10);              //10 histories of beamlet_1
beam_0.append_beamlet(beamlet_1, 20, coordinate0); //with coordinate transformation.

So, then you can generate histories beamlet-by-beamlet

    //generate histories from a beam model
    for(size_t id = 0 ; id < beam.total_beamlets() ; ++id){
       auto beamlet = beam[id];
       auto beamlet_dist = std::get<0> (beamlet);
       auto histories = std::get<0> (beamlet);
       while(histories--){
           auto history = beamlet_dist();
       }
    }

or history by history

    const size_t total_histories = beam.total_histories();
    auto  h =0; 
    while(h_id++ < total_histories){
        auto history = beam(h);
    }
⚠️ **GitHub.com Fallback** ⚠️