How to implement a separate energy (or dose) scoring grid?

[Jlambeck]

I have a geometry that does not fit well as an XYZ geometry (lot of curved pieces using ConeStacks) but would like to get a 3D dose map based on geometric location, not just region number. Essentially, I would like to create a 3d grid-like "overlay structure" that records energy or dose deposited in each small cube of the grid, no matter what region of the overall geometry it was deposited in. Is there a way to do this?

Can I create an XYZGeometry in the location of my main geometry but not give it a media then use a 3ddose ausgab scoring object?

Or make the XYZ over the main geometry with a random media (air/vacuum/etc) with a Union prioritizing my main geometry? Will that keep track of the XYZ region number where energy is deposited or does the Union mean the code 'forgets' the XYZGeometry is there if it's overwritten by another geometry?

The simulations are being done in egs_chamber with egs_2021 if that matters. (I cannot change/update the version as the copy on our computer cluster is controlled by someone else)

[ftessier]

This is an interesting problem that hinges on the mechanics of the electron transport algorithm. While it might seem straightforward to modify the source code and use a separate grid geometry to record energy deposition events in grid cells as they occur, this approach would generally yield incorrect results.

For example, consider a simple transport region the size of two cubes side by side, and a separate scoring grid composed of those two cubes. When an electron is transported through the actual region, the algorithm only provides the total energy deposited within the entire region. There's no way to determine how that energy is partitioned between the two sides of the region—unless, that is, we transport the electron in two separate, real cubes! And we are back to square (or cube) one.

Even with a region large enough to contain a few electron steps, there would still be a bias due to the unknown curvilinear path between the endpoints of the straight electron displacement step.

Note the Geant4 implements this "overlay" idea with Parallel Geometries, see also doi:10.1109/NSSMIC.2008.4774535

Here are three potential avenues that come to mind for EGSnrc, ordered by implementation difficulty:

1. Create an XYZ geometry grid

• Replicate your original geometry within each scoring grid cell
• This creates n×N regions (n = original regions, N = grid cells)
• Energy deposition for each grid cell must be summed across its n regions
• Performance impact may be minimal if physical electron step size limit already commensurate with cell size

2. Use an octree geometry

• Currently needs technical fixes but conceptually sound
• Automatically subdivides space into variable-sized boxes
• Creates larger regions where medium is constant
• Adds refinement where medium changes occur
• May lose some geometric detail (pixelation of geometry) but suitable for heatmap visualization, for example.

3. Modify source code to handle dual geometries

• Similar to Parallel Geometries in Geant4.
• Add howfar() and hownear() calls to check both geometries before each particle step
• Consider the nearest boundary to select the transport geometry
• Score energy deposition in both geometries (based on isWhere() location in the non-selected geometry)
• Most flexible solution, and in fact relatively straightforward to implement now that I am thinking about it!!

[mainegra]

This is already possible using an EGS_XYZGeometry geometry!

There is an undocumented option for passing a geometry to an EGS_XYZGeometry to be voxelized.

Here is an example on how to do this for a geometry named seed:

:start geometry:
    library = egs_ndgeometry
    type    = EGS_XYZGeometry
    name    = seed_voxelized
    x-slabs = -0.3 0.00333 200
    y-slabs = -0.3 0.00333 200
    z-slabs = -0.3 0.00333 200
    voxelize geometry = seed
    #delete geometry = yes  # frees memory used by original geometry
:stop geometry:

The voxel midpoints of geometry seed_voxelized are mapped onto the original geometry seed using the geometry method isWhere(), and each voxel's medium is assigned based on the region where its midpoint falls.

The main challenge arises when objects within the geometry vary significantly in size. Small objects may be missed entirely, and finer details could be lost. Adjusting the grid resolution is necessary to mitigate these issues.

For instance if one voxelizes the geometry hdr1000_with_seed using 1.2 mm voxels, one ends up losing the seed:

    :start geometry:
        library = egs_ndgeometry
        type    = EGS_XYZGeometry
        name    = hdr1000_voxelized
        x-slabs = -12. 0.12 200
        y-slabs = -12. 0.12 200
        z-slabs = -12. 0.12 200
        voxelize geometry = hdr1000_with_seed
        #delete geometry = yes
    :stop geometry:

Original:

Voxelized (misses the seed):

In those cases, one could voxelize the geometry into different geometries and then combine them.

[mchamberland]

Whoa! Nice feature!

[Jlambeck]

This looks like exactly what I needed! Thank you!