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.

[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).

[blakewalters]

Update Oct. 2, 2024

Further simulations of the 225 kVp reflection anode (W thickness=0.1 mm, angle=4.92 degrees) from the previous update have confirmed consistent differences on the order of 1% between BEAMnrc and egs++ photon spectra, regardless of whether egs++ is used with UBS or DBS. This implies some slight difference between the two codes in the way the anode geometry is defined (maybe traceable to differences in boundary tolerance used?) and indicates that, at least for now, if we want to compare reflection spectra, the way forward is to compare egs++_DBS against egs++_UBS. The assumption that egs++_UBS is correct seems fair given that photon splitting is handled inside EGSnrc in this case.

The plots below show the fractional difference between photon spectra for the reflection anode generated using egs++_DBS and egs++_UBS. Note that, instead of analyzing phase space data as I did previously, I've scored spectra using the egs_fluence ausgab object. This facilitates running more histories to improve the statistics of the comparison and eliminates the requirement for storing and analyzing phase space data. The first plot shows the fractional difference with atomic relaxations turned off (i.e., smartCompton used), while the second includes atomic relaxations (smartCompton not used) and electron impact ionization:

With relaxation off, spectra agree to within ~0.05%, while with atomic relaxation on there is a significant difference of 0.2% at the photon peak at ~62 keV.

For reference, the plot below shows photon spectra with relaxation on with and without electron impact ionization (eii). Note that, while eii has a significant effect on the spectrum, any differences between egs++_UBS and egs++_DBS are too small to see at this scale:

I then turned eii off and, once again, egs++_DBS agrees with egs++_UBS to within ~0.05%:

It seems, then, that eii results in additional fluorescent photons that aren't handled correctly in egs++_DBS. Alternatively, it may magnify the effects of a bug in the handling of these photons whether eii is on or not.

The adventure continues! Any and all ideas about what might be going on with eii are welcome.

[rtownson]

@blakewalters In the EII sampling routine it copies the particle properties from the parent electron like this (line 6174 in egsnrc.mortran):

$TRANSFER PROPERTIES TO (np) FROM (np_save);

And this is copied to relaxation secondaries. For the fluorescent photon should the weight be adjusted? I'm not quite sure the solution but certainly this is the spot to look. Probably at the latch and what it should be for fluorescence v.s. auger electrons.

[blakewalters]

Thanks, @rtownson! I'll start there.

[blakewalters]

Update Oct 9, 2024

Using the egs_fluence ausgab object together with egs_beam_source, I've been able to use egs++ to calculate BEAMnrc photon spectra down to ~0.05% uncertainty. As a result, we can more precisely compare BEAMnrc and egs++ results.

The plot below shows the fractional difference between 225 kVp reflection spectra with relaxation effects on and eii on where now egs++ spectra with DBS (egs++_DBS) and UBS (egs++_UBS) are compared against the BEAMnrc spectrum with DBS (ssd=10 cm, r=10 cm, nsplit=10000).

The plot clearly shows the ~1% difference between egs++ and BEAMnrc consistent with the previous update and which we surmise is due to slight differences in anode geometry.

Note, however, that egs++_UBS is ~0.2% higher than BEAMnrc at the peak at ~62 keV (this shows up as a ~0.2% dip in fractional difference relative to BEAMnrc). This is similar to the fractional difference between egs++_UBS and egs++_DBS at this energy, and we may cautiously advance that the treatment of fluorescent photons resulting from eii events in egs++_DBS is consistent with that in BEAMnrc with DBS (or BEAMnrc_DBS, if that's easier). This makes sense, given that egs++_DBS is heavily based on the DBS implementation in BEAMnrc.

The plot below shows a comparison of 225 kVp spectra using a transmission target, where no differences between BEAMnrc and egs++ geometries are expected:

In this case, egs++_DBS and BEAMnrc_DBS show agreement to < 0.05% with the exception of the lowest energy bins, where uncertainties are > 1%. The spectrum calculated using egs++_UBS, however, is ~0.15% higher at the ~62 keV peak (showing up as a ~0.15% dip in the above plot). This is consistent with the ~0.2% higher peak seen in the egs++_UBS spectrum from the reflection anode (first plot). We also note a ~0.05% difference between egs++_UBS and BEAMnrc at ~190 keV, although this difference is likely within statistical uncertainty.

We can hesitatingly conclude, then, that treatment of fluorescent photons resulting from eii is the same in egs++_DBS and BEAMnrc_DBS. The question now is: Why are they both ~0.2% low at the K-shell peak(s)?

[blakewalters]

@mainegra, @ftessier, @rtownson:

Continuing on with testing of the DBS implementation in EGSRadiativeSplitting, I've been mulling over a couple of questions:

1. How do we implement DBS so that it will work with transformed sources? Right now, the implementation assumes the source is at Z=0, with the SSD and splitting radius defined relative to this. It would be easy to add more inputs to the AO to allow the user to translate/rotate the splitting field, but would it be more convenient and useful to somehow be able to "attach" a DBS splitting field to a particular source? That way, when a source is transformed, the relevant DBS splitting parameters are automatically transformed along with it. I guess an option would be to implement an "egs++_simulation_source" (something we might want to do anyway), in which case the DBS AO would just be part of the source simulation that can then be transformed in any way the user wishes in the second-stage simulation. Come to think of it, that might be the best way to accomplish this. Thoughts?

2. What about BCSE? I'm assuming we ultimately want to implement BCSE along with DBS, but this will require hooking into MORTRAN macros inside of EGSnrc. It's currently fairly straightforward to access subroutines inside EGSnrc, but I don't know if macros can be modified, short of making changes to interface coding (e.g., egs_c_interface2.macros). Let me know if any of you know another way!

[mainegra]

• Adding enough info during input to allow for arbitrary splitting field definition would be my choice
• BCSE is implemented in egs_brachy, perhaps a good place to start?

[rtownson]

The DBS should not be tied to the source or the way it is transformed, it should be completely independent IMO.

Ideally, it should let you set filter where it takes place by geometry, include regions and include media like in egs_isotropic_source. Then use a generic shape (not just a circle?) with a transform to denote the field, like the target shape for egs_collimated_source. This might be tricky to check, so if a circular cone is easier that's fine with me!

[ftessier]

I agree with @mainegra and @rtownson, except that I feel there should be an option to transform the target shape with the source. That is, when dealing with an egs_transformed_source, one could opt to apply the transformation to the target shape as well. Seems simple enough to implement, a flag and corresponding conditional to multiply by T, the transformation. Perhaps this ought to be a generic source option, useful for egs_collimated_source as well, for example.

[rtownson]

I agree with @mainegra and @rtownson, except that I feel there should be an option to transform the target shape with the source. That is, when dealing with an egs_transformed_source, one could opt to apply the transformation to the target shape as well. Seems simple enough to implement, a flag and corresponding conditional to multiply by T, the transformation. Perhaps this ought to be a generic source option, useful for egs_collimated_source as well, for example.

I thought this was already essentially the case. Since egs_transformed_source transforms whatever particle you get out of the original source, which already had its direction towards the target. Thus its direction toward the target gets transformed to point toward a new target, as per the rotation in the transformation. Or are you imagining something else?

[ftessier]

Yes, of course you are right! What I had in mind is that the DBS target shape would be defined within the source, not within the ausgab object. Or perhaps one could define a completely decoupled DBS field in the ausgab object? I can think of scenarios where either makes sense.

[rtownson]

Ah I see. I was certainly imagining a completely decoupled ausgab object. It feels unintuitive to place inputs for an ausgab object inside of a source. What if you had a collection of sources? You'd have to just have one target, in the source collection, but a user might put it in each of their sub sources, etc..

I can also imagine wanting several DBS ausgab objects. What if you have a simulation of multiple x-ray tubes?

[blakewalters]

@ftessier and @rtownson, this is why I was thinking that a shared library source might be the best way to take care of this: DBS is defined as an ausgab object in the source simulation and then that source can be duplicated/transformed in the 2nd stage simulation without having to worry about defining multiple DBS AO's. Kind of like what already happens when we use BEAMnrc to simulate an X-ray tube and then use this as a source in egs++. Of course, there may be issues with optimizing DBS parameters across multiple transformed sources.

[ftessier]

Now that I think about it, would it makes sense to implement DBS as a source feature with a target shape, instead of an ausgab object? I know that DBS is implemented as a region-based VRT, so to speak, and does not fit in the source... But perhaps there is a way to adjust our thinking and make it a source feature? 🤔

[blakewalters]

@ftessier: After today's chat with @mainegra and @rtownson, we've arrived at the following steps for egs++ DBS implementation:

1. Within the AO, implement the option to transform (rotate/translate) the DBS splitting cone.

2. Also implement an option to read in the (active) source transform so that the the DBS splitting cone can also either: a) move with a dynamic source or b) be applied to the current source in a source collection. This transform would be applied on top of any transform input in the AO itself (i.e., step 1).

Priority, of course, is benchmarking the DBS routine for accuracy, and the hope is that all of us can undertake more extensive testing of the routine at some point between implementing step 1 and step 2.

There is also an issue with the UBS implementation in egs++ resulting in a ~0.15% discrepancy in the photon spectrum at the K-shell peak when electron impact ionization (EII) is on (see above). Although peripherally related to DBS, we would like to solve this as well.

@mainegra and @rtownson, please let me know if I'm missing anything essential from our discussion here.

[blakewalters]

Update 12/03/2024:

Testing with IBRDST=1

Continuing testing of the egs++ DBS implementation, I was interested in seeing its accuracy when the full Koch-Motz (KM) distribution is used for bremsstrahlung angular sampling (IBRDST=1). This is done by setting "brems angular sampling=KM" in the transport parameter input block. Note that comparisons above only the leading term of the KM distribution (IBRDST=0), or "brems angular sampling=simple" in the transport parameter inputs.

The plot below shows the fractional difference in the 225 kVp transmission photon spectra between the baseline case using the "simple" brems angular distribution and when the full KM expression is used. BEAMnrc was used to generate both spectra. Worth noting here is that in current egs++ DBS implementation (as in BEAMnrc DBS) smartBrems is not used when IBRDST=1. This differs from Iwan's beampp implementation of DBS, in which smartBrems does have conditions for handling IBRDST=1.

Use of the full KM distribution results in differences in the energy spectrum ranging from -0.2% at 120 keV to 0.8% at 225 keV. Given that egs++ with DBS agrees with BEAMnrc with DBS to within <~0.05%, then, potential issues arising when IBRDST=1 is used with the egs++ DBS routine should be visible in a comparison with BEAMnrc (also using IBRDS=1).

The plot below shows the fractional difference in transmission photon spectra generated using egs++ with DBS and BEAMnrc with DBS where the full KM distribution is used for brems angular sampling in both simulations (black). The plot also shows the fractional difference when egs++ with UBS is used (red line).

This shows that, even using the full KM distribution, egs++ with DBS agrees with BEAMnrc with DBS to within ~0.05%. We also note that the ~0.2% difference at the K-shell peak (~60 keV) between BEAMnrc and egs++ with UBS when electron impact ionization (EII) is on persists with IBRDST=1.

A side issue: Comparison of egs++ with UBS and BEAMnrc with UBS against BEAMnrc_DBS:

In making these IBRDST=1 comparisons and attempting to locate the source of the discrepancy between egs++_UBS and BEAMnrc_DBS (and egs++_DBS) at ~60 keV, I also introduced into the comparison runs with BEAMnrc using UBS (nsplit=100) and found some "interesting" results overall.

Both figures below show the same plot of fractional difference between the photon spectra generated with egs++_UBS and BEAMnrc_UBS and our baseline BEAMnrc_DBS case at two different scales. Simulations with UBS were run with charged particle Russian Roulette on and off. In all cases, I reverted to the simple brems angular distribution (IBRDST=0).

The full-scale plot (upper figure) shows significant differences between BEAMnrc with UBS and Russian Roulette on (green line) and BEAMnrc with DBS at energies < 50 keV, with differences of ~40% at energies <= 20 MeV. These large differences disappear when Russian Roulette is turned off (magenta line). The close-up (lower figure) shows that the ~0.2% difference between egs++ with UBS and BEAMnrc with DBS at the K-shell persists when Russian Roulette is turned off (red and black lines). Moreover, we don't see this difference with BEAMnrc_UBS when Russian Roulette is off (magenta line), a further indication of an issue with the UBS implementation in egs++ (with EII on).

Therefore, while not distracting ourselves too much from the original egs++ DBS implementation, these tests have given rise to a couple of UBS-related questions that need to be addressed:

1. What is the source of the large discrepancy at low energies with BEAMnrc UBS when Russian Roulette is on? How does the implementation of Russian Roulette in BEAMnrc differ from that in egs++ (and EGSnrc)? Is this somehow an artifact of how I'm using BEAMnrc as a simulation source in an egs++ simulation when generating BEAMnrc spectra?
2. What is causing the discrepancy at ~60 keV in the egs++ implementation of UBS (with EII on)?

[ftessier]

1. Within the AO, implement the option to transform (rotate/translate) the DBS splitting cone.

2. Also implement an option to read in the (active) source transform so that the the DBS splitting cone can also either: a) move with a dynamic source or b) be applied to the current source in a source collection. This transform would be applied on top of any transform input in the AO itself (i.e., step 1).

Why are we considering a cone, instead of a general target shape for DBS?

[blakewalters]

@ftessier, only because the smartBrems and smartCompton routines in DBS pre-calculate the number of split particles to generate based on the photon angular distributions for the respective interactions and a splitting field with a circular base (i.e. a cone). I guess we could use a rejection technique to make the shape of the field more general, and ultimately we may want to try this, but it doesn't seem like a priority at this point in the implementation.

[ikawrakow]

I guess we could use a rejection technique to make the shape of the field more general, and ultimately we may want to try this, but it doesn't seem like a priority at this point in the implementation.

It is obvious, right? And so elegant - just do UBS and Russian Roulette away the photons not going towards the FOI. What was Kawrakow thinking?

It was 20 years ago, but if IIRC, the speed gain from estimating the probability to emit bremsstrahlung or Compton-scatter into the FOI was very significant. The probability estimate used in the BEAMnrc implementation is just that - an estimate. It overestimates, and the overestimation can be quite significant in some cases. When I was implementing the treatment head simulation code for ViewRay I was able to derive a better estimate (compared to what is in BEAMnrc) and that sped it up by a factor of 2, IIRC. If you want to give the user the ability to define arbitrary FOIs, that will come at very significant efficiency cost.

[blakewalters]

@ikawrakow, you really are at large (but hopefully not at loose ends)! Yes, I suspected it would be a nontrivial matter to give the user the option to define an arbitrarily-shaped splitting field. Thanks for confirming this! And that offering the flexibility of defining a sawtooth splitting field or one that says "Bernie 2028" is really not worth the hit to efficiency.

[blakewalters]

Update 01/08/2025

Recently, the ability to apply an affine transform to the DBS splitting cone (defined by the splitting ssd and radius, fs, of the splitting field) has been implemented. This allows the user to user to orient the splitting geometry to match that of the e- beam and bremsstrahlung target.

To put this functionality through its paces, I evaluated the difference between the 225 kVp photon transmission spectra generated with the transforms illustrated in the figure below (T1-T3) and that with no transform applied to the splitting geometry (T0). In each transformed case, the transform applied to the DBS splitting cone matches an equivalent transform of the incident e- beam/target geometry. In all cases, DBS splitting parameters are: nsplit=1000, ssd=10 cm, fs=10 cm.

In T1, the target geometry/splitting cone has been translated +10 cm in both X- and Y-directions; in T2, the geometry has been rotated 45-degrees counterclockwise about the Y-axis; and in T3, the geometry has first been rotated 180-degrees about the Y-axis then translated -10cm in the Z-direction.

Fractional differences between the spectrum with no transform and spectra with transforms T1-T3 are shown in the figure below:

Note that I've exaggerated the -ve portion of the vertical axis to show that differences in low energy (<20 keV) components of the spectra are within uncertainty, although the photon fluence in this energy range is a factor of 2-3 lower than elsewhere in the spectra.

With the exception of low energy photons, this plot shows agreement to within ~0.05% between results with transforms applied to the splitting cone and the untransformed case, thus giving us confidence that transforming the splitting cone does not alter the results. Even so, some questions remain:

1. In order to obtain approximately the same uncertainty with T2 (45-degree rotation of the geometry about the Y-axis) as in the untransformed case, I had to increase the no. of incident histories by a factor of 4. Why?
2. Despite the fact that spectra < 20 keV agree within uncertainty, why are the differences so large, and why are they consistently -ve (i.e. higher than in the untransformed case)?

I make no claim that these cases represent an exhaustive test of the effects of transforming the splitting geometry, and suggests for further testing are welcome!