Implement range rejection AO

HEN_HOUSE/egs++/Makefile

@@ -67,7 +67,7 @@ shape_libs = egs_circle egs_ellipse egs_extended_shape egs_gaussian_shape \
              egs_circle_perpendicular
 
 aobject_libs = egs_track_scoring egs_dose_scoring egs_radiative_splitting \
-	       egs_phsp_scoring egs_fluence_scoring
+	       egs_phsp_scoring egs_fluence_scoring egs_range_rejection
 
 all_libs = $(geometry_libs) $(source_libs) $(shape_libs) $(aobject_libs)
 lib_objects = $(addprefix $(DSO1), $(addsuffix .$(obje), $(all_libs)))

HEN_HOUSE/egs++/ausgab_objects/egs_range_rejection/Makefile

@@ -0,0 +1,48 @@
+
+###############################################################################
+#
+#  EGSnrc egs++ makefile to build range rejection ausgab object
+#  Copyright (C) 2023 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:          Ernesto Mainegra-Hing, 2023
+#
+#  Contributors:
+#
+###############################################################################
+
+
+include $(EGS_CONFIG)
+include $(SPEC_DIR)egspp.spec
+include $(SPEC_DIR)egspp_$(my_machine).conf
+
+DEFS = $(DEF1) -DBUILD_RANGE_REJECTION_DLL
+
+library = egs_range_rejection
+lib_files = egs_range_rejection
+my_deps = $(common_ausgab_deps)
+extra_dep = $(addprefix $(DSOLIBS), $(my_deps))
+
+include $(SPEC_DIR)egspp_libs.spec
+
+$(make_depend)
+
+test:
+	@echo "common_h2: $(common_h2)"
+	@echo "extra_dep: $(extra_dep)"

HEN_HOUSE/egs++/ausgab_objects/egs_range_rejection/egs_range_rejection.cpp

@@ -0,0 +1,120 @@
+/*
+###############################################################################
+#
+#  EGSnrc egs++ range rejection asugab object
+#  Copyright (C) 2018 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:          Ernesto Mainegra-Hing, 2018
+#
+#  Contributors:
+#
+###############################################################################
+#
+#  A general range rejection tool.
+#
+#  TODO:
+#
+#  - Add Russian roulette (RR)
+#
+###############################################################################
+*/
+
+
+/*! \file egs_range_rejection.cpp
+ *  \brief A range rejection ausgab object: implementation
+ *  \EM
+ */
+
+#include <fstream>
+#include <string>
+#include <cstdlib>
+
+#include "egs_range_rejection.h"
+#include "egs_input.h"
+#include "egs_functions.h"
+
+EGS_RangeRejection::EGS_RangeRejection(const string &Name,
+        EGS_ObjectFactory *f) :
+    algorithm(rr), E_max(0) {
+    otype = "EGS_RangeRejection";
+}
+
+EGS_RangeRejection::EGS_RangeRejection(const string &Name,
+        EGS_ObjectFactory *f, const EGS_Float &E) :
+    algorithm(rr), E_max(E) {
+    otype = "EGS_RangeRejection";
+}
+
+EGS_RangeRejection::~EGS_RangeRejection() {
+}
+
+void EGS_RangeRejection::setApplication(EGS_Application *App) {
+    EGS_AusgabObject::setApplication(App);
+    if (!app) {
+        return;
+    }
+
+    char buf[32];
+
+    EGS_Float ecut = app->getEcut();
+
+    description = "\n===========================================\n";
+    description +=  "range rejection Object (";
+    description += name;
+    description += ")\n";
+    description += "===========================================\n";
+    if (E_max >= ecut) {
+        description +="\n - Range rejection below ";
+        sprintf(buf,"%g MeV\n\n",E_max);
+        description += buf;
+        // Set EGSnrc internal range rejection number .
+        app->setRangeRejection(E_max);
+    }
+    else {
+        description +="\n - NO range rejection";
+    }
+    description += "\n===========================================\n\n";
+}
+
+//*********************************************************************
+// Process input for this ausgab object
+//
+//**********************************************************************
+extern "C" {
+
+    EGS_RANGE_REJECTION_EXPORT EGS_AusgabObject *createAusgabObject(EGS_Input *input,
+            EGS_ObjectFactory *f) {
+        const static char *func = "createAusgabObject(range_rejection)";
+        if (!input) {
+            egsWarning("%s: null input?\n",func);
+            return 0;
+        }
+
+        EGS_Float E_rr = 0.0;
+        int err = input->getInput("maximum energy",E_rr);
+
+        //=================================================
+
+        /* Setup range rejection object with input parameters */
+        EGS_RangeRejection *result = new EGS_RangeRejection("",f,E_rr);
+        result->setName(input);
+        return result;
+    }
+}

HEN_HOUSE/egs++/ausgab_objects/egs_range_rejection/egs_range_rejection.h

@@ -0,0 +1,125 @@
+/*
+###############################################################################
+#
+#  EGSnrc egs++ range rejection ausgab object headers
+#  Copyright (C) 2023 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:          Ernesto Mainegra-Hing, 2023
+#
+#  Contributors:
+#
+###############################################################################
+#
+#  A general range rejection tool.
+#
+#  TODO:
+#
+#  - Add Russian roulette (RR)
+#
+###############################################################################
+*/
+
+
+/*! \file egs_range_rejection.h
+ *  \brief A range rejection ausgab object: header
+ *  \EM
+ */
+
+#ifndef EGS_RANGE_REJECTION_
+#define EGS_RANGE_REJECTION_
+
+#include "egs_ausgab_object.h"
+#include "egs_application.h"
+#include "egs_scoring.h"
+#include "egs_base_geometry.h"
+
+#ifdef WIN32
+
+    #ifdef BUILD_RANGE_REJECTION_DLL
+        #define EGS_RANGE_REJECTION_EXPORT __declspec(dllexport)
+    #else
+        #define EGS_RANGE_REJECTION_EXPORT __declspec(dllimport)
+    #endif
+    #define EGS_RANGE_REJECTION_LOCAL
+
+#else
+
+    #ifdef HAVE_VISIBILITY
+        #define EGS_RANGE_REJECTION_EXPORT __attribute__ ((visibility ("default")))
+        #define EGS_RANGE_REJECTION_LOCAL  __attribute__ ((visibility ("hidden")))
+    #else
+        #define EGS_RANGE_REJECTION_EXPORT
+        #define EGS_RANGE_REJECTION_LOCAL
+    #endif
+
+#endif
+
+/*! \brief A range rejection object: header
+
+\ingroup AusgabObjects
+
+This ausgab object can be used to increase the efficiency of radiative
+events. This ausgab object is specified via:
+\verbatim
+:start ausgab object:
+    library        = egs_range_rejection
+    name           = some_name
+    maximum energy = E # Energy below which to reject charged particles [MeV]
+:stop ausgab object:
+\endverbatim
+
+TODO:
+ - Add Russian roulette (RR)
+
+*/
+
+class EGS_RANGE_REJECTION_EXPORT EGS_RangeRejection : public EGS_AusgabObject {
+
+public:
+
+    /*! Rejection algorithm type */
+    enum AlgorithmType {
+        rr,  // EGSnrc intrinsic range rejection (approximation)
+        rrRR // Range rejection with Russin roulette (true VRT)
+    };
+
+    EGS_RangeRejection(const string &Name="", EGS_ObjectFactory *f = 0);
+    EGS_RangeRejection(const string &Name="", EGS_ObjectFactory *f = 0, const EGS_Float &E = 0);
+
+    ~EGS_RangeRejection();
+
+    void setApplication(EGS_Application *App);
+
+    int processEvent(EGS_Application::AusgabCall iarg) {
+        return 0;
+    };
+
+    int processEvent(EGS_Application::AusgabCall iarg, int ir) {
+        return 0;
+    };
+
+protected:
+    /*! Maximum energy below which to reject charged particles. Rest mass included. */
+    EGS_Float E_max;
+    AlgorithmType algorithm;
+
+};
+
+#endif

HEN_HOUSE/egs++/egs_advanced_application.cpp

@@ -1182,9 +1182,9 @@ void EGS_AdvancedApplication::resetRNGState() {
     }
 }
 
-//************************************************************
-// Utility functions for use with ausgab dose scoring object
-//************************************************************
+//****************************************************
+// Utility functions for use with dose scoring object
+//****************************************************
 // Returns density for medium ind
 EGS_Float EGS_AdvancedApplication::getMediumRho(int ind) {
     // handle the negative medium index for vacuum
@@ -1213,11 +1213,33 @@ EGS_Float EGS_AdvancedApplication::getPcut() {
 EGS_Float EGS_AdvancedApplication::getRM() {
     return the_useful->rm;
 }
-// Turns ON/OFF radiative splitting
+
+//*********************************************************
+// Utility functions for radiative splitting ausgab object
+//*********************************************************
+
+//! Turns ON/OFF radiative splitting
 void EGS_AdvancedApplication::setRadiativeSplitting(const EGS_Float &nsplit) {
     the_egsvr->nbr_split = nsplit;
 }
 
+//*****************************************************
+// Utility functions for range rejection ausgab object
+//*****************************************************
+void setRangeRejection(const EGS_Float &E) {
+
+    the_egsvr->e_max_rr = E;
+
+    if ( E > 0.0 ){
+       the_egsvr->i_do_rr  = 1;
+    }
+    else{
+       the_egsvr->i_do_rr  = 0;
+    }
+
+};
+
+
 //************************************************************
 // Utility functions for fluence scoring objects
 //************************************************************

HEN_HOUSE/egs++/egs_advanced_application.h

@@ -203,9 +203,9 @@ class APP_EXPORT EGS_AdvancedApplication : public EGS_Application {
     */
     void setEIIData(EGS_I32 len);
 
-    //************************************************************
-    // Utility functions for use with ausgab dose scoring objects
-    //************************************************************
+    //*****************************************************
+    // Utility functions for use with dose scoring objects
+    //*****************************************************
     EGS_Float getMediumRho(int ind);
     EGS_Float getEdep();
     void setEdep(EGS_Float edep);
@@ -229,9 +229,18 @@ class APP_EXPORT EGS_AdvancedApplication : public EGS_Application {
 
     /* Needed by some sources */
     EGS_Float getRM();
-    /* Turn ON/OFF radiative splitting */
+
+    //*********************************************************
+    // Utility functions for radiative splitting ausgab object
+    //*********************************************************
+    /*! Turn ON/OFF radiative splitting */
     void setRadiativeSplitting(const EGS_Float &nsplit);
 
+    //*****************************************************
+    // Utility functions for range rejection ausgab object
+    //*****************************************************
+    void setRangeRejection(const EGS_Float &E) {};
+
 protected:
 
     int              nmed;      //!< number of media