egs_isotropic_source.h

Go to the documentation of this file.

/*
###############################################################################
#
#  EGSnrc egs++ isotropic source headers
#  Copyright (C) 2015 National Research Council Canada
#
#  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:          Iwan Kawrakow, 2005
#
#  Contributors:    Long Zhang
#                   Frederic Tessier
#                   Reid Townson
#                   Ernesto Mainegra-Hing
#                   Hugo Bouchard
#                   Hubert Ho
#                   Marc Chamberland
#                   Reid Townson
#
###############################################################################
*/
 
 
#ifndef EGS_ISOTROPIC_SOURCE_
#define EGS_ISOTROPIC_SOURCE_
 
#include "egs_vector.h"
#include "egs_base_source.h"
#include "egs_rndm.h"
#include "egs_shapes.h"
#include "egs_base_geometry.h"
#include "egs_math.h"
 
 
#ifdef WIN32
 
    #ifdef BUILD_ISOTROPIC_SOURCE_DLL
        #define EGS_ISOTROPIC_SOURCE_EXPORT __declspec(dllexport)
    #else
        #define EGS_ISOTROPIC_SOURCE_EXPORT __declspec(dllimport)
    #endif
    #define EGS_ISOTROPIC_SOURCE_LOCAL
 
#else
 
    #ifdef HAVE_VISIBILITY
        #define EGS_ISOTROPIC_SOURCE_EXPORT __attribute__ ((visibility ("default")))
        #define EGS_ISOTROPIC_SOURCE_LOCAL  __attribute__ ((visibility ("hidden")))
    #else
        #define EGS_ISOTROPIC_SOURCE_EXPORT
        #define EGS_ISOTROPIC_SOURCE_LOCAL
    #endif
 
#endif
 
class EGS_ISOTROPIC_SOURCE_EXPORT EGS_IsotropicSource :
    public EGS_BaseSimpleSource {
 
public:
 
    enum GeometryConfinement {
        IncludeAll      = 0,
        ExcludeAll      = 1,
        IncludeSelected = 2,
        ExcludeSelected = 3
    };
 
    EGS_IsotropicSource(int Q, EGS_BaseSpectrum *Spec, EGS_BaseShape *Shape,
                        EGS_BaseGeometry *geometry,
                        const string &Name="", EGS_ObjectFactory *f=0) :
        EGS_BaseSimpleSource(Q,Spec,Name,f), shape(Shape),
        geom(geometry), regions(0), min_theta(85.), max_theta(95.),
        buf_1(1), buf_2(-1), min_phi(0), max_phi(2*M_PI),
        nrs(0), gc(IncludeAll), media(0), nms(0), gcm(IncludeSelected) {
        setUp();
    };
 
    EGS_IsotropicSource(EGS_Input *, EGS_ObjectFactory *f=0);
    ~EGS_IsotropicSource() {
        EGS_Object::deleteObject(shape);
        if (geom) {
            if (!geom->deref()) {
                delete geom;
            }
        }
        if (nrs > 0 && regions) {
            delete [] regions;
        }
        if (nms > 0 && media) {
            delete [] media;
        }
    };
 
    void getPositionDirection(EGS_RandomGenerator *rndm,
                              EGS_Vector &x, EGS_Vector &u, EGS_Float &wt) {
        bool ok = true;
        bool ok2 = true;
        int ntry = 0;
        do {
            x = shape->getRandomPoint(rndm);
            if (geom) {
                // Filter by region
                if (gc == IncludeAll) {
                    ok = geom->isInside(x);
                }
                else if (gc == ExcludeAll) {
                    ok = !geom->isInside(x);
                }
                else if (gc == IncludeSelected) {
                    ok = false;
                    int ireg = geom->isWhere(x);
                    for (int j=0; j<nrs; ++j) {
                        if (ireg == regions[j]) {
                            ok = true;
                            break;
                        }
                    }
                }
                else if (gc == ExcludeSelected) {
                    ok = true;
                    int ireg = geom->isWhere(x);
                    for (int j=0; j<nrs; ++j) {
                        if (ireg == regions[j]) {
                            ok = false;
                            break;
                        }
                    }
                }
 
                // Filter by medium
                if (ok && nms > 0) {
                    if (gcm == IncludeSelected) {
                        ok2 = false;
                        int ireg = geom->isWhere(x);
                        int currentMedium = geom->medium(ireg);
                        for (int j=0; j<nms; ++j) {
                            if (currentMedium == media[j]) {
                                ok2 = true;
                                break;
                            }
                        }
                    }
                    else if (gcm == ExcludeSelected) {
                        ok2 = true;
                        int ireg = geom->isWhere(x);
                        for (int j=0; j<nrs; ++j) {
                            if (ireg == media[j]) {
                                ok2 = false;
                                break;
                            }
                        }
                    }
                }
 
                // Combine the filtering by region and medium
                if (ok2 == false) {
                    ok = false;
                }
            }
            ntry++;
            if (ntry > 100000) {
                egsFatal("\nEGS_IsotropicSource::getPositionDirection:\n"
                         " For target shape %s, which is of type %s:\n"
                         " Failed to create a particle from the source after 100 000 attempts\n"
                         " Please ensure your .egsinp source definition -> source"
                         " -> shape encompasses the source region.\n",
                         shape->getObjectName().c_str(),
                         shape->getObjectType().c_str());
            }
        }
        while (!ok);
        u.z = buf_1 - rndm->getUniform()*(buf_1 - buf_2);
        EGS_Float sinz = 1-u.z*u.z;
        if (sinz > epsilon) {
            sinz = sqrt(sinz);
            EGS_Float cphi, sphi;
            EGS_Float phi = min_phi +(max_phi - min_phi)*rndm->getUniform();
            cphi = cos(phi);
            sphi = sin(phi);
            u.x = sinz*cphi;
            u.y = sinz*sphi;
        }
        else {
            u.x = 0;
            u.y = 0;
        }
        wt = 1;
    };
 
    EGS_Float getFluence() const {
        return count;
    };
 
    bool storeFluenceState(ostream &) const {
        return true;
    };
 
    bool setFluenceState(istream &) {
        return true;
    };
 
    bool isValid() const {
        return (s != 0 && shape != 0);
    };
 
protected:
 
    EGS_BaseShape    *shape;  
    EGS_BaseGeometry *geom;
    int              *regions, *media;
 
    void setUp();
 
    EGS_Float min_theta, max_theta;
    EGS_Float buf_1, buf_2; 
    EGS_Float min_phi, max_phi;
 
    int                 nrs, nms;
    GeometryConfinement gc, gcm;
};
 
#endif