egs_ensdf.h

Go to the documentation of this file.

/*
###############################################################################
#
#  EGSnrc egs++ ensdf 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:
#
###############################################################################
*/


#ifndef EGS_ENSDF_
#define EGS_ENSDF_

#include "egs_libconfig.h"
#include "egs_functions.h"
#include "egs_math.h"
#include "egs_alias_table.h"
#include "egs_atomic_relaxations.h"

#include <iostream>
#include <fstream>
#include <algorithm>
#include <vector>
#include <map>

using namespace std;

template <class T> class Branch {
public:

    Branch() {}

    ~Branch() {
        for (typename vector<T *>::iterator it  = branchLeaves.begin();
                it!=branchLeaves.end(); it++) {
            (*it)->removeBranch();
        }
        branchLeaves.clear();
    }

    void addLeaf(T *leaf) {
        branchLeaves.push_back(leaf);
    }

    void removeLeaf(T *leaf) {
        branchLeaves.erase(std::remove(branchLeaves.begin(),
                                       branchLeaves.end(),
                                       leaf), branchLeaves.end());
    }

    vector<T *> getLeaves() const {
        return branchLeaves;
    }

    // A new == operator for this class
    bool operator==(const Branch<T> &rhs) const {
        for (typename vector<T *>::const_iterator it = branchLeaves.begin();
                it!=branchLeaves.end(); it++) {

            bool foundLeaf = false;
            for (typename vector<T *>::const_iterator irhs =
                        rhs.branchLeaves.begin();
                    irhs!=rhs.branchLeaves.end(); irhs++) {

                if (*irhs != 0 && *it != 0) {
                    if (*irhs == *it) {
                        foundLeaf = true;
                    }
                }
            }

            if (!foundLeaf) {
                return false;
            }
        }
        return true;
    }

protected:
    vector<T *> branchLeaves;
};

template <class T> class Leaf {
public:

    Leaf(T *existingBranch) {
        branch = existingBranch;
        if (branch) {
            branch->addLeaf(this);
        }
    }

    ~Leaf() {
        if (branch) {
            branch->removeLeaf(this);
        }
        branch = 0;
    }

    virtual T *getBranch() const {
        return branch;
    }

    void removeBranch() {
        branch = 0;
    }

    // A new == operator for this class
    bool operator== (const T &rhs) const {
        if (branch==0 && rhs.branch==0) {
            return true;
        }
        else if ((branch==0) && rhs.branch!=0) {
            return false;
        }
        else if ((branch!=0) && rhs.branch==0) {
            return false;
        }
        else if (branch!=0 && rhs.branch!=0) {
            return *branch == *(rhs.branch);
        }
    }

private:
    T *branch;
};

// The Record class
class Record {
public:
    Record();
    Record(vector<string> ensdf);
    virtual ~Record();
    vector<string> getRecords() const;

protected:
    double recordToDouble(int startPos, int endPos);
    string recordToString(int startPos, int endPos);
    double getTag(string searchString, string notAfter);
    double parseHalfLife(int startPos, int endPos);
    double parseStdUncertainty(string value, string stdUncertainty);
    string getStringAfter(string searchString, size_t len);

    // All the lines corresponding to this record type
    vector<string> lines;
};

// Comment Record
class CommentRecord : public Record {
public:
    CommentRecord(vector<string> ensdf);
    string getComment();

private:
    string comment;
    void processEnsdf();
};

// Parent Record
class ParentRecord : public Record, public Branch<Leaf<ParentRecord> > {
public:
    ParentRecord(vector<string> ensdf);
    double getHalfLife() const;
    double getQ() const;

protected:
    double  halfLife,
            Q;

private:
    void processEnsdf();
};

class ParentRecordLeaf : public Leaf<ParentRecord> {
public:
    ParentRecordLeaf(ParentRecord *myRecord);
    virtual ParentRecord *getParentRecord() const;
};

// Normalization Record
class NormalizationRecord : public Record, public
    Branch<Leaf<NormalizationRecord> >, public ParentRecordLeaf {
public:
    NormalizationRecord(vector<string> ensdf, ParentRecord *parent);
    double getRelativeMultiplier() const;
    double getTransitionMultiplier() const;
    double getBranchMultiplier() const;
    double getBetaMultiplier() const;
    EGS_AtomicRelaxations *getRelaxations() const;
    double getBindingEnergy(int shell) const;
    int getNShell() const;
    void relax(int shell,
               EGS_Float ecut, EGS_Float pcut,
               EGS_RandomGenerator *rndm, double &edep,
               EGS_SimpleContainer<EGS_RelaxationParticle> &particles);

protected:
    double normalizeRelative;
    double normalizeTransition;
    double normalizeBeta;
    double normalizeBranch;

private:
    void processEnsdf();
    EGS_AtomicRelaxations *relaxations;
    int nshell, Z;
};

class NormalizationRecordLeaf : public Leaf<NormalizationRecord> {
public:
    NormalizationRecordLeaf(NormalizationRecord *myRecord);
    virtual NormalizationRecord *getNormalizationRecord() const;
};

// Level Record
class EGS_EXPORT LevelRecord : public Record, public Branch<Leaf<LevelRecord> > {
public:
    LevelRecord();
    LevelRecord(vector<string> ensdf);
    void resetDisintegrationIntensity();
    void cumulDisintegrationIntensity(double disintIntensity);
    double getDisintegrationIntensity() const;
    void setLevelCanDecay(bool canDecay);
    bool levelCanDecay() const;
    double getEnergy() const;
    double getHalfLife() const;

protected:
    double disintegrationIntensity;
    double energy;
    double halfLife;
    bool canDecay;

private:
    void processEnsdf();
};

class LevelRecordLeaf : public Leaf<LevelRecord> {
public:
    LevelRecordLeaf(LevelRecord *myRecord);
    virtual LevelRecord *getLevelRecord() const;
};

// Generic beta record
class EGS_EXPORT BetaRecordLeaf : public Record, public ParentRecordLeaf, public
    NormalizationRecordLeaf, public LevelRecordLeaf {
public:
    BetaRecordLeaf(vector<string> ensdf, ParentRecord *myParent,
                   NormalizationRecord *myNormalization, LevelRecord *myLevel);

    virtual double getFinalEnergy() const = 0;
    virtual double getBetaIntensity() const = 0;
    virtual double getPositronIntensity() const {
        return 0;
    };
    virtual double getECIntensity() const {
        return 0;
    };
    virtual void relax(int shell,
                       EGS_Float ecut, EGS_Float pcut,
                       EGS_RandomGenerator *rndm, double &edep,
                       EGS_SimpleContainer<EGS_RelaxationParticle> &particles) {};
    virtual void setBetaIntensity(double newIntensity)  = 0;
    int getCharge() const;
    void incrNumSampled();
    EGS_I64 getNumSampled() const;
    unsigned short int getZ() const;
    unsigned short int getAtomicWeight() const;
    unsigned short int getForbidden() const;
    void setSpectrum(EGS_AliasTable *bspec);
    EGS_AliasTable *getSpectrum() const;
    vector<double> ecShellIntensity;

protected:
    EGS_I64 numSampled;
    double finalEnergy;
    double betaIntensity;
    int q;
    unsigned short int Z;
    unsigned short int A;
    unsigned short int forbidden;
    EGS_AliasTable *spectrum;
};

// Beta- record
class EGS_EXPORT BetaMinusRecord : public BetaRecordLeaf {
public:
    BetaMinusRecord(vector<string> ensdf, ParentRecord *myParent,
                    NormalizationRecord *myNormalization, LevelRecord *myLevel);

    double getFinalEnergy() const;
    double getBetaIntensity() const;
    double getBetaIntensityUnc() const;
    void setBetaIntensity(double newIntensity);

private:
    void processEnsdf();
    double betaIntensityUnc;
};

// Beta+ Record (and Electron Capture)
class EGS_EXPORT BetaPlusRecord : public BetaRecordLeaf {
public:
    BetaPlusRecord(vector<string> ensdf, ParentRecord *myParent,
                   NormalizationRecord *myNormalization, LevelRecord *myLevel);

    double getFinalEnergy() const;
    double getBetaIntensity() const;
    double getPositronIntensity() const;
    double getPositronIntensityUnc() const;
    double getECIntensityUnc() const;
    void setBetaIntensity(double newIntensity);
    void setPositronIntensity(double newIntensity);
    void relax(int shell,
               EGS_Float ecut, EGS_Float pcut,
               EGS_RandomGenerator *rndm, double &edep,
               EGS_SimpleContainer<EGS_RelaxationParticle> &particles);

protected:
    double  ecIntensity,
            positronIntensity,
            ecIntensityUnc,
            positronIntensityUnc;

private:
    void processEnsdf();
};

// Gamma record
class EGS_EXPORT GammaRecord : public Record, public ParentRecordLeaf,
    public NormalizationRecordLeaf, public LevelRecordLeaf {
public:
    GammaRecord(vector<string> ensdf, ParentRecord *myParent,
                NormalizationRecord *myNormalization,
                LevelRecord *myLevel);
    GammaRecord(GammaRecord *gamma);

    double getDecayEnergy() const;
    double getTransitionIntensity() const;
    double getGammaIntensity() const;
    double getGammaIntensityUnc() const;
    double getICIntensity() const;
    double getICIntensityUnc() const;
    double getIPIntensity() const;
    double getIPIntensityUnc() const;
    void setTransitionIntensity(double newIntensity);
    void setGammaIntensity(double newIntensity);
    void setICIntensity(double newIntensity);
    double getMultiTransitionProb() const;
    void setMultiTransitionProb(double newIntensity);
    int getCharge() const;
    LevelRecord *getFinalLevel() const;
    void setFinalLevel(LevelRecord *newLevel);
    void incrGammaSampled();
    void incrICSampled();
    void incrIPSampled();
    EGS_I64 getGammaSampled() const;
    EGS_I64 getICSampled() const;
    EGS_I64 getIPSampled() const;
    vector<double> icIntensity;
    double getBindingEnergy(int shell) const;
    void relax(int shell,
               EGS_Float ecut, EGS_Float pcut,
               EGS_RandomGenerator *rndm, double &edep,
               EGS_SimpleContainer<EGS_RelaxationParticle> &particles);

protected:
    EGS_I64 numGammaSampled, numICSampled, numIPSampled;
    double  decayEnergy;
    double  transitionIntensity,
            multipleTransitionProb,
            gammaIntensity,
            gammaIntensityUnc,
            icCoeff,
            icCoeffUnc,
            ipCoeff,
            ipCoeffUnc;
    int q;
    LevelRecord *finalLevel;

private:
    void processEnsdf();
};

// Alpha record
class EGS_EXPORT AlphaRecord : public Record, public ParentRecordLeaf, public
    NormalizationRecordLeaf, public LevelRecordLeaf {
public:
    AlphaRecord(vector<string> ensdf, ParentRecord *myParent,
                NormalizationRecord *myNormalization, LevelRecord *myLevel);

    double getFinalEnergy() const;
    double getAlphaIntensity() const;
    double getAlphaIntensityUnc() const;
    int getCharge() const;
    void setAlphaIntensity(double newIntensity);
    void incrNumSampled();
    EGS_I64 getNumSampled() const;

protected:
    EGS_I64 numSampled;
    double  finalEnergy,
            alphaIntensity,
            alphaIntensityUnc;
    int q;

private:
    void processEnsdf();
};

class EGS_EXPORT EGS_Ensdf {

public:

    EGS_Ensdf(const string nuclide, const string ensdf_filename="",
              const string relaxType="eadl", const bool allowMultiTrans=false, int verbosity=1);

    ~EGS_Ensdf();

    vector<Record * > getRecords() const;
    vector<BetaRecordLeaf *> getBetaRecords() const;
    vector<ParentRecord * > getParentRecords() const;
    vector<LevelRecord * > getLevelRecords() const;
    vector<AlphaRecord * > getAlphaRecords() const;
    vector<GammaRecord * > getGammaRecords() const;
    vector<GammaRecord * > getMetastableGammaRecords() const;
    vector<GammaRecord * > getUncorrelatedGammaRecords() const;
    vector<double > getXRayIntensities() const;
    vector<double > getXRayEnergies() const;
    vector<double > getAugerIntensities() const;
    vector<double > getAugerEnergies() const;

    string radionuclide;
    int verbose;
    string relaxationType;
    unsigned short int Z;
    double decayDiscrepancy;
    bool allowMultiTransition;

    void normalizeIntensities();

protected:

    unsigned short int findAtomicWeight(string element);
    void parseEnsdf(vector<string> ensdf);
    void buildRecords();

    void getEmissionsFromComments();

    ifstream ensdf_file;
    unsigned short int A;

    vector<Record * > myRecords;
    vector<CommentRecord * > myCommentRecords;
    vector<ParentRecord * > myParentRecords;
    vector<NormalizationRecord * > myNormalizationRecords;
    vector<LevelRecord * > myLevelRecords;
    vector<BetaRecordLeaf *> myBetaRecords;
    vector<BetaMinusRecord * > myBetaMinusRecords;
    vector<BetaPlusRecord * > myBetaPlusRecords;
    vector<AlphaRecord * > myAlphaRecords;
    vector<GammaRecord * > myGammaRecords;
    vector<GammaRecord * > myMetastableGammaRecords;
    vector<GammaRecord * > myUncorrelatedGammaRecords;

private:

    vector<vector<string> > recordStack;
    vector<string> commentLines;
    vector<double>  xrayEnergies,
           xrayIntensities,
           augerEnergies,
           augerIntensities;
    ParentRecord *previousParent;
};




#endif