Merge #216: implement a radionuclide source

HEN_HOUSE/doc/src/pirs898-egs++/Doxyfile

@@ -223,6 +223,7 @@ ALIASES                = "internwarning=\par Warning:\n This class is not part o
                          "FT=\author Frederic Tessier, NRC" \
                          "EM=\author Ernesto Mainegra-Hing, NRC" \
                          "BW=\author Blake Walters, NRC" \
+                         "RT=\author Reid Townson, NRC" \
                          "GG=\author Georgi Gerganov"
 
 # This tag can be used to specify a number of word-keyword mappings (TCL only).

HEN_HOUSE/egs++/Makefile

@@ -41,7 +41,8 @@ egspp_files = egs_input egs_base_geometry egs_library egs_transformations \
              egs_projectors egs_alias_table egs_object_factory egs_spectra \
              egs_base_source egs_functions egs_application egs_run_control \
              egs_scoring egs_interpolator egs_atomic_relaxations \
-             egs_ausgab_object egs_particle_track egs_fortran_geometry
+             egs_ausgab_object egs_particle_track egs_fortran_geometry \
+             egs_ensdf
 
 egspp_objects = $(addprefix $(DSO1), $(addsuffix .$(obje), $(egspp_files)))
 config1h = $(IEGS1)$(DSEP)egs_config1.h egs_libconfig.h egs_functions.h
@@ -54,7 +55,8 @@ geometry_libs = egs_planes egs_cd_geometry egs_gtransformed egs_nd_geometry \
 
 source_libs = egs_collimated_source egs_isotropic_source egs_parallel_beam \
               egs_point_source egs_source_collection egs_transformed_source \
-              egs_beam_source egs_phsp_source egs_angular_spread iaea_phsp_source
+              egs_beam_source egs_phsp_source egs_angular_spread \
+              iaea_phsp_source egs_radionuclide_source
 
 shape_libs = egs_circle egs_ellipse egs_extended_shape egs_gaussian_shape \
              egs_line_shape egs_polygon_shape egs_rectangle egs_shape_collection \

HEN_HOUSE/egs++/ausgab_objects/egs_dose_scoring/egs_dose_scoring.h

@@ -195,6 +195,34 @@ class EGS_DOSE_SCORING_EXPORT EGS_DoseScoring : public EGS_AusgabObject {
         return 0;
     };
 
+    int processEvent(EGS_Application::AusgabCall iarg, int ir) {
+
+        if (ir == -1) {
+            ir = app->top_p.ir;
+        }
+
+        int imed = ir>=0 ? app->getMedium(ir):-1;
+        EGS_Float edep = app->getEdep();
+
+        /**** energy deposition in a medium ***/
+        if (iarg <= 4 && imed >= 0 && edep > 0 && doseM) {
+            doseM->score(imed, edep*app->top_p.wt);
+        }
+
+        /*** Check if scoring in current region ***/
+        if (dose) {
+            if (d_reg_index[ir]<0) {
+                return 0;
+            }
+        }
+
+        /**** energy deposition in current region ***/
+        if (iarg <= 4 && ir >= 0 && edep > 0 && dose) {
+            dose->score(d_reg_index[ir], edep*app->top_p.wt);
+        }
+        return 0;
+    };
+
     bool needsCall(EGS_Application::AusgabCall iarg) const {
         if (iarg <= 4) {
             return true;

HEN_HOUSE/egs++/egs_advanced_application.cpp

@@ -1114,7 +1114,18 @@ EGS_Float EGS_AdvancedApplication::getMediumRho(int ind) {
 EGS_Float EGS_AdvancedApplication::getEdep() {
     return the_epcont->edep;
 }
-//************************************************************
+// Sets edep
+void EGS_AdvancedApplication::setEdep(EGS_Float edep) {
+    the_epcont->edep = edep;
+}
+// Get the ecut
+EGS_Float EGS_AdvancedApplication::getEcut() {
+    return the_bounds->ecut;
+}
+// Get the pcut
+EGS_Float EGS_AdvancedApplication::getPcut() {
+    return the_bounds->pcut;
+}
 // Returns rest mass
 EGS_Float EGS_AdvancedApplication::getRM() {
     return the_useful->rm;

HEN_HOUSE/egs++/egs_advanced_application.h

@@ -208,6 +208,9 @@ class APP_EXPORT EGS_AdvancedApplication : public EGS_Application {
     //************************************************************
     EGS_Float getMediumRho(int ind);
     EGS_Float getEdep();
+    void setEdep(EGS_Float edep);
+    EGS_Float getEcut();
+    EGS_Float getPcut();
     /* Needed by some sources */
     EGS_Float getRM();
 

HEN_HOUSE/egs++/egs_application.cpp

@@ -104,11 +104,17 @@ EGS_EXPORT void EGS_Application::setActiveApplication(EGS_Application *a) {
     active_egs_application = a;
 };
 
-int EGS_Application::userScoring(int iarg) {
+int EGS_Application::userScoring(int iarg, int ir) {
     if (a_objects) {
         int early_return = 0;
         for (int j=0; j<a_objects[iarg].size(); ++j) {
-            int res = a_objects[iarg][j]->processEvent((AusgabCall)iarg);
+            int res;
+            if (ir > -1) {
+                res = a_objects[iarg][j]->processEvent((AusgabCall)iarg, ir);
+            }
+            else {
+                res = a_objects[iarg][j]->processEvent((AusgabCall)iarg);
+            }
             if (res < 0) {
                 return res;
             }
@@ -120,7 +126,12 @@ int EGS_Application::userScoring(int iarg) {
             return early_return;
         }
     }
-    return ausgab(iarg);
+    if (ir > -1) {
+        return 0;
+    }
+    else {
+        return ausgab(iarg);
+    }
 }
 
 void EGS_Application::checkEnvironmentVar(int &argc, char **argv,

HEN_HOUSE/egs++/egs_application.h

@@ -676,7 +676,7 @@ class EGS_EXPORT EGS_Application {
       registered with the application and then proceeds to call ausgab().
 
     */
-    int userScoring(int iarg);
+    int userScoring(int iarg, int ir=-1);
 
     /*! \brief User scoring function.
 
@@ -1101,6 +1101,13 @@ class EGS_EXPORT EGS_Application {
     virtual EGS_Float getEdep() {
         return 0.0;
     };
+    virtual void setEdep(EGS_Float edep) {};
+    virtual EGS_Float getEcut() {
+        return 0.0;
+    };
+    virtual EGS_Float getPcut() {
+        return 0.0;
+    };
     virtual EGS_Float getRM() {
         return -1.0;
     };

HEN_HOUSE/egs++/egs_atomic_relaxations.cpp

@@ -23,7 +23,7 @@
 #
 #  Author:          Iwan Kawrakow, 2008
 #
-#  Contributors:
+#  Contributors:    Reid Townson
 #
 ###############################################################################
 */
@@ -231,7 +231,7 @@ class EGS_LOCAL EGS_RelaxImplementation {
 public:
 
     EGS_RelaxImplementation(const char *data_path) : elements(0), nz(0) {
-        data_file = egsJoinPath(data_path,"relax.data");
+        data_file = egsJoinPath(data_path,"relax_onebyte.data");
     };
 
     ~EGS_RelaxImplementation() {
@@ -278,6 +278,10 @@ class EGS_LOCAL EGS_RelaxImplementation {
         return elements[Z-1]->shells[shell].be;
     };
 
+    int getNShell(int Z) {
+        return elements[Z-1]->nshell;
+    };
+
     void setBindingEnergy(int Z, int shell, EGS_Float new_be) {
         checkData(Z,shell);
         elements[Z-1]->shells[shell].be = new_be;
@@ -526,6 +530,10 @@ EGS_Float EGS_AtomicRelaxations::getBindingEnergy(int Z, int shell) {
     return p->bindingEnergy(Z,shell);
 }
 
+int EGS_AtomicRelaxations::getNShell(int Z) {
+    return p->getNShell(Z);
+}
+
 void EGS_AtomicRelaxations::setBindingEnergy(int Z, int shell,
         EGS_Float new_be) {
     p->setBindingEnergy(Z,shell,new_be);

HEN_HOUSE/egs++/egs_atomic_relaxations.h

@@ -23,7 +23,7 @@
 #
 #  Author:          Iwan Kawrakow, 2008
 #
-#  Contributors:
+#  Contributors:    Reid Townson
 #
 ###############################################################################
 */
@@ -106,6 +106,9 @@ class EGS_EXPORT EGS_AtomicRelaxations {
     /*! Returns the binding energy of shell \a shell in element \a Z */
     EGS_Float getBindingEnergy(int Z, int shell);
 
+    /*! Returns the number of shells in element \a Z */
+    int      getNShell(int Z);
+
     /*! Set the binding energy of shell \a shell in element \a Z to
         \a new_be. This method can be used to overwrite the default
         data loaded with the relaxation probabilities when using

HEN_HOUSE/egs++/egs_ausgab_object.h

@@ -88,6 +88,7 @@ class EGS_EXPORT EGS_AusgabObject : public EGS_Object {
      *
      */
     virtual int processEvent(EGS_Application::AusgabCall iarg) = 0;
+    virtual int processEvent(EGS_Application::AusgabCall iarg, int ir) {};
 
     /*! \brief Is the ausgab call \a iarg relevant for this object?
      *

HEN_HOUSE/egs++/egs_base_source.h

@@ -169,6 +169,28 @@ class EGS_EXPORT EGS_BaseSource : public EGS_Object {
      */
     virtual EGS_Float getFluence() const = 0;
 
+    /*! \brief Get the time of emission for the most recently sampled particle
+     *
+     * This method is only reimplemented by EGS_RadionuclideSource. It
+     * returns the emission time of the particle that was most recently sampled.
+     */
+    virtual double getTime() const {};
+
+    /*! \brief Get the shower index for radionuclide emissions
+     *
+     * This method is only reimplemented by EGS_RadionuclideSource. It
+     * gets the index of the most recent shower.
+     */
+    virtual EGS_I64 getShowerIndex() const {};
+
+    /*! \brief Prints out the sampled emissions for radionuclide spectra
+     *
+     * This method is only reimplemented by EGS_RadionuclideSource. It
+     * prints the actual sampled intensity of each type of emission from
+     * the radionuclide spectra.
+     */
+    virtual void printSampledEmissions() {};
+
     /*!  \brief Store the source state into the stream \a data_out.
      *
      * Every source should reimplement this method to store
@@ -309,6 +331,7 @@ class EGS_EXPORT EGS_BaseSource : public EGS_Object {
  *  \link EGS_DoubleGaussianSpectrum a double-Gaussian spectrum\endlink,
  *  \link EGS_UniformSpectrum a uniform spectrum\endlink and
  *  \link EGS_TabulatedSpectrum a tabulated spectrum\endlink.
+ *  \link EGS_RadionuclideSpectrum a radionuclide spectrum\endlink.
  *
  */
 class EGS_EXPORT EGS_BaseSpectrum {
@@ -350,6 +373,60 @@ class EGS_EXPORT EGS_BaseSpectrum {
         return e;
     };
 
+    /*! \brief Get the charge for the most recently sampled particle
+     *
+     * This method is only reimplemented by EGS_RadionuclideSpectrum. It
+     * returns the charge of the particle that was most recently sampled
+     * using sampleEnergy().
+     */
+    virtual int getCharge() const {};
+
+    /*! \brief Get the time of emission for the most recently sampled particle
+     *
+     * This method is only reimplemented by EGS_RadionuclideSpectrum. It
+     * returns the emission time of the particle that was most recently sampled
+     * using sampleEnergy().
+     */
+    virtual double getTime() const {};
+
+    /*! \brief Get the shower index for radionuclide emissions
+     *
+     * This method is only reimplemented by EGS_RadionuclideSpectrum. It
+     * gets the index of the most recent shower produced using sampleEnergy().
+     */
+    virtual EGS_I64 getShowerIndex() const {};
+
+    /*! \brief Get the spectrum weight for radionuclide spectra
+     *
+     * This method is only reimplemented by EGS_RadionuclideSpectrum. It
+     * gets the weight of the spectrum to balance emissions from multiple
+     * spectra.
+     */
+    virtual EGS_Float getSpectrumWeight() const {};
+
+    /*! \brief Set the spectrum weight for radionuclide spectra
+     *
+     * This method is only reimplemented by EGS_RadionuclideSpectrum. It
+     * sets the weight of the spectrum to balance emissions from multiple
+     * spectra. This allows a source to normalize the spectrum weights.
+     */
+    virtual void setSpectrumWeight(EGS_Float newWeight) {};
+
+    /*! \brief Prints out the sampled emissions for radionuclide spectra
+     *
+     * This method is only reimplemented by EGS_RadionuclideSpectrum. It
+     * prints the actual sampled intensity of each type of emission from
+     * the radionuclide spectra.
+     */
+    virtual void printSampledEmissions() {};
+
+    /*! \brief Get energy that should be deposited locally from relaxations/alphas.
+     *
+     * This method is only reimplemented by EGS_RadionuclideSpectrum. It
+     * gets the energy deposited locally during spectrum generation.
+     */
+    virtual EGS_Float getEdep() const {};
+
     /*! \brief Get the maximum energy of this spectrum.
      *
      * This pure virtual method must be reimplemented by derived classes