egs_spectra.cpp

Go to the documentation of this file.

/*
###############################################################################
#
#  EGSnrc egs++ spectra
#  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:    Frederic Tessier
#                   Reid Townson
#                   Ernesto Mainegra-Hing
#
###############################################################################
*/
 
 
#include "egs_base_source.h"
#include "egs_rndm.h"
#include "egs_alias_table.h"
#include "egs_input.h"
#include "egs_math.h"
 
#include <cstdio>
#include "egs_math.h"
#include <fstream>
#ifndef NO_SSTREAM
    #include <sstream>
    #define S_STREAM std::istringstream
#else
    #include <strstream>
    #define S_STREAM std::istrstream
#endif
#include <limits>
 
 
using namespace std;
 
 
 
class EGS_EXPORT EGS_MonoEnergy : public EGS_BaseSpectrum {
 
public:
 
    EGS_MonoEnergy(EGS_Float energy) : EGS_BaseSpectrum(), E(energy) {
        char buf[1024];
        sprintf(buf,"monoenergetic %g MeV",E);
        type = buf;
    };
    ~EGS_MonoEnergy() {};
    EGS_Float expectedAverage() const {
        return E;
    };
    EGS_Float maxEnergy() const {
        return E;
    };
 
protected:
 
    EGS_Float sample(EGS_RandomGenerator *) {
        return E;
    };
 
    EGS_Float E; 
 
};
 
static inline EGS_Float getGaussianRN(EGS_RandomGenerator *rndm) {
    static bool have_x = false;
    static EGS_Float the_x;
    if (have_x) {
        have_x = false;
        return the_x;
    }
    EGS_Float r = sqrt(-2*log(1-rndm->getUniform()));
    EGS_Float cphi, sphi;
    rndm->getAzimuth(cphi,sphi);
    the_x = r*sphi;
    have_x = true;
    return r*cphi;
};
 
class EGS_EXPORT EGS_GaussianSpectrum : public EGS_BaseSpectrum {
 
public:
 
    EGS_GaussianSpectrum(EGS_Float mean_energy, EGS_Float Sigma) :
        EGS_BaseSpectrum(), Eo(mean_energy), sigma(Sigma) {
        if (Eo <= 0) egsFatal("EGS_GaussianSpectrum: attempt to construct "
                                  "a spectrum with a negative mean energy (%g)\n",Eo);
        if (sigma < 0) {
            sigma = -sigma*0.4246609;    // i.e. assume
        }
        // the user has specified FWHM
        char buf[1024];
        sprintf(buf,"Gaussian spectrum with Eo = %g and sigma = %g",Eo,sigma);
        type = buf;
        if (Eo - 5*sigma > 0) {
            Emax = Eo + 5*sigma;
        }
        else {
            Emax = 2*Eo;
        }
    };
    ~EGS_GaussianSpectrum() {};
    EGS_Float expectedAverage() const {
        return Eo;
    };
    EGS_Float maxEnergy() const {
        return Emax;
    };
 
    EGS_Float sample(EGS_RandomGenerator *rndm) {
        EGS_Float E;
        do {
            E = Eo + sigma*getGaussianRN(rndm);
        }
        while (E <= 0 || E > Emax);
        return E;
    };
 
protected:
 
    EGS_Float  Eo;    
    EGS_Float  sigma; 
    EGS_Float  Emax;
};
 
class EGS_EXPORT EGS_DoubleGaussianSpectrum : public EGS_BaseSpectrum {
 
public:
 
    EGS_DoubleGaussianSpectrum(EGS_Float mean_energy, EGS_Float sig_left,
                               EGS_Float sig_right) : EGS_BaseSpectrum(), Eo(mean_energy),
        sleft(sig_left), sright(sig_right) {
        if (sleft < 0) {
            sleft = -sleft*0.4246609;
        }
        if (sright< 0) {
            sright= -sright*0.4246609;
        }
        Emax = Eo + 4*sright;
        if (Eo - 4*sleft < 0) egsWarning("EGS_DoubleGaussianSpectrum: "
                                             "for Eo=%g, sigma=%g there will be negative energy sampled\n");
        p = sleft/(sleft + sright);
        char buf[1024];
        sprintf(buf,"Double Gaussian spectrum with Eo = %g sig(left) = %g"
                " and sig(right) = %g",Eo,sleft,sright);
        type = buf;
    };
    ~EGS_DoubleGaussianSpectrum() {};
    EGS_Float maxEnergy() const {
        return Emax;
    };
    EGS_Float expectedAverage() const {
        return Eo + sqrt(2/M_PI)*(sright-sleft);
    };
 
protected:
 
    EGS_Float Eo;     
    EGS_Float Emax;   
    EGS_Float sleft;  
    EGS_Float sright; 
    EGS_Float p;      
    EGS_Float sample(EGS_RandomGenerator *rndm) {
        EGS_Float E;
        if (rndm->getUniform() < p) {
            do {
                E = Eo-sleft*fabs(getGaussianRN(rndm));
            }
            while (E <= 0);
        }
        else {
            do {
                E = Eo+sright*fabs(getGaussianRN(rndm));
            }
            while (E > Emax);
        }
        return E;
    };
 
};
 
class EGS_EXPORT EGS_UniformSpectrum : public EGS_BaseSpectrum {
 
public:
 
    EGS_UniformSpectrum(EGS_Float emin, EGS_Float emax) : EGS_BaseSpectrum(),
        Emin(emin), Emax(emax), de(emax - emin) {
        char buf[1024];
        sprintf(buf,"Uniform spectrum with Emin = %g Emax = %g",Emin,Emax);
        type = buf;
    };
    ~EGS_UniformSpectrum() {};
    EGS_Float maxEnergy() const {
        return Emax;
    };
    EGS_Float expectedAverage() const {
        return (Emin+Emax)/2;
    };
 
protected:
 
    EGS_Float sample(EGS_RandomGenerator *rndm) {
        return Emin + rndm->getUniform()*de;
    };
 
    EGS_Float Emin;  
    EGS_Float Emax;  
    EGS_Float de;    
};
 
class EGS_EXPORT EGS_TabulatedSpectrum : public EGS_BaseSpectrum {
 
public:
 
    EGS_TabulatedSpectrum(int N, const EGS_Float *x, const EGS_Float *f,
                          int Type = 1, const char *fname = 0) : EGS_BaseSpectrum(),
        table(new EGS_AliasTable(N,x,f,Type)) {
        setType(Type,fname);
    };
 
    ~EGS_TabulatedSpectrum() {
        delete table;
    };
    EGS_Float maxEnergy() const {
        return table->getMaximum();
    };
    EGS_Float expectedAverage() const {
        return table->getAverage();
    };
 
protected:
 
    EGS_AliasTable *table; 
    void setType(int Type,const char *fname) {
        if (Type == 0) {
            type = "tabulated line spectrum";
        }
        else if (Type == 1) {
            type = "tabulated histogram spectrum";
        }
        else {
            type = "tabulated spectrum";
        }
        if (fname) {
            type += " defined in ";
            type += fname;
        }
        else {
            type += " defined inline";
        }
    };
 
    EGS_Float sample(EGS_RandomGenerator *rndm) {
        return table->sample(rndm);
    };
 
};
 
//
// The following function is used to skip potentially present
// commas as separators in spectrum files. With it we can write
//   stream >> input1 >> skipsep >> input2;
// and it will work for white space and for comma separated input.
//
istream &skipsep(istream &in) {
    char c;
    for (EGS_I64 loopCount=0; loopCount<=loopMax; ++loopCount) {
        if (loopCount == loopMax) {
            egsFatal("skipsep: Too many iterations were required! Input may be invalid, or consider increasing loopMax.");
            return in;
        }
        in.get(c);
        if (in.eof() || in.fail() || !in.good()) {
            break;
        }
        if (c == ',') {
            break;
        }
        if (!isspace(c)) {
            in.putback(c);
            break;
        }
    }
    return in;
}
 
static char spec_msg1[] = "EGS_BaseSpectrum::createSpectrum:";
 
EGS_BaseSpectrum *EGS_BaseSpectrum::createSpectrum(EGS_Input *input) {
    if (!input) {
        egsWarning("%s got null input?\n",spec_msg1);
        return 0;
    }
    EGS_Input *inp = input;
    bool delete_it = false;
    if (!input->isA("spectrum")) {
        inp = input->takeInputItem("spectrum");
        if (!inp) {
            egsWarning("%s no 'spectrum' input!\n",spec_msg1);
            return 0;
        }
        delete_it = true;
    }
    string stype;
    int err = inp->getInput("type",stype);
    if (err) {
        egsWarning("%s wrong/missing 'type' input\n",spec_msg1);
        if (delete_it) {
            delete inp;
        }
        return 0;
    }
    EGS_BaseSpectrum *spec = 0;
    if (inp->compare(stype,"monoenergetic")) {
        EGS_Float Eo;
        err = inp->getInput("energy",Eo);
        if (err) egsWarning("%s wrong/missing 'energy' input for a "
                                "monoenergetic spectrum\n",spec_msg1);
        else {
            spec = new EGS_MonoEnergy(Eo);
        }
    }
    else if (inp->compare(stype,"Gaussian")) {
        EGS_Float Eo, sig, fwhm;
        int err1 = inp->getInput("mean energy",Eo);
        int err2 = inp->getInput("sigma",sig);
        int err3 = inp->getInput("fwhm",fwhm);
        if (err1 || (err2 && err3)) {
            if (err1) egsWarning("%s wrong/missing 'mean energy' input"
                                     " for a Gaussian spectrum\n",spec_msg1);
            else egsWarning("%s wrong/missing 'sigma' and 'FWHM' input for a "
                                "Gaussian spectrum\n",spec_msg1);
 
        }
        else {
            if (Eo <= 0) egsWarning("%s mean energy must be positive but your"
                                        " input was %g\n",Eo);
            else {
                if (!err2) {
                    if (sig <= 0) egsWarning("%s sigma must be positive"
                                                 " but your input was %g\n",spec_msg1,sig);
                    else {
                        spec = new EGS_GaussianSpectrum(Eo,sig);
                    }
                }
                else {
                    if (fwhm <= 0)  egsWarning("%s FWHM must be positive"
                                                   " but your input was %g\n",spec_msg1,fwhm);
                    else {
                        spec = new EGS_GaussianSpectrum(Eo,-fwhm);
                    }
                }
            }
        }
    }
    else if (inp->compare(stype,"Double Gaussian")) {
        EGS_Float Eo;
        vector<EGS_Float> sig, fwhm;
        int err1 = inp->getInput("mean energy",Eo);
        int err2 = inp->getInput("sigma",sig);
        int err3 = inp->getInput("fwhm",fwhm);
        if (!err1 && Eo <= 0) {
            err1 = 1;
        }
        if (!err2 && sig.size() != 2) {
            err2 = 1;
        }
        if (!err2 && (sig[0] <= 0 || sig[1] <= 0)) {
            err2 = 1;
        }
        if (!err3 && fwhm.size() != 2) {
            err3 = 1;
        }
        if (!err3 && (fwhm[0] <= 0 || fwhm[1] <= 0)) {
            err3 = 1;
        }
        if (err1 || (err2 && err3)) {
            if (err1) egsWarning("%s wrong/missing 'mean energy' input"
                                     " for a Double Gaussian spectrum\n",spec_msg1);
            if (err2 && err3) egsWarning("%s wrong/missing 'sigma' and 'FWHM'"
                                             " input for a Double Gaussian spectrum\n",spec_msg1);
        }
        else {
            if (!err2 && !err3) egsWarning("%s found 'sigma' and 'FWHM' "
                                               "input, using 'sigma'\n",spec_msg1);
            if (!err2) {
                spec = new EGS_DoubleGaussianSpectrum(Eo,sig[0],sig[1]);
            }
            else {
                spec = new EGS_DoubleGaussianSpectrum(Eo,-fwhm[0],-fwhm[1]);
            }
        }
    }
    else if (inp->compare(stype,"uniform")) {
        vector<EGS_Float> range;
        EGS_Float Emin, Emax;
        int err1 = inp->getInput("range",range);
        int err2 = inp->getInput("minimum energy",Emin);
        int err3 = inp->getInput("maximum energy",Emax);
        if (!err2 && !err3 && Emin > Emax) {
            egsWarning("%s Emin (%g) is greater than Emax (%g)?\n",
                       spec_msg1,Emin,Emax);
            err2 = 1;
            err3 = 1;
        }
        if (err1 && err2 && err3) egsWarning("%s wrong/missing 'range' and"
                                                 " 'minimum/maximum energy' input\n",spec_msg1);
        else {
            if (!err2 && !err3) {
                spec = new EGS_UniformSpectrum(Emin,Emax);
            }
            else {
                if (range[0] < range[1]) {
                    spec = new EGS_UniformSpectrum(range[0],range[1]);
                }
                else {
                    spec = new EGS_UniformSpectrum(range[1],range[0]);
                }
            }
        }
    }
    else if (inp->compare(stype,"tabulated spectrum")) {
        string spec_file;
        err = inp->getInput("spectrum file",spec_file);
        // Expands FIRST environment variable found in spec_file
        spec_file = egsExpandPath(spec_file);
        if (!err) {
            ifstream sdata(spec_file.c_str());
            if (!sdata) egsWarning("%s failed to open spectrum file %s\n",
                                       spec_msg1,spec_file.c_str());
            else {
                char title[1024];
                sdata.getline(title,1023);
                if (sdata.eof() || sdata.fail() || !sdata.good()) {
                    egsWarning("%s error while reading title of spectrum file"
                               "%s\n",spec_msg1,spec_file.c_str());
                    if (delete_it) {
                        delete inp;
                    }
                    return 0;
                }
                if (sdata.eof() || sdata.fail() || !sdata.good()) {
                    egsWarning("%s error while reading spectrum type and "
                               "number of bins in spectrum file %s\n",
                               spec_msg1,spec_file.c_str());
                    if (delete_it) {
                        delete inp;
                    }
                    return 0;
                }
                EGS_Float dum;
                int nbin, mode;
                sdata >> nbin >> skipsep >> dum >> skipsep >> mode;
                if (sdata.eof() || sdata.fail() || !sdata.good()) {
                    egsWarning("%s error while reading spectrum type and "
                               "number of bins in spectrum file %s\n",
                               spec_msg1,spec_file.c_str());
                    if (delete_it) {
                        delete inp;
                    }
                    return 0;
                }
                if (nbin < 2) {
                    egsWarning("%s nbin in a spectrum must be at least 2\n"
                               "  you have %d in the spectrum file %s\n",
                               spec_msg1,nbin,spec_file.c_str());
                    if (delete_it) {
                        delete inp;
                    }
                    return 0;
                }
                if (mode < 0 || mode > 3) {
                    egsWarning("%s unknown spectrum type %d in spectrum file"
                               " %s\n",spec_msg1,mode,spec_file.c_str());
                    if (delete_it) {
                        delete inp;
                    }
                    return 0;
                }
                EGS_Float *en_array, *f_array;
                int ibin;
                f_array = new EGS_Float [nbin];
                if (mode == 0 || mode == 1) {
                    en_array = new EGS_Float [nbin+1];
                    en_array[0] = dum;
                    ibin=1;
                }
                else {
                    en_array = new EGS_Float [nbin];
                    ibin=0;
                }
                for (int j=0; j<nbin; j++) {
                    sdata >> en_array[ibin++] >> skipsep >> f_array[j];
                    if (sdata.eof() || sdata.fail() || !sdata.good()) {
                        egsWarning("%s error on line %d in spectrum file %s\n",
                                   spec_msg1,j+2,spec_file.c_str());
                        if (delete_it) {
                            delete inp;
                        }
                        delete [] en_array;
                        delete [] f_array;
                        return 0;
                    }
                    if (mode != 2 && ibin > 1) {
                        if (en_array[ibin-1] <= en_array[ibin-2]) {
                            egsWarning("%s energies must be in increasing "
                                       "order.\n   This is not the case for input on "
                                       "lines %d,%d in spectrum file %s\n",
                                       spec_msg1,j+2,j+1,spec_file.c_str());
                            if (delete_it) {
                                delete inp;
                            }
                            return 0;
                        }
                    }
                    if (mode == 0) {
                        f_array[j]/=(en_array[ibin-1]-en_array[ibin-2]);
                    }
                }
                int itype = 1;
                if (mode == 2) {
                    itype = 0;
                }
                else if (mode == 3) {
                    itype = 2;
                }
                int nb = itype == 1 ? nbin+1 : nbin;
                spec = new EGS_TabulatedSpectrum(nb,en_array,f_array,itype,
                                                 spec_file.c_str());
                delete [] en_array;
                delete [] f_array;
            }
        }
        else {
            vector<EGS_Float> eners, probs;
            int itype=1, mode;
            int err1 = inp->getInput("energies",eners);
            int err2 = inp->getInput("probabilities",probs);
            int err3 = inp->getInput("spectrum mode",mode);      // according to EGSnrc convention
            if (err3) {
                err3 = inp->getInput("spectrum type",mode);      // deprecated
            }
            if (err3) {
                egsWarning("%s wrong/missing 'spectrum mode' input\n",spec_msg1);
                if (delete_it) {
                    delete inp;
                }
                return 0;
            }
            else {
                if (mode < 0 || mode > 3) {
                    egsWarning("%s unknown spectrum 'mode' %d"
                               " %s\n",spec_msg1,mode,spec_file.c_str());
                    if (delete_it) {
                        delete inp;
                    }
                    return 0;
                }
                if (mode == 2) {
                    itype = 0;
                }
                else if (mode == 3) {
                    itype = 2;
                }
            }
            if (err1 || err2) {
                if (err1) egsWarning("%s wrong/missing 'energies' input\n",
                                         spec_msg1);
                if (err2) egsWarning("%s wrong/missing 'probabilities' "
                                         "input\n",spec_msg1);
            }
            else {
                if (itype == 1 && probs.size() != eners.size()-1)
                    egsWarning("%s for spectrum type 1 the number of energies"
                               " must be the number of probabilities + 1\n",
                               spec_msg1);
                else if ((itype == 0 || itype == 2) &&
                         probs.size() != eners.size())
                    egsWarning("%s for spectrum types 0 and 2 the number of "
                               "energies must be equal to the number of probabilities\n",
                               spec_msg1);
                else {
                    int nbin = eners.size();
                    int nbin1 = itype == 1 ? nbin-1 : nbin;
                    EGS_Float *x = new EGS_Float [nbin],
                    *f = new EGS_Float [nbin1];
                    int ibin = 0;
                    if (itype == 1) {
                        ibin = 1;
                        x[0] = eners[0];
                    }
                    for (int j=0; j<nbin1; j++) {
                        x[ibin] = eners[ibin];
                        ibin++;
                        f[j] = mode == 0 ? probs[j]/(eners[ibin-1]-eners[ibin-2]) : probs[j];
                        if (itype != 0 && ibin > 1) {
                            if (x[ibin-1] <= x[ibin-2]) {
                                egsWarning("%s energies must be given in "
                                           "increasing order\n  This is not the case"
                                           " for inputs %d and %d (%g,%g) %d\n",
                                           spec_msg1,ibin-2,ibin-1,x[ibin-2],x[ibin-1],
                                           j);
                                if (delete_it) {
                                    delete inp;
                                }
                                delete [] x;
                                delete [] f;
                                return 0;
                            }
                        }
                    }
                    spec = new EGS_TabulatedSpectrum(nbin,x,f,itype,0);
                    delete [] x;
                    delete [] f;
                }
            }
        }
    }
    else {
        egsWarning("%s unknown spectrum type %s\n",spec_msg1,stype.c_str());
    }
    if (delete_it) {
        delete inp;
    }
    return spec;
}