EGSnrc
EGSnrc copied to clipboard
nrc-cnrc
Inconsistent Results between EGSnrc v4 (rev. 2.2.3) and EGSnrc v2021
Discussed in https://github.com/nrc-cnrc/EGSnrc/discussions/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
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!
@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.
@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?
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.
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.
@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.
@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...
No, @agsponer-metas: I haven't yet installed the EGSnrc_V4 (or equivalent) on my system to see why the slab seems to be placed at the wrong position.
Any updates on this in the new year? Please let me know if I can help you with any additional information.
@agsponer-metas, I'm unable to reproduce this on the oldest release available on the EGSnrc wiki (EGSnrc 2015), so I'll track down EGSnrcV4, install on another machine, and see if I can reproduce it there. Stay tuned, then.
@ftessier, I was able to find a dropbox link you'd made to EGSnrcV4 with self-extracting installers, but do you know where I might find the standalone tarred/compressed directories for this distribution together with an installation script? I'd like to try installing on Mac, and, as I remember, this was the only way to do it.
Hi @agsponer-metas: I ran your cool_cap_test with EGSnrcV4 and still wasn't able to reproduce the missing airgap above the slab. I ran with BEAMnrc rev 1.78 (same as you) and EGSnrc rev 1.16. Although the EGSnrc rev is possibly not the same as yours, I'm not concerned because the geometry should all be taken care of in BEAMnrc in any case.
Can you send me any .egslog files from your EGSnrcV4 runs? Also, if you could zip up the contents of the $HEN_HOUSE/omega/beamnrc directory (including subdirectories) and somehow attach it, that would be helpful as well.
Hi @blakewalters Thanks for looking into the issue, I have attached the simulation files of a fresh run (BEAM_cool_cap_test.zip) as well as the beamnrc directory (beamnrc.zip).
Hi @agsponer-metas: Thanks for sending the files! They certainly helped sort out what I think is the source of the discrepancy between old and new results: In some of the CMs, SLABS among them, the Z-distance between the start of the CM geometry (i.e. Z-position of the first slab) and the end of the previous CM is tested against a minimum air gap thickness, $MIN_GAP. If this distance is found to be < $MIN_GAP, then the Z-position of the first slab is moved up to coincide with the end of the previous CM. So, you eliminate the air gap AND make that first slab (the only slab in your cool_cap_test case) thicker by $MIN_GAP. Now, $MIN_GAP is a macro defined in $OMEGA_HOME/beamnrc/beamnrc_user_macros.mortran and has been set to 0.01 cm since 2006 or earlier. The difference is that, at a more recent point in EGStory, we moved from single precision to double precision variables, including for those defining CM geometries. As a result, in EGSnrc from 2004-2006, roundoff error results in (Z of first slab) - (Z of back of previous CM) being < $MIN_GAP, while in the current version of the system this does not happen.
In summary, then, in the current version of the system, this Z-distance test is less likely (i.e. likely will never) fail, the geometry simulated is identical to that input and results are more accurate. Just to make sure that test never fails, though, I would suggest going into $OMEGA_HOME/beamnrc/beamnrc_user_macros.mortran and setting $MIN_GAP to a smaller number, say 1e-5, so this no longer becomes a concern.
I'm also going to suggest that we set $MIN_GAP to 1e-5 in the distribution. Given that we're using double precision now, there's no point in potentially introducing unnecessary inaccuracies by not modeling the actual Z-dimension of the first layer in a CM!