Add dbs to radiative splitting

[blakewalters]

This is a draft PR for the DBS option in the EgsRadiativeSplitting ausgab object and will be migrated to a regular PR once testing is complete, some of the more obvious bugs are ironed out, and it's been rebased on the current develop branch.

EGS developers are encouraged to download and try out the branch and provide comments and corrections as needed. Note also that posting of test results in the comment feed below is encouraged.

[blakewalters]

Plots showing results of some initial tests in which I compared spectra using DBS in egs++, BEAMnrc and beampp in 2 different applications:

1. EX16MVp

2. A 450 kVp X-ray tube (spectra scored in vacuum)

i) With bound Compton (i.e. not using SmartCompton in DBS):

ii) Without bound Compton (i.e. using SmartCompton in DBS):

From 2.ii, it's clear that using SmartCompton in the egs++ DBS algorithm has introduced some bias. Currently, the algorithm uses the same SmartCompton implementation as beampp. I plan to replace this implementation with the one found in BEAMnrc and see if the difference persists.

[ftessier]

This is awesome @blakewalters. I move to merge this in develop soon, so that we have a lead testing time in the field before the next release! Let's review this on the short term!

[blakewalters]

Update 13/05/2023: After implementing BEAMnrc's SmartCompton algorithm (translated into C++) in the egs_radiative_splitting DBS routine, agreement between BEAMnrc and egs++ in the 450 kVp tube case shown above is inching ever closer:

In fact, at this point we might say there's no significant difference.

Continued testing is needed, though, given that, during one of these 450 kVp validation runs, I noted a phat photon crossing the scoring plane with r < DBS splitting radius.

[blakewalters]

Update 07/07/2023:

After finding/fixing the bug in the SmartCompton method which was causing the odd high weight photon to be scored inside the splitting field (see update from 13/05/2023 above), 450 kVp X-ray photon spectra calculated using the DBS routine implemented in egs_radiative_splitting both with and without SmartCompton (i.e. using the standard EGSnrc Compton routine) look to be in good agreement as shown in the figure below.

This is borne out in a plot of the fractional difference between egs++ and BEAMnrc spectra. However, this plot also highlights significant differences (>5% at the extreme ends of the spectrum) between BEAMnrc and egs++ results (Note: In a previous edit of this comment, I'd erroneously labeled the vertical axis (BEAMnrc-egs++)/egs++ when, in fact, fractional differences are given relative to BEAMnrc results).

At this point, I can't say much about the source of these differences except that they're most likely not coming from egs++ DBS's handling of Compton events!

I think it might be worth repeating the egs++ simulations with UBS just to confirm that the difference isn't due to any slight differences in geometry between BEAMnrc and egs++.

[ftessier]

Wonderful @blakewalters!

[blakewalters]

Update 09/22/2023:

Implementation of beampp's SmartCompton algorithm in egs++ (as opposed to the BEAMnrc SmartCompton algorithm from the previous update) results in similar differences with BEAMnrc:

Erratum: The first edit of this post from 09/15/2023 showed agreement between BEAMnrc and egs++ with the regular Compton routine used in place of Smart Compton. This was incorrect and, moreover, contradicts the results above from 07/07/2023.

A comparison of photon fluence vs radial position for BEAMnrc (SmartCompton) vs egs++ with both SmartCompton and regular Compton shows fairly close agreement:

While photon energy fluence is visibly significantly lower for the egs++ cases:

Indicating that, perhaps, the issue is actually "upstream" of Compton interactions and may be in the energy distribution of brems photons...

[mainegra]

@blakewalters great job! Is the above the current status?

[blakewalters]

That's where it's at, @mainegra. My next debugging step will be to see if there's something going on with the energies of the split brem photons. These energies are determined in a separate function after the interaction has been split and the photons generated.

[blakewalters]

Update 31/05/2024:

There've been a couple of major fixes in the DBS splitting routine over the past couple of months. Most significant of these stems from the realization that nonphat photons about to undergo Compton, photoelectric, pair or Rayleigh events should be subject to Russian Roulette--and made phat if they survive--prior to interacting. This is a key element of DBS, and I can't believe I didn't notice this omission in my egs++ implementation before! In addition, however, there was a bug in the function, doUphi21, determining angles for Rayleigh scattered particles as well as a C++-->Fortran indexing issue in the getBremsEnergies() function calculating energies for brem photons generated using doSmartBrems, among other bugs.

I've also since moved to doing BEAMnrc/egs++ comparisons using transmission photon spectra instead of reflection spectra to reduce the number of potential discrepancies between BEAMnrc-simulated and egs++-simulated source/target geometries.

Here's the result of a comparison of 450 kVp W spectra with relaxation effects turned on. In this case, bound Compton effects are also simulated and, thus, smartCompton is not used:

Spectra in 10 cm vaccum (within r=10 cm splitting field):

Fractional difference between BEAMnrc and egs++:

In this case, I'm pretty confident of the match between BEAMnrc and egs++ spectra, although this could be verified to < 0.5% by running more histories.

Below is the result for the same case with relation effects off. Here, bound Compton is turned off and, by default, smartCompton routines are used in both BEAMnrc and egs++:

Spectra in 10 cm vacuum (within r=10 cm splitting field):

Fractional Difference between BEAMnrc and egs++

Although we might say the two spectra agree within uncertainty, there is a drift in the egs++ spectrum from lower to higher values with increasing energy (see fractional differences) that's a little troubling. This needs to be verified by running more histories, but if it persists, then it is almost certainly due to issues in the egs++ smartCompton method. Currently, this method uses the beampp implementation of smartCompton, which differs slightly from the BEAMnrc implementation.

Oh yes, for those interested in helping out with the testing/debugging of egs++ DBS, a good place to start might be the BEAMnrc (EXslabs) and tutor7pp input files used in this latest round of tests. Examples of these are attached here: 450kVp_BEAMnrc_split_test_relax_off.txt 450kVp_tutor7pp_split_test_relax_off.txt

[ftessier]

📌 Run astyle for code styling, and remove end-of-line whitespace.

[ftessier]

Awesome, thanks @blakewalters!

[blakewalters]

Continuing with the tests from the last update, here's the fractional difference between BEAMnrc and egs++ photon spectra with relaxation effects off when smartCompton is not used in egs++ DBS:

which shows good agreement and seems to confirm that the bias seen in the previous comparison with relaxation off is due to the current smartCompton implementation in egs++ DBS.

As mentioned above, the smartCompton implementation in egs++ is currently based on the one in beampp, and so there are a couple of possible next steps:

1. Comb through the egs++ smartCompton implementation for bugs
2. Implement smartCompton following BEAMnrc's logic (I've tried this before with previous iterations of egs++ DBS. However, with recent fixes to the algorithm, it bears repeating...again!)

[blakewalters]

After (re)implementing the BEAMnrc smartCompton algorithm in egs++ DBS (and noting some errors I'd made in my original implementation of this along the way), we get good agreement between BEAMnrc and egs++ spectra with relaxation effects off:

We have to decide, then, if we're going to go with our original vision of having more than one flavour of egs++ DBS (one modeled after BEAMnrc, one modeled after beampp, another modeled after?). If so, then we'll have to revisit and debug the beampp smartCompton algorithm. For now, though, I think egs++ DBS is ready for some more intense testing!

[blakewalters]

Update 07/04/2024:

Results of recent tests comparing BEAMnrc/DBS and egs++/DBS in the case of a 225 kVp e- beam incident on a reflection anode (1 mm W, Cu backing, anode angle = 4.92 deg.). Note that I've used a pencil beam of e- (parallel beam incident at a point) to try to mitigate any differences due to source geometry. For the BEAMnrc simulation, this meant using source 13 (e- beam pointed in the -X direction at the face of the anode) with rectangular dimensions (0.0002 x 0.0002 cm, i.e. really small).

Here are the photon spectra at SSD=10 (splitting SSD) and r=10 cm (splitting radius). The DBS splitting number in all cases is 10,000. I've also included results for egs++ with UBS (splitting no. = 100) and egs++/DBS with smartBrems turned off:

With the exception of the peak at ~60 kV, where egs++/DBS looks a little low, a visual comparison indicates good agreement. The fractional difference between BEAMnrc and egs++ results, however, tells a slightly different story:

Here we see that egs++/DBS is up to 2.5% greater than BEAMnrc at low energies, egs++/DBS with smartBrems is up to ~1% high at energies > 50 kV. egs++/UBS and egs++/DBS without smartBrems are up to ~0.5% high for energies > 50 kV, although better stats are required to determine the significance of this difference. In the case of egs++/DBS with smartBrems, though, it seems clear that some bias is introduced.

Next steps:

1. Simulate to better uncertainty to confirm differences
2. Comb through the smartBrems routine once again to check for possible differences between egs++ and BEAMnrc implementation and, possibly, reimplement it in egs++ to match that in BEAMnrc

[blakewalters]

Update 08/02/2024:

After running the 225 kVp reflection anode results down to ~0.1% uncertainty, a plot of the difference between BEAMnrc/DBS and egs++/DBS photon fluence spectra confirms the differences of ~1% seen in the last update. The "good" news is that, after running more histories in the egs++/UBS (100x splitting) case, the difference between BEAMnrc and this case is confirmed as well:

Moreover, egs++/DBS and egs++/UBS results are within uncertainty, leading one to suspect that differences between egs++ and BEAMnrc are due to differences in modeling the X-ray tube geometry (I'm assuming, of course, that my double checking that all MC transport parameters are identical doesn't require triple checking!).

A plot of the difference in photon fluence along the anode-cathode axis (incident e- in the -X direction) might indicate a slightly smaller anode angle in the BEAMnrc simulation:

On the other hand, this may also indicate slight differences in the effective thickness of the W target layer.

More investigation is required to track potential geometric differences down. However, it now seems they're adding an unnecessary layer of complication in comparisons between BEAMnrc/DBS and egs++/DBS when a reflection anode is used. At this point, then, I suggest going ahead with a comparison at 225 kVp using a transmission anode (i.e a simple W slab).