How can I model an x-ray system in EGSnrc?

[bdnkth]

Hi,

I'm pretty new to the world of EGSnrc. For our upcoming xray project we want to run simulations in order to design our radiation protection.

I've started implementing our machine in the egs world which so far has worked pretty well. However, when it comes to simulating the radiation, the calculated values seem way to low. I've tried to fix this but so far I was not able to. So hopefully someone here can help me.

I'm currently using an implementation which is based on the egs_chamber example. I've stumpled across the BEAMnrc package which seems to be better for xray but is not as user friendly as the egs_chamber stuff. Can I do my simulations with the egs_chamber or do I have to switch to BEAMnrc? Our xray source can handle up to 70 kV so that's the maximum which I have to simulate. Currently I'm using our measured spectrum of the xray source as the radiation source in egsnrc.

[rtownson]

Either egs_chamber or BEAMnrc can be used for this, and both have some advantages or disadvantages. Since you're already working with egs_chamber, you might as well stay there. They use the same transport physics and data!

Can you provide some numbers - how do the experimental numbers compare? Note that MC simulations are best used for relative comparisons, and all results are reported "per source fluence" or "per initial particle". Calculating absolute dose with simulations is generally not as reliable, since in the real world you don't know exactly how many source particles there were.

[bdnkth]

Good to hear that egs_chamber also works.

1 PB700ICRU       11.340   3430.0000 3.0951e-03 +/- 0.240  % 1.2747e-17 +/- 0.240  %

These are the results for a lead plate I've placed approximately 5cm below the radiation source. For our understanding the dose area product is really small.

[rtownson]

What are you comparing to? Like I said, this is "dose per source particle", so it will always be a very small number.

[ftessier]

@bdnkth, the expected dose result would be on the order of the energy deposited in the lead by a single representative source photon, divided by the mass of the lead plate (in Joule/kg). Here your results show that the energy deposited on average is 3.0951e-03 +/- 0.240 % MeV, so about 3 keV. Does this make sense, theoretically, according to the cross section of ~70 keV photons incident on a thin lead plate?

[ftessier]

What is the thickness of your lead plate?

[ftessier]

Presumably, the average energy of the photons is also much less than 70 keV, if that is the maximum x-ray spectrum energy. What does the spectrum look like?

[ftessier]

Also, if you are using an isotropic source, the results are normalized by the total number of emitted photons, in 4*pi. To get a result normalized by fluence within a given aperture, you can use a collimated source.

[ftessier]

You may also want to look at particle tracks in egs_view to get a sense, visually, of whether or not the lead shield is stopping photons, and to what extent.

[bdnkth]

The lead plate is defined as box size = 64.7 45.7 1 and the source is above the plate in z-dimension. That the deposited energy is from a single source photon explains why I was expecting a higher dosis.

I'm currently using a pencil_beam source, I'll try with something different.

I've attached the used spectrum

[ftessier]

Sounds good. A pencil beam is fine, I was just worried that you were using an isotropic source, which would decreased the dose per photon on account of the particles that are emitted away from the plate.

Results are normalized by particle (or fluence), because we want the results to converge. We typically adjust the number of histories to reach the sought level of uncertainty, and we don't want the dose value to change all the time, depending on N. The history number is in fact the only unit of "time" in these simulations, so you can view the results as rates (per particle emitted).

[bdnkth]

Okay, that makes sense. Since the deposited energy is per photon, is there a way to calculate the number of particles hitting the plate? I guess there is, probably I can't see the forest for the trees

[bdnkth]

Any updates on this? I'm still struggling to find out the number of photons hitting my object

[ftessier]

The geometry lies on the +z axis, but the particles are emitted along -z: direction = 0 0 -1.

[ftessier]

Oh I see you are emitting towards the thin steel plate, which is in -z. The issue is that you have to set the calculation geometry to machine, instead of below_source:

    :start calculation geometry:
        geometry name = machine
        cavity regions = 8
        cavity mass = 100000
    :stop calculation geometry:

The geometry name is the name of the entire composite geometry, and then you select the scoring regions (the cavity regions) within that geometry.

[ftessier]

Note that you cannot have strictly overlapping planes in egs++, from two disjoint geometries. These will show up as speckle pattern, showing when particles hit or miss the second plane (owing to the limit of floating-point precision). For example, in you geometry:

[ftessier]

Another note: the union is not a very efficient geometry, because at every step the code must calculate the distance to all geometries in the union, since in general there is not constraint on where they are located with respect to each other.

[ftessier]

Here is a more efficient way to build the main box:

    :start geometry:
        name = conveyor_channel
        library = egs_ndgeometry
        type = EGS_XYZGeometry
        x-planes = -35.15 -34.95 34.95 35.15
        y-planes = -64 -63.7 63.7 64
        z-planes = 0 0.25 0.75 2 29.7 30
        :start media input:
            media = STEEL700ICRU AIR700ICRU
            set medium = 4 1
            set medium = 22 1
            set medium = 28 1
            set medium = 31 1
            set medium = 34 1
        :stop media input:
    :stop geometry:
    simulation geometry = conveyor_channel

(clipped to show internal structure)

[ftessier]

@bdnkth here is a better way to contruct your geometry, without unions. The idea is to build the shield, the conveyor, and the radiation box, and then "glue" them together as three layers along z using a CD geometry (there is a sligh complication because the radiation box has to "punch" through the conveyor):

:start geometry definition:

    ## big air box
    :start geometry:
        name     = air_box
        library  = egs_box
        box size = 200  200 400
        :start media input:
            media = AIR700ICRU
        :stop media input:
    :stop geometry:

    ############################################################################v
    ### conveyor
    ############################################################################

    ### conveyor channel
    :start geometry:
        name     = conveyor
        library  = egs_ndgeometry
        type     = EGS_XYZGeometry
        x-planes = -35.15 -34.95 34.95 35.15
        y-planes = -64 -63.7 63.7 64
        z-planes = -40 0.25 0.75 2 29.7 40
        :start media input:
            media = STEEL700ICRU AIR700ICRU
            set medium = 4 1
            set medium = 22 1
            set medium = 28 1
            set medium = 31 1
            set medium = 34 1
        :stop media input:
    :stop geometry:

    ### hole in top
    :start geometry:
        name     = hole
        library  = egs_ndgeometry
        type     = EGS_XYZGeometry
        x-planes = -32.2 32.2
        y-planes = -22.7 22.7
        z-planes = 29 39
        :start media input:
            media = AIR700ICRU
        :stop media input:
    :stop geometry:

    ### conveyor with hole
    :start geometry:
        name    = hole_in_conveyor
        library = egs_genvelope
        base geometry = conveyor
        inscribed geometries = hole
    :stop geometry:

    ### conveyor in air
    :start geometry:
        name    = converyor_box_in_air
        library = egs_genvelope
        base geometry = air_box
        inscribed geometries = hole_in_conveyor
    :stop geometry:


    ############################################################################
    ### radiation box
    ############################################################################

    ### radiation box
    :start geometry:
        name     = radiation_box
        library  = egs_ndgeometry
        type     = EGS_XYZGeometry
        x-planes = -32.5 -32.2 32.2 32.5
        y-planes = -23 -22.7 22.7 23
        z-planes = 0 85.2 85.5
        :start media input:
            media = STEEL700ICRU AIR700ICRU
            set medium = 4 1
        :stop media input:
    :stop geometry:

    ### radiation box in air
    :start geometry:
        name    = radiation_box_in_air
        library = egs_genvelope
        base geometry = air_box
        inscribed geometries = radiation_box
    :stop geometry:


    ############################################################################
    ### shield
    ############################################################################

    ### shield
    :start geometry:
        name     = shield
        library  = egs_ndgeometry
        type     = EGS_XYZGeometry
        x-planes = -35.15 35.15
        y-planes = -15 15
        z-planes = -1 10
        :start media input:
            media = PB700ICRU
        :stop media input:
    :stop geometry:

    ### shield in air
    :start geometry:
        name    = shield_in_air
        library = egs_genvelope
        base geometry = air_box
        inscribed geometries = shield
    :stop geometry:


    ############################################################################
    ### machine
    ############################################################################

    ### machine planes
    :start geometry:
        name      = machine_planes
        library   = egs_planes
        type      = EGS_Zplanes
        positions = -100 0 30 100
    :stop geometry:

    ### machine
    :start geometry:
        name    = machine
        library = egs_cdgeometry
        base geometry = machine_planes
        set geometry = 0 shield_in_air
        set geometry = 1 converyor_box_in_air
        set geometry = 2 radiation_box_in_air
    :stop geometry:

    # ### simulation geometry
    simulation geometry = machine

:stop geometry definition:

Here is a snapshot (with energy set to 300 keV to get some tracks, and cavity regions set to 1):

[ftessier]

Running 1000000 histories

  Batch             CPU time        Result   Uncertainty(%)
==========================================================
      1                0.56     1.3614e-14           0.07
      2                1.11    1.36192e-14           0.05
      3                1.66    1.36209e-14           0.04
      4                2.21    1.36188e-14           0.04
      5                2.76    1.36255e-14           0.03
      6                3.30    1.36255e-14           0.03
      7                3.85    1.36273e-14           0.03
      8                4.40    1.36265e-14           0.03
      9                4.95    1.36259e-14           0.02
     10                5.49     1.3626e-14           0.02


Finished simulation

Total cpu time for this run:            5.50 (sec.) 0.0015(hours)
Histories per hour:                     6.54545e+08   
Number of random numbers used:          150854217     
Number of electron CH steps:            5.50715e+06   
Number of all electron steps:           5.52284e+06   


 last case = 1000000 fluence = 33333.3

Geometry                        Cavity dose      
-----------------------------------------------
machine                   1.3626e-14 +/- 0.023  % 

[bdnkth]

Tank you for that lenthly response. I think I'll need some time to process all that information. Originally I've started with the xyzgeometries but I've struggled with their syntax and switched to the boxes since they are really straight forward to work with.

Can you check if this file connected all your code in the right way? Since I do not get the same result as you

[rtownson]

When you changed the position input to a transformation did you include the entire transformation input block in your shape input?
https://nrc-cnrc.github.io/EGSnrc/doc/pirs898/classEGS__AffineTransform.html

e.g.

        :start shape:
            type = box
            box size = 3 3 1
            :start transformation:
                translation =  0 -21.125 41.5 # in cm
            :stop transformation:
        :stop shape:

[bdnkth]

Yes I did. I also tried it without having the transformation, to see if that's the problem, but somehow the problem prevailed. I've attached the file, to show it
currentStatus.txt
.

[rtownson]

Are you running the simulation using egs_app? The inputs you have for scoring, under :start scoring options: are only used by egs_chamber. Your input file gives results for me when I use egs_chamber.

[bdnkth]

I guess that's the problem. I was laking some code from the egs_chamber.cpp in my project cpp file which I now added. But then make failes with
Mortran compiling EGSnrc sources ...
gfortran -O2 -mtune=native -fPIC -c -o egsnrc_linux64.o egsnrc_linux64.F
egsnrc_linux64.F:8:9:

8 |       d %C80
  |         1

Error: Symbol ‘d’ at (1) has no IMPLICIT type
make: *** [/mnt/c/Users/Kath/EGSnrc/HEN_HOUSE/makefiles/cpp_makefile:70: egsnrc_linux64.o] Error 1

[rtownson]

I'm not sure why you're seeing this, there must be a typo or bug in the code you have added. Why are you writing your own code, what are you trying to do that egs_chamber can't?