Question about doses varying between using original source and derived phasespace as source of particles

[nexustowers]

I am using egs_chamber from EGSnrc 2025 to create a simple model of a linac head with a phantom underneath. The head model consists of a target, primary collimator etc etc and the primary source of particles is an electron beam hitting the target creating bremstrahlung. The surrounding medium is air.

The head assembly is surrounded by a spherical shell of air which I am using to create a IAEA phase space of particles exiting the shell using the from-regions / to-regions method . I am interested in head leakage, hence the spherical phase space geometry instead of a plane under the jaws. I am tracking all particle types.

The phantom is underneath the jaws and outside the phase space shell. I am calculating the dose to a dosebox on the cax inside the phantom at some depth.

I run the simulation with ncase, say, 1e8 electrons, and this creates a phasespace with original particles close to (but not exactly) 1e8 and, say, 1e6 total particles in the phasespace. I note the dose to the dosebox.

I run the simulation again, this time with the phasespace as particle source and set ncase in my .egsinp file to 1e6. I understand this should use all the particles in the phasespace file once and normalize the results to the number of original particles used to create the phase space file ie 1e8. I am not doing any recycling or variance reduction.

Since using the phase space is essentially recreating the original scenario and the normalization should be consistent, I am expecting the dose reported to the dosebox from the phasespace simulation to be the same as the first simulation (within uncertainties, which in the case I have calculated is about 15%) but it is not - it is about an order of magnitude different with the phase space simulation dose being the higher. Can anyone please explain any conceptual reasons why they're not the same?

[rtownson]

Are you sure that you're actually capturing the entire source in the phase-space that you're scoring? If some of the particles managed to get around the phase-space scoring "to regions" and "from regions" pairs, then the normalization will end up different.

Just to be sure, I've tested this myself. An enclosed sphere for me scored to a phase-space, which is then used as the source, gives the same dose. However, if I use a source only partially covered by the phase-space scoring regions, the dose ends up different.

[nexustowers]

Hi Reid, thanks for your reply. Apologies for my delay in replying but I needed to park this somewhere for a while to try to clear my thoughts a bit.

In quick answer, yes, I'm pretty sure my phase space surface enclosed the whole of the source sphere.

Anyway, I tried reducing the scenario to a very simple geometry to test the concept. I got rid of the sphere and the complications and used a 6MV photon parallel beam source aimed at a phantom with a dosebox inside. I put a wide but thin (in the beam direction) box just above the phantom to score the phasespace. Running the 2 sims, one with the parallel beam only and the other with the phase space still gave different doses (this time by a factor of about 1000). The phase space scoring box is easily of enough extent to score all the parallel beam. The number of particles in the phase space which is created is almost identical to the number of particles in the parallel beam, which is as expected.

Here is my script if you'd be willing to have a look. It should just run

I feel like I must be doing something fundamentally very wrong, but just can't pin down where it is..

thank you

#######################################################################
#######################################################################
#  
#       Simple test egs_chamber script to investigate equivalence of dose calculated
#       from parallel beam source and from derived phasespace

#       Parallel beam (6MV photon) with circular cross-section impinges onto phantom which
#       contains dose box at 1cm deep.  A phase space is scored in the planar_phsp_box object
#       situated approx 3cm above the phantom and can also be used as source to generate the dose

 
#       (1)   run with "parallel_beam_photons" source selected and phase space ausgab object
#               turned on to derive dose due to the parallel beam and also generate the phase space
#       (2)   run with "parallel_beam_phasespace" source selected and phase space ausgab object
#               commented out to generate dose due to phasespace


#######################################################################
# Media (pegsless mode)
#######################################################################

:start media definition:

    AP = 0.005
    UP = 13
    AE = 0.521
    UE = 13.511


    :start air:
        density correction file = air_dry_nearsealevel
    :stop air:

    :start air2:
        density correction file = air_dry_nearsealevel
    :stop air2:

    :start iso_water:
        density correction file = water_liquid
    :stop iso_water:

    :start icru_tissue:
        density correction file = tissue_soft_icrufour-component
    :stop icru_tissue:

 
:stop media definition:



######################################################################
# RNG and run control
######################################################################

:start rng definition:

    initial seeds = 91 2556
    type = ranmar
    high resolution = yes

:stop rng definition:



:start run control:
                                
    ncase = 5000
    nbatch = 10
    calculation = first
    geometry error limit = 1
                         
:stop run control:



############################################################################
# Source Definition
############################################################################

:start source definition:

   :start source:
        name = parallel_beam_photons
        library = egs_parallel_beam
        :start shape:
            type = cylinder
            radius = 10
            height = 0.1
            axis = 0 0 1
            midpoint = 0 0 100
        :stop shape:
        direction = 0 0 -1
        charge = 0
        :start spectrum:
            type = monoenergetic
            energy = 6
        :stop spectrum:
    :stop source:


   :start source:
        name = parallel_beam_phasespace
        library = iaea_phsp_source
        iaea phase space file = zbunker_head_leakage_phasespace_10X.1
        particle type = all
    :stop source:


    simulation source = parallel_beam_photons            ###### Select source  ########
#   simulation source = parallel_beam_phasespace

:stop source definition:



###############################################################################
# Build the geometry   NB isocentre is at the origin of the coordinate system
###############################################################################

:start geometry definition:

    :start geometry:                                    
        name = wider_world
        library = egs_box
        box size = 2000 1200 1000
            :start media input:
                media = air
            :stop media input:
        :stop geometry:


    # Create iso_box - 50x50x30cm phantom at 100cm fsd with dosebox at 1cm deep

    :start geometry:                           
        name = iso_box
        library = egs_ndgeometry
            type = EGS_XYZGeometry
            x-planes = -25 -5 5 25
            y-planes = -25 -5 5 25
            z-planes = -30 -1.10 -0.9 0
            :start media input:
                media = iso_water icru_tissue
                set medium = 13 1
            :stop media input:
            set label = iso_box_label 13
    :stop geometry:



    # Create box just above iso_box phantom for planar phase space generation

    :start geometry:                           
        name = planar_phsp_box
        library = egs_ndgeometry
            type = EGS_XYZGeometry
            x-planes = -30 30
            y-planes = -30 30
            z-planes = 3 4         
            :start media input:
                media = air2
            :stop media input:
    :stop geometry:



    # create simulation geometry
 
    :start geometry:
        name = test_geometry
        library = egs_gunion
        geometries = planar_phsp_box iso_box wider_world
    :stop geometry:


    simulation geometry = test_geometry


:stop geometry definition:


#################################################################
# Scoring options - NB these are egs_chamber specific
#################################################################

    :start scoring options:

            :start calculation geometry:
                geometry name = test_geometry
                cavity regions = iso_box_label
                cavity mass = 20
                cavity geometry = iso_box
                enhance regions = iso_box_label
                enhancement = -64
            :stop calculation geometry:

    :stop scoring options:


############################################################################
### AUSGAB OBJECTS
############################################################################

:start ausgab object definition:


   :start ausgab object:                   #####  Select to create phase space, otherwise comment out  ######
        name = zbunker_head_leakage_phasespace_10X
        library = egs_phsp_scoring
        output format = IAEA
        phase space geometry = planar_phsp_box
        score particles on = entry and exit
        particle type = all
    :stop ausgab object:


    :start ausgab object:
        name = tracks
        library = egs_track_scoring
        score photons = on
        score electrons = on
        score positrons = on
        start scoring = 1
        stop scoring = 1000
        buffer size = 5000
    :stop ausgab object:



:stop ausgab object definition:



#################################################################
# Variance reduction - NB these are egs_chamber specific
#################################################################

:start variance reduction:            #### No variance reduction is used #####

    cs enhancement = 0                                    # turn on for primary beam doses
    TmpPhsp = 0                                                # turn on/off IPSS
    :start range rejection:                                    # turn on russian roulette
        rejection = 0                                               # 0 is off, or 1/128 is the survival probability
        Esave = 0.521 
        cavity geometry = iso_box                         # this is a dummy name
        rejection range medium = iso_water          # not doing range rejection so this doesn't mattter
    :stop range rejection:

:stop variance reduction:



############################################################################
### END OF FILE
############################################################################

[rtownson]

This is because a parallel beam uses the surface area of the source shape as a fluence normalization factor. E.g. in egs_parallel_beam, the getFluence function is:

    EGS_Float getFluence() const {
        return count/shape->area();
    };

And for a cylinder shape:

    EGS_Float area() const {
        return 2*M_PI*R*(R+h);
    };

This value of fluence is used to normalize the energy/dose deposited. However, the fluence is not what is used as the "number of primary particles" for the scored phase-space. Therefore the source normalization is lost when you go to use the phase-space. I'll have to give some thought to whether we should make a change. For now, you can recover the original normalization manually, now that you know how it is calculated.

[nexustowers]

Thank you Reid and yes, when I account for the surface area of the cylindrical source then the parallel beam / phase space source particle derived doses become consistent with each other.

I knew there are normalization considerations with collimated sources and phase space sources but didn't know about parallel beam sources. I also now see this question was previously addressed in question #1282 on Apr 5 2025.

Thank you again