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]

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!

[blakewalters]

Update 02/26/2025

In the 3 most recent commits to this PR, I've implemented a charged particle splitting option in egs++ DBS. Inputs for this option existed prior to this, but were not enabled.

Charged particle splitting can be used when one is interested in recovering some of the contaminant electrons which are otherwise subject to Russian Roulette when DBS is used in simulating photon beams. This then allows the user to examine quantities such as the spectrum and fluence of these electrons, contaminant dose, etc.

In BEAMnrc DBS, e-/e+ splitting is done at a horizontal region boundary, the "splitting plane", within a component module (CM--always FLATFILT). Charged particles crossing this region boundary are split nsplit (i.e. the bremsstrahlung splitting no.) times with their weight reduced by a factor of 1/nsplit, with an option to redistribute split particles evenly about a circle of radius sqrt(x^2+y^2), where x and y are the X- and Y-positions of the particle at the splitting plane prior to splitting. This latter option is known as "radial redistribution" of split particles, and can improve spatial sampling of charged particles for beams with radial symmetry. In addition, the user is also asked to input the Z-position of a Russian Roulette plane (zrr) below which Russian Roulette is no longer played with charged particles and Russian Roulette IS played with photons about to undergo photoelectric, pair production or Compton events. These events are then split to generate nsplit (or nsplit*2 in the case of pair production) low-weight charged particles. In general, zrr is < the Z-position of the splitting plane. For more information on e-/e+ splitting, see the BEAMnrc Manual.

The implementation of charged particle splitting in egs++ DBS is similar to that in BEAMnrc DBS. Inputs are also similar, with the difference that the splitting plane is not restricted to a region boundary within a CM but can be any region or group of regions in the geometry. On entering this(these) region(s), phat charged particles are immediately split and possibly radially redistributed about the Z-axis. This offers greater flexibility, but may also prove a pitfall, in that it is difficult to imagine a situation where the splitting plane should not be perpendicular to the beam axis and encompass the entire beam field (at that plane).

Also, in egs++ the charged particle splitting plane is subject to any transform applied to the DBS splitting cone (see above). Thus, radial redistribution of split particles about the transformed plane would actually place particles outside their original (untransformed) region. To avoid this possibility, if the radial redistribution option is used, we always redistribute split particles evenly about the Z-axis in the original (untransformed) geometry. This mens that split particles may not end up being evenly distributed about the beam axis (which, presumably, is the same as the transformed axis of the splitting cone), and so we advise against radial redistribution of split charged particles in cases where you are also applying a transform to the DBS splitting cone.

Onward, then, to testing of the e-/e+ splitting option, in which I repeat many of the tests above for the 225 kVp transmission beam but now examine both photon and electron spectra at the splitting ssd. In the case of electron spectra, I've increased the field radius considered to 20 cm (splitting radius=10 cm) in an attempt to improve the statistics of the e- spectrum (and recognizing that the fluence of contaminant e- is not restricted to the DBS splitting field).

More to come, stay tuned!

[blakewalters]

Update 02/27/2025

...continued from previous update.

As mentioned above, tests of e-/e+ splitting use the same 225 kVp photon transmission setup used in earlier tests, with a point source of e- incident on an infinite W slab (target) of thickness 0.01 cm. Below the target is vacuum.

The first comparisons between BEAMnrc and egs++ DBS with the e-/e+ splitting option involved placing the splitting plane at the bottom of the W target (at Z=0.01 cm) and the Russian Roulette plane (zrr) at 1 cm (i.e. in vacuum below the target). Since interactions resulting in charged particles (pair production, photoelectric, Compton) are only split below zrr, this ensures that all charged particles reaching the splitting plane are phat and, thus, serves as a test of splitting (and radial redistribution) at this plane.

The figure below shows the e- spectra 10 cm below the target (i.e. at the DBS ssd) in a circle of radius 20 cm:

In BEAMnrc and egs++, charged particle splitting at the splitting plane has resulted in improved statistics in the e- spectra, and superficial agreement between BEAMnrc and egs++ e- spectra can be seen. A plot of particle weights (not shown) indicates that all e- reaching the scoring plane have weight=1/1000, as expected.

The plot below shows the fractional difference between the BEAMnrc and egs++ e- spectra in this case:

This shows agreement between BEAMnrc and egs++ e- spectra within an uncertainty of 2-10%.

The figure below shows the fractional difference between BEAMnrc and egs++ photon spectra (scored in a field with r=10 cm at the ssd) for this e-/e+ splitting case:

Similar to earlier comparisons without e-/e+ splitting, BEAMnrc and egs++ photon spectra agree to within ~0.05%

In order to make a more definitive statement about the agreement between BEAMnrc and egs++ e-spectra, I wanted to generate e- spectra with uncertainties < 0.1%. I found that with zrr=1 cm, however, this required > 100x10^9 histories, which presented a challenge given the computing resources available. With this in mind, I varied the value of zrr in an attempt to optimize e- statistics, ultimately arriving at zrr=0.005 cm (i.e., halfway through the target), which allowed me to generate spectra having < 0.1% uncertainty at their peak (kinetic energy=40 keV) using 40x10^9 histories.

The plot below shows the fractional difference between BEAMnrc and egs++ e- spectra with this new value of zrr:

From this plot, we can infer that BEAMnrc and egs++ e- spectra agree to within < 0.5% over most of the energy range, certainly a stronger statement regarding agreement than we were able to make above.

Comparing the photon spectra with this value of zrr (and no. of histories), though, we obtain the following:

Which shows that, while agreeing within the uncertainty of the calculation, the egs++ photon spectrum is consistently higher than the BEAMnrc spectrum by ~0.01%. For reference, the plot also shows the comparison between photon spectra using identical simulation parameters but without e-/e+ splitting. This latter comparison exhibits something closer to the expected variations around zero, leading us to conclude that use of e-/e+ splitting introduces a systematic positive bias in the photon spectrum.

The plot below shows the fractional difference between photon spectra with e-/e+ splitting on and off for both BEAMnrc and egs++:

The plot clearly indicates that e-/e+ splitting introduces a systematic increase at the ~0.01% level in the egs++ photon results, while this is not the case for BEAMnrc.

Although this is a very small-scale (higher-order?) issue, it may end up having a larger effect at different energies and/or source configurations, and so it seems that some further debugging of the e-/e+ splitting algorithm in egs++ DBS is warranted. Stay tuned...again!