egs_radionuclide_source.h

Go to the documentation of this file.

/*
###############################################################################
#
#  EGSnrc egs++ radionuclide 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:          Reid Townson, 2016
#
#  Contributors:    Martin Martinov
#                   Patrick Saull
#
###############################################################################
*/
 
 
#ifndef EGS_RADIONUCLIDE_SOURCE_
#define EGS_RADIONUCLIDE_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"
#include "egs_application.h"
#include "egs_ensdf.cpp"
 
#include <algorithm>
#include <complex>
 
 
#ifdef WIN32
 
    #ifdef BUILD_RADIONUCLIDE_SOURCE_DLL
        #define EGS_RADIONUCLIDE_SOURCE_EXPORT __declspec(dllexport)
    #else
        #define EGS_RADIONUCLIDE_SOURCE_EXPORT __declspec(dllimport)
    #endif
    #define EGS_RADIONUCLIDE_SOURCE_LOCAL
 
#else
 
    #ifdef HAVE_VISIBILITY
        #define EGS_RADIONUCLIDE_SOURCE_EXPORT __attribute__ ((visibility ("default")))
        #define EGS_RADIONUCLIDE_SOURCE_LOCAL  __attribute__ ((visibility ("hidden")))
    #else
        #define EGS_RADIONUCLIDE_SOURCE_EXPORT
        #define EGS_RADIONUCLIDE_SOURCE_LOCAL
    #endif
 
#endif
 
class EGS_EXPORT EGS_RadionuclideBetaSpectrum {
 
public:
    EGS_RadionuclideBetaSpectrum(EGS_Ensdf *decays, const string outputBetaSpectra) {
 
        EGS_Application *app = EGS_Application::activeApplication();
        rm = app->getRM();
 
        vector<BetaRecordLeaf *> myBetas = decays->getBetaRecords();
 
        for (vector<BetaRecordLeaf *>::iterator beta = myBetas.begin();
                beta != myBetas.end(); beta++) {
 
            // Skip electron capture records
            if ((*beta)->getCharge() == 1 &&
                    (*beta)->getPositronIntensity() == 0) {
                continue;
            }
 
            unsigned short int daughterZ = (*beta)->getZ();
 
            egsInformation("EGS_RadionuclideBetaSpectrum: "
                           "Energy, Z, A, forbidden: %f %d %d %d\n",
                           (*beta)->getFinalEnergy(), daughterZ,
                           (*beta)->getAtomicWeight(), (*beta)->getForbidden()
                          );
 
            const int nbin=1000;
            EGS_Float *e = new EGS_Float [nbin];
            EGS_Float *spec = new EGS_Float [nbin];
            EGS_Float *spec_y = new EGS_Float [nbin];
 
            double de, s_y, factor, se_y;
 
            ncomps=1; // if we increase this, then we must fill the remainder
            area[0]=1.0;
            rel[0]=1.0;
 
            emax = (*beta)->getFinalEnergy();
            zzz[0] = (double)daughterZ;
            rmass = (*beta)->getAtomicWeight();
 
            // These are some special cases where fudge factors are used!
            // Specify lamda[0]=4 for special shape factors
 
            // For Cl-36 (ref: nuc. phys. 99a,  625,(67))
            // Only the beta- spectrum
            if (daughterZ == 18 && (*beta)->getCharge() == -1) {
                lamda[0] = 4;
            }
            // For I-129 (ref: phys. rev. 95, 458, 54))
            // The beta- spectrum with 151 keV endpoint
            else if (daughterZ == 54 && emax < 0.154 && emax > 0.150) {
                lamda[0] = 4;
            }
            // For Cs-137 (ref: nuc. phys. 112a, 156, (68))
            // The beta- spectrum with 1175 keV endpoint
            else if (daughterZ == 56 && emax > 1.173 && emax < 1.177) {
                lamda[0] = 4;
            }
            // For Tl-204 (ref: can. j. phys., 45, 2621, (67))
            // There is only 1 beta- spectrum
            else if (daughterZ == 82) {
                lamda[0] = 4;
            }
            // For Bi-210 (ref: nuc. phys., 31, 293, (62))
            // There is only 1 beta- spectrum
            else if (daughterZ == 84) {
                lamda[0] = 4;
            }
            else {
                lamda[0] = (*beta)->getForbidden();
            }
 
            // For positrons from zzz negative (just how the spectrum code
            // was designed)
            if ((*beta)->getCharge() == 1) {
                zzz[0] *= -1;
            }
 
            etop[0]=emax;
 
            // prbs july 9, 2007 moved here from before src loop.
            // also, now tabulate based on
            // endpoint, using roughly nbin bins across spectrum. Do this by
            // rounding
            // endpoint E0 up to nearest 100 keV, and dividing by NBIN to get
            // binwidth.
 
            de=((int)(etop[0]*10.0+1)/10.)/nbin; // round up to nearest 100kev;
            // /=     NBIN
            //cout << "Binwidth " << de << endl;
 
            for (int ib=0; ib<nbin; ib++) {
                e[ib]=de+ib*de;
//                 egsInformation("%.12f, %.12f\n", e[ib], etop[0]);
            }
 
            s_y=0.0;
            se_y=0.0;
            for (int ib=0; ib<nbin; ib++) {
 
                if (e[ib]<=emax) {
                    sp(e[ib],spec_y[ib],factor);
                }
                else {
                    spec_y[ib]=0.0;
                }
 
                s_y=s_y+spec_y[ib];
                se_y=se_y+spec_y[ib]*e[ib];
            }
 
            for (int ib=0; ib<nbin; ib++) {
                spec[ib]=1/de*(spec_y[ib]/s_y);
//                 cout << e[ib] << " " << spec[ib] << endl;
            }
 
            EGS_AliasTable *bspec = new EGS_AliasTable(nbin,e,spec,1);
            (*beta)->setSpectrum(bspec);
 
            // Write the spectrum to a file
            if (outputBetaSpectra == "yes") {
 
                string fname(app->getOutputFile());
                fname += "_" + decays->radionuclide + "_" + std::to_string(emax);
                if (!egsIsAbsolutePath(fname)) {
                    fname = egsJoinPath(app->getAppDir(),fname);
                }
                fname += ".spec";
 
 
                egsInformation("EGS_RadionuclideBetaSpectrum: Outputting beta spectrum to file: %s\n", fname.c_str());
 
                ofstream specStream;
                specStream.open(fname.c_str());
                for (int ib=0; ib<nbin; ib++) {
                    spec[ib]=1/de*(spec_y[ib]/s_y);
                    specStream << e[ib] << " " << spec[ib] << endl;
                }
                specStream.close();
            }
        }
    }
 
protected:
 
    complex<double> cgamma(complex<double> z) {
 
        static const int g=7;
        static const double pi = 3.1415926535897932384626433832795028841972;
        static const double p[g+2] = {0.99999999999980993, 676.5203681218851,
                                      -1259.1392167224028,
                                      771.32342877765313,
                                      -176.61502916214059,
                                      12.507343278686905,
                                      -0.13857109526572012,
                                      9.9843695780195716e-6,
                                      1.5056327351493116e-7
                                     };
 
        if (real(z)<0.5) {
            return pi / (sin(pi*z)*cgamma(double(1.0)-z));
        }
 
        z -= 1.0;
        complex<double> x=p[0];
        for (int i=1; i<g+2; i++) {
            x += p[i]/(z+complex<double>(i,0));
        }
        complex<double> t = z + (g + double(0.5));
 
        return double(sqrt(2.*pi)) * pow(t,z+double(0.5)) * exp(-t) * x;
    }
 
    complex<double> clgamma(complex<double> z) {
        complex<double> u, v, h, p, r;
 
        static const double pi = 3.1415926535897932384626433832795028841972;
        static const double c1 = 9.189385332046727e-1;
        static const double c2 = 1.144729885849400;
        static const double c[10] = {8.333333333333333e-2,
                                     -2.777777777777777e-3,
                                     7.936507936507936e-4,
                                     -5.952380952380952e-4,
                                     8.417508417508417e-4,
                                     -1.917526917526917e-3,
                                     6.410256410256410e-3,
                                     -2.955065359477124e-2,
                                     1.796443723688305e-1,
                                     -1.392432216905901
                                    };
 
        static const double hf = 0.5;
 
        double x = real(z);
        double y = imag(z);
        h = 0;
 
        if (y == 0 && -abs(x) == int(x)) {
            return 0;
        }
        else {
            double ya = abs(y);
            if (x < 0) {
                u = double(1.) - complex<double>(x, ya);
            }
            else {
                u = complex<double>(x, ya);
            }
 
            h = 0;
            double ur = real(u);
            double ui, a;
            if (ur < 7.) {
                ui = imag(u);
                a = atan2(ui,ur);
                h = u;
                for (int i=1; i<=6-int(ur); i++) {
                    ur = ur + 1;
                    u = complex<double>(ur, ui);
                    h = h * u;
                    a = a + atan2(ui, ur);
                }
                h = complex<double>(hf * log(pow(real(h),2) + pow(imag(h),2)),
                                    a);
 
                u = double(1.) + u;
            }
 
            r = double(1.) / pow(u,2);
            p = r * c[9];
 
            for (int i=8; i>=1; i--) {
                p = r * (c[i] + p);
            }
 
            h = c1 + (u-hf)*log(u) - u + (c[0]+p) / u - h;
 
            if (x < 0.) {
                ur = double(int(x)) - 1.;
                ui = pi * (x-ur);
                x = pi * ya;
                double t = exp(-x-x);
                a = sin(ui);
                t = x + hf * log(t*pow(a,2)+pow(hf*(1.-t),2));
                a = atan2(cos(ui)*tanh(x),a) - ur*pi;
                h = c2 - complex<double>(t,a) - h;
            }
            if (y < 0) {
                h = conj(h);
            }
        }
 
        return h;
    }
 
    void slfact(double p, double z, double radf, double xl[4]) {
 
        double ff[4];
        double dfac[4]= {1.0, 3.0, 15.0, 105.0};
        double pi,c137,az,w,rad,pr,y,x1,gk,bb,cc,dd,x2;
 
        complex<double> aa;
 
        pi  = acos(-1.0);
        c137= 137.036; // 1/ fine structure constant
        az  = z/c137;
        w   = sqrt(p*p+1.0);
        rad = radf/386.159;
        pr  = p*rad;
        y   = az*w/p;
 
        for (int k=1; k<=4; k++) {
 
            gk = sqrt(k*k-az*az);
            x1 = pow((pow(p,k-1)/dfac[k-1]),2);
 
            aa = clgamma(complex<double>(gk,y));
            double aa_real = real(aa);
 
            bb=lgamma((double)k);
            cc=lgamma(2.0*k +1.0);
            dd=lgamma(2.0*gk+1.0);
 
            ff[k-1] =
                pow(2.0*pr, 2.0*(gk-k)) *
                exp(pi*y+2.0*(aa_real+cc-bb-dd)) *
                (k+gk)/(2.0*k);
 
            x2 =
                1.0 -
                az*pr*(2.0*w*(2.0*k+1.0)/(p*k*(2.0*gk+1.0)) -
                       2.0*p*gk/(w*k*(2.0*gk+1.0))) -
                2.0*k*pr*pr/((2.0*k+1.0)*(k+gk));
            xl[k-1] = x1*ff[k-1]/ff[0]*x2;
        }
 
        return;
    }
 
    void bsp(double e, double &bspec, double &factor) {
 
        // *****************************************************************
        //     Calculates n(e) (unnormalized) for one spectral component,
        //     specified by:zz,emax,lamda, where
        //
        //     lamda = 1  first  forbidden
        //     lamda = 2  second forbidden
        //     lamda = 3  third  forbidden
        //     lamda = 0  otherwise
        //     lamda = 4  for a few nuclides whose experimental shape don't
        //                fit theory. fudge factors are applied to the allowed
        //                shape.
        //     zz         charge of the daughter nucleus
        //     emax       maximum beta energy in mev.
        //     w and tsq   in mc**2 units
        //     e and emax  in mev   units.
        //     v is the screening correction for atomic electrons
        //     for the thomas-fermi model of the atom
        //     v = 1.13 * (alpha)**2  * z**(4/3)
        //
        //
 
        double pi,c137,zab,v,z,x,w,psq,p,y,qsq,g,cab,f,radf;
 
        double xl[4];
        complex<double> c;
        complex<double> a;
 
        bspec=0.0;
        if (e>emax) {
            return;
        }
 
        pi  =acos(-1.);
        c137=137.036;           // 1/ fine structure constant
 
        zab = abs(zz);
        v   = 1.13*pow(zab,1.333)/pow(c137,2); // Screening correction
        v   = copysign(v,zz);
        z   = zab/c137;
        x   = sqrt(1.0-z*z);        // s parameter
        w   = 1.0+(e/rm)-v;    // Total energy of b particle
        if (w<1.0000001) {
            bspec = 0.;
            return;
        }
        //if(w<1.00001) w=1.00001;
        psq = w*w-double(1.0);
        p   = sqrt(psq);          // Momemtum of beta particle
        y   = z*w/p;              // eta = alpha * z * e / p
        y   = copysign(y,zz);
        qsq = 3.83*pow(emax-e,2);
 
        if (e <= 1.0e-5) {
            g=0.0;              // Low energy approximation
            if (zz>=0.0) {
                g = qsq*2.0*pi*pow(z,(2.0*x-1.0));
            }
        }
        else {
            a   = complex<double>(x,y);
            c   = cgamma(a);
            cab = abs(c);
            f   = pow(psq,x-1.0)*exp(pi*y)*pow(cab,2);
            g   = f*p*w*qsq;
        }
 
        factor  = 1.0; // Necessary to calculate kurie plot (not done)
        bspec   = g;
        if (lam == 0) {
            return;
        }
 
        radf = 1.2*pow(rmass,0.333); // Nuclear radius
        slfact(p,zz,radf,xl);
 
        if (lam==1) {
            bspec = g*(qsq*xl[0]+9.0*xl[1]);
            return;
        }
        else if (lam==2) {
            bspec = g*(pow(qsq,2)*xl[0]+30.0*qsq*xl[1]+225.0*xl[2]);
            return;
        }
        else if (lam==3) {
            bspec = g*(pow(qsq,3.0)*xl[0]+63.0*pow(qsq,2)*xl[1]+
                       1575.0*qsq*xl[2] + 11025.0*xl[3]);
            return;
        }
        else { // lam==4
 
            // Fudge factors for nuclides whose experimental spectra don't
            // seem to fit theory.
 
            //     for cl36 (ref: nuc. phys. 99a,  625,(67))
            if (zab == 18.0) {
                bspec = bspec*(qsq*xl[0]+20.07*xl[1]);
            }
 
            //     for i129 (ref: phys. rev. 95, 458, 54))
            if (zab == 54.) {
                bspec = bspec*(psq+10.0*qsq);
            }
 
            //     for cs-ba137 (ref: nuc. phys. 112a, 156, (68))
            if (zab == 56.) {
                bspec = bspec*(qsq*xl[0]+0.045*xl[1]);
            }
 
            //     for tl204 (ref: can. j. phys., 45, 2621, (67))
            if (zab == 82.) {
                bspec = bspec*(1.0-1.677*e+ 2.77*e*e);
            }
 
            //     for bi210 (ref: nuc. phys., 31, 293, (62))
            if (zab == 84.) {
                bspec = bspec*(1.78-2.35*e+e*e);
            }
 
            return;
        }
    }
 
    // Sums weighted, normalized spectral components to give total spectrum
    void sp(double e, double &spec, double &factor) {
 
        double bspec;
 
        spec=0.0;
        for (int icomp=0; icomp<ncomps; icomp++) {
            zz  = zzz[icomp];
            emax= etop[icomp];
            lam = lamda[icomp];
            bsp(e,bspec,factor);
            spec= spec+bspec*rel[icomp]/area[icomp];
        }
    }
 
private:
    EGS_Float rm;
    double zz,emax,rmass;
    double zzz[9],etop[9],rel[9],area[9],lamda[9];
    int lam, ncomps;
};
 
class EGS_EXPORT EGS_RadionuclideSpectrum {
 
public:
    EGS_RadionuclideSpectrum(const string nuclide, const string ensdf_file,
                             const EGS_Float relativeActivity, const string relaxType, const string outputBetaSpectra, const bool scoreAlphasLocally, const bool allowMultiTransition);
 
    ~EGS_RadionuclideSpectrum() {
        if (decays) {
            delete decays;
        }
        if (betaSpectra) {
            delete betaSpectra;
        }
    };
 
    EGS_Float maxEnergy() const {
        return Emax;
    };
 
    int getCharge() const {
        return currentQ;
    }
 
    double getTime() const {
        return currentTime;
    }
 
    EGS_I64 getShowerIndex() const {
        return ishower;
    }
 
    EGS_Float getEdep() const {
        return edep;
    }
 
    EGS_Float getSpectrumWeight() const {
        return spectrumWeight;
    }
 
    void setSpectrumWeight(EGS_Float newWeight) {
        spectrumWeight = newWeight;
    }
 
    unsigned int getEmissionType() const {
        return emissionType;
    }
 
    void printSampledEmissions();
 
    EGS_Ensdf *getRadionuclideEnsdf() {
        return decays;
    }
 
    bool storeState(ostream &data) const {
        return egsStoreI64(data,ishower);
    }
 
    bool setState(istream &data) {
        return egsGetI64(data,ishower);
    }
 
    void resetCounter() {
        currentLevel = 0;
        currentTime = 0;
        ishower = -1;
        totalGammaEnergy = 0;
    }
 
    EGS_Float sample(EGS_RandomGenerator *rndm);
 
private:
 
    EGS_Ensdf                   *decays;
    vector<BetaRecordLeaf *>    myBetas;
    vector<AlphaRecord *>       myAlphas;
    vector<GammaRecord *>       myGammas,
           myMetastableGammas,
           myUncorrelatedGammas;
    vector<LevelRecord *>       myLevels;
    vector<double>              xrayIntensities,
           xrayEnergies,
           augerIntensities,
           augerEnergies;
    vector<EGS_I64>             numSampledXRay,
           numSampledAuger;
    vector<const LevelRecord *> multiTransitions;
    EGS_SimpleContainer<EGS_RelaxationParticle> relaxParticles;
    const LevelRecord           *currentLevel;
    int                         currentQ;
    unsigned int                emissionType;
    EGS_Float                   currentTime,
                                Emax,
                                spectrumWeight,
                                totalGammaEnergy,
                                edep;
    EGS_I64                     ishower;
    string                      relaxationType;
    bool                        scoreAlphasLocal;
 
    EGS_RadionuclideBetaSpectrum *betaSpectra;
    EGS_Application             *app;
};
 
class EGS_RADIONUCLIDE_SOURCE_EXPORT EGS_RadionuclideSource :
    public EGS_BaseSource {
 
public:
 
    EGS_RadionuclideSource(EGS_Input *, EGS_ObjectFactory *f=0);
 
    ~EGS_RadionuclideSource() {
        if (baseSource)
            if (!baseSource->deref()) {
                delete baseSource;
            }
 
        for (vector<EGS_RadionuclideSpectrum * >::iterator it =
                    decays.begin();
                it!=decays.end(); it++) {
            delete *it;
            *it=0;
        }
        decays.clear();
    };
 
    EGS_I64 getNextParticle(EGS_RandomGenerator *rndm,
                            int &q, int &latch, EGS_Float &E, EGS_Float &wt,
                            EGS_Vector &x, EGS_Vector &u);
 
    EGS_Float getEmax() const {
        return Emax;
    };
 
    EGS_Float getFluence() const {
        return (ishower+1)*(baseSource->getFluence()/sCount); 
    };
 
    double getTime() const {
        return time;
    };
 
    double getLastTime() const {
        return lastDisintTime;
    };
 
    double getExperimentTime() const {
        return experimentTime;
    };
 
    EGS_I64 getShowerIndex() const {
        return ishower;
    };
 
    unsigned int getEmissionType() const {
        return emissionType;
    }
 
    void printSampledEmissions() {
        egsInformation("\n======================================================\n");
        egsInformation("Start of source emissions statistics:\n");
        for (unsigned int i=0; i<decays.size(); ++i) {
            decays[i]->printSampledEmissions();
        }
        egsInformation("End of source emissions statistics\n");
        egsInformation("======================================================\n\n");
    };
 
    bool isValid() const {
        return baseSource;
    };
 
    bool storeState(ostream &data_out) const;
 
    bool addState(istream &data);
 
    void resetCounter();
 
    bool setState(istream &data);
 
    EGS_RadionuclideSpectrum *createSpectrum(EGS_Input *input);
 
    vector<EGS_Ensdf *> getRadionuclideEnsdf() {
        vector<EGS_Ensdf *> decayEnsdf;
        for (auto dec: decays) {
            decayEnsdf.push_back(dec->getRadionuclideEnsdf());
        }
        return decayEnsdf;
    }
 
private:
    EGS_Application *app;
 
    EGS_I64             count; 
    EGS_Float           Emax; 
 
    void setUp();
 
    string sName; 
    EGS_I64 sCount; 
    EGS_BaseSource *baseSource; 
 
    vector<int>         q_allowed; 
    vector<EGS_RadionuclideSpectrum *> decays; 
    EGS_Float           activity; 
 
    bool                q_allowAll; 
    bool                disintegrationOccurred; 
    EGS_Float           time, 
                        experimentTime, 
                        lastDisintTime; 
    EGS_I64             ishower; 
    EGS_Vector          xOfDisintegration; 
 
    unsigned int emissionType;
};
 
#endif