egs++ cavity: estimating half-value layer (HVL) #939

walleludvig · 2022-11-17T04:41:34Z

walleludvig
Nov 17, 2022

Hi all,

I seek help to understand how I can implement the HVL scoring option in egs cavity package appropriately as I keep getting this wrong.

I wish to estimate the HVL of an X-ray beam penetrating different absorber thicknesses where the absorber is made of PMMA, and the air kerma rate is scored at the same reference plane (at a distance behind the absorber) for each geometry. I wish to do this using a phase space source input, however for troubleshooting problems with the cavity .egsinp file I am using a collimated point source. The most recent attempted input file looks as follows:

###############################################################################
#
#  PMMA HVL Scoring
#
###############################################################################

:start run control:

    #############################################
    # Define the run control
    #############################################

  ncase = 2E6
:stop run control:

###############################################################################

:start geometry definition:

#############################################
#   Define the air rectlinear cylinder to place phantom and detector within
#############################################
    :start geometry:
        library = egs_ndgeometry
        type = EGS_XYZGeometry
        name = surrounding_air
        x-planes = -50 50
        y-planes = -50 50
        z-planes = -0.5 50
        :start media input:
            media = air
        :stop media input:
    :stop geometry


############### Define the phantoms

    :start geometry:
        library = egs_ndgeometry
        type = EGS_XYZGeometry
        name = p_14
        x-planes = -15 15
        y-planes = -15 15
        z-planes = 0.5 14.5
        :start media input:
            media = pmma
        :stop media input:
    :stop geometry:
    ...
    :start geometry:
        library = egs_ndgeometry
        type = EGS_XYZGeometry
        name = p_15
        x-planes = -15 15
        y-planes = -15 15
        z-planes = 0.5 15.5
        :start media input:
            media = pmma
        :stop media input:
    :stop geometry:
    ...
    :start geometry:
        library = egs_ndgeometry
        type = EGS_XYZGeometry
        name = p_16
        x-planes = -15 15
        y-planes = -15 15
        z-planes = 0.5 16.5
        :start media input:
            media = pmma
        :stop media input:
    :stop geometry:
    ...
    :start geometry:
        library = egs_ndgeometry
        type = EGS_XYZGeometry
        name = p_17
        x-planes = -15 15
        y-planes = -15 15
        z-planes = 0.5 17.5
        :start media input:
            media = pmma
        :stop media input:
    :stop geometry:
    ...
    :start geometry:
        library = egs_ndgeometry
        type = EGS_XYZGeometry
        name = p_18
        x-planes = -15 15
        y-planes = -15 15
        z-planes = 0.5 18.5
        :start media input:
            media = pmma
        :stop media input:
    :stop geometry:

############### Define the scoring plane

    :start geometry:
        library = egs_ndgeometry
        type = EGS_XYZGeometry
        name = air_scoring
        x-planes = -21.5 21.5
        y-planes = -21.5 21.5
        z-planes = 30.5 33.7
        :start media input:
            media = air
        :stop media input:

        set label = scoring_plane 0
    :stop geometry:


############### Reference Geometry

    :start geometry:
        library = egs_genvelope
        name = p_0in
        base geometry = surrounding_air
        inscribed geometries = air_scoring
    :stop geometry:

############### Place phantoms and scoring plane in situ

    :start geometry:
        library = egs_genvelope
        name = p_14in
        base geometry = surrounding_air
        inscribed geometries = p_14 air_scoring
    :stop geometry:
    ...
    :start geometry:
        library = egs_genvelope
        name = p_15in
        base geometry = surrounding_air
        inscribed geometries = p_15 air_scoring
    :stop geometry:
    ...
    :start geometry:
        library = egs_genvelope
        name = p_16in
        base geometry = surrounding_air
        inscribed geometries = p_16 air_scoring
    :stop geometry:
    ...
    :start geometry:
        library = egs_genvelope
        name = p_17in
        base geometry = surrounding_air
        inscribed geometries = p_17 air_scoring
    :stop geometry:
    ...
    :start geometry:
        library = egs_genvelope
        name = p_18in
        base geometry = surrounding_air
        inscribed geometries = p_18 air_scoring
    :stop geometry:

    simulation geometry = p_18in

:stop geometry definition:

#############################################
#   Define the source input
#############################################

:start source definition:

    :start source:
        name = spectrum
        library = egs_phsp_source
        phase space file = /Users/ludvigwalle/EGSnrc/egs_home/BEAM_XSample_mod/XSample_5pRe_120SSD.egsphsp1
        particle type = photons
        cutout = -15 15 -15 15
        ## recycle photons = 100
    :stop source:

    :start source:
        library = egs_collimated_source
        name = coll_source
        :start source shape:
            type = point
            position = 0 0 -29.5
        :stop source shape:
            :start target shape:
                library   = egs_rectangle
                rectangle = -15 -15 15 15
            :stop target shape:
            distance = 30
            charge = 0
        :start spectrum:
            type = monoenergetic
            energy = 0.1
        :stop spectrum:
    :stop source:

simulation source = coll_source

:stop source definition:

#############################################
#   Define scoring
#############################################

:start ausgab object definition: # Output options
    :start ausgab object: # Track scoring (view in egs_view)
        name    = my_tracks
        library = egs_track_scoring
    :stop ausgab object:
:stop ausgab object definition:

############### HVL scoring

:start scoring options:

    calculation type = HVL

    :start calculation geometry:
        geometry name = p_0in
        cavity regions = scoring_plane
        cavity mass = 7.1285
    :stop calculation geometry:
    ...
    :start calculation geometry:
        geometry name = p_14in
        cavity regions = scoring_plane
        cavity mass = 7.1285
    :stop calculation geometry:
    ...
    :start calculation geometry:
        geometry name = p_15in
        cavity regions = scoring_plane
        cavity mass = 7.1285
    :stop calculation geometry:
    ...
    :start calculation geometry:
        geometry name = p_16in
        cavity regions = scoring_plane
        cavity mass = 7.1285
    :stop calculation geometry:
    ...
    :start calculation geometry:
        geometry name = p_17in
        cavity regions = scoring_plane
        cavity mass = 7.1285
    :stop calculation geometry:
    ...
    :start calculation geometry:
        geometry name = p_18in
        cavity regions = scoring_plane
        cavity mass = 7.1285
    :stop calculation geometry:

    correlated geometries = p_0in p_14in
    correlated geometries = p_0in p_15in
    correlated geometries = p_0in p_16in
    correlated geometries = p_0in p_17in
    correlated geometries = p_0in p_18in

    :start HVL scoring:
        scoring circle = 0 0 30.5 0.5
        scoring plane normal = 0 0 1
        muen file = /Users/ludvigwalle/EGSnrc/egs_home/cavity/muen_rho_pmma.txt
        #absorber thicknesses = 0 14 15 16 17 18
        scatter = no
    :stop HVL scoring:


############### Fluence scoring

    :start fluence scoring:
        minimum energy = 0.01
        maximum energy = 0.2
        number of bins = 200
        scale          = linear
    :stop fluence scoring:

:stop scoring options:

############### Pegsless mode media definition

:start media definition:

    ae = 0.512
    ue = 2.511
    ap = 0.001
    up = 2.00


    :start pmma:
        density correction file     = /Users/ludvigwalle/EGSnrc/HEN_HOUSE/pegs4/density_corrections/compounds/polymethylmethacrylate__lucite___perspex___plexiglas_.density
        bremsstrahlung correction= NRC
    :stop pmma:

    :start air:
        density correction file     = /Users/ludvigwalle/EGSnrc/HEN_HOUSE/pegs4/density_corrections/compounds/air_dry_nearsealevel.density
        bremsstrahlung correction= NRC
    :stop air:

:stop media definition:

My problem is that I get non-sensible results from the simulation with > +/- 100% uncertainty at times (sometimes inf). Please see example of an output copy pasted below:

Running 2000000 histories

  Batch             CPU time        Result   Uncertainty(%)
==========================================================
      1                9.64    7.49635e-14          16.18
      2               19.27    9.35801e-14           9.75
      3               28.90    9.73008e-14           7.82
      4               38.51    9.48924e-14           7.06
      5               48.11      9.662e-14           6.17
      6               57.72    9.83053e-14           5.62
      7               67.32    9.58092e-14           5.23
      8               76.94    9.32565e-14           4.97
      9               86.55    9.32503e-14           4.64
     10               96.19    9.21308e-14           4.41


Finished simulation

Total cpu time for this run:            96.19 (sec.) 0.0267(hours)
Histories per hour:                     7.48519e+07   
Number of random numbers used:          4411095479    
Number of electron CH steps:            2.0351e+08    
Number of all electron steps:           2.64603e+08   


 last case = 2000000 fluence = 2222.22


***** This is a HVL calculation.
      scatter is NOT included!
      KERMA is scored in a circle with radius 0.5
      with midpoint (0,0,30.5) in a plane with normal (0,0,1)

Geometry                           KERMA         
-----------------------------------------------
p_0in                     8.364475e-13 +/- 4.40540  % 
p_14in                    2.190993e+254 +/- inf      % 
p_15in                    5.035280e+253 +/- inf      % 
p_16in                    1.074350e+254 +/- inf      % 
p_17in                    1.489751e+254 +/- inf      % 
p_18in                    5.057166e+253 +/- inf      % 


Geometry 1           Geometry 2              KERMA ratio
p_0in                p_14in                   0.00000  +/- inf    
p_0in                p_15in                   0.00000  +/- inf    
p_0in                p_16in                   0.00000  +/- inf    
p_0in                p_17in                   0.00000  +/- inf    
p_0in                p_18in                   0.00000  +/- inf

The full output is attached here for reference:
example_output.txt.

I have read the ins -and outs of the documentation pages (https://nrc-cnrc.github.io/EGSnrc/doc/pirs898/cavity.html#cavity_options), however there is obviously something I am misunderstanding.

Any guidance and help indicating where I am going wrong would be greatly appreciated!

Answered by rtownson

Nov 17, 2022

Just an idea from first glance - have you tried scoring plane normal = 0 0 -1? Since your source is coming in from -z toward +z, that makes more sense to me.

View full answer

rtownson · 2022-11-17T14:11:55Z

rtownson
Nov 17, 2022
Maintainer

Just an idea from first glance - have you tried scoring plane normal = 0 0 -1? Since your source is coming in from -z toward +z, that makes more sense to me.

2 replies

walleludvig Nov 17, 2022
Author

Thanks for your suggestion @rtownson. I had assumed that choosing between scoring plane normal = 0 0 -1 and scoring plane normal = 0 0 1 would be arbitrary in this case because it's scoring the photon fluence crossing the scoring plane, but I'm probably wrong? - Either or, I have made the suggested change and re-ran but issue was not resolved.

I have been wondering however, does my choice of simulation geometry matter for the scoring of air kerma? i.e., should it matter whether I choose simulation geometry = <some geometry with an absorber> or simulation geometry = <geometry without absorber> other than for viewing in egs_view?

rtownson Nov 18, 2022
Maintainer

I ran this input file with my own muen file (for air, actually), and it gave normal results. So I guess there is an issue with the muen file. The choice of simulation geometry only effects egs_view for cavity and egs_chamber input files.

 last case = 2000000 fluence = 2222.22


***** This is a HVL calculation.
      scatter is NOT included!
      KERMA is scored in a circle with radius 0.5
      with midpoint (0,0,30.5) in a plane with normal (0,0,1)

Geometry                           KERMA         
-----------------------------------------------
p_0in                     8.205240e-13 +/- 4.63120  % 
p_14in                    1.614223e-13 +/- 9.85765  % 
p_15in                    1.563332e-13 +/- 9.72312  % 
p_16in                    1.510532e-13 +/- 9.68375  % 
p_17in                    1.191270e-13 +/- 11.02534 % 
p_18in                    1.004478e-13 +/- 13.55423 % 


Geometry 1           Geometry 2              KERMA ratio
p_0in                p_14in                   5.08309  +/- 0.55313
p_0in                p_15in                   5.24856  +/- 0.56486
p_0in                p_16in                   5.43202  +/- 0.58270
p_0in                p_17in                   6.88781  +/- 0.82327
p_0in                p_18in                   8.16866  +/- 1.16961

The muen format should look like this:

399
   0.0010000   3.6039495
   0.0010140   3.5195360
   0.0010280   3.4377441
   0.0010420   3.3588822
   0.0010570   3.2775164
...etc for 399 items, plus a blank new line at the end...

walleludvig · 2022-11-19T02:16:37Z

walleludvig
Nov 19, 2022
Author

Hi @rtownson,

Thanks very much for your input. I have looked into the problem with the muen file and believe the problem either arose from:

case sensitivity of scientific formatting in the numbers. i.e., E*muen numbers were specified in format 2.484955e-02 instead of 2.484955E-02
OR:
some energies listed in the file was missing the pre-decimal 0. i.e., 0.07 was written as .07.

Running with a corrected muen file I get results with expected trends in air kerma rates.

Again thank you kindly for looking into and resolving this.

1 reply

rtownson Jan 3, 2023
Maintainer

Could you send along your muen file that was not working? I tried both of your suggestions of the problem and they worked for me.

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

egs++ cavity: estimating half-value layer (HVL) #939

{{title}}

Replies: 2 comments 3 replies

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

Select a reply

egs++ cavity: estimating half-value layer (HVL) #939

walleludvig Nov 17, 2022

Replies: 2 comments · 3 replies

rtownson Nov 17, 2022 Maintainer

walleludvig Nov 17, 2022 Author

rtownson Nov 18, 2022 Maintainer

walleludvig Nov 19, 2022 Author

rtownson Jan 3, 2023 Maintainer

walleludvig
Nov 17, 2022

Replies: 2 comments 3 replies

rtownson
Nov 17, 2022
Maintainer

walleludvig Nov 17, 2022
Author

rtownson Nov 18, 2022
Maintainer

walleludvig
Nov 19, 2022
Author

rtownson Jan 3, 2023
Maintainer