nrc-cnrc / EGSnrc

Toolkit for Monte Carlo simulation of ionizing radiation — Trousse d'outils logiciels pour la simulation Monte Carlo du rayonnement ionisant
http://nrc-cnrc.github.io/EGSnrc
GNU Affero General Public License v3.0
242 stars 146 forks source link

egs_chamber dose output for cylinders inscribed in water phantom #1220

Open erjft opened 3 days ago

erjft commented 3 days ago

Describe the bug Hi I have a 10 cm thick cylinder of 0.25 cm radius that is cut up into 1 mm sections along the z axis, and inscribed in a square water phantom of 15x15x15 cm^3. I want to score the dose in all the cylindrical volumes from an electron/photon beam.

The ausgab output, when I choose the same regions, is different from the scoring options output. For scoring options this would mean repeating :start calculation geometry: definition 100 times which may be painful.

:start geometry definition:
# WATER PHANTOM #
    :start geometry:
        name        = surrounding_water_box
        library     = egs_box
        box size    = 15 15 15
        :start media input:
            media = H2O521ICRU
        :stop media input:
        :start transformation:
            translation = 0 0 7.5
        :stop transformation:
    :stop geometry:

    # INFINITE CYLINDER OF RADIUS 2.5 mm
    :start geometry:
        library  = egs_cylinders
        type     = EGS_ZCylinders
        name     = sensitive_cylinder
        midpoint = 0 0 0
        radii    = 0.25 # cm
    :stop geometry:

    # PLANES FOR 0 to 100 mm CYLINDERS step of 1 mm
    :start geometry:
        library   = egs_planes
        type      = EGS_Zplanes
        name      = cylind_planes
        positions = 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1\
                      1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2\
                      2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3\
                      3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 4\
                      4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 5\
                      5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9 6\
                      6.1 6.2 6.3 6.4 6.5 6.6 6.7 6.8 6.9 7\
                      7.1 7.2 7.3 7.4 7.5 7.6 7.7 7.8 7.9 8\
                      8.1 8.2 8.3 8.4 8.5 8.6 8.7 8.8 8.9 9\
                      9.1 9.2 9.3 9.4 9.5 9.6 9.7 9.8 9.9 10
    :stop geometry:

    # SPLIT INFINITE CYLINDER WITH PLANES
    :start geometry:
        library    = egs_ndgeometry
        name       = phantom_cylinders
        dimensions = cylind_planes sensitive_cylinder
        :start media input:
            media = H2O521ICRU
        :stop media input:
    :stop geometry:

    # UNION INSTEAD OF ENVELOPE BUT HAVE TRIED BOTH  WAYS AND HAS SAME RESULT (FAR BELOW) #
    :start geometry:
        name        = water_phantom_not_moved
        library = egs_gunion
        geometries = phantom_cylinders surrounding_water_box
    :stop geometry:
:stop geometry definition:
.... later in the file

:start ausgab object definition:
  # WATER AGAIN
    :start ausgab object:
        library = egs_dose_scoring
        name = doseScoring
        volume = 1
        dose start region = 2     # first cylinder along z axis is region number 2
        dose stop region  = 101 # last cylinder along z axis is region number 101
   :stop ausgab object:
:stop ausgab object definition:

Then the output is:

======================================================
Dose Scoring Object(doseScoring)
======================================================
=> last case = 100000 fluence = 100000

==> Summary of region dosimetry (per particle)

 ir           medium  rho (g/cm^3)  Volume (cm^3)   Edep (MeV)                   D (Gy)                   
----------------------------------------------------------------------------------------------------------
  2  H2O521ICRU        1.00000000       1.000000    3.0056e-03  +/-  0.249  %    4.8149e-13  +/-  0.249  %
  3  H2O521ICRU        1.00000000       1.000000    4.9013e-02  +/-  2.048  %    7.8519e-12  +/-  2.048  %
  4  H2O521ICRU        1.00000000       1.000000    9.8427e-03  +/-  0.174  %    1.5768e-12  +/-  0.174  %
  5  H2O521ICRU        1.00000000       1.000000    9.7217e-02  +/-  1.574  %    1.5574e-11  +/-  1.574  %
  6  H2O521ICRU        1.00000000       1.000000    5.5793e-04  +/-  0.430  %    8.9380e-14  +/-  0.430  %
  7  H2O521ICRU        1.00000000       1.000000    6.5169e-03  +/-  4.937  %    1.0440e-12  +/-  4.937  %
  8  H2O521ICRU        1.00000000       1.000000    1.4542e-02  +/-  0.142  %    2.3296e-12  +/-  0.142  %
  9  H2O521ICRU        1.00000000       1.000000    3.1015e-01  +/-  0.780  %    4.9686e-11  +/-  0.780  %
 10  H2O521ICRU        1.00000000       1.000000    5.0571e-04  +/-  0.568  %    8.1015e-14  +/-  0.568  %
 11  H2O521ICRU        1.00000000       1.000000    1.8754e-02  +/-  2.610  %    3.0044e-12  +/-  2.610  %
 12  H2O521ICRU        1.00000000       1.000000    2.2330e+00  +/-  0.193  %    3.5772e-10  +/-  0.193  %
 13  H2O521ICRU        1.00000000       1.000000    1.4073e-03  +/-  5.535  %    2.2546e-13  +/-  5.535  %
 14  H2O521ICRU        1.00000000       1.000000    3.5947e-04  +/-  0.605  %    5.7587e-14  +/-  0.605  %
 15  H2O521ICRU        1.00000000       1.000000    2.2663e-02  +/-  2.426  %    3.6307e-12  +/-  2.426  %
 16  H2O521ICRU        1.00000000       1.000000    2.1591e-10  +/-  100.000%    3.4589e-20  +/-  100.000%
 17  H2O521ICRU        1.00000000       1.000000    2.5338e-04  +/-  10.333 %    4.0591e-14  +/-  10.333 %
 18  H2O521ICRU        1.00000000       1.000000    4.1515e-07  +/-  17.103 %    6.6507e-17  +/-  17.103 %
 19  H2O521ICRU        1.00000000       1.000000    1.2453e-02  +/-  2.094  %    1.9950e-12  +/-  2.094  %
 20  H2O521ICRU        1.00000000       1.000000    1.0782e+01  +/-  0.040  %    1.7273e-09  +/-  0.040  %
 21  H2O521ICRU        1.00000000       1.000000    0.0000e+00  +/-  100.000%    0.0000e+00  +/-  100.000%
 22  H2O521ICRU        1.00000000       1.000000    0.0000e+00  +/-  100.000%    0.0000e+00  +/-  100.000%
 23  H2O521ICRU        1.00000000       1.000000    0.0000e+00  +/-  100.000%    0.0000e+00  +/-  100.000%
 24  H2O521ICRU        1.00000000       1.000000    0.0000e+00  +/-  100.000%    0.0000e+00  +/-  100.000%
 25  H2O521ICRU        1.00000000       1.000000    0.0000e+00  +/-  100.000%    0.0000e+00  +/-  100.000%
 26  H2O521ICRU        1.00000000       1.000000    0.0000e+00  +/-  100.000%    0.0000e+00  +/-  100.000%
 27  H2O521ICRU        1.00000000       1.000000    0.0000e+00  +/-  100.000%    0.0000e+00  +/-  100.000%
 28  H2O521ICRU        1.00000000       1.000000    1.9939e-02  +/-  1.855  %    3.1943e-12  +/-  1.855  %
 29  H2O521ICRU        1.00000000       1.000000    1.0904e-10  +/-  51.395 %    1.7468e-20  +/-  51.395 %
 30  H2O521ICRU        1.00000000       1.000000    1.8126e-04  +/-  0.791  %    2.9038e-14  +/-  0.791  %
 31  H2O521ICRU        1.00000000       1.000000    7.4943e-10  +/-  100.000%    1.2006e-19  +/-  100.000%
 32  H2O521ICRU        1.00000000       1.000000    2.6512e-03  +/-  3.989  %    4.2472e-13  +/-  3.989  %
 33  H2O521ICRU        1.00000000       1.000000    1.0550e-06  +/-  20.718 %    1.6901e-16  +/-  20.718 %
 34  H2O521ICRU        1.00000000       1.000000    6.9709e-05  +/-  4.289  %    1.1167e-14  +/-  4.289  %
 35  H2O521ICRU        1.00000000       1.000000    0.0000e+00  +/-  100.000%    0.0000e+00  +/-  100.000%
 **Same message from 36-100 too**
101  H2O521ICRU        1.00000000       1.000000    0.0000e+00  +/-  100.000%    0.0000e+00  +/-  100.000%

When I use the scoring options with start calculation geometry it gives good results. I don't have the command line output but can attach a csv - it's just a depth dose for photons of 6 MV.

:start scoring options:
##### EACH CYLINDER SPECIFIED INDIVIDUALLY - so 100 of these definitions, see below
:start calculation geometry:
    geometry name  = water_phantom_not_moved
    cavity regions = 2 # repeated from 2 to 101
    cavity mass    = 1
:stop calculation geometry:
:stop scoring options:

To Reproduce

  1. Input 6 MeV electron beam (0.2 mm sigma x and y) with 0.5 mm Tungsten target or without target - happens in both cases
  2. Setup geometry in the way I have with output from ausgab, include enclosing geometry definition, MC, rng and source parts too
  3. Can make a direct comparison if you have the ausgab and scoring options together in the same egsinp file

Expected behavior Expect results to be consistent for dose between ausgab and scoring options versions of obtaining the dose in the regions

Operating system

EGSnrc version Most recent - installed from Github on 18/11/24

Notes

Each application may output certain scoring information like these quantities, in the output of the run or in data files. Each application is different with how it does this, so check the docs for the app you're using (e.g. https://nrc-cnrc.github.io/EGSnrc/doc/pirs898/egs_cbct.html for egs_cbct). Historically users actually edited the applications to score their quantities of interest, but this requires a fairly advanced understanding.

Having said that, all the egs++ apps can include the inputs for the 'ausgab' scoring options. You will want to look at the dose scoring one (https://nrc-cnrc.github.io/EGSnrc/doc/pirs898/classEGS__DoseScoring.html). This is limited in that it just gives you dose or energy deposited in regions or materials that you list. There is no kerma scoring here, but we will be adding an egs_kerma app soon. To get dose vs y-axis, you could only do this reliably if you're modelling an XYZ voxel geometry. If you are, then you can get it to output a 3ddose file, and then use a tool like statdose or CERR to extract the profile. I can explain how to do that further if you want.

rtownson commented 3 days ago

I can't reproduce the problem. Here is my input file, it gives identical dose in region 2 between ausgab dose scoring and a calculation geometry:

:start geometry definition:
# WATER PHANTOM #
    :start geometry:
        name        = surrounding_water_box
        library     = egs_box
        box size    = 15 15 15
        :start media input:
            media = H2O521ICRU
        :stop media input:
        :start transformation:
            translation = 0 0 7.5
        :stop transformation:
    :stop geometry:

    # INFINITE CYLINDER OF RADIUS 2.5 mm
    :start geometry:
        library  = egs_cylinders
        type     = EGS_ZCylinders
        name     = sensitive_cylinder
        midpoint = 0 0 0
        radii    = 0.25 # cm
    :stop geometry:

    # PLANES FOR 0 to 100 mm CYLINDERS step of 1 mm
    :start geometry:
        library   = egs_planes
        type      = EGS_Zplanes
        name      = cylind_planes
        positions = 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1\
                      1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2\
                      2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3\
                      3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 4\
                      4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 5\
                      5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9 6\
                      6.1 6.2 6.3 6.4 6.5 6.6 6.7 6.8 6.9 7\
                      7.1 7.2 7.3 7.4 7.5 7.6 7.7 7.8 7.9 8\
                      8.1 8.2 8.3 8.4 8.5 8.6 8.7 8.8 8.9 9\
                      9.1 9.2 9.3 9.4 9.5 9.6 9.7 9.8 9.9 10
    :stop geometry:

    # SPLIT INFINITE CYLINDER WITH PLANES
    :start geometry:
        library    = egs_ndgeometry
        name       = phantom_cylinders
        dimensions = cylind_planes sensitive_cylinder
        :start media input:
            media = H2O521ICRU
        :stop media input:
    :stop geometry:

    # UNION INSTEAD OF ENVELOPE BUT HAVE TRIED BOTH  WAYS AND HAS SAME RESULT (FAR BELOW) #
    :start geometry:
        name        = water_phantom_not_moved
        library = egs_gunion
        geometries = phantom_cylinders surrounding_water_box
    :stop geometry:

    simulation geometry = water_phantom_not_moved
:stop geometry definition:

:start scoring options:
    ##### EACH CYLINDER SPECIFIED INDIVIDUALLY - so 100 of these definitions, see below
    :start calculation geometry:
        geometry name  = water_phantom_not_moved
        cavity regions = 2 # repeated from 2 to 101
        cavity mass    = 1
    :stop calculation geometry:
:stop scoring options:

:start run control:
    #ncase = 1e9
    ncase = 1e5
    nbatch = 10
:stop run control:

:start ausgab object definition:
  # WATER AGAIN
    :start ausgab object:
        library = egs_dose_scoring
        name = doseScoring
        volume = 1
        dose start region = 2     # first cylinder along z axis is region number 2
        dose stop region  = 101 # last cylinder along z axis is region number 101
   :stop ausgab object:
:stop ausgab object definition:

:start source definition:
    :start source:
        library = egs_parallel_beam
        name = my_source
        :start shape:
            library = egs_circle
            radius   = 1
        :stop shape:
        direction = 0 0 1
        charge = 0
        :start spectrum:
            type = monoenergetic
            energy = 6
        :stop spectrum:
    :stop source:
:stop source definition:
erjft commented 2 days ago

Hi Reid, thanks for the reply. I had a look into the problem today, and I was able to reproduce the problem when I had an ausgab definition, but no calculation geometry definition for a given geometry - it could be a definition that references any region in that geometry. Otherwise it works fine when you have both the ausgab and the geometry definition as you have above. I tested even putting the region number of the surrounding water envelope and scoring the ausgab of the inscribed regions and that works fine but the cavity dose in the surrounding water envelope is then 0. I also realised having two calculation geometries and one ausgab seemed to mess with the resulting values. So in summary I need one calculation geometry and one ausgab together referencing the same region to give nice results.

:start geometry definition:
    :start geometry:
        name        = world_box
        library     = egs_box
        box size    = 50 50 70 
        :start media input:
            media = vacuum
        :stop media input:
        :start transformation:
            translation = 0 0 0
        :stop transformation:
    :stop geometry:
# WATER PHANTOM #
    :start geometry:
        name        = surrounding_water_box
        library     = egs_box
        box size    = 15 15 15
        :start media input:
            media = H2O521ICRU
        :stop media input:
        :start transformation:
            translation = 0 0 7.5
        :stop transformation:
    :stop geometry:

    # INFINITE CYLINDER OF RADIUS 2.5 mm
    :start geometry:
        library  = egs_cylinders
        type     = EGS_ZCylinders
        name     = sensitive_cylinder
        midpoint = 0 0 0
        radii    = 0.25 # cm
    :stop geometry:

    # PLANES FOR 0 to 100 mm CYLINDERS step of 1 mm
    :start geometry:
        library   = egs_planes
        type      = EGS_Zplanes
        name      = cylind_planes
        positions = 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1\
                      1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2\
                      2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3\
                      3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 4\
                      4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 5\
                      5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9 6\
                      6.1 6.2 6.3 6.4 6.5 6.6 6.7 6.8 6.9 7\
                      7.1 7.2 7.3 7.4 7.5 7.6 7.7 7.8 7.9 8\
                      8.1 8.2 8.3 8.4 8.5 8.6 8.7 8.8 8.9 9\
                      9.1 9.2 9.3 9.4 9.5 9.6 9.7 9.8 9.9 10
    :stop geometry:

    # SPLIT INFINITE CYLINDER WITH PLANES
    :start geometry:
        library    = egs_ndgeometry
        name       = phantom_cylinders
        dimensions = cylind_planes sensitive_cylinder
        :start media input:
            media = H2O521ICRU
        :stop media input:
    :stop geometry:

    # UNION INSTEAD OF ENVELOPE BUT HAVE TRIED BOTH  WAYS AND HAS SAME RESULT (FAR BELOW) #
    :start geometry:
        name        = water_phantom_not_moved
        library = egs_gunion
        geometries = phantom_cylinders surrounding_water_box
    :stop geometry:

      ########### EXAMPLE DISK
    :start geometry:
    library = egs_cones
    type    = EGS_ConeStack
    name    = example_disk
    axis    = 0 0 20 0 0 1 
    :start layer:
        thickness    = 0.5
        top radii    = 3
        bottom radii = 3
        media        = CU521ICRU 
    :stop layer:
    :stop geometry:
    ########### EXAMPLE DISK

     :start geometry:
        name        = full_simulation_geometry
        library     = egs_genvelope
        base geometry = world_box

        inscribed geometries = water_phantom_not_moved example_disk

        :start transformation:
            translation = 0 0 0
        :stop transformation:
    :stop geometry:

    simulation geometry = full_simulation_geometry

:stop geometry definition:

:start scoring options:
    :start calculation geometry:
        geometry name  = example_disk
        cavity regions = 0 
        cavity mass    = 1
    :stop calculation geometry:

    # UNCOMMENT FOR AUSGAB TO WORK
    #:start calculation geometry:
    #    geometry name  = water_phantom_not_moved
    #    cavity regions = 2
    #    cavity mass    = 1
    #:stop calculation geometry:
:stop scoring options:

:start run control:
    #ncase = 1e9
    ncase = 1e5
    nbatch = 10
:stop run control:

:start ausgab object definition:
  # WATER AGAIN
    :start ausgab object:
        library = egs_dose_scoring
        name = doseScoring
        volume = 1
        dose start region = 2     # first cylinder along z axis is region number 2
        dose stop region  = 101 # last cylinder along z axis is region number 101
   :stop ausgab object:
:stop ausgab object definition:

:start source definition:
    :start source:
        library = egs_parallel_beam
        name = my_source
        :start shape:
            library = egs_circle
            radius   = 1
        :stop shape:
        direction = 0 0 1
        charge = -1
        :start spectrum:
            type = monoenergetic
            energy = 6
        :stop spectrum:
    :stop source:
:stop source definition:
rtownson commented 2 days ago

I tried with just an ausgab object in egs_app, and just a calculation geometry in egs_chamber, and the results agree within statistical uncertainty. Of course the results won't be identical, because the two apps take different logical paths which means the same random number seeds don't result it the same simulation. Do your results disagree outside the uncertainties?