Changing egs++ geometry getMass to getVolume

The egs_base_geometry getMass virtual function is exchanged with
a getVolume function that returns region volume instead of mass.
This change was performed because typically, initGeometry() is
invoked before any media data is read in egs_applications, thus
getMass' dependence on medium data could be unreliable.  The
getMass function was thus replaced with a getVolume function,
which only relies on internal gemoetry definitions to return a
value.
The geometries in egs_autoenvelope, egs_nd_geometry, egs_rz,
and egs_spheres are affected.  The ausgab dose scoring object
which outputs a 3ddose file was also changed to reflect the
change from getMass to getVolume.  This modification of
egs_dose_scoring also now includes relativeRho in the mass
calculation.

HEN_HOUSE/egs++/ausgab_objects/egs_dose_scoring/egs_dose_scoring.cpp

@@ -132,7 +132,7 @@ void EGS_DoseScoring::setApplication(EGS_Application *App) {
     // Get the number of media in the input file
     nmedia = app->getnMedia();
     // determine maximum medium name length
-    char buf[32];
+    char buf[64];
     int count = 0;
     int imed=0;
     max_medl = 6;  // length of the string "medium" is the minimum
@@ -261,7 +261,6 @@ void EGS_DoseScoring::setApplication(EGS_Application *App) {
     description += buf;
     description += "--------------------------------------\n";
     for (imed=0; imed < nmedia; imed++) {
-        sprintf(buf,"%-*s",max_medl,app->getMediumName(imed));
         description += buf;
         description += "  ";
         sprintf(buf,"%11.8f",app->getMediumRho(imed));
@@ -451,7 +450,7 @@ void EGS_DoseScoring::outputDoseFile(const EGS_Float &normD) {
         //divide dose by mass and output
         for (int i=0; i<nx*ny*nz; i++) {
             doseF->currentResult(i,r,dr);
-            EGS_Float mass = dose_geom->getMass(i); //local reg.
+            EGS_Float mass = dose_geom->getVolume(i)*getRealRho(i); //local reg.
             dose=r*normD/mass;
             df_out << dose << " ";
         }

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

@@ -164,7 +164,7 @@ If one of the geometries in the simulation is an EGS_XYZGeometry, then the user
 has the option to output the dose in the voxels of this geometry to a file.  Currently,
 the only available output format is the .3ddose format, which is familiar to users
 of DOSXYZnrc (see DOSXYZnrc users manual for more details). In this case,
-masses of each voxel are available through the member function getMass of
+masses of each voxel are available through the member function getVolume of
 EGS_XYZGeometry.  An example input to obtain a .3ddose file for an EGS_XYZGeometry is given below:
 
 \verbatim
@@ -298,6 +298,11 @@ class EGS_DOSE_SCORING_EXPORT EGS_DoseScoring : public EGS_AusgabObject {
         }
         return (int)log10((float)imax);
     };
+
+	EGS_Float getRealRho(int ireg) {
+		int med = dose_geom->medium(ireg);
+		return dose_geom->getRelativeRho(ireg)*app->getMediumRho(med);
+	}
 
     void setVol(const vector<EGS_Float> volin) {
         vol_list=volin;

HEN_HOUSE/egs++/egs_base_geometry.h

@@ -231,11 +231,12 @@ class EGS_EXPORT EGS_BaseGeometry {
      */
     virtual EGS_Float hownear(int ireg, const EGS_Vector &x) = 0;
 
-    /*! \brief Calculates the volume*relative rho (rhor) of region ireg.
+    /*! \brief Calculates the volume of region ireg.
 
-      Currently only implemented in EGS_XYZGeometry
+      Currently only implemented EGS_XYZGeometry, EGS_cSpheres,
+	  EGS_cSphericalShell, EGS_AEnvelope, and EGS_RZGeometry
     */
-    virtual EGS_Float getMass(int ireg) {
+    virtual EGS_Float getVolume(int ireg) {
         return 1.0;
     }
 

HEN_HOUSE/egs++/geometry/egs_autoenvelope/egs_autoenvelope.cpp

@@ -59,7 +59,7 @@
 
 string EGS_AENVELOPE_LOCAL EGS_AEnvelope::type = "EGS_AEnvelope";
 string EGS_AENVELOPE_LOCAL EGS_ASwitchedEnvelope::type = "EGS_ASwitchedEnvelope";
-/* only geometries that support getting regional volume/mass can be used as phantoms */
+/* only geometries that support getting regional volume can be used as phantoms */
 const string EGS_AENVELOPE_LOCAL EGS_AEnvelope::allowed_base_geom_types[] = {"EGS_cSpheres", "EGS_cSphericalShell", "EGS_XYZGeometry", "EGS_RZ"};
 
 static char EGS_AENVELOPE_LOCAL geom_class_msg[] = "createGeometry(AEnvelope): %s\n";
@@ -82,7 +82,7 @@ EGS_AEnvelope::EGS_AEnvelope(EGS_BaseGeometry *base_geom,
 
         string msg(
             "EGS_AEnvelope:: Volume correction is not available for geometry type '%s (%s)'. "
-            "Geometry types must implement getMass.  Valid choices are:\n\t"
+            "Geometry types must implement getVolume.  Valid choices are:\n\t"
         );
 
         int end = (int)(sizeof(allowed_base_geom_types)/sizeof(string));
@@ -99,12 +99,6 @@ EGS_AEnvelope::EGS_AEnvelope(EGS_BaseGeometry *base_geom,
     is_convex = base_geom->isConvex();
     has_rho_scaling = getHasRhoScaling();
 
-    // initialize uncorrected/corrected masses
-    for (int ir = 0; ir < nregbase; ir++) {
-        uncorrected_mass.push_back(base_geom->getMass(ir));
-        corrected_mass.push_back(base_geom->getMass(ir));
-    }
-
     ninscribed = inscribed.size();
     if (ninscribed == 0) {
         egsFatal("EGS_AEnvelope: no inscribed geometries!\n");
@@ -156,12 +150,6 @@ EGS_AEnvelope::EGS_AEnvelope(EGS_BaseGeometry *base_geom,
         copy(it->second.begin(), it->second.end(), back_inserter(geoms_in_region[it->first]));
     }
 
-    // set mass correction based on volume correction
-    for (int ir = 0; ir < nregbase; ir++) {
-        double volume_ratio = vc_results.corrected_volumes[ir]/vc_results.uncorrected_volumes[ir];
-        corrected_mass[ir] = uncorrected_mass[ir]*volume_ratio;
-    }
-
     if (getNRegWithInscribed() == 0) {
         egsFatal("EGS_AEnvelope:: Failed to find any regions with inscribed geometries\n");
     }
@@ -618,24 +606,24 @@ int EGS_AEnvelope::getMaxStep() const {
     return nstep + inscribed_geoms.size();
 };
 
-EGS_Float EGS_AEnvelope::getMass(int ireg) {
+EGS_Float EGS_AEnvelope::getVolume(int ireg) {
 
     if (ireg < 0) {
         return -1;
     }
     else if (ireg < nregbase) {
-        return corrected_mass[ireg];
+        return vc_results.corrected_volumes[ireg];
     }
 
     volcor::GeomRegPairT local = global_reg_to_local[ireg];
 
-    return local.first->getMass(local.second);
+    return local.first->getVolume(local.second);
 };
 
-EGS_Float EGS_AEnvelope::getMassCorrectionRatio(int ireg) {
+EGS_Float EGS_AEnvelope::getCorrectionRatio(int ireg) {
 
     if (0 <= ireg && ireg < nregbase) {
-        return corrected_mass[ireg]/uncorrected_mass[ireg];
+        return vc_results.corrected_volumes[ireg]/vc_results.uncorrected_volumes[ireg];
     }
 
     return 1;
@@ -653,15 +641,15 @@ void EGS_AEnvelope::printInfo() const {
     if (debug_info) {
 
         for (int ir=0; ir < nregbase; ir++) {
-            if (fabs(uncorrected_mass[ir] - corrected_mass[ir]) > 1E-8) {
-                egsInformation("    mass of region %d was corrected from %.5E g to %.5E g\n",
-                               ir, uncorrected_mass[ir], corrected_mass[ir]);
+            if (fabs(vc_results.uncorrected_volumes[ir] - vc_results.corrected_volumes[ir]) > 1E-8) {
+                egsInformation("    volume of region %d was corrected from %.5E g to %.5E g\n",
+                               ir, vc_results.uncorrected_volumes[ir], vc_results.corrected_volumes[ir]);
             }
         }
 
         for (int ir=0; ir < nregbase; ir++) {
             if (geoms_in_region[ir].size() > 0) {
-                egsInformation("    region %d has %d insribed geometries\n", ir, geoms_in_region[ir].size());
+                egsInformation("    region %d has %d inscribed geometries\n", ir, geoms_in_region[ir].size());
             }
         }
     }

HEN_HOUSE/egs++/geometry/egs_autoenvelope/egs_autoenvelope.h

@@ -402,9 +402,9 @@ class EGS_AENVELOPE_EXPORT EGS_AEnvelope : public EGS_BaseGeometry {
 
     int getMaxStep() const;
 
-    virtual EGS_Float getMass(int ireg);
+    virtual EGS_Float getVolume(int ireg);
 
-    virtual EGS_Float getMassCorrectionRatio(int ireg);
+    virtual EGS_Float getCorrectionRatio(int ireg);
 
     virtual const string &getType() const {
         return type;
@@ -462,10 +462,6 @@ class EGS_AENVELOPE_EXPORT EGS_AEnvelope : public EGS_BaseGeometry {
     //keep track of which geometries are present in which base geometry regions
     vector<EGS_BaseGeometry *> *geoms_in_region;
 
-    // base geometry mass correction by region
-    vector<EGS_Float> uncorrected_mass;
-    vector<EGS_Float> corrected_mass;
-
     static string type;    //!< Geometry type
 
     const static string allowed_base_geom_types[];

HEN_HOUSE/egs++/geometry/egs_autoenvelope/volcor.h

@@ -542,7 +542,7 @@ vector<EGS_Float> applyFileVolumeCorrections(VCOptions *opts, vector<RegVolume>
 vector<EGS_Float> getUncorrectedVolumes(EGS_BaseGeometry *base) {
     vector<EGS_Float> uncorrected;
     for (int ir=0; ir < base->regions(); ir++) {
-        EGS_Float vol = base->getMass(ir)/base->getRelativeRho(ir);
+        EGS_Float vol = base->getVolume(ir);
         uncorrected.push_back(vol);
     }
     return uncorrected;

HEN_HOUSE/egs++/geometry/egs_nd_geometry/egs_nd_geometry.cpp

@@ -42,7 +42,6 @@
 #include "egs_input.h"
 #include "egs_functions.h"
 #include "egs_transformations.h"
-#include "egs_application.h"
 
 #ifdef HAS_GZSTREAM
     #include "gzstream.h"
@@ -172,8 +171,6 @@ EGS_XYZGeometry::EGS_XYZGeometry(EGS_PlanesX *Xp, EGS_PlanesY *Yp,
     xpos = xp->getPositions();
     ypos = yp->getPositions();
     zpos = zp->getPositions();
-
-    setApplication(EGS_Application::activeApplication());
 }
 
 EGS_XYZGeometry::~EGS_XYZGeometry() {

HEN_HOUSE/egs++/geometry/egs_nd_geometry/egs_nd_geometry.h

@@ -588,7 +588,6 @@ class EGS_NDG_EXPORT EGS_NDGeometry : public EGS_BaseGeometry {
 
 #ifdef EXPLICIT_XYZ
 #include "../egs_planes/egs_planes.h"
-#include "egs_application.h"
 
 /*! \brief An XYZ-geometry
 
@@ -1119,51 +1118,12 @@ class EGS_NDG_EXPORT EGS_XYZGeometry : public EGS_BaseGeometry {
     };
 
 
-    EGS_Float getMass(int ireg) {
-        if (ireg >= 0) {
-            int iz = ireg/nxy;
-            int ir = ireg - iz*nxy;
-            int iy = ir/nx;
-            int ix = ir - iy*nx;
-            EGS_Float vol = (xpos[ix+1]-xpos[ix])*(ypos[iy+1]-ypos[iy])*
-                            (zpos[iz+1]-zpos[iz]);
-
-            EGS_Float dens;
-            if (has_rho_scaling) {
-                dens = getRelativeRho(ireg);
-            }
-            else {
-                // We should only get here for the XYZ geometry
-
-                // Calculate the voxel midpoint
-                EGS_Vector tp;
-                EGS_Float minx,maxx,miny,maxy,minz,maxz;
-                minx=getBound(0,ix);
-                maxx=getBound(0,ix+1);
-                miny=getBound(1,iy);
-                maxy=getBound(1,iy+1);
-                minz=getBound(2,iz);
-                maxz=getBound(2,iz+1);
-                tp.x=(minx+maxx)/2.;
-                tp.y=(miny+maxy)/2.;
-                tp.z=(minz+maxz)/2.;
-
-                // Get the global region number from the current voxel midpoint
-                int g_reg = app->isWhere(tp);
-
-                // Get the medium index for this region
-                int imed = app->getMedium(g_reg);
-
-                // Get the density
-                dens = app->getMediumRho(imed);
-            }
-            EGS_Float mass = vol*dens;
-
-            return mass;
-        }
-        else {
-            return 1.0;
-        }
+    EGS_Float getVolume(int ireg) {
+		int iz = ireg/nxy;
+		int ir = ireg - iz*nxy;
+		int iy = ir/nx;
+		int ix = ir - iy*nx;
+		return (xpos[ix+1]-xpos[ix])*(ypos[iy+1]-ypos[iy])*(zpos[iz+1]-zpos[iz]);
     }
 
     EGS_Float getBound(int idir, int ind) {
@@ -1274,10 +1234,6 @@ class EGS_NDG_EXPORT EGS_XYZGeometry : public EGS_BaseGeometry {
         zp->getLabelRegions(str, regs);
     }
 
-    void setApplication(EGS_Application *App) {
-        app = App;
-    };
-
 protected:
 
     EGS_PlanesX      *xp;

HEN_HOUSE/egs++/geometry/egs_rz/egs_rz.cpp

@@ -59,8 +59,7 @@ EGS_RZGeometry::EGS_RZGeometry(vector<EGS_BaseGeometry *> geoms,
         radii.insert(radii.begin(), 0.);
     }
 
-    vector<EGS_Float> mass;
-    int ir = 0;
+    vector<EGS_Float> vol;
     for (size_t r=0; r < radii.size()-1; r++) {
 
         EGS_Float rmin = radii[r];
@@ -71,10 +70,7 @@ EGS_RZGeometry::EGS_RZGeometry(vector<EGS_BaseGeometry *> geoms,
         for (size_t plane = 0; plane < zbounds.size()-1; plane++) {
             EGS_Float zmin = zbounds[plane];
             EGS_Float zmax = zbounds[plane+1];
-            EGS_Float rho = getRelativeRho(ir);
-            EGS_Float vol = (zmax-zmin)*area;
-            reg_mass.push_back(rho*vol);
-            ir++;
+            reg_vol.push_back((zmax-zmin)*area);
         }
     }
 
@@ -101,11 +97,11 @@ int EGS_RZGeometry::getNRegDir(int dir) {
     return 0;
 }
 
-EGS_Float EGS_RZGeometry::getMass(int ireg) {
+EGS_Float EGS_RZGeometry::getVolume(int ireg) {
     if (ireg < 0 || ireg >= nreg) {
         return -1;
     }
-    return reg_mass[ireg];
+    return reg_vol[ireg];
 }
 
 

HEN_HOUSE/egs++/geometry/egs_rz/egs_rz.h

@@ -118,7 +118,7 @@ class EGS_RZ_EXPORT EGS_RZGeometry: public EGS_NDGeometry {
 
     vector<EGS_Float> radii; /*! cylinder radii. Note radii[0] is always set to 0 */
     vector<EGS_Float> zbounds;
-    vector<EGS_Float> reg_mass; /* calculated mass of each region */
+    vector<EGS_Float> reg_vol; /* calculated vol of each region */
     static string RZType;
 
 public:
@@ -153,7 +153,7 @@ class EGS_RZ_EXPORT EGS_RZGeometry: public EGS_NDGeometry {
     int getNRegDir(int dir);
 
     /*! \brief get mass of a given region  */
-    EGS_Float getMass(int ireg);
+    EGS_Float getVolume(int ireg);
 
 };
 

HEN_HOUSE/egs++/geometry/egs_spheres/egs_spheres.cpp

@@ -70,8 +70,7 @@ EGS_cSpheres::EGS_cSpheres(int ns, const EGS_Float *radius,
         rbounds.push_back(radius[ireg]);
         EGS_Float router = rbounds[ireg];
         EGS_Float rinner = ireg > 0 ? rbounds[ireg-1] : 0;
-        EGS_Float volume = (4./3.)*M_PI*(router*router*router - rinner*rinner*rinner);
-        mass.push_back(getRelativeRho(ireg) * volume);
+        vol.push_back((4./3.)*M_PI*(router*router*router - rinner*rinner*rinner));
     }
 }
 
@@ -347,9 +346,9 @@ int EGS_cSpheres::getNRegDir(int dir) {
 }
 
 
-EGS_Float EGS_cSpheres::getMass(int ireg) {
+EGS_Float EGS_cSpheres::getVolume(int ireg) {
     if (ireg >= 0 && ireg < nreg) {
-        return mass[ireg];
+        return vol[ireg];
     }
     return 1;
 }
@@ -382,8 +381,7 @@ EGS_cSphericalShell::EGS_cSphericalShell(int ns, const EGS_Float *radius,
     for (int ireg=0; ireg < nreg; ireg++) {
         EGS_Float rinner = R[ireg];
         EGS_Float router = R[ireg + 1];
-        EGS_Float volume = (4./3.)*M_PI*(router*router*router - rinner*rinner*rinner);
-        mass.push_back(getRelativeRho(ireg) * volume);
+        vol.push_back((4./3.)*M_PI*(router*router*router - rinner*rinner*rinner));
     }
 }
 
@@ -660,9 +658,9 @@ int EGS_cSphericalShell::getNRegDir(int dir) {
 }
 
 
-EGS_Float EGS_cSphericalShell::getMass(int ireg) {
+EGS_Float EGS_cSphericalShell::getVolume(int ireg) {
     if (ireg >= 0 && ireg < nreg) {
-        return mass[ireg];
+        return vol[ireg];
     }
     return 1;
 }