Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

egs_chamber dose output for cylinders inscribed in water phantom #1220

Open
erjft opened this issue Nov 19, 2024 · 5 comments
Open

egs_chamber dose output for cylinders inscribed in water phantom #1220

erjft opened this issue Nov 19, 2024 · 5 comments
Labels
bug

Comments

@erjft
Copy link

erjft commented Nov 19, 2024 •
edited
Loading

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

  • macOS Sequoia 15.1

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.

  • I tried with loads of particles - but it doesn't make sense the doses anyway in the scored regions, it should be a depth dose so build up and come down
  • Has ausgab output got a maximum number of regions it can score?
  • Why do you think it misses region 21-27?
  • I updated my MacOS to 15.1 Sequoia and the egs_view doesn't work now but can still run my code so can't show images (doing make in $HEN_HOUSE/egs++/view had issues, so I had to edit the Makefile_osx, so instead of the default it is now
    CC = /opt/homebrew/opt/gcc@14/bin/gcc-14 # was /Library/Developer/CommandLineTools/usr/bin/gcc
    CXX = /opt/homebrew/opt/gcc@14/bin/g++-14 # was /Library/Developer/CommandLineTools/usr/bin/g++
    LINK = /opt/homebrew/opt/gcc@14/bin/g++-14 #was /Library/Developer/CommandLineTools/usr/bin/g++
    then I had to create aliases to the files in my ~/.zshrc:
    alias egs_view=$HEN_HOUSE/bin/osx/egs_view.app/Contents/MacOS/egs_view
    alias egs_gui=$HEN_HOUSE/bin/osx/egs_gui.app/Contents/MacOS/egs_gui
    alias egs_inprz=$HEN_HOUSE/bin/osx/egs_inprz.app/Contents/MacOS/egs_inprz
    So fixed a different problem here)
@erjft erjft added the bug label Nov 19, 2024
@rtownson
Copy link
Collaborator

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
Copy link
Author

erjft commented Nov 20, 2024 •
edited
Loading

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
Copy link
Collaborator

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?

@erjft
Copy link
Author

erjft commented Nov 27, 2024

Thanks for the message. Yeah my results are slightly different between the ausgab and the scoring options value for dose, not by much at all and when I do a ratio of the dose D1/D2 for ausgab and D1/D2 for the scoring options value its practically identical for the same random numbers, in the same run. I'm not particularly worried about it. I realised that an issue was occurring because when I put a label on a water phantom, then moved it - the unmoved water phantom was being scored instead of the moved one, so my results were a little odd. Instead I fixed this by just waiting until I had my complete model, making sure to instantiate my water phantom first in the largest world geometry (so the region numbers were smallest), and finding these numbers in the world geometry, creating a label there, then getting output by referencing this world geometry and the corresponding label

@rtownson
Copy link
Collaborator

Well the label should certainly work when you put in directly in the lowest level geometry, not in the world geometry. That's the whole point of labels! So I'm not sure what the issue was there. To extract the local region numbers to use in the label, did you ensure you visualized only that specific geometry in egs_view, using the simulation geometry drop box in the GUI?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
bug
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
2 participants
@rtownson @erjft