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
#
###############################################################################
*/


#ifndef EGS_BASE_SOURCE_
#define EGS_BASE_SOURCE_

#include "egs_vector.h"
#include "egs_object_factory.h"
#include "egs_functions.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) {};

    EGS_BaseSource(EGS_Input *input, EGS_ObjectFactory *f = 0) :
        EGS_Object(input,f) {};
    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) { };

    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 getMu() {
        return -1;
    };

    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);

protected:

    string description;

};

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;

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