Inconsistent Results between EGSnrc v4 (rev. 2.2.3) and EGSnrc v2021

[ftessier]

Discussed in #772

^{Originally posted by agsponer-metas September 29, 2021}

I'm new to EGSnrc and have been recently tasked with "resurrecting" an old BEAMnrc simulation of a treatment head attached to an electron accelerator. At the time the old simulations were performed (~2003-2006), a cluster of office type machines was used to run the simulations. In order to retire this old cluster and because we want to run new simulations which modify the existing geometry of the treatment head, it was decided to migrate the simulations to new hardware and the new version of EGSnrc (v2021). The old simulations were performed using EGSnrc v4 rev. 2.2.3. (which I will call v4 for brevity) and BEAMnrc Rev. 1.78.

To verify that the results obtained by EGSnrc v2021 are unchangend compared to the old results (EGSnrc v4), I ran the simulations using EGSnrc v2021 and the same code that was used for the older runs (same .egsinp files, BEAM accelerator modules, PEGS4 files). For electron beams this worked fine and I was able to reproduce the exact same results as EGSnrc v4 (inside the statistical uncertainty).

However, for photon beams (which use a gold conversion target) a statistically significant difference in the energy spectrum of primary and secondary electrons (which reach the scoring plane after the first part of the treatment head) was observed, see EGSnrc_v4(2.2.3)_vs_v2021_for_entire_treatment_head.pdf in the attachments. In the histogram, the particle count has been normalized to the incident electrons originating from the accelerator "source" and to the area of the scoring plane (33 x 33 cm^2).
The rest energy of the electron has been subtracted for the electrons/positrons and the error bars are assumed to be sqrt(N), where N is the amount of particles in a bin. In the v2021 version of EGSnrc, around 20% more primary electrons reach the scoring plane. This is significant for us, for example when we want to characterize the electron contamination in the photon beam.

In order to try and narrow down if this difference stems from a single component module or is the result from interactions along the entire length of the treatment head, I ran simulations using subsets of the treatment head geometry component modules. If I just consider the first two component modules (the target in its holder targ_hol and the cone behind it tar_cone), the results from EGSnrc v4 and v2021 seem to match. However, if the next CM (a copper cooling SLAB cool_cap) is added the results start to diverge (see "EGSnrc_v4(2.2.3)_vs_v2021.pdf").

For this situation, I provided a minimum reproducible example (BEAM_cool_cap_test in the attachments). Looking at the .egslst output files (also attached) for this third component module, I see a difference in the geometry:

v4 (rev 2.2.3) ("cool_cap_test_w1_EGSnrc_v4.egslst"):
 ........
 coolcap geometry parameters:
 -----------------------------
 Distance of front of CM from reference plane =        2.66000 cm
 Half-width of outer boundary of CM =         2.00000 cm

 slab #    Z front    thickness
            face               
            (cm)        (cm)   
    1       2.660      0.050
 ........

v2021 ("cool_cap_test_w1_EGSnrc_v2021.egslst"):
 ........
 coolcap geometry parameters:
 -----------------------------
 Distance of front of CM from reference plane =        2.66000 cm
 Half-width of outer boundary of CM =         2.00000 cm

 slab #    Z front    thickness
            face               
            (cm)        (cm)   
 airgap     2.660      0.010
    1       2.670      0.050
 ........

cool_cap_test.egsinp:
........
*********** start of CM SLABS with identifier cool_cap  ***********
2.0, RMAX
Target holder cooling circuit cap
1, NSLABS
2.67, ZMIN
0.05, 0, 0, 0, 2, 0
CU
........

For some reason, EGSnrc v4 inserts the cool_cap SLAB at a position of 2.66 cm, even though Z_MIN is specified as 2.67 in the .egsinp file (see above). This means that the SLAB is flush with the CONS3R module in front of it, which extends to 2.66 cm. EGSnrc v2021, however, correctly inserts an air gap of 0.01 cm and puts the SLAB at a position of 2.67 cm.

The interesting thing is now, if I increase the thickness of the copper cool_cap SLAB by 0.01 cm in EGSnrc v2021 (see below and the energy spectrum in EGSnrc_v4(2.2.3)_vs_EGSnrc_v2021+0.01_cm_Cu.pdf), I can reproduce the same energy spectrum as EGSnrc v4!

cool_cap_test_0.01cm_additionnal_SLAB.egsinp:
........
*********** start of CM SLABS with identifier cool_cap  ***********
2.0, RMAX
Target holder cooling circuit cap
1, NSLABS
2.67, ZMIN
0.06, 0, 0, 0, 2, 0
CU
........

Does this mean that EGSnrc v4 (rev. 2.2.3) erroneously inserts 0.01 cm of additional material? Is this a bug in EGSnrc v4 which has been fixed at some point between v4 and v2021?

Adding the 0.01 cm of copper does not fix the issue for the simulation of the entire treatment head, as there seems to be a similar error occurring in a second CM later on (which I could not pinpoint yet). I will have to go from component module to component module to find out which of them introduces further errors. Compared to the cool_cap SLAB, however, I have not been able to find any further discrepancies in the .egslst log files.

As I have around 15 simulations (with different electron beam energies and flattening filters) to go through, I would be glad for any aid in explaining these differences between EGSnrc v4 (rev. 2.2.3) and v2021. Was there maybe some change to the geometry generation algorithm? Would this imply that the results using the old version of EGSnrc are incorrect?

I have not yet played around with the Monte Carlo transport parameters (where some default values between the EGSnrc versions might have changed which I do not explicitly define in the .egsinp file), but from the electron beam simulations (without the conversion target) I did not yet find any transport parameters which lead to significant differences.

I hope that the information I provided is complete and clear enough to describe my issue.

Thanks in advance for your reply,
Andreas

attachments.zip

[ftessier]

I see only cosmetic changes in the SLAB component module's history on git, back to 2015. @agsponer-metas if you can run a test with EGSnrc v2016 quickly, that would prove useful. We have CVS history back to V4 2.2.3, but it is harder to untangle that... So knowing if the change happened before or after v2016 (when EGSnrc was ported to git) would help me steer my efforts in the right direction!

[agsponer-metas]

@ftessier, I tried using the v2016 tag, but could not get this version of EGSnrc to compile... (see compile_error_v2016.txt ). However, I was able to compile and run the simulation using the v2015 version (is there a significant difference between v2015 and v2016?).

With v2015 I still get the same result as v2021 (an excess of electrons, see plot). Similar to before, an airgap of 0.01 cm is inserted before the SLAB (see .egslst).

Therefore, I assume the relevant changes in the code happened some time before EGSnrc was ported to git / 2015.

[ftessier]

@agsponer-metas, thank you for checking, this is helpful, I will dig up older changes to the SLABS component module, unless @blakewalters can recall any such change from memory alone. If I understand correctly, the SLAB CM is now behaving as expected (i.e., with the extra air added), right?

[agsponer-metas]

Hi @ftessier
Thanks for looking into this issue! Yes, the SLAB CM is now (v2021) behaving as expected. What I'm trying to figure out is if this means that the old results obtained using EGSnrc v4 rev 2.2.3 are incorrect, or if the differences could be caused by some other parameters which have changed between EGS versions.

[agsponer-metas]

On a related note, I found the second CM in the simulation of the treatment head where I get differences between the two EGSnrc versions. As mentioned previously, even after adding the 0.01 cm of material to the SLAB mentioned above the results for the simulation using all CMs do not agree between the EGS versions.

This second CM is the flattening filter, which is constructed using FLATFILT. I ran a simulation using just this CM (and all other CMs removed, see filter_only_treathead_phot_8_5mev.egsinp.txt). In the resulting energy spectrum (total_energy.pdf) I see a very large difference for the primary electrons, but identical results for all other particle types. What intrigues me is that the peak around 8 MeV is identical for both EGSnrc versions.To me, if there was some additional absorption (which leads to EGSnrc v4 having fewer primary electrons at lower energies as in the plot), I would expect that this would also have an effect on the maximum electron energy.

In contrast to the SLAB module, I can not find any differences in the geometry (.egslst) generated by the two EGSnrc versions (see EGSnrc_v4 .egslt and EGSnrc_v2021 .egslst.

I will try splitting the FLATFILT in half and run the simulations again to see if I can identifiy a single layer which is responsible for this.

[ftessier]

@agsponer-metas thank you for the detailed investigation, this proves extremely valuable to track down this issue. I will first diff through your .egslst files; some default values for the transport parameters have changed over the years.

[agsponer-metas]

@ftessier Any updates on this?

[ftessier]

No yet from my end.

[blakewalters]

@ftessier and @agsponer-metas: I've only had time to follow this from a distance up til now, but should have time to look at it in more detail this week! Just off the bat, I can't think of any changes to SLABS that would explain the (lack of) insertion of the 0.01 cm airgap between EGSnrc V4 and the 2021 release. In fact, I can't think of any substantive changes to that piece of coding at all since it was first released in 1995! Perhaps BEAMnrc has been changed somehow, but, even there, I can't think of one that was introduced that would be relevant to this issue. In any case, I'll take a closer look this week...