Enhancement of the EGSnrc code egs_chamber for fast fluence calculations of charged particles

Purpose Simulation of absorbed dose deposition in a detector is one of the key tasks of Monte Carlo (MC) dosimetry methodology. Recent publications (Hartmann and Zink, 2018; Hartmann and Zink, 2019; Hartmann et al., 2021) have shown that knowledge of the charged particle fluence differential in energy contributing to absorbed dose is useful to provide enhanced insight on how response depends on detector properties. While some EGSnrc MC codes provide output of charged particle spectra, they are often restricted in setup options or limited in calculation efficiency. For detector simulations, a promising approach is to upgrade the EGSnrc code egs_chamber which so far does not offer charged particle calculations. Methods Since the user code cavity offers charged particle fluence calculation, the underlying algorithm was embedded in egs_chamber. The modified code was tested against two EGSnrc applications and DOSXYZnrc which was modified accordingly by one of the authors. Furthermore, the gain in efficiency achieved by photon cross section enhancement was determined quantitatively. Results Electron and positron fluence spectra and restricted cema calculated by egs_chamber agreed well with the compared applications thus demonstrating the feasibility of the new code. Additionally, variance reduction techniques are now applicable also for fluence calculations. Depending on the simulation setup, considerable gains in efficiency were obtained by photon cross section enhancement. Conclusion The enhanced egs_chamber code represents a valuable tool to investigate the response of detectors with respect to absorbed dose and fluence distribution and the perturbation caused by the detector in a reasonable computation time. By using intermediate phase space scoring, egs_chamber offers parallel calculation of charged particle fluence spectra for different detector configurations in one single run.


Introduction
Publication of the international IAEA Code of Practice TRS-398 [4] in 2000 turned out to be the basis for a considerable quality improvement in dosimetry for external radiotherapy [5]. Since that, a number of further developments have taken place in the field of dosimetry for external radiotherapy. An important advancement is due to Monte Carlo (MC) simulation of radiation transport which has become a widely used technique in the accurate calculation of dosimetric quantities for all beam types. Such quantities are stopping-power ratios, dose conversion factors, or perturbation correction factors required both for reference ionization chamber dosimetry and for dosimetry in non-reference conditions. A comprehensive review was recently published [6]. MC calculated data have indeed superseded many of the approximations used to determine the data in TRS-398 and are therefore explicitly considered in a forthcoming update of TRS-398. Additionally, MC calculated beam quality correction factors are already established in other dosimetry protocols such as the AAPM's addendum to the TG-51 report [7] and the German DIN report 6800-2:2020-08 [8].
A key role is now taken by the MC calculated dose conversion factor which is defined as the ratio between the absorbed dose to water at a point of interest and the absorbed dose in the sensitive volume of a detector used for dosimetry. With this method, the factorization of the dose conversion factor into the stopping power ratio and one or more perturbation correction factors and hence the associated limitation to Bragg-Gray conditions can be avoided [1,9]. A particularly appropriate tool to produce data for this approach is the MC code egs_chamber [10]. In order to reduce calculation time, a series of variance reduction techniques (VRTs) have been implemented in this user code. Additional to the VRTs photon splitting and Russian Roulette [11] which are already available in cavity, the photon cross section enhancement (XCSE) [12,10] is implemented. In XCSE, free electrons, positrons or scattered photons are induced along a photon's trajectory corresponding to increased reaction cross sections, while the original photon stays on its original track. To preserve the conservation of energy the generated particles have a reduced statistical weight. In addition, the methods of intermediate phasespace scoring (IPSS) and correlated sampling (CS) are included which allow the simulation of more than one geometry at the same time.
Another advantage of using MC simulations is to enable a closer look at the role of the charged particles, i.e. both primary and secondary electrons and positrons. There are a series of papers in which the fluence of charged particles differential in energy was recently considered in great detail [13,14,1,15]. It particularly turned out that understanding how the charged particle fluence differential in energy is influenced by radiation conditions and by detector properties can considerably contribute to better understand the response characteristics of this detector [2].
Recent papers particularly address the dosimetrical quantity restricted cema (acronym of: converted energy per mass) which can serve as bridge between absorbed dose and the charged particle fluence [16,3]. In these papers, a restricted cema based formalism for the determination of absorbed dose was suggested. For this purpose, but also for many other dosimetry related studies on the role of charged particles as listed above, there is an increasing use of detailed spectral fluence data. However, only a few MC codes provide these data. Examples are the codes FLURZnrc and cavity, both included in the EGSnrc system. Both codes have some limitations in application: the code FLURZnrc uses a cylindrical coordinate system which is not always appropriate to simulate a desired detector geometry while the code cavity is comparably inefficient for detector simulations. The code egs_chamber, on the other hand, better meets these requirements, however, it does not yet provide fast calculation of spectral charged particle fluence. The aim of this paper is to describe an enhancement of this code to become a general and at the same time fast MC calculation tool to provide data on absorbed dose and spectral fluence for various real detector and irradiation conditions. Results for calculated charged particle fluence, restricted cema as well as on the improved calculation efficiency are presented and compared with the results of existing MC codes for fluence calculations. The quantity fluence, U, is given by U ¼ dN =da, where dN is the number of particles incident on a sphere of cross-sectional area da. Thus, it refers to a mathematical point in a medium. However, it can be well approximated by the average fluence in a small volume which is assigned to that point. This approximation is most appropriate for MC simulations and therefore used in this work. The MC calculation of fluence within a volume is frequently based on a theorem of Kellerer [17], where the mean fluence is given by the sum of the track lengths of the particles, ds, per unit volume dV : with unit cm À2 . The quantity of interest in this work is the mean fluence differential in energy, U E , with unit cm À2 Á MeV À1 . The method based on Eq. (1) is used, for example, in the code FLURZnrc. A compilation of further methods appropriate for the implementation in MC codes is found in [15]. Since one of the methods was particularly useful for this work, it is also shortly described in the appendix.

Restricted cema
Based on the definition of ICRU Report 90 [18], the restricted converted energy per mass (cema) is calculated as: where U E is the distribution of the total electron and positron fluence with respect to energy E at a point in a medium, L D =q is the corresponding mass linear energy transfer and TE D denotes the track end term [19,18,20] representing the sum of the kinetic energy of all electrons and positrons in the scoring volume with energies below the selected energy threshold value D. Restricted stopping power S D;tot =q which is sometimes used in the definition of restricted cema is not exactly equal to linear energy transfer, however, at a value of D ¼ 10 keV, the difference is negligible [20].

Determination of efficiency
The metric to quantify the efficiency e of a Monte Carlo simulation can be determined as using the total CPU hours T needed on a given computer to achieve a certain estimated statistical uncertainty r. Since this is a general formalism, many quantities such as absorbed dose, fluence and restricted cema can be evaluated with that relationship. Any approach which decreases the time corresponding to a certain uncertainty or vice versa but does not injure the original result is called variance reduction technique.

Fluence calculations
In this work, four different MC codes of the EGSnrc system were applied for the calculation of the spectral fluence of charged particles and compared: The first two codes are generally available, whereas the two modified codes represent in-house developments that would be ready for further distribution on request. Since cavity is not restricted to RZ-geometries and provides virtually all investigated quantities out of the box we consider it as the archetype of this work. The code DOSXYZnrc was originally developed for dose calculations in a rectilinear voxel system [22]. DOSXYZnrc was modified by one of the authors in order to offer more versatile applications, for instance for dose calculations in other volumes such as cylinders, spheres and also for many common detector types, or to permit fluence calculations in such volumes. Additional information on the method for fluence scoring is provided in [1,15].
The modifications made to the user code egs_chamber in this work are briefly described next. Since cavity provides the full capability of fluence scoring and uncertainty estimation, the algorithm extracted from the corresponding section ausgab has been adopted. With this approach the built-in variance reduction techniques from egs_chamber such as photon cross section enhancement, Russian Roulette, intermediate phase space scoring and a region-based ECUT are conserved out of the box. However, the remaining VRT onegeom 2 and correlated sampling are not provided.

Calculation of restricted cema
Using the data of the spectral electron and positron fluence obtained as described in the section above, the computation of restricted cema is performed by where U E;i is the spectral fluence including liberated secondary charged particles, S D =q ð Þ i is the restricted total mass stopping power at the energy of each bin center i beginning at the first bin above D (i.e. i D ), dW is the energy corresponding to the bin width, and TE D is the track end term representing the kinetic energy (E kin ) of all charged particles with energies below the threshold value for particle tracking. The restricted total stopping power and restricted linear energy transfer are equal in a good approximation since their difference for D ¼ 10 keV is negligible in our study [20,23]. Based on the assumption that all energy contributions E dep below D are deposited locally, the track end term could be therefore directly calculated [3] as: This formula was implemented in the modified versions of DOSXYZnrc and egs_chamber. However, this formula is usually not implemented in standard codes. Under the condition that a MC code provides the data of spectral fluence the restricted cema without track end term can be computed. If the absorbed dose is accessible too, the track end term can then be obtained due to the close agreement between absorbed dose D and restricted cema [16,18] by: For all EGSnrc user codes, the calculation of the relative standard uncertainty of restricted cema was performed using the history-by-history statistical estimator method [24] as applied for the uncertainty determination of absorbed dose.

Test setups and comparison
Results obtained with the modified user code egs_chamber were compared with those from FLURZnrc (which was specifically designed for fluence computation), from a modified version of DOSXYZnrc, and from cavity. Results were obtained for fluence spectra of electrons and positrons, for restricted cema, and for data on efficiency. Fig. 1 shows the investigated test setups. Two external radiation setups were simulated for a collimated point source geometry at SSD = 95 cm. As external sources the 6 MV Mohan spectrum [25] and the Co-60 Mora spectrum [26] were used. The third setup consisted of a point source emitting the spectrum of an Ir-192 microSelectron v2 HDR source immersed at the center of the water phantom. Results were scored in a small water disc of 0.1 cm radius and height or in the sensitive volume of a fully modelled PTW 31014 ionization chamber. Table 1 shows the chosen EGSnrc Monte Carlo parameters and simulation setups.
However, normalization of charged particle fluence computation is not performed identically in the investigated MC codes. One reason are different definitions of the incident primary photon spectrum: For a collimated point source FLURZnrc and DOSXYZnrc use the primary photon fluence per beam area (unit cm 2 ) at the SSD. Because the beam area is cancelled out, the secondary charged particle fluence is then given in units of MeV À1 . In cavity and egs_chamber the incident photon spectrum is referred to steradian and independent on the SSD. For comparison, all external beam radiation spectra were therefore normalized per incident photon fluence per MeV À1 independent of the SSD but normalized to the field size at reference depth. The immersed point source is normalized per initial particle in every user code thus the secondary charged particle spectra have units MeV À1 cm À2 . Regarding comparable examples as shown in Fig. 1 and no use of VRTs, all investigated user codes yield similar calculation efficiencies. Therefore all efficiency enhancements can be compared to the original cavity code. The variance reduction technique from egs_chamber used in this work was photon cross section enhancement. To investigate the resulting efficiency gain, the detectors were embedded in water shells of varying sizes with radii r in which XCSE was applied, as shown in Fig. 1. For the ionization chamber the water shells had a fixed length extending outside the active volume. For water disc cases both radii and height were extended by the same amount. Each simulation had one fixed XCSE value that was assigned to all detector parts including an extended stem and the surrounding shell.
It is noteworthy that the investigation of VR-techniques can become quite tedious due to the interference of the Figure 1. Schemes of the investigated beam setups. In all cases the XCSE water region in which XCSE is switched on in the egs_chamber simulations is marked by red cylinders with radius r. The water disc XCSE shell (shown in (a)) also varies the height in beam direction by the same. In case of an ionization chamber (shown in (b) and (c)) the height of the cylinder is fixed to 1.5 cm outstanding the active volume whereas the radius varies. Example (a): External point source collimated to a 4 Â 4 cm 2 square field size at 100 cm distance to the detector. The detector (here: blue water disc) is placed at 5 cm depth in a 30 Â 30 Â 30 cm 3 cubic water phantom. Example (b): Same as (a) but with different detector such as PTW 31014 ionization chamber. Example c): Immersed point source located at the center of a 30 Â 30 Â 30 cm 3 water phantom. The detector's reference point (here: PTW 31014 ionization chamber with its sensitive volume in white) was at the same spot translated one centimeter away from the source. An analog immersed-source setup replacing the ionization chamber by a water voxel and XCSE region is not explicitly shown. Drawings are not to scale. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) different techniques [31] or the dependency on radiation setups because several VRT parameters have to be considered together in every beam setup [10]. For this reason we investigated the efficiency gain due to XCSE only. Additionally, it is found in [31] that XCSE offers the greatest efficiency enhancement when the VRTs are considered individually. The aim at this point is to compare the charged particle fluence scoring techniques and the corresponding efficiency gain resulting from XCSE in egs_chamber systematically in a given scenario.

Charged particle fluence spectra and restricted cema computation
As initial verification -and based on the same initial random numbers -the enhanced egs_chamber as well as the original egs_chamber and cavity user codes yield exactly the same results regarding all investigated quantities if no additional VRTs from egs_chamber are used. This is even the case for charged particle fluence data although not provided by the original egs_chamber user code. Fig. 2 shows the electron spectra differential in energy and the corresponding restricted cema in the volume of interest obtained using the user code cavity and the corresponding relative deviations from egs_chamber, FLURZnrc and DOSXYZnrc. For the sake of clarity the spectra of secondary positrons differential in energy are not shown since they are at most responsible for only about 1 % of the response in a 6 MV Bremsstrahlung photon beam. The analyzed volume refers to the sensitive, air-filled cavity of the fully modelled ionization chamber, type PTW 31014.
All datasets belonging to one simulation setup agree with each other within their statistical uncertainties. Numerical deviations between cavity and the other user codes exceed 5 % only at high energies. For the line spectra of Co-60 and Ir-192 the corresponding Compton edges at about Absorbed dose, restricted cema, track end terms and charged particle fluence differential in energy to a water disc (radius and height 1 mm) or ionization chamber (PTW 31014). Detailed simulation setups are shown in Fig. 1. # histories/ statistical uncertainty   Table 2 shows the results concerning absorbed dose, restricted cema and track end term for external irradiation with Co-60 photons as in Fig. 1(a). The normalization was performed independent of the SSD for a beam area of 4 Â 4 cm 2 at 100 cm distance to the source. Again the results agree within their statistical uncertainties.
For radiation qualities of 6 MV Bremsstrahlung and higher energies the contributions from positrons to restricted cema exceed 1 % and are therefore not negligible. Additionally, the relative contributions from TE D are energy dependent and vary from 6% for 6 MV Bremsstrahlung up to 14 % for Ir-192. Fig. 3 shows the relative efficiency gain of egs_chamber simulations with various XCSE factors and shell thicknesses in the calculation of absorbed dose or restricted cema without track end term. The track end term was omitted in the efficiency evaluation because the restricted cema was calculated in the user code based on the charged particle fluence scoring algorithm which ignores particle energies below D. For these calculations a PTW 31014 ionization chamber was placed in a water phantom and externally irradiated with Co-60 photons as shown in Fig. 1(b). In this example the efficiencies first increase with increasing enhancement of the cross sections. This behavior reaches a maximum after which the computation time to trace the growing number of secondary charged particles initiated by the higher cross sections begins to increasingly reduce the efficiency.

Efficiency evaluation
The photon cross section enhancement shell thickness also has an impact on the simulation efficiency. Similar to an increased XCSE value, a thicker shell around the detector leads to more photon interactions causing an increase in charged particles scored in the detection volume. And again more interactions require more computation resources. Both with enhancement of the value of the cross sections and enhancement of the volume in which XCSE is applied, the relative efficiency gain depends on the computed quantity. In the example shown in Fig. 3 the maximum achievable relative efficiency gains concerning restricted cema without TE D and absorbed dose calculations Table 2 Absorbed dose D, restricted cema C D (as calculated by Eq. (4) including TE D as calculated by Eq. (5) and track end term TE D per incident photon fluence (i.e. in Gy cm À2 ) calculated for the irradiation of a water disc with Co-60 photons in the external beam setup from Fig. 1 a). The relative standard uncertainty of the final digit is shown in parentheses for each value.
are different and the efficiency peaks are achieved at different values of XCSE. Table 3 shows the relative efficiency gain for all investigated beam setups at the optimum XCSE value and shell thickness. Several observations can be made: First, the relative efficiency gain depends on the radiation quality and increases with decreasing particle energies. Second, the optimum cross section enhancement value increases with decreasing particle energies whereas the optimum shell thickness shows an inverse behavior. Third, the efficiency improvement potential is larger in ionization chamber simulations in every investigated setup. Fourth, the optimum XCSE conditions for absorbed dose and restricted cema without TE D calculations are equal for water voxels. The parameters differ only slightly for the ionization chamber simulations where the setups referring to restricted cema without TE D tend to smaller enhancement values.

Discussion
As stated in the EGSnrc C++ class library [21] the track length based fluence scoring algorithm as used in FLURZnrc assumes the stopping powers to be constant along a charged particle step. This might affect spectra at regions where the progression of the stopping powers is steep. Our results show that this applies only to the first bin above D (i.e. 10 keV) in water disc simulations, which yields about 5 % more than the other algorithms in certain energy bins. One workaround is the reduction of the MC parameter ESTEPE so that the condensed history steps are small enough for an accurate scoring with constant stopping powers. This miscalculation on charged particle fluence scoring in FLURZnrc does not impair the resulting restricted cema since its impact is considerably smaller than the statistical uncertainty of restricted cema calculation. Thus the effect is negligible in our study and may be accounted for in the overall uncertainty budget.
The investigated VR-technique photon cross section enhancement is well established since several years in various EGSnrc user codes [11,28,12]. With the introduction of XCSE in egs_chamber, the application to charged particle fluence computations was explicitly mentioned in [10]. Furthermore the EGSnrc standard statistical uncertainty history-by-history estimator is also explicitly dedicated to arbitrary quantities such as fluence [24]. XCSE as used in egs_chamber offers a dramatic efficiency gain especially for small scoring regions compared to the irradiated volume [10]. Our results confirm this especially regarding an internal irradiation as shown in Fig. 1(c).
The relative efficiency gains between absorbed dose and restricted cema without TE D differ for ionization chamber simulations. Fluence scoring in cavity is done using track-length estimation by precomputed stopping powers from the EGSnrc database. Restricted cema can then be obtained as described in Section 3.2. Absorbed dose scoring tallies up the energy deposited due to ionization (continuous energy losses) and discrete interactions. As the energy of the charged particles changes, the scoring can be more or less efficient and hence the dose scoring may lead to a different statistical uncertainty than restricted cema computation for the same number of traced particles. A similar mechanism could occur in the calculation of kerma using photon fluence and mass-energy absorption coefficients in comparison to using the kerma approximation of depositing the energy of secondary charged particles on the spot.
By using other variance reduction techniques such as intermediate phase space scoring in the enhanced egs_chamber code the computation of several detector compositions needed for perturbation correction factors may be obtained in one single simulation run. Fig. 4 depicts one practical example of the combined simulation of the components related to the dose conversion factor for an ionization chamber split into stopping power ratio ( Fig. 4(a) to (b)) and several perturbation correction factors (Fig. 4(b) to Table 3 Relative efficiency gain in the calculation of absorbed dose and restricted cema without track end term with egs_chamber at the optimum values of photon cross section enhancement conditions compared to no use of VRTs. The parameter r depicts the surrounding XCSE shell where r ¼ 0 means an enhancement of the cross sections of only the components of the ionization chamber. The calculations were performed in the external beam setup for 6 MV Bremsstrahlung and Co-60 photons as shown in Fig. 1(a) and (b), and in the immersedsource setup for Ir-192 photons as shown in Fig. 1(c (c) and (c) to (d)). Applying this example to the setups as shown in Fig. 1 with external or immersed source and adjusted XCSE conditions for each detector, a further relative efficiency enhancement of 1.52, 1.57 and 1.61 for external 6 MV Bremsstrahlung, external Co-60 and immersed source Ir-192 photon beams could be observed respectively. If the phase space scoring volume is large enough and contains all detector components in every position, several detector positions inside a phantom can be simulated simultaneously. An additional enhancement would be the implementation of correlated sampling [28,21] for charged particle fluence scoring so that spectral resolved fluence or restricted cema quotients can be calculated.

Conclusions
This work shows a straightforward enhancement of egs_chamber to calculate charged particle fluence using methods already present in the EGSnrc framework and in particular in cavity. The EGSnrc applications are built modular so that there is not much interference with the inserted code. The comparison with several EGSnrc user codes show good agreement which proves that the implemented charged particle fluence scoring in egs_chamber does not impair the original code. Furthermore, in our calculations, the different radiation setups (external or immersed source) and the use of variance reduction techniques showed no influence on the charged particle fluence spectra.
Using VRT provided by EGSnrc such as XCSE an improved simulation efficiency could be observed. There are even more techniques to raise the efficiency implemented in egs_chamber such as a region-based ECUT, intermedi-ate phase space scoring (IPSS), Correlated Sampling and onegeom. The latter allows several simulation geometries to be equal while using different scoring regions parallel. Region-based ECUT, XCSE and IPSS work out of the box whereas Correlated Sampling and onegeom require additional modifications.
With the given approach it is possible to achieve spectral information even for sophisticated demands within acceptable computing time. Additionally the restricted cema as a link between absorbed dose and charged particle fluence can also be calculated. Since the restricted cema methodology has the advantage of being applicable in setups which do not fulfill Bragg-Gray conditions [16] the enhanced egs_chamber code is well suited for further investigations on detector developments and non-reference conditions. tron path length is obtained either in the single scattering mode, i.e. between two scattering events, or in the condensed history mode, i.e. between two hard (inelastic) interaction sites. In both cases, the relation between a given path length to the energy loss over this path and vice versa must be known. For this relation, the 'continuous slowing down approximation' is assumed where the energy loss results from sub-threshold processes only. Thus, the energy loss DE of an electron or positron is formally linked to its corresponding path length s through: where L D is the linear energy transfer, also called restricted stopping power. It follows that an electron or positron path can be characterized by three parameters of kinetic energy: The energy loss DE associated to the path length s, the energy E A at the path start, and E B ¼ E A À DE at its stop. The position of this parameter with respect to a linear energy grid 0, E 1 , E 2 , . . . E max is shown in Fig. 5. The method of fluence calculation is now essentially based on two steps: Step 1: Each energy deposition (equivalent with the energy loss DE) is split into a number of single subsections adjusted to the energy binning which will be used to present the electron or positron fluence differential in energy. Using the example of Fig. 5, the first subsection is E i -E B , the intermediate subsections are identical to the bin width dW , and the last subsection is E A -E i+j . The quotient of each subsection by DE is termed a.
Step 2: By means of the quotient a i;j where the index j refers to a specific energy loss DE during the MC simulation, the spectral fluence U E;i at bin i within the scoring volume V can now directly be derived.
The two methods used in this work are: 1. track length based: U E;i ¼ P M j¼1 1 V Á ai;j dW Á s j where j refers to one particle track out of the total number of histories M, and s j is the path length associated to the energy loss DE j ; this method is implemented in FLURZnrc.

energy deposition based:
DEj L D;i where L D;i is the linear energy transfer of the material in the scoring volume at bin i. In the modified version of DOSXYZnrc, L D;i is taken at the energy of bin center, whereas in cavity, it is -according to the documentation -at the lower energy of each bin. Actually, however, it is also at bin center.