tutor2pp.cpp

The tutor2pp application calculates the deposited, transmitted and reflected energy fractions for any source of electrons or photons that can be defined with the egspp source package, incident on any geometry that can be constructed using the egspp geometry package. The corresponding original tutorial mortran user code tutor2, which calculates the deposited, transmitted and reflected energy loss fractions for 20 MeV electrons perpendicularly incident on a semi-infinite 1 mm plate of Tantalum, only provides a very small fraction of the functionality offered by tutor2pp.

The complete source code of the application is found at the end of this documentation page.

To obtain the required functionality we implement a Tutor2_Application class derived from EGS_SimpleApplication

class APP_EXPORT Tutor2_Application : public EGS_SimpleApplication {
 
    EGS_ScoringArray *edep;   // our scoring object
    int              nreg;    // number of regions in the geometry.
 
public:

The meaning of the private data members edep and nreg will become clear in the following. The constructor of the Tutor2_Application class takes the command line arguments given to the application as arguments and passes them to the EGS_SimpleApplication constructor,

    Tutor2_Application(int argc, char **argv) :
        EGS_SimpleApplication(argc,argv) {

which reads the input file and initializes the source, the geometry, the cross section data and the random number generator, opens I/O units necessary for the mortran back-end, creates a temporary working directory for the run and checks that the cross section data covers the energy range defined by the source. This corresponds to steps 0, 1, 2, 3, 4 and 6 from the recipe for writing EGSnrc mortran user codes (see PIRS–701). We then obtain the number of regions from the geometry

        nreg = g->regions();

and initialize a EGS_ScoringArray object that will be used to collect the energy deposited in the various regions

        edep = new EGS_ScoringArray(nreg+2);
    };

This corresponds to step 5 from the recipe for writing EGSnrc mortran user codes. We initialize the scoring array to have 2 more regions than the geometry so that we can collect the energy reflected (defined as the energy of particles exiting the geometry and moving backwards and scored in region 0 of edep) and the energy transmitted (defined as the energy of particles exiting the geometry and moving forward and scored in region nreg+1 of edep) in addition to the energy deposited in the nreg regions of the geometry (scored in regions 1 ... nreg of edep). This is all that is needed in the Tutor2_Application constructor.

To demonstrate good coding habits we provide a destructor for the Tutor2_Application that deallocates the memory needed by the scoring array

    ~Tutor2_Application() {

The next step is the implementation of the actual scoring which is accomplished by re-implementing the ausgab function to collect energy deposition. Energy deposition is indicated by values of its argument less or equal to 4

    int ausgab(int iarg) {
        if (iarg <= 4) {

We obtain the index of the top particle on the stack and its region index from the the_stack, which is a pointer to the Fortran common block STACK

            int np = the_stack->np - 1;

Note that the mortran back-end uses Fortran style indexing and therefore we have to use the_stack->np-1 for the particle index. The convention for region numbers used by the geometry package is that the outside region has index -1 and inside regions are between 0 and nreg-1. However, this is translated to the convention typically adopted by EGSnrc mortran user codes that the outside region is region 1 and inside regions are 2 to nreg+1 by EGS_SimpleApplication. We have to therefore subtract 1 from the region index and change it to nreg+1 if it is 0 (particle is outside) and the particle is movinge forward to obtain the edep region for collecting the energy:

We then simply use the score() function of the scoring array object to score the energy

This is all about ausgab. We automatically obtain history-by-history statistical analysis, provided we inform our scoring object each time the simulation of a new history begins. This is accomplished by reimplementing the startHistory function, which is called from within the shower loop before transporting a new particle from the source:

The final step in the implementation of the Tutor2_Application class is to report the simulation results. For this purpose we define a separate function reportResults, which will be called from our main program. In the implementation of this function we simply use the reportResults() method of the scoring array. This method divides the accumulated results by the number of statistically independent events. Hence, in order to get fractions of the energy imparted into the phantom, we have to multiply by this number (available in the last_case data member) and divide by the total energy that entered the phantom (available in the Etot data member). Our reportResults tfunction herefore looks like this

and completes the implementation of the Tutor2_Application class.

The main program of the tutor2pp application is very simple, if fact we use the defined macros APP_MAIN or APP_LIB, from EGS_SimpleApplication. APP_MAIN creates the application object, runs a simulation, outputs results and finishes. APP_LIB just creates and returns the application object.

That's all. With 35 lines of code (excluding comments and blank lines) we have implemented a complete EGSnrc C++ application, which provides much more functionality than the original tutor2 mortran application (97 lines of code excluding comments and blank lines).

We now need a Makefile for the tutor2pp application. We include the general EGSnrc config file, the general egspp config file and the config file describing our C++ compiler:

include $(EGS_CONFIG)
include $(SPEC_DIR)egspp1.spec
include $(SPEC_DIR)egspp_$(my_machine).conf

We then set the name of our user code

USER_CODE = tutor2pp

We will be using the generic Makefile for C++ applications provided with the distribution (see below). This generic Makefile offers the possibility to add extra files with mortran macros and mortran code when building the EGSnrc mortran functions and subroutines by setting the EGSPP_USER_MACROS make variable. In our case we don't use this capability and therefore leave EGSPP_USER_MACROS empty:

EGSPP_USER_MACROS =

We also need to specify the base class from which our application class was derived:

EGS_BASE_APPLICATION = egs_simple_application

and add dependences for the tutor2pp.cpp file. In our case the only additional dependence is the egs_scoring.h file, which we include in order to get access to the definition of the EGS_ScoringArray class:

other_dep_user_code = $(ABS_EGSPP)egs_scoring.h

The make variable ABS_EGSPP is defined in egspp1.spec to be the absolute path to the egspp sources ($HEN_HOUSE/egs++/) Finally we simply include the generic Makefile for EGSnrc C++ applications, which provides rules for compiling the various sources and building the application:

include $(HEN_HOUSE)makefiles$(DSEP)cpp_makefile

We now can build our application by simply going to the tutor2pp directory and typing make. The resulting executable will be called tutor2pp (tutor2pp.exe on Windows) and will be put into the bin/$my_machine subdirectory of $EGS_HOME.

To run the application, we simply type

tutor2pp -i test1 -p tutor_data

with test1 being the name of the input file without extension (test1.egsinp is an example input file provided with the distribution) and tutor_data the name of the PEGS data file for the media of the geometry defined in test1. Alternatively, we can use the egs_gui graphical user interface. We can also set the number of histories to be run on the command line with the -n option:

tutor2pp -i test1 -p tutor_data -n 1000000

(all applications derived from EGS_SimpleApplication automatically understand the -n option).

Here is the complete source code of the tutor2pp application:

/*
###############################################################################
#
#  EGSnrc egs++ tutor2pp application
#  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:    Ernesto Mainegra-Hing
#                   Reid Townson
#
###############################################################################
#
#  A simple EGSnrc application using the C++ interface. It implements the
#  functionality of the original tutor2 tutorial code written in mortran
#  except that now, due to the use of the general geometry and source packages
#  included in egs++, any geometry or any source can be used, not just a
#  monoenergetic 20 MeV source of electrons incident on a 1 mm plate of
#  tantalum as un tutor2.
#
###############################################################################
*/
 
 
#include "egs_simple_application.h"
#include "egs_scoring.h"
#include "egs_interface2.h"
#include "egs_functions.h"
 
class APP_EXPORT Tutor2_Application : public EGS_SimpleApplication {
 
    EGS_ScoringArray *edep;   // our scoring object
    int              nreg;    // number of regions in the geometry.
 
public:
 
    Tutor2_Application(int argc, char **argv) :
        EGS_SimpleApplication(argc,argv) {
        nreg = g->regions();
        edep = new EGS_ScoringArray(nreg+2);
    };
 
    ~Tutor2_Application() {
        delete edep;
    };
 
    int ausgab(int iarg) {
        if (iarg <= 4) {
            int np = the_stack->np - 1;
            int ir = the_stack->ir[np]-1;
            if (ir == 0 && the_stack->w[np] > 0) {
                ir = nreg+1;
            }
            edep->score(ir,the_epcont->edep*the_stack->wt[np]);
        }
 
        return 0;
    };
 
    void startHistory(EGS_I64 icase) {
        edep->setHistory(icase);
    };
 
 
    void reportResults() {
        double norm = ((double)last_case)/Etot;
        egsInformation(" last case = %d Etot = %g\n",(int)last_case,Etot);
        edep->reportResults(norm,
                            "Reflected/deposited/transmitted energy fraction",false,
                            "  %d  %9.5f +/- %9.5f %c\n");
    };
 
};
 
#ifdef BUILD_APP_LIB
APP_LIB(Tutor2_Application);
#else
APP_MAIN(Tutor2_Application);
#endif