egs_focal_spot_source.h

Go to the documentation of this file.

/*
###############################################################################
#
#  EGSnrc egs++ focal spot source headers
#  Copyright (C) 2025 Marvin Apel
#
#  This file is part of EGSnrc.
#
#  EGSnrc is free software: you can redistribute it and/or modify it under
#  the terms of the GNU Affero General Public License as published by the
#  Free Software Foundation, either version 3 of the License, or (at your
#  option) any later version.
#
#  EGSnrc is distributed in the hope that it will be useful, but WITHOUT ANY
#  WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
#  FOR A PARTICULAR PURPOSE.  See the GNU Affero General Public License for
#  more details.
#
#  You should have received a copy of the GNU Affero General Public License
#  along with EGSnrc. If not, see <http://www.gnu.org/licenses/>.
#
###############################################################################
#
#  Author:          Marvin Apel, 2025
#
#  Contributors:
#
###############################################################################
*/
 
 
#ifndef EGS_FOCAL_SPOT_
#define EGS_FOCAL_SPOT_
 
#include "egs_base_source.h"
#include "egs_vector.h"
#include "egs_rndm.h"
#include "egs_transformations.h"
 
 
#ifdef WIN32
 
    #ifdef BUILD_FOCAL_SPOT_DLL
        #define EGS_FOCAL_SPOT_EXPORT __declspec(dllexport)
    #else
        #define EGS_FOCAL_SPOT_EXPORT __declspec(dllimport)
    #endif
    #define EGS_FOCAL_SPOT_LOCAL
 
#else
 
    #ifdef HAVE_VISIBILITY
        #define EGS_FOCAL_SPOT_EXPORT __attribute__ ((visibility ("default")))
        #define EGS_FOCAL_SPOT_LOCAL  __attribute__ ((visibility ("hidden")))
    #else
        #define EGS_FOCAL_SPOT_EXPORT
        #define EGS_FOCAL_SPOT_LOCAL
    #endif
 
#endif
 
class EGS_FOCAL_SPOT_EXPORT EGS_FocalSpot : public EGS_BaseSimpleSource {
    bool valid;               
    bool is_deviating=false;  
    bool is_rotated=false;    
    // The beam is always pointing in the + z-direction
    EGS_Float  z_pos;               
    EGS_Float  sigma_x_space;       
    EGS_Float  sigma_y_space;       
    EGS_Float  space_cutoff_x;      
    EGS_Float  space_cutoff_y;      
    EGS_Float  sigma_x_angle=0;     
    EGS_Float  sigma_y_angle=0;     
    EGS_Float  x_translation=0;     
    EGS_Float  y_translation=0;     
    EGS_Float  x_rotation=0;        
    EGS_Float  y_rotation=0;        
    EGS_Float  z_point_of_rotation; 
 
public:
 
    EGS_FocalSpot(EGS_Input *, EGS_ObjectFactory *f=0);
    ~EGS_FocalSpot() {};
 
    void getPositionDirection(EGS_RandomGenerator *rndm, EGS_Vector &x, EGS_Vector &u, EGS_Float &wt) {
        wt = 1;
        // 1. Sample position
        do {
            double r = sqrt(-2*log(1-rndm->getUniform())); // store value characteristic for Gaussian Sampling
            double phi = PI2*rndm->getUniform();             // store value for phi
            x.x = sigma_x_space * r * cos(phi);
            x.y = sigma_y_space * r * sin(phi);
            // Only Accept points within cutoff limits
        }
        while (pow(x.x/space_cutoff_x,2) + pow(x.y/space_cutoff_y,2) > 1);
        x.x += x_translation; // save computation time by only applying translation after finding a valid point
        x.y += y_translation; // otherwise it will be done in the while condition above repeatedly
        x.z = z_pos;
 
        // 2. Sample direction
        switch (angle_mode) {
        case 1: { // Deviation in x and y
            double phi = PI2 * rndm->getUniform(); // store value of rotation arc phi !
            do {
                u.z = sqrt(-2*pow((sigma_x_angle*sigma_y_angle), 2)/
                           (pow(sigma_x_angle*sin(phi),2) + pow(sigma_y_angle*cos(phi),2)) * (log(1- rndm->getUniform())));
                u.z = cos(u.z);
            }
            while (u.z < 0);
 
            double r = sqrt((1-u.z)*(1+u.z));
            u.y = r * sin(phi);
            u.x = r * cos(phi);
 
            break;
        }
 
        case 2: // Deviation only in x, deviation is planar only along x-direction ! (u.y = 0 !!!)
            do {
                u.z = cos(sigma_x_angle * sqrt(-2*log(rndm->getUniform())));
            }
            while (u.z < 0);
            u.x = sqrt((1-u.z)*(1+u.z));
            u.x = (rndm->getUniform() < 0.5) ? u.x : -1*u.x;
            u.y = 0;
            break;
 
        case 3: // Deviation only in y, deviation is planar only along y-direction ! (u.x = 0 !!!)
            do {
                u.z = cos(sigma_y_angle * sqrt(-2*log(rndm->getUniform())));
            }
            while (u.z < 0);
            u.y = sqrt((1-u.z)*(1+u.z));
            u.y = (rndm->getUniform() < 0.5) ? u.y : -1*u.y;
            u.x = 0;
            break;
 
        default: // No Deviation
            u.x = 0;
            u.y = 0;
            u.z = 1;
        }
 
        // 3. Rotate X, Y, Z, U, V, W and project them into the z-pos plane
        if (is_rotated) {
            rotation.transform(x);   // position: rotation + pivot translation
            rotation.rotate(u);      // direction: rotation only, no translation
 
            // project back onto z_pos plane (ray-plane intersection)
            double t = (z_pos - x.z) / u.z;
            x.x += t * u.x;
            x.y += t * u.y;
            x.z = z_pos;
        }
    };
 
    EGS_Float getFluence() const {
        return count;
    };
 
    bool storeFluenceState(ostream &) const {
        return true;
    };
 
    bool setFluenceState(istream &) {
        return true;
    };
 
    bool isValid() const {
        return (valid && s != 0);
    };
 
protected:
    static constexpr double PI  = 3.141592653589793;
    static constexpr double PI2 = 6.283185307179586;
    static constexpr double DEGREE_TO_RAD = 0.017453292519943295;
    int angle_mode;
    EGS_AffineTransform rotation;
    void setUp();
 
};
 
#endif