egs_application.h

Go to the documentation of this file.

/*
###############################################################################
#
#  EGSnrc egs++ application 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:    Frederic Tessier
#                   Ernesto Mainegra-Hing
#                   Blake Walters
#                   Reid Townson
#                   Alexandre Demelo
#                   Hannah Gallop
#
###############################################################################
*/
 
 
#ifndef EGS_APPLICATION_
#define EGS_APPLICATION_
 
#include "egs_base_geometry.h"
#include "egs_base_source.h"
#include "egs_simple_container.h"
#include "egs_interpolator.h"
#include "egs_input_struct.h"
#include "egs_run_control.h"
#include "egs_scoring.h"
#include "egs_rndm.h"
 
#include <memory>
#include <string>
#include <iostream>
#include <dirent.h>
using namespace std;
 
class EGS_Input;
class EGS_BaseSource;
class EGS_RandomGenerator;
class EGS_RunControl;
class EGS_GeometryHistory;
class EGS_AusgabObject;
class EGS_Interpolator;
//template <class T> class EGS_SimpleContainer;
 
inline void addmcBlock(shared_ptr<EGS_InputStruct> blockPtr) {
    shared_ptr<EGS_BlockInput> mcBlock = blockPtr->addBlockInput("MC transport parameter");
    mcBlock->addSingleInput("Global ECUT", false, "Global electron transport cutoff");
    mcBlock->addSingleInput("Global PCUT", false, "Global photon transport cutoff");
    mcBlock->addSingleInput("Global SMAX", false, "Global maximum step-size restriction for e-transport");
    mcBlock->addSingleInput("ESTEPE", false, "The maximum fractional energy loss in an electron step due to continuous energy loss. Defaults to 0.25");
    mcBlock->addSingleInput("XIMAX", false, "Maximum first GS moment per step. Defaults to 0.5, maximum value is 1.");
    mcBlock->addSingleInput("Boundary crossing algorithm", false, "Defaults to exact", {"exact", "PRESTA-I"});
    mcBlock->addSingleInput("Skin depth for BCA", false, "The distance from a boundary to turn on electron single scattering, in mean free paths of the current medium. Set to a very large value to do a single scattering simulation. Defaults to 3 for exact boundary crossing");
    mcBlock->addSingleInput("Electron-step algorithm", false, "Defaults to PRESTA-II", {"PRESTA_II", "PRESTA_I"});
    mcBlock->addSingleInput("Spin effects", false, "Includes relativistic spin effects in electron scattering (On is highly recommended). Defaults to On", {"On", "Off"});
    mcBlock->addSingleInput("Brems angular sampling", false, "Defaults to KM", {"KM", "Simple"});
    mcBlock->addSingleInput("Brems cross sections", false, "Defaults to BH", {"BH", "NIST"});
    mcBlock->addSingleInput("Pair angular sampling", false, "Defaults to Simple", {"Simple", "Off", "KM"});
    mcBlock->addSingleInput("Pair cross sections", false, "Defaults to BH. Use NRC at low energies to consider asymmetry in positron-electron energy distribution.", {"BH", "NRC"});
    mcBlock->addSingleInput("Triplet production", false, "Defaults to On", {"On", "Off"});
    mcBlock->addSingleInput("Electron Impact Ionization", false, "For more accurate kV x-ray simulations, turn this on. Defaults to Off", {"On", "Off", "casnati", "kolbenstvedt", "gryzinski"});
    mcBlock->addSingleInput("Bound Compton scattering", false, "Defaults to norej", {"On", "Off", "Simple", "norej"});
    mcBlock->addSingleInput("Radiative Compton corrections", false, "Defaults to Off", {"On", "Off"});
    mcBlock->addSingleInput("Rayleigh scattering", false, "Defaults to Off", {"On", "Off", "custom"});
    mcBlock->addSingleInput("Photoelectron angular sampling", false, "Defaults to On", {"On", "Off"});
    mcBlock->addSingleInput("Atomic relaxations", false, "Defaults to on", {"On", "Off"});
    mcBlock->addSingleInput("Photon cross sections", false, "Defaults to xcom, can also be user-supplied", {"si", "epdl", "xcom"});
    mcBlock->addSingleInput("Photon cross-sections output", false, "Defaults to Off", {"On", "Off"});
    mcBlock->addSingleInput("Compton cross sections", false, "User-supplied, Defaults to comp-xsections");
    mcBlock->addSingleInput("Photonuclear attenuation", false, "Defaults to Off", {"On", "Off"});
    mcBlock->addSingleInput("Photonuclear cross sections", false, "Defaults to default, or is user-supplied", {"default"});
}
 
inline void addMediaDefBlock(shared_ptr<EGS_InputStruct> blockPtr) {
    shared_ptr<EGS_BlockInput> mediaBlockInput = blockPtr->addBlockInput("media definition");
    mediaBlockInput->addSingleInput("ae", false, "Lowest energy for electron production (kinetic+0.511) in MeV. Minimum and most accurate value is 0.521.");
    mediaBlockInput->addSingleInput("ap", false, "Lowest energy for photon production (kinetic) in MeV. Minimum and most accurate value is 0.001.");
    mediaBlockInput->addSingleInput("ue", false, "Maximum energy for electrons (kinetic+0.511) in MeV. Set to just above the maximum source energy.");
    mediaBlockInput->addSingleInput("up", false, "Maximum energy for photons (kinetic) in MeV. Set to just above the maximum source energy.");
 
    mediaBlockInput->addSingleInput("material data file", false, "The full file path to a file containing material data. See HEN_HOUSE/pegs4/data/material.dat for an example containing all the materials from 521ICRU.pegs4dat and 700ICRU.pegs4dat.");
 
    shared_ptr<EGS_BlockInput> mediumBlock = mediaBlockInput->addBlockInput("myMediumName");
    mediumBlock->addSingleInput("elements", false, "A list of elements, separated by commas");
    auto atomsPtr = mediumBlock->addSingleInput("number of atoms", false, "For each of the elements, the number of atoms in a molecule of the medium, separated by commas. If the medium is a mixture, use mass fractions instead.");
    auto fractionsPtr = mediumBlock->addSingleInput("mass fractions", false, "For each of the elements, the mass fraction of that element in the mixture, comma separated.");
    mediumBlock->addSingleInput("rho", false, "The bulk density of the medium in g/cm^3");
    mediumBlock->addSingleInput("sterncid", false, "Material ID for the STERNHEIMER-SELTZER-BERGER (SSB) lookup table");
    mediumBlock->addSingleInput("stopping powers", false, "This parameter allows various stopping powers to be calculated, rather than just the restricted total stopping powers. This allows for simulations that model various types of CSDA calculations. Defaults to restricted total.", {"restricted total", "unrestricted collision", "unrestricted collision and radiative", "unrestricted collision and restricted radiative", "restricted collision and unrestricted radiative", "unrestricted radiative"});
    mediumBlock->addSingleInput("bremsstrahlung correction", false, "Corrections to apply to calculated bremsstrahlung cross-sections. Defaults to KM", {"KM", "NRC", "None"});
    mediumBlock->addSingleInput("gas pressure", false, "The pressure in atm of this medium, assuming it is a gas. This is only applied if a density correction file is not used. To apply a gas pressure in that case, generate a new density correction file with the appropriately scaled density.");
    mediumBlock->addSingleInput("e- stopping power output file", false, "An output filename to write electron stopping powers to.");
 
    atomsPtr->addDependency(fractionsPtr, "", true);
    fractionsPtr->addDependency(atomsPtr, "", true);
}
 
inline string addmcExample() {
    string example = {
        R"(
:start MC transport parameter:
    Global ECUT                    = 0.521          # Global electron transport cutoff
    Global PCUT                    = 0.001          # Global photon transport cutoff
    Global SMAX                    = 1e10           # Global maximum step-size restriction for e- transport
    ESTEPE                         = 0.25           # Defaults to 0.25
    XIMAX                          = 0.5            # Defaults to 0.5, max. value is 1.
    Boundary crossing algorithm    = exact          # exact (default), PRESTA-I
    Skin depth for BCA             = 3              # Default value is 3 for exact boundary crossing
    Electron-step algorithm        = PRESTA-II      # PRESTA-II (default),PRESTA-I
    Spin effects                   = On             # On (default), Off
    Brems angular sampling         = KM             # Simple, KM (default)
    Brems cross sections           = BH             # BH (default), NIST, NRC
    Pair angular sampling          = Simple         # Off, Simple (default), KM
    Pair cross sections            = BH             # NRC, BH (default)
    Triplet production             = Off            # On, Off (default)
    Electron Impact Ionization     = Off            # On, Off (default), casnati, kolbenstvedt, gryzinski
    Bound Compton scattering       = norej          # On, Off, Simple, norej (default)
    Radiative Compton corrections  = Off            # On,  Off (default)
    Rayleigh scattering            = Off            # On, Off (default), custom
    Photoelectron angular sampling = On             # On (default), Off
    Atomic relaxations             = On             # On (default), Off
    Photon cross sections          = xcom           # si, epdl, xcom (default) or user-supplied
    Photon cross-sections output   = Off            # Off (default), On
    Compton cross sections         = comp_xsections # User-supplied
    Photonuclear attenuation       = Off            # Off (default), On
    Photonuclear cross sections    = default        # default (default), user-supplied
:stop MC transport parameter:
)"};
    return example;
}
 
inline string addMediaExample() {
    string example = {
    // Not completed yet, need to fill in numbers for example
        R"(
:start media definition:
    ae  = 0.521                 # lowest  energy for electron production (kinetic+0.511)
    ap  = 0.01                  # lowest  energy for photon production   (kinetic)
    ue  = 50.511                # maximum energy for electrons (kinetic+0.511)
    up  = 50                    # maximum energy for photons (kinetic) # maximum energy for photons (kinetic)
 
    # Here is an example for defining water
    :start water:
        density correction file = water_liquid
    :stop water:
 
:stop media definition:
)"};
    return example;
}
 
struct EGS_Particle {
    int        q;      
    int        latch;  
    int        ir;     
    EGS_Float  E;      
    EGS_Float  wt;     
    EGS_Vector x;      
    EGS_Vector u;      
};
 
class EGS_EXPORT EGS_Application {
 
public:
 
    EGS_Application(int argc, char **argv);
 
    virtual ~EGS_Application();
 
    virtual int initSimulation();
 
    virtual void setSimulationChunk(EGS_I64 nstart, EGS_I64 nrun, int npar, int nchunk);
 
    virtual int runSimulation();
 
    virtual int finishSimulation();
 
    virtual void describeSimulation();
 
    virtual int simulateSingleShower();
 
    virtual void getCurrentResult(double &sum, double &sum2, double &norm,
                                  double &count) {
        sum = 0;
        sum2 = 0;
        norm = 1;
        count = 0;
    };
 
    virtual void analyzeResults() {};
 
    virtual void outputResults() {};
 
    virtual int combineResults();
 
    virtual int combinePartialResults();
 
    int howManyJobsDone();
 
 
    virtual int  outputData();
 
    virtual int  readData();
 
    EGS_Input *getInput() {
        return input;
    };
 
    const string &getAppName() const {
        return app_name;
    };
 
    const string &getEgsHome() const {
        return egs_home;
    };
 
    const string &getHenHouse() const {
        return hen_house;
    };
 
    const string &getOutputFile() const {
        return output_file;
    };
 
    const string &getFinalOutputFile() const {
        return final_output_file;
    };
 
    string constructIOFileName(const char *extension, bool with_run_dir) const;
 
    const string &getAppDir() const {
        return app_dir;
    };
 
    const string &getRunDir() const {
        return run_dir;
    };
 
    const string &getWorkDir() const {
        return run_dir;
    };
 
    enum AusgabCall {
        BeforeTransport = 0,     
        EgsCut = 1,              
        PegsCut = 2,             
        UserDiscard = 3,         
        ExtraEnergy = 4,         
        AfterTransport = 5,      
        BeforeBrems = 6,         
        AfterBrems = 7,          
        BeforeMoller = 8,        
        AfterMoller = 9,         
        BeforeBhabha = 10,       
        AfterBhabha = 11,        
        BeforeAnnihFlight = 12,  
        AfterAnnihFlight = 13,   
        BeforeAnnihRest = 28,    
        AfterAnnihRest = 14,     
        BeforePair = 15,         
        AfterPair = 16,          
        BeforeCompton = 17,      
        AfterCompton = 18,       
        BeforePhoto = 19,        
        AfterPhoto = 20,         
        EnteringUphi = 21,       
        LeavingUphi = 22,        
        BeforeRayleigh = 23,     
        AfterRayleigh = 24,      
        FluorescentEvent = 25,   
        CosterKronigEvent = 26,  
        AugerEvent = 27,         
        BeforePhotoNuc = 29,     
        AfterPhotoNuc = 30,      
        BeforeEII = 31,          
        AfterEII = 32,           
        AfterSubPhoton = 33,     
        AfterSubElectron = 34,   
        UnknownCall = 35         
    };
 
    virtual void setAusgabCall(AusgabCall call, bool on_or_off) {
        ausgab_flag[call] = on_or_off;
    };
 
    int getNparallel() const {
        return n_parallel;
    };
 
    int getIparallel() const {
        return i_parallel;
    };
 
    int getFirstParallel() const {
        return first_parallel;
    };
 
    inline int howfar(int ireg, const EGS_Vector &x, const EGS_Vector &u,
                      EGS_Float &t, int *newmed) {
 
        geometry->resetErrorFlag();
        EGS_Float twant = t;
        int inew = geometry->howfar(ireg,x,u,t,newmed);
        storeGeometryStep(ireg,inew,x,u,twant,t);
        if (geometry->getLastError()) {
            reportGeometryError();
        }
        return inew;
 
        //return geometry->howfar(ireg,x,u,t,newmed);
    };
 
    inline EGS_Float hownear(int ireg,const EGS_Vector &x) {
        return geometry->hownear(ireg,x);
    };
 
    inline int getMedium(int ireg) {
        return geometry->medium(ireg);
    };
 
    bool isRealRegion(int ireg) {
        return geometry->isRealRegion(ireg);
    }
    int  isWhere(EGS_Vector &r) {
        return geometry->isWhere(r);
    }
 
    int getGlobalRegionOffset(const string geomName) {
        return geometry->getGlobalRegionOffset(geomName);
    }
 
    void getNumberRegions(const string &str, vector<int> &regs) {
        geometry->getNumberRegions(str, regs);
    }
 
    void getLabelRegions(const string &str, vector<int> &regs) {
        geometry->getLabelRegions(str, regs);
    }
 
    EGS_Float getTimeIndex() {
        return source->getTimeIndex();
    }
 
    void setTimeIndex(EGS_Float temp_time) {
        source->setTimeIndex(temp_time);
    }
 
    int userScoring(int iarg, int ir=-1);
 
    virtual int ausgab(int) {
        return 0;
    };
 
    virtual void startNewParticle() { };
 
    virtual void enterNewRegion() { };
 
    virtual void fillRandomArray(int n, EGS_Float *rns);
 
    static EGS_Application *activeApplication();
 
    static void setActiveApplication(EGS_Application *);
 
    virtual EGS_I64 randomNumbersUsed() const;
 
    virtual void getElectronSteps(double &ch_steps, double &all_steps) const {
        ch_steps = 0;
        all_steps = 0;
    };
 
    virtual int addState(istream &data);
 
    virtual void resetCounter();
 
    virtual void describeUserCode() const {};
 
    virtual void appInformation(const char *);
 
    virtual void appWarning(const char *);
 
    virtual void appFatal(const char *);
 
    void checkDeviceFull(FILE *);
 
    static bool getArgument(int &argc, char **argv,
                            const char *name1, const char *name2, string &arg);
 
    static void checkEnvironmentVar(int &argc, char **argv, const char *env,
                                    const char *n1, const char *n2, string &var);
 
protected:
 
    virtual int initGeometry();
 
    virtual int initSource();
 
    virtual int initCrossSections() {
        return 0;
    };
 
    virtual int initScoring() {
        return 0;
    };
 
    virtual int initRunControl();
 
    virtual int initRNG();
 
    virtual int initEGSnrcBackEnd() {
        return 0;
    };
 
    void initAusgabObjects();
 
    void addAusgabObject(EGS_AusgabObject *o);
 
    virtual int startNewShower();
 
    virtual int finishShower() {
        return 0;
    };
 
    virtual int shower() {
        return 0;
    };
 
    virtual void finishRun() { };
 
    void storeGeometryStep(int ireg, int inew, const EGS_Vector &x,
                           const EGS_Vector &u, EGS_Float twant, EGS_Float t);
 
    void reportGeometryError();
 
    EGS_Input           *input;         
    EGS_BaseGeometry    *geometry;      
    EGS_BaseSource      *source;        
    EGS_RandomGenerator *rndm;          
    EGS_RunControl      *run;           
 
    bool    ausgab_flag[UnknownCall]; 
 
    string  app_name;           
    string  egs_home;           
    string  hen_house;          
    string  app_dir;            
    string  run_dir;            
    string  egs_config;         
    string  input_file;         
    string  output_file;        
    string  final_output_file;  
    string  pegs_file;          
    string  abs_pegs_file;      
 
    int     n_parallel,  
            i_parallel,  
            first_parallel; 
    bool    batch_run;   
    bool    simple_run;  
    bool    uniform_run; 
    bool    is_pegsless; 
 
    EGS_Particle p;      
    EGS_I64      current_case; 
    EGS_I64      last_case;    
 
    ostream *data_out;
 
    istream *data_in;
 
    int app_index;
 
    EGS_SimpleContainer<EGS_AusgabObject *> a_objects_list;
 
    EGS_SimpleContainer<EGS_AusgabObject *> *a_objects;
 
    EGS_GeometryHistory *ghistory;
 
private:
 
    static int n_apps; 
 
public:
 
    EGS_Particle top_p;  
    int          Np;     
    //************************************************************
    // Utility functions for use with ausgab dose scoring objects
    //************************************************************
    EGS_Float getFluence() {
        return source->getFluence();
    };
    int       getnRegions() {
        return geometry->regions();
    };
    int       getnMedia() {
        return geometry->nMedia();
    };
    const char *getMediumName(int ind) {
        return geometry->getMediumName(ind);
    };
    virtual EGS_Float getMediumRho(int ind) {
        (void)ind;
        return -1.0;
    };
    virtual EGS_Float getEdep() {
        return 0.0;
    };
    virtual void setEdep(EGS_Float edep) {
        (void)edep;
    };
    virtual EGS_Float getEcut() {
        return 0.0;
    };
    virtual EGS_Float getPcut() {
        return 0.0;
    };
    virtual EGS_Float getRM() {
        return -1.0;
    };
 
    //************************************************
    // For use with ausgab radiative splitting objects
    //************************************************
    virtual void setRadiativeSplitting(const EGS_Float &nsplit) {
        (void)nsplit;
    };
    virtual void setRussianRoulette(const EGS_Float &iSwitchRR) {
        (void)iSwitchRR;
    };
    virtual void splitTopParticleIsotropically(const EGS_Float &fsplit) {
        (void)fsplit;
    }
 
    //************************************************************
    // Utility functions for use with ausgab fluence scoring objects
    //************************************************************
    virtual EGS_Float getTVSTEP() {
        return 0.0;
    };
 
    virtual EGS_Interpolator *getDEDX(const int &imed, const int &iq) {
        (void)imed;
        (void)iq;
        return 0;
    };
 
    string sourceType() {
        return source->getObjectType();
    }
 
    int sourceCharge() {
        return source->getCharge();
    }
 
    int sourceEmax() {
        return source->getEmax();
    }
 
    virtual void setLatch(const int &ip, const int &latch) {
        (void)ip;
        (void)latch;
    };
 
    virtual void incLatch(const int &ip, const int &increment) {
        (void)ip;
        (void)increment;
    };
 
    virtual int getNp() {
        return 0;
    };
 
    virtual int getNpOld() {
        return 0;
    };
 
    //************************************************************
    // Utility function for ausgab phase space scoring objects
    //************************************************************
    virtual void setLatch(int latch) {
        (void)latch;
    };
 
    static unique_ptr<EGS_InputStruct> inputStructure;
 
    bool containsDynamic() {
        bool hasDynamic = false;
        geometry->containsDynamic(hasDynamic);
        if (!hasDynamic) {
            source->containsDynamic(hasDynamic);
        }
        return hasDynamic;
    }
};
 
#define APP_MAIN(app_name) \
    int main(int argc, char **argv) { \
        app_name app(argc,argv); \
        int err = app.initSimulation(); \
        if( err ) return err; \
        err = app.runSimulation(); \
        if( err < 0 ) return err; \
        return app.finishSimulation(); \
    }
 
#define APP_SIMPLE_MAIN(app_name) \
    int main(int argc, char **argv) { \
        app_name app(argc,argv); \
        app.run(); \
        app.reportResults(); \
        app.finish(); \
        return 0; \
    }
 
#define APP_LIB(app_name) \
    extern "C" {\
        APP_EXPORT EGS_Application* createApplication(int argc, char **argv) {\
            return new app_name(argc, argv);\
        }\
        APP_EXPORT shared_ptr<EGS_InputStruct> getAppInputs() {\
            shared_ptr<EGS_InputStruct> appInputStruct = make_shared<EGS_InputStruct>();\
            addmcBlock(appInputStruct);\
            addRngDefinitionBlock(appInputStruct);\
            addRunControlBlock(appInputStruct);\
            addMediaDefBlock(appInputStruct);\
            return appInputStruct;\
        }\
        APP_EXPORT string getmcExample() {\
            return addmcExample();\
        }\
        APP_EXPORT string getRunControlExample() {\
            return addRunControlExample();\
        }\
        APP_EXPORT string getRngDefinitionExample() {\
            return addRngDefinitionExample();\
        }\
        APP_EXPORT string getMediaExample() {\
            return addMediaExample();\
        }\
    }\
 
#endif