Geometry error in CD geometry with conestack

[rtownson]

It appears that one of the changes we've made (it was probably me) to fix other geometry errors has introduced a new one for CD geometries, that arises in this very simple input file. It only seems to occur immediately with a phase-space source, for some reason. I haven't narrowed down the exact commit where this occurs yet, but it does disappear in v2017 and appear in v2018.

Edit: I have narrowed it down, see comment below.

Here is the original reddit thread where this was reported: https://www.reddit.com/r/EGSnrc/comments/kv1b0o/geometry_errors_in_egs_chamber_using_iaea_phase/

##################################################################################
##################################################################################
#
#   Chamber in air at 0,0,100
#   Varian 6X Truebeam (v2) phase Space source is at 26.7cm
#   Also uses EGS collimated source for comparison
#
##################################################################################
##################################################################################


#################################################################
# MC transport parameters
#################################################################

:start MC transport parameter:
    Global Ecut = 0.521
    Global Pcut = 0.01
:stop MC transport parameter:



#################################################################
# Random number generator seeds
#################################################################

:start rng definition:
    initial seeds = 91 2556
:stop rng definition:



###############################################################################
# Media definitions
###############################################################################

:start media definition:

    ae = 0.521
    ap = 0.01
    ue = 50.511
    up = 50

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

    :start polystyrene:
        density correction file = polystyrene
    :stop polystyrene:

    :start c552:
        density correction file = c-552air-equivalentplastic
    :stop c552:


:stop media definition:


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

:start source definition:

    :start source:
        name = phsps
        library = iaea_phsp_source
        iaea phase space file = /home/rwtownson/Downloads/Varian_TrueBeam6MV_02
        particle type = all
        cutout = -1.33 1.33 -1.33 1.33    #gives 10x10 @100cm
        recycle photons = 0
        recycle electrons = 0
    :stop source:


    :start source:
        name = photons
        library = egs_collimated_source
        charge = 0

        :start source shape:
            type = point
            position = 0 0 200
        :stop source shape:

        :start target shape:
            library = egs_rectangle
            rectangle = -5 -5 5 5
        :stop target shape:

        :start spectrum:
           type = monoenergetic
           energy = 6
        :stop spectrum:
    :stop source:


    :start source:                  # translate/rotate collimated source to mimic phasespace source position and beam direction
        name = trans_photons
        library = egs_transformed_source
        source name = photons
        :start transformation:
            translation = 0 0 200
            rotation vector = 0 0 -1
        :stop transformation:
    :stop source:


   #simulation source = trans_photons        # choose one or the other
   simulation source = phsps

:stop source definition:


############################################################################
# Build the geometry
############################################################################

:start geometry definition:

    # create chamber tip as sphere #
    :start geometry:
        name = chamber_tip
        library = egs_spheres
        midpoint = 0.0 0.0 100.0
        radii = 0.05 0.3 0.34                    # units cm
        :start media input:
            media = c552 dryair
            set medium = 0 0
            set medium = 1 1
            set medium = 2 0
        :stop media input:
    :stop geometry:


    # create chamber body as conestack
    :start geometry:
        name = chamber_body
        library = egs_cones
        type = EGS_ConeStack
        axis = 0.5 0.0 100.0 -1 0 0           # ie x-axis, first module starts at 0.5, note extra 5mm to provide overlap

        # unsheathed electrode
        :start layer:
            thickness = 1.1
            top radii =    0.05 0.3 0.34
            bottom radii = 0.05 0.3 0.34
            media = c552 dryair c552
        :stop layer:

        # sheathed electrode
        :start layer:
            thickness = 0.23
            top radii =    0.06 0.3 0.34
            bottom radii = 0.06 0.3 0.34
            media = c552 dryair c552
        :stop layer:

        # chamber body start
        :start layer:
            thickness = 0.5
            top radii =    0.34
            bottom radii = 0.34
            media = polystyrene
        :stop layer:

        # next chamber body piece
        :start layer:
            thickness = 1.3
            top radii = 0.34
            bottom radii = 0.34
            media = polystyrene
        :stop layer:

    :stop geometry:


    # define chopping planes for final chamber composite geometry
    :start geometry:
        name = cd_planes
        library = egs_planes
        type = EGS_Xplanes
        positions = -4.0 0.0 4.0
    :stop geometry:


    # create the composite geometry chamber from tip and body
    :start geometry:
        name = chamber
        library = egs_cdgeometry
        base geometry = cd_planes
        set geometry = 0 chamber_body
        set geometry = 1 chamber_tip
    :stop geometry:


    # create air column
    :start geometry:
        name = air_column
        library = egs_box
        box size = 50 50 400
        :start media input:
            media = dryair
        :stop media input:
    :stop geometry:


    # inscribe chamber into air column
    :start geometry:
        name = chamber_in_air
        library = egs_genvelope
        base geometry = air_column
        inscribed geometries = chamber
    :stop geometry:


    simulation geometry = chamber_in_air


:stop geometry definition:



################################################################################
# Run Control
################################################################################

:start run control:

    ncase = 1e6
    nbatch = 10
    geometry error limit = 20

:stop run control:



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

:start variance reduction:

#   photon splitting = 100                          # commented out, not used
    cs enhancement = 0                              # turned off
#   :start range rejection:                         # commented out, not used
        rejection = 256
        Esave = 0.521
        cavity geometry = air_column
        rejection range medium = dryair
    :stop range rejection:

:stop variance reduction:



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

:start scoring options:

    :start calculation geometry:
        geometry name = chamber_in_air
        cavity regions = 2 5 14                        # chamber cavity regions
        cavity mass = 0.1
        cavity geometry = chamber                  # for RR, not used
        enhance regions = 2 5                          # for xcse, not used
        enhancement = -128                             # use -ve to set all regions, not used

    :stop calculation geometry:

:stop scoring options:



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

:start ausgab object definition:

    :start ausgab object:
        name = tracks
        library = egs_track_scoring
        score photons = yes
        score electrons = yes
        score positrons = yes
        start scoring = 0
        stop scoring = 5000
        buffer size = 20000
    :stop ausgab object:

:stop ausgab object definition:

[ojalaj]

What kind of error (message) you actually get?

[ojalaj]

Ah, OK, I found it on Reddit.

[rtownson]

Apparently adding edep_local to struct EGS_Epcont in egs_interface2.h was the cause of this. Commenting out one or more of the EGS_Epcont parameters in egs_interface2.h makes the geometry errors go away! @mainegra @ftessier any idea what's going on here?? It seems that the number of parameters in EGS_Epcont plays some role, but the parameters from the interface match what you find in egsnrc.macros. Is there some step in the interface process that has been missed?

The troublesome commit: https://github.com/nrc-cnrc/EGSnrc/commit/85d4b7d6a6344400a4cb5f6e74df255a9ebda2d6

[ftessier]

Does that change not require a modification on the mortran back end as well, so the structs remain aligned? Oh wait, it is there in egsnrc.macros.

[ftessier]

Woah, there is an unrelated bug with the order of the variables Epcont->vstep and Epcont->tvstep, compared to the Mortran side! Fortunately, these variables are never used on the egs++, but I will fix that in passing.

[rtownson]

I saw that also!

[ftessier]

I can't figure it out. The structs seems correct (modulo the inversion of vstep and tvstep). So my only possible conclusion at this point in time is that the code is correct, and by making it correct it is revealing an existing bug elsewhere? That's a stretch I know. Other things to consider:

• are you testing from a fresh clone?
• is it necessary to reconfigure the code to get the interface linked on the Mortran side?
• are there compiler issues or flags that cause a misalignment of the structs between g++ and fortran? unlikely.
• are both fortran and egs++ compiled in double presision?

[rtownson]

It was not a fresh clone in my case, but it was for the users who found the bug. I don't believe it's necessary to reconfigure, all you need to do is recompile the application for the interface to be updated. I'm using the default compiler flags except with C++11, but the other users were using the defaults only.

I think you could be right that this part of the code is fine and just highlighting a different bug.. I'll have to try going back to yet older versions or (sob) actually stepping through the geometry error to see what is happening.