Setting medium dose = yes
to get some statistics as part of the output
:start ausgab object:
library = egs_dose_scoring
name = BONE
medium dose = yes
region dose = no
:start output dose file:
geometry name = phantom4 # must name an EGS_XYZGeometry
file type = 3ddose
:stop output dose file:
:stop ausgab object:
With original phantom but only a single region for simplicity
:start geometry:
name = phantom4
library = egs_ndgeometry
type = EGS_XYZGeometry
x-slabs = -30 60 1
y-slabs = -30 60 1
z-slabs = 0 1000 1
:start media input:
media = ICRPBONE521ICRU #H2O521ICRU
:stop media input:
:stop geometry:
By running for both media individually (and manually combining the results for convenience),
==> Summary of media dosimetry (per fluence)
medium Edep/[MeV*cm2] D/[Gy*cm2]
-------------------------------------------------------------------
H2O521ICRU 1.8255e+02 +/- 0.104 % 2.9245e-08 +/- 0.104 %
ICRPBONE521ICRU 1.9001e+02 +/- 0.102 % 1.6454e-08 +/- 0.102 %
-------------------------------------------------------------------
With a big phantom,
:start geometry:
name = phantom4
library = egs_ndgeometry
type = EGS_XYZGeometry
x-slabs = -200 400 1
y-slabs = -200 400 1
z-slabs = 0 1000 1
:start media input:
media = ICRPBONE521ICRU #H2O521ICRU
:stop media input:
:stop geometry:
==> Summary of media dosimetry (per fluence)
medium Edep/[MeV*cm2] D/[Gy*cm2]
-------------------------------------------------------------------
H2O521ICRU 1.9119e+02 +/- 0.102 % 3.0629e-08 +/- 0.102 %
ICRPBONE521ICRU 1.9153e+02 +/- 0.102 % 1.6586e-08 +/- 0.102 %
-------------------------------------------------------------------
So the phantom is only leaking a bit of energy, and a quick look at the dose shows that it behaves as expected (D_bone = D_water / 1.85). The value written to the 3ddose file is wrong and assumes rho = 1.0. That can be easily confirmed by hooking into the (dose scoring routine)[https://github.com/nrc-cnrc/EGSnrc/blob/master/HEN_HOUSE/egs%2B%2B/ausgab_objects/egs_dose_scoring/egs_dose_scoring.cpp#L446] with
std::cout << "Edep: " << r << ", mass: " << mass << std::endl;
and indeed, mass
is 1.0 even in the ICRPBONE
case. Seems like a good old fashioned bug to me, but I'm not that familiar with the egs++ side of things.