egs_base_source.h

Go to the documentation of this file.

/*
###############################################################################
#
#  EGSnrc egs++ base 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:    Reid Townson
#                   Ernesto Mainegra-Hing
#                   Blake Walters
#                   Alexandre Demelo
#
###############################################################################
*/
 
 
#ifndef EGS_BASE_SOURCE_
#define EGS_BASE_SOURCE_
 
#include "egs_vector.h"
#include "egs_object_factory.h"
#include "egs_functions.h"
#include "egs_ensdf.h"
#include "egs_input_struct.h"
 
#include <string>
#include <iostream>
#include "egs_math.h"
 
using namespace std;
 
class EGS_Input;
class EGS_RandomGenerator;
 
static shared_ptr<EGS_BlockInput> srcBlockInput = make_shared<EGS_BlockInput>("source");
inline void setBaseSourceInputs(bool isSimpleSource = true, bool includeSpectrumBlock = true) {
    srcBlockInput->addSingleInput("library", true, "The type of source, loaded by shared library in egs++/dso.");
    srcBlockInput->addSingleInput("name", true, "The user-declared unique name of this source. This is the name you may refer to elsewhere in the input file");
 
    if (isSimpleSource) {
        includeSpectrumBlock = true;
        srcBlockInput->addSingleInput("charge", true, "The type of particle to emit from the source, as defined by the charge. Use 0 for photons, -1 for electrons and 1 for positrons.", {"0", "1", "-1"});
    }
    if (includeSpectrumBlock) {
        shared_ptr<EGS_BlockInput> specBlock = srcBlockInput->addBlockInput("spectrum");
        auto typePtr = specBlock->addSingleInput("type", true, "The type of energy distribution for the spectrum.", {"monoenergetic", "Gaussian", "Double Gaussian", "uniform", "tabulated spectrum", "radionuclide"});
 
        // Monoenergetic
        specBlock->addSingleInput("energy", false, "The kinetic energy of the source particles in MeV.")->addDependency(typePtr, "monoenergetic");
 
        // Gaussian & double Gaussian
        specBlock->addSingleInput("mean energy", false, "The mean kinetic energy of the source particles in MeV.")->addDependency(typePtr, "Gaussian");
        auto sigmaPtr = specBlock->addSingleInput("sigma", false, "The sigma of the spectrum. For a double Gaussian, input two values.");
        sigmaPtr->addDependency(typePtr, "Gaussian");
        sigmaPtr->addDependency(typePtr, "Double Gaussian");
        auto fwhmPtr = specBlock->addSingleInput("fwhm", false, "The full-width-at-half-maximum of the spectrum. For a double Gaussian, input two values.");
        fwhmPtr->addDependency(typePtr, "Gaussian");
        fwhmPtr->addDependency(typePtr, "Double Gaussian");
        fwhmPtr->addDependency(sigmaPtr, "", true);
        sigmaPtr->addDependency(fwhmPtr, "", true);
 
        // Uniform
        auto rangePtr = specBlock->addSingleInput("range", false, "The minimum and maximum energy for the spectrum, in MeV.");
        rangePtr->addDependency(typePtr, "uniform");
        auto minEPtr = specBlock->addSingleInput("minimum energy", false, "The minimum energy for the spectrum, in MeV.");
        minEPtr->addDependency(typePtr, "uniform");
        auto maxEPtr = specBlock->addSingleInput("maximum energy", false, "The maximum energy for the spectrum, in MeV.");
        maxEPtr->addDependency(typePtr, "uniform");
        minEPtr->addDependency(rangePtr, "", true);
        maxEPtr->addDependency(rangePtr, "", true);
        rangePtr->addDependency(minEPtr, "", true);
        rangePtr->addDependency(maxEPtr, "", true);
 
        // Tabulated
        auto specFilePtr = specBlock->addSingleInput("spectrum file", false, "The full file path to the spectrum file. See documentation for the format of the file.");
        specFilePtr->addDependency(typePtr, "tabulated spectrum");
        auto modePtr = specBlock->addSingleInput("spectrum mode", false, "The mode number that denotes how to create the spectrum. Use 0 for histogram counts/bin, 1 for counts/MeV, 2 for a line spectrum and 3 for an interpolated spectrum.");
        modePtr->addDependency(typePtr, "tabulated spectrum");
        auto energiesPtr = specBlock->addSingleInput("energies", false, "A list of energies for the spectrum, in MeV. When applicable, this is the upper edge of the bin.");
        energiesPtr->addDependency(typePtr, "tabulated spectrum");
        auto probsPtr = specBlock->addSingleInput("probabilities", false, "A list of probabilities for the spectrum. Does not need to be normalized.");
        probsPtr->addDependency(typePtr, "tabulated spectrum");
        modePtr->addDependency(specFilePtr, "", true);
        energiesPtr->addDependency(specFilePtr, "", true);
        probsPtr->addDependency(specFilePtr, "", true);
 
        // Radionuclide
        specBlock->addSingleInput("nuclide", true, "The name of the nuclide to model, e.g. Co-60 or Tc-99m. If the 'ensdf file' input is not specified, then the ENSDF file will be searched for as $HEN_HOUSE/spectra/lnhb/ensdf/nuclide.txt, where nuclide is the text you input. Note that a radionuclide spectrum is ONLY compatible with a radionuclide source.")->addDependency(typePtr, "radionuclide");
        specBlock->addSingleInput("ensdf file", false, "The full path to the ENSDF file to use.")->addDependency(typePtr, "radionuclide");
        specBlock->addSingleInput("relative activity", false, "If multiple radionuclide spectra are specified for a single radionuclide source, this is the relative weight of this spectrum. Defaults to 1.")->addDependency(typePtr, "radionuclide");
        specBlock->addSingleInput("atomic relaxations", false, "The model to use for atomic relaxations resulting from radionuclide decay. Defaults to EADL.", {"eadl", "ensdf", "off"})->addDependency(typePtr, "radionuclide");
        specBlock->addSingleInput("output beta spectra", false, "Whether or not to output as files the beta spectra that are used for sampling beta decay energies. Defaults to No.", {"yes", "no"})->addDependency(typePtr, "radionuclide");
        specBlock->addSingleInput("alpha scoring", false, "The model to use for scoring alpha particles during radionuclide decay. Defaults to Discard.", {"local", "discard"})->addDependency(typePtr, "radionuclide");
        specBlock->addSingleInput("extra transition approximation", false, "Whether or not to use the option that automatically balances transition intensities. Defaults to Off.", {"on","off"})->addDependency(typePtr, "radionuclide");
    }
}
 
class EGS_EXPORT EGS_BaseSource : public EGS_Object {
 
public:
 
    EGS_BaseSource(const string &Name="", EGS_ObjectFactory *f = 0) :
        EGS_Object(Name,f), time_index(-1) {};
 
    EGS_BaseSource(EGS_Input *input, EGS_ObjectFactory *f = 0) :
        EGS_Object(input,f), time_index(-1)  {};
    virtual ~EGS_BaseSource() {};
 
    const char *getSourceDescription() const {
        return description.c_str();
    };
 
    virtual EGS_I64 getNextParticle(EGS_RandomGenerator *rndm,
                                    int &q, int &latch,                // charge and latch
                                    EGS_Float &E, EGS_Float &wt,       // energy and weight
                                    EGS_Vector &x, EGS_Vector &u) = 0; // position and direction
 
    virtual void setSimulationChunk(EGS_I64 nstart, EGS_I64 nrun, int npar, int nchunk) {
        (void)nstart;
        (void)nrun;
        (void)npar;
        (void)nchunk;
    };
 
    virtual int getCharge() const  {
        return -99;
    };
    virtual EGS_Float getEmax() const = 0;
 
    virtual EGS_Float getFluence() const = 0;
 
    /* A virtual function which can be implemented in derived classes
    * to return a fractional monitor unit associated with each source
    * particle.  Currently only makes sense for IAEA_PhspSource and
    * EGS_BeamSource.
    */
    //virtual EGS_Float getTimeIndex() {
    //return -1;
    //};
 
    //virtual void setTimeIndex(EGS_Float temp_time) {};
 
    virtual void printSampledEmissions() {};
 
    virtual vector<EGS_Ensdf *> getRadionuclideEnsdf() {
        return vector<EGS_Ensdf *>();
    };
 
    virtual bool storeState(ostream &data_out) const {
        (void)data_out;
        return true;
    };
 
    virtual bool setState(istream &data_in) {
        (void)data_in;
        return true;
    };
 
    virtual bool addState(istream &data_in) {
        (void)data_in;
        return true;
    };
 
    virtual void resetCounter() {};
 
    static EGS_BaseSource *createSource(EGS_Input *);
 
    static EGS_BaseSource *getSource(const string &Name);
 
    static void addKnownSource(EGS_BaseSource *o);
 
    static void addKnownTypeId(const char *name);
 
    /* Centralize the time index parameter so that it can be saved to and
     * accessed from a single point. In almost any instance where a time index
     * parameter is created, it is saved in the source object using
     * setTimeIndex. When other objects would like to access the time index
     * (most relevant example being the dynamic geometry checking if the source
     * has provided a time index before setting its own), they can call
     * getTimeIndex. In many cases, this method is indirectly called via the
     * getTimeIndex in the application class, which returns the results of
     * getTimeIndex call on the simulation source. Note there are two cases
     * which behave slightly differently and may need some modifications. While
     * the dynamicSource time index implementation was completely absorbed into
     * the basesource, this was not done for the beam source and the iaea_phsp
     * source, as they had slightly more involved time index implementations
     * that seemed best left alone. They do not have setTimeIndex methods, \
     * and may not return the right thing if we set the time using the geometry,
     * as calling get may try to get the beam or iaea_phsp time index and not
     * the basesource index we set with the geometry. */
    EGS_Float getTimeIndex() {
        return time_index;
 
    };
 
    void setTimeIndex(EGS_Float temp_time) {
        time_index=temp_time;
    };
 
    /* This method is essentially used to determine whether the simulation
     * source contains a dynamic source. The only
     * non-empty implementations of this function are in composite sources
     * (where it simply calls containsDynamic on its components),
     * where it will update the boolean reference to true and
     * call containsDynamic on its base geometry, and sources that
     * can contain time indices (dynamic, phsp, beam sources).
     * This function was conceived to be used in the
     * view/viewcontrol (to determine whether time index objects are visible or
     * hidden), and track scoring */
    virtual void containsDynamic(bool &hasdynamic) {
        (void)hasdynamic;
    };
 
protected:
 
    string description;
 
    EGS_Float time_index;
 
};
 
class EGS_EXPORT EGS_BaseSpectrum {
 
public:
 
    EGS_BaseSpectrum() : count(0), sum_E(0), sum_E2(0),
        type("Unknown spectrum") {};
 
    virtual ~EGS_BaseSpectrum() {};
 
    const string &getType() const {
        return type;
    };
 
    inline EGS_Float sampleEnergy(EGS_RandomGenerator *rndm) {
        EGS_Float e = sample(rndm);
        count++;
        sum_E += e;
        sum_E2 += e*e;
        return e;
    };
 
    virtual EGS_Float maxEnergy() const = 0;
 
    virtual EGS_Float expectedAverage() const = 0;
 
    virtual bool storeState(ostream &data_out) const {
        if (!egsStoreI64(data_out,count)) {
            return false;
        }
        data_out << " " << sum_E << " " << sum_E2 << endl;
        if (!data_out.good() || data_out.fail()) {
            return false;
        }
        return true;
    };
 
    virtual bool setState(istream &data_in) {
        if (!egsGetI64(data_in,count)) {
            return false;
        }
        data_in >> sum_E >> sum_E2;
        if (data_in.eof() || !data_in.good() || data_in.fail()) {
            return false;
        }
        return true;
    };
 
    virtual bool addState(istream &data_in) {
        EGS_I64 count_save = count;
        double sum_E_save = sum_E, sum_E2_save = sum_E2;
        if (!setState(data_in)) {
            return false;
        }
        count += count_save;
        sum_E += sum_E_save;
        sum_E2 += sum_E2_save;
        return true;
    };
 
    virtual void resetCounter() {
        count = 0;
        sum_E = 0;
        sum_E2 = 0;
    };
 
    static EGS_BaseSpectrum *createSpectrum(EGS_Input *inp);
 
    void getSampledAverage(EGS_Float &e, EGS_Float &de) const {
        if (count > 1) {
            e = sum_E/count;
            de = sum_E2/count;
            de -= e*e;
            if (de > 0) {
                de = sqrt(de/(count-1));
            }
        }
    };
 
    void reportAverageEnergy() const {
        egsInformation("expected average energy: %g\n",expectedAverage());
        EGS_Float e=0,de=0;
        getSampledAverage(e,de);
        egsInformation("sampled  average energy: %g +/- %g\n",e,de);
    };
 
protected:
 
    virtual EGS_Float sample(EGS_RandomGenerator *rndm) = 0;
 
    EGS_I64 count;
 
    double  sum_E;
 
    double  sum_E2;
 
    string  type;
 
};
 
class EGS_EXPORT EGS_BaseSimpleSource : public EGS_BaseSource {
 
public:
 
    EGS_BaseSimpleSource(int Q, EGS_BaseSpectrum *Spec,
                         const string &Name="", EGS_ObjectFactory *f=0) :
        EGS_BaseSource(Name,f), q(Q), s(Spec), count(0) { };
 
    EGS_BaseSimpleSource(EGS_Input *input, EGS_ObjectFactory *f=0);
 
    ~EGS_BaseSimpleSource() {
        if (s) {
            delete s;
        }
    };
 
    virtual bool isValid() const {
        return (s != 0);
    };
 
    virtual EGS_I64 getNextParticle(EGS_RandomGenerator *rndm,
                                    int &Q, int &latch, EGS_Float &E, EGS_Float &wt,
                                    EGS_Vector &x, EGS_Vector &u) {
        Q = q;
        E = s->sampleEnergy(rndm);
        getPositionDirection(rndm,x,u,wt);
        setLatch(latch);
        return ++count;
    };
 
    virtual void getPositionDirection(EGS_RandomGenerator *rndm,
                                      EGS_Vector &x, EGS_Vector &u, EGS_Float &wt) = 0;
 
    virtual EGS_Float getEmax() const {
        return s->maxEnergy();
    };
 
    int getCharge() const {
        return q;
    };
 
    virtual bool storeFluenceState(ostream &data_out) const {
        (void)data_out;
        return true;
    };
 
    virtual bool storeState(ostream &data_out) const {
        if (!egsStoreI64(data_out,count)) {
            return false;
        }
        if (!s->storeState(data_out)) {
            return false;
        }
        if (!storeFluenceState(data_out)) {
            return false;
        }
        return true;
    };
 
    virtual bool addState(istream &data) {
        EGS_I64 count_save = count;
        if (!egsGetI64(data,count)) {
            return false;
        }
        if (!s->addState(data)) {
            return false;
        }
        if (!addFluenceData(data)) {
            return false;
        }
        count += count_save;
        return true;
    };
 
    virtual void resetCounter() {
        count = 0;
        s->resetCounter();
        resetFluenceCounter();
    };
 
    virtual bool addFluenceData(istream &data) {
        (void)data;
        return true;
    }
 
    virtual void resetFluenceCounter() { };
 
    virtual bool setFluenceState(istream &data) {
        (void)data;
        return true;
    };
 
    virtual bool setState(istream &data) {
        if (!egsGetI64(data,count)) {
            return false;
        }
        if (!s->setState(data)) {
            return false;
        }
        if (!setFluenceState(data)) {
            return false;
        }
        return true;
    };
 
protected:
 
    virtual void setLatch(int &latch) {
        latch = 0;
    };
 
    int              q;
 
    EGS_BaseSpectrum *s;
 
    string           type;
 
    EGS_I64          count;
 
};
 
template <class T>
EGS_BaseSource *createSourceTemplate(EGS_Input *input,
                                     EGS_ObjectFactory *f, const char *name) {
    EGS_BaseSource::addKnownTypeId(typeid(T).name());
    if (!input) {
        egsWarning("createSource(%s): null input?\n",name);
        return 0;
    }
    T *res = new T(input,f);
    if (!res->isValid())  {
        egsWarning("createSource(%s): the input is not "
                   "sufficient to create a valid source\n",name);
        delete res;
        return 0;
    }
    return res;
};
 
#endif