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 <string>
#include <iostream>
#include "egs_math.h"
 
using namespace std;
 
class EGS_Input;
class EGS_RandomGenerator;
 
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) { };
 
    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 {
        return true;
    };
 
    virtual bool setState(istream &data_in) {
        return true;
    };
 
    virtual bool addState(istream &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) { };
 
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 {
        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) {
        return true;
    }
 
    virtual void resetFluenceCounter() { };
 
    virtual bool setFluenceState(istream &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