Pencil beam kernel-based dose calculations on CT data for a mixed neutron-gamma fission field applying tissue correction factors

Objective. For fast neutron therapy with mixed neutron and gamma radiation at the fission neutron therapy facility MEDAPP at the research reactor FRM II in Garching, no clinical dose calculation software was available in the past. Here, we present a customized solution for research purposes to overcome this lack of three-dimensional dose calculation. Approach. The applied dose calculation method is based on two sets of decomposed pencil beam kernels for neutron and gamma radiation. The decomposition was performed using measured output factors and simulated depth dose curves and beam profiles in water as reference medium. While measurements were performed by applying the two-chamber dosimetry method, simulated data was generated using the Monte Carlo code MCNP. For the calculation of neutron dose deposition on CT data, tissue-specific correction factors were generated for soft tissue, bone, and lung tissue for the MEDAPP neutron spectrum. The pencil beam calculations were evaluated with reference to Monte Carlo calculations regarding accuracy and time efficiency. Main results. In water, dose distributions calculated using the pencil beam approach reproduced the input from Monte Carlo simulations. For heterogeneous media, an assessment of the tissue-specific correction factors with reference to Monte Carlo simulations for different tissue configurations showed promising results. Especially for scenarios where no lung tissue is present, the dose calculation could be highly improved by the applied correction method. Significance. With the presented approach, time-efficient dose calculations on CT data and treatment plan evaluations for research purposes are now available for MEDAPP.


Introduction
For the application of external radiotherapy with fast neutrons, only two facilities are actively treating cancer patients.At the clinical neutron therapy system (CNTS) of the University of Washing in Seattle, USA, and at the Tomsk Polytechnic University in Tomsk, Russia, fast neutrons with high energies are generated from protons with E = 50.5MeV respectively deuterium ions with E = 13.6 MeV incident on beryllium targets (Moffitt et al 2018, Specht et al 2015and Polytechnic University Tomsk 2023).The high-energy neutrons produced by the CNTS have depth-dose characteristics similar to 6 MV x-rays (Moffitt et al 2018) which offers additional skin sparing in contrast to lower neutron energies.The planning (Moffitt et al 2020 andViscariello et al 2021) and delivery of intensity modulated neutron therapy (IMNT) was introduced with the CNTS in October 2022 (Sandison et al 2023 andMoffitt et al 2023).
The experience collected over the last decades with fast neutron therapy (FNT) as a therapeutic modality with high linear energy transfer (LET) radiation is summarized in several general and facility-specific review articles for example by Jones and Wambersie (2007), Specht et al (2015), Jones (2020), and Gordon et al (2022).All mentioned publications indicate a high biological effectiveness and a low enhancement by cell oxygenation for fast neutrons as well as a low dependency of the biological effect on the cell cycle.Adequate dose calculation and advanced beam collimation are crucial in order to minimize dose to organs at risk (Jones 2020).
At the research reactor FRM II in Garching near Munich, the medical applications facility MEDAPP for FNT with fast fission neutrons is currently undergoing a general upgrade program in order to re-establish FNT in Garching for specific palliative treatment indications.The upgrade consists of the development and installation of an improved multi leaf collimator (MLC) introduced in 2019 (Sommer et al 2023) and an ongoing dosimetric revision of the mixed neutron-gamma radiation field.
Up to now, a total of 841 patients were treated with fast fission neutrons in reactor-based FNT between 1985 and 2015 in Garching at FRM II and its predecessor research reactor.At MEDAPP, the mixed neutron-gamma radiation field is produced in two uranium converter plates providing an effective area source of 150 × 150 mm 2 .The neutron and gamma dose rates at 20 mm depth of a water phantom are  = -D 0.52 Gy min n 1 and  = g -D 0.20 Gy min 1 , respectively (Wagner et al 2008) and the half maximum depth of the neutron dose rate in water is between 50 and 60 mm (Wagner et al 2008 andSommer et al 2023).Neutrons have a mean energy of ¯= E 1.9 MeV n (Breitkreutz et al 2008).Further information on MEDAPP and the spectral characterization of the neutron and gamma components of the beam were discussed by Wagner et al (2008), Breitkreutz et al (2008), Garny et al (2011a), andJungwirth et al (2012).
The depth dose characteristics of the MEDAPP fission beam with the absence of a build-up region and the steep dose fall-off only allow the treatment of superficial lesions and target volumes close to the surface.At MEDAPP, the main treatment indications were skin lesions followed by head and neck carcinomas, lymph nodes and sarkomas (Specht et al 2015).As discussed by Specht et al (2015), the treatments led to an improved quality of live for individual patients due to good tumor responses.
As highlighted by Jones (2020), accurate dose calculation is of highest importance.While first steps towards a treatment planning software were performed in the past (Garny et al 2009 andGarny et al 2011b), no software for dose calculation on patient CT data was available for MEDAPP.In the work presented here, the decomposition of pencil beam kernels (PBK) according to Bortfeld et al (1993) for three-dimensional therapeutic dose calculations for the mixed neutron-gamma radiation field is presented as the next step to overcome the shortage in dose calculation software.
The similarities between neutrons and photons in the characteristics of depth dose curves (DDC) and the exponential attenuation in water motivate the application of dose calculation algorithms originally developed for photons also for neutron radiation (Söderberg et al 2003 andKalet et al 2013).Two sets of PBKs were generated, one for neutron radiation and one for gamma radiation.The beam parameters necessary for the generation of PBKs and their validation are output factors (OF), percentage depth dose (PDD) curves and beam profiles.While OF were measured in water, Monte Carlo (MC) simulations were run to include the build-up region close to the phantom surface and to distinguish between the direct and the scatter components of the neutron and the gamma dose.Dose distributions from MC calculations in water were validated with reference to measured data and dose distributions calculated using the PBK approach were validated against MC calculations.For the purpose of dose calculation using PBKs, the implementation of the PBK algorithm in the open source research treatment planning software matRad was used (Wieser et al 2017).A correction method for tissue heterogeneities suggested by Söderberg et al (2003) was applied.

Measurement of dose distributions in water
Dose measurements from the initial dosimetric characterization of the new MLC from 2020 (Sommer et al 2023) were used for verification.The determination of the dose deposition in water was performed applying the two chamber method suggested for dosimetry in mixed neutron-gamma radiation fields in ICRU Report 26 (ICRU 1976), AAPM report no.7 (AAPM 1980), and DIN 6802-6 (German Institute for Standardization 2013).The ionization chambers were manufactured by PTW-Freiburg.While the TM33053 chamber is made of tissue equivalent (TE) A150 plastic, the TM33054 chamber is made of magnesium.They provide high respectively low sensitivities to neutron radiation and approximately similar sensitivities to gamma radiation.During the measurements, the chambers were flushed with methane based TE gas and argon, respectively.Even though recommended in AAPM report no.7 and DIN 6802-6, no displacement correction factors for the effective point of measurement of the chambers were applied.Adopting multiplicative correction factors provided by both dosimetry guidelines for dosimetry in a mixed neutron-gamma field would be very challenging and is out of the scope of this work.
Dose measurements were carried out in a cuboidal water phantom of size 520 × 635 × 635 mm 3 manufactured by PTW-Freiburg.OF were measured on the central beam axis in a reference depth d ref = 50 mm in the phantom for all square fields of side lengths 34 mm, 58 mm, 82 mm, 116 mm, 146 mm, and 176 mm at the MLC exit.Due to limitation in beam time at FRM II, neutron and gamma DDC and beam profiles were only measured for measured fields with side length 82 mm and 116 mm.The beam profiles were measured in 30 mm and 50 mm depth.Due to the finite extension of the ionization chambers and the 20 mm PMMA tank wall, the measurements started in 30 mm depth.
2.2.Monte Carlo simulation of dose deposition in water and tissue Dose distributions were simulated using the Monte Carlo code MCNP version 6.2 (Werner 2017).As spectral input for the MC calculations, the results from Breitkreutz et al (2008) were used for neutrons and a modified version of the gamma spectrum discussed by Jungwirth et al (2012) was used.Two separate simulations-one for the neutron component and one for the gamma component-were run for each square field.In order to save computational resources, the transport of neutrons, gammas, and secondary charged particles was not followed through the beamline repetitively.Instead, a simplified geometry of the MEDAPP beamline was used to generate approximate particle velocity vector distributions at the MLC exit with the spectra known at patient position.The beamline was voided and except for components directly in the beam path, radiation transport was switched off so that only the effect of the geometric limitations and the finite extension of the fission source was incorporated.Therefore, no scattering effects from the beamline or from interactions in the MLC were included.
For the calculation of the dose, a water cuboid with the same outer dimensions as the water phantom was positioned at patient position downstream from the MLC exit.Tank walls were not included in the simulation.Particle generation and transport were switched on for neutrons, photons, electrons, positrons, protons, and all heavier ions available in MCNP.A cutoff energy E n,cut = 0 MeV was used for neutrons and for all other particles the cutoff energy was set to E cut = 1 keV.
Total and particle-specific energy deposition on the central beam axis for the depth dose curves were tallied using +F6-tallies and F6-heating tallies, respectively, so that correct summation of energy transfer from primary to secondary particles was guaranteed (Werner 2017).Along the central beam axis, the tallies were defined as cylindrical cells of 2 mm height with a diameter of d cyl = 15 mm and the radius vector oriented perpendicular to the beam axis.To tally the contribution of scattered particles towards the central axis for the decomposition (Bortfeld et al 1993), each tally was flagged to record the energy deposition inside the respective cylinder and the contribution from particles entering the cylinder through the lateral area.The MCNP importance inside the cylinders was set to 7 for variance reduction.
To tally the lateral beam profiles, TMESH-tallies with a grid resolution of 2.5 mm and a total number of approximately 5 × 10 5 voxels were used.TMESH-tallies are equivalent to F6-tallies but are defined as a threedimensional grid.Mesh tallies were used only for neutrons in order to save computational resources.
In order to check the correction method for different tissue types suggested by Söderberg et al (2003), simulations were run where the phantom dimensions and the tallies on the central axis were kept unchanged but water as transport medium was replaced by soft tissue, bone, and lung tissue or by different sequences of tissue slabs.One slab phantom was modeled to be comparable to the one used by Söderberg et al but with the difference that only three tissue types were used in this work and the separation of soft tissue into adipose and muscle was not done.The elemental tissue compositions were defined as tabulated in ICRU Report 46 (ICRU 1992) for adults and tissue densities for the calculation of correction factors were set according to the data for material compositions provided by the National Institute of Standards and Technology (NIST 2022).Element cross sections for the material definitions were taken from ENDF/B-VII.1 (Chadwick et al 2011) and only for hydrogen, the older ENDF/B-VI.6:Xwas used.For each slab phantom, one simulation with a field of side length 82 mm was run.The importance in the central cells was reduced to 3 and no TMESH-tallies were used.
The simulations were run in parallel mode on a virtual Windows 10 machine which runs on a cluster and is attributed with 48 Xeon Platinum 8160 CPUs with 2.1 GHz and two logical processors each and 130 GB of RAM.For water as transport medium, 10 8 starting neutrons and gammas were sampled to ensure accurate particle velocity distributions at the MLC exit.Simulations took 2300 h of computation time for neutrons and 1600 h for gammas.For the simulations of dose deposition in tissue with 10 9 starting particles, the maximum simulation time was 260 h for bone.The reduced runtime was mainly due to the reduced importance.

Dose calculations using pencil beam kernels
For the dose calculation in voxelized CT data, the method described by Bortfeld et al (1993) in combination with a ray tracing algorithm according to Siddon (1985) is implemented in matRad (Wieser et al 2017) and was used in the work presented here.Basic scripts and instructions for the generation of PBKs for matRad are provided on the matRad GitHub web page by Bangert and Ensminger (2020) and were adopted here.
Reference positions: Since the original formulation of the pencil beam decomposition is based on a point source approximation, the extension of the converter plates made it necessary to use a virtual point source position for the dose calculation.The source to axis distance (SAD) of this virtual point source was approximated to be infinite by setting it to SAD virt = 50 m.The SAD virt value was defined in a heuristic manner in order to allow both the calculation in a fan-line coordinate system as described by Bortfeld et al (1993) and implemented in matRad (see Wieser et al 2017) and to account for the extended area source at MEDAPP.
For the description of the field sizes at reference position 550 mm downstream from the MLC exit, the actual geometrical distance between the converter plates and the reference position SAD geom = 578 cm for MEDAPP was used.Side lengths of the fields at reference position in air were calculated to be 38 mm, 64 mm, 91 mm, 128 mm, 161 mm, and 195 mm using a point source approximation for the square field sizes measured at the MLC exit.
Pencil beam kernel decomposition: The calculation of the dose deposition D irr at a lateral position (x p , y p ) for irregular radiation fields F(x, y) along a pencil beam in depth d introduced by Bortfeld et al (1993) is performed according to: Here, the components D i (d) are the contributions of the direct and scattered contribution for the smallest square field with side length 34 mm for i = 1 and i = 2, respectively, and the scattered contribution for a large square field with side length 146 mm for i = 3.In addition to the analytical approach suggested by Bortfeld et al (1993) to describe D i (d), a second exponential term was introduced in the following equation to allow higher flexibility in the fitting process of D i (d) to DDCs calculated from measured OFs and MC simulations: Here, Σ indicates the total macroscopic cross section.For gamma radiation, Σ is replaced by μ.The weights w i in equation (1) were obtained by differentiation and normalization of the initial weights W i from the fitting process of D i (d) and the products of w i and D i (d) is referred to as pencil beam kernels.The source dependent primary energy fluence matrix Ψ(x, y) was calculated according to the geometric approach described by Breitkreutz et al (2008) taking into account the effective area of the fission source and the geometrical limitation of the radiation field by the MLC.The calculation was done for a square field with 176 mm side length at the MLC exit.The convolution of equation (1) with a Gaussian function originally suggested to account for the spatial extension of the source was omitted.
Implementation of dose calculation: PBK calculated by the procedure described above were generated in a way to be compatible with the open source matRad version Alan v2.1.0.The following adjustments were made for the calculation of mixed neutron-gamma dose distributions for the MEDAPP beamline.First, the second exponential term used in the convolution of equation (2) was added to the dose calculation function.Second, the convolution with the Gaussian was switched off.Finally, the dose calculation was performed in two steps in order to calculate both the dose contribution by neutrons and by photons.
Corrections for tissue heterogeneities: For photons, the radiological depth of each point in the voxelized geometry is calculated in matRad using a lookup table to convert Hounsfield units (HU) to a voxel-specific relative electron density to scale the linear attenuation coefficients along the path through the voxels (Wieser et al 2017 andSiddon 1985).For neutrons a similar approach can be used due to the similar exponential depthdependence of the neutron attenuation.Söderberg et al (2003) suggest that for the case of neutron propagation through tissue the linear attenuation coefficient should be replaced by the total macroscopic cross section Σ t and scaled relative to the cross section of water by a tissue correction factor c t ≔ Σ t /Σ w in order to calculate the radiological depth.The total macroscopic cross section is calculated from the total atomic cross section σ t by multiplication with the atom density N. Average atomic cross sections sF t , for the neutron fluence Φ n (E) present at MEDAPP were calculated from energy-dependent cross section data for water, soft tissue, lung tissue, and bone according to the following equation: The final correction factor for lung tissue was gained by an iterative minimization the sum of squares of residuals between the PBK result and MC calculation.Söderberg et al also suggest that in order to calculate the dose to tissue from the dose calculated using the PBK algorithm the result is multiplied by a KERMA correction factor k t .It is defined as the ratio of the average KERMA factor for the tissue t over the average KERMA factor for water.Averages KERMA factors for the MEDAPP neutron spectrum were calculated similar to equation (3) for water, soft tissue, lung tissue, and bone.Data for the KERMA factors were taken from ICRU Report 63 (ICRU 2000).In order to apply the correction factors in dose calculations on CT data, Hounsfield unit (HU) intervals were used for the categorization of individual voxels as suggested by Schneider et al (1996) or DeMarco et al (1998).In this approach, CT voxels are categorized as air, lung tissue, soft tissue, and bone using HU intervals from −1000 to −950, from −950 to −170, from −170 to 280, and from 280 to 4000, respectively.

Pencil beam kernels for neutron and gamma components of the fission beam
In figure 1(a), the direct and scattered contribution D i (d) of secondary ions on the central beam axis from neutron interactions simulated for the small square field and the scattered contribution for the larger square field are shown.Analytical fits to the simulated data according to equation (2) are plotted in black.The DDC were normalized to the total dose deposition for the respective field size in 50 mm depth.For neutron interaction in water, the macroscopic cross section indicated by Σ in equation (2) fitted to the direct depth dose curve of the smallest field was obtained to be Σ w,fit = 0.023 mm −1 and calculated using equation (3) to be Σ w,calc = 0.026 mm −1 .The deviations between the fitted and the calculated value is −12%.β-values fitted according to equation (2) are listed for direct and scattered component in table 1. Dose curves for the different components were also simulated and fitted for the summed primary and secondary gamma radiation.For the gamma radiation (not shown here), the build-up effect is more pronounced.An attenuation coefficient  Field-size dependent weights W i (r) for the three components of the neutron kernel are shown in figure 1(b).While the weights for the direct component are rather independent of the field size, the scattered component for the small and large field show strong dependence on the field size.The fitted weights for the scattered contribution are dominated by the large field contribution for four out of six fields.A negative value is reported for the small field scatter contribution for field sizes larger or equal to 128 mm.This result is interpreted as compensation for an overestimation of the large field scatter contribution in the fitting process (see Bortfeld et al 1993).

Evaluation of calculated dose deposition in water
For neutron radiation, figures 2(a) and (b) show comparisons between calculated PDDs from MCNP (green) and PBK (purple) calculations and between measured PDDs (blue) and PDDs from MCNP calculations in water for fields with side lengths of 82 mm and 116 mm, respectively.All PDDs were normalized to d ref = 50 mm.The black solid curves associated to the right ordinate indicate the percentage deviation of the PBK results with respect to the MC dose calculations.For MC dose calculation as reference, relative errors between ±10% are achieved within the first 100 mm for all field sizes.For larger depths, an increasing overestimation of the dose deposition by the PBK calculations is reported with reference to the MC calculations.In general, the build-up region calculated by MC simulations is well reproduced by PBK calculation with a low relative error.Deviations between MC calculations and depth dose measurements for the two available field sizes are shown as dotted black curves in figure 2. Relative errors below ±15% are calculated for the first 150 mm.Larger errors are reported for larger depth but no recognizable trend is present.
In figures 3(a) and (b), a comparison between calculated and measured beam profiles in 50 mm depth is shown for field sizes of 82 mm and 116 mm, respectively.The color coding of dose curves and error plots is identical to figure 2. The comparison of beam profiles calculated with PBKs and MC in two depth at 30 mm and 50 mm shows good agreement in the plateau region with maximum relative errors of about ±10%.Only for the smallest field of size 34 mm, larger deviations of up to ±20% are reported for the central region of the beam  profile.In the fall-off and low dose regions, an overestimation of the dose deposition with large relative errors by the PBK algorithm is reported in comparison to the MC calculations.The comparisons between the MC simulations and measurements given by the black dotted lines in figure 3 show better agreement also in the falloff and low dose region of the beam profile.
Figure 4 shows a comparison between calculated and measured DDC for the gamma component of the mixed radiation field in water for fields of side lengths 82 mm and 116 mm with normalization to 50 mm depth.For the comparison of the PBK calculation and MC simulations indicated by the black solid lines, large relative deviations up to 25% are reported for the build-up region in the first few millimeters for all field sizes.For larger depths, relative errors mostly within ±5% with outliers up to ±7% were calculated.From a direct comparison of the MC and the PBK calculations it can be seen that dose deposition from secondary gammas from neutron capture reactions in the fall-off region of the depth dose curve are poorly represented.Especially for the large field of side length 116 mm this effect is visible for depths between 60 and 150 mm in figure 4(b).
As indicated by the dotted black curve, the comparison between measured data and data from MC simulations shows increasing relative errors with increasing depth.A clear overestimation of the dose deposition by the MC calculation up to 16% is visible.
A comparison of normalized beam profiles with PBK calculations with measured data for the two available field sizes at 30 mm and 50 mm depths shows maximum deviations in the plateau regions between 5% and 10% and again larger deviations in the fall-off and low dose regions.

Evaluation of calculated dose deposition in heterogeneous slab geometry
Correction factors for the radiological path length were calculated to be c l ≈ 0.25, c st ≈ 0.95, and c b ≈ 0.91 for lung (l), soft tissue (st), and bone (b), respectively.Best agreement for optimizing c l ä [0.25, 0.45] with reference to MC calculations with lung tissue as transport medium was achieved for c l = 0.36.After the optimization for lung tissue, relative errors of ±10% were reported for the first 100 mm in the comparison between MC and PBK calculations for all tissue types.KERMA correction factors were calculated to be k l = 0.96, k st = 0.93, and k b = 0.38.
In figure 5(a) calculated neutron DDC for alternating slabs of soft tissue and bone are shown.Material boundaries are indicated as black vertical lines and soft tissue and bone slabs are colored in light and dark gray, respectively.In figure 5(b), an additional layer of lung tissue of thickness 106 mm indicated by an intermediate gray was inserted as fourth layer.This last slab phantom was modeled to be comparable to the one used by Söderberg et al (2003) but with the difference that only three tissue types were used in the work presented here.In both plots, results from MC calculations as simulated for the different tissue compositions are indicated in red.Results from the PBK algorithm with both corrections for attenuation and dose to tissue conversion applied are shown in blue.Relative deviations between the MC and the PBK calculation are plotted as solid black line.All curves were normalized to 20 mm depth to avoid a normalization with respect to the boundary region between soft tissue and lung in figure 5(b).
Figure 5(a) shows that the results from the PBK calculations can be highly improved by the application of the correction procedure for soft tissue and bone.Nonetheless, high deviations of 10%-20% from the MC reference calculation are reported for the first 200 mm and even higher deviations are present for larger depths.Large deviations are reported for the lung slab in figures 5(b) where the PBK algorithm overestimates the dose to lung tissue by up to 55%.

Evaluation of calculated dose deposition on CT data
In the following, a qualitative assessment of the performance of the PBK algorithm on CT data is provided.As an example, the planning CT of an actual treatment at MEDAPP was used.For the calculation, a radiation field of 58 × 58 mm 2 was defined and a voxel size of 1.5 × 1.5 × 3 mm 3 was used where the CT slice thickness was 3 mm.The original planning target volume (PTV) and the spinal cord are indicated in figure 6 in purple and dark red, respectively.The axial view shows the CT slice that incorporated the matRad isocenter lying within the PTV.The dose distribution is normalized to 1.0 Gy in the isocenter and isodose lines for the dose distribution overlaid on the CT data were set to 0.25 Gy, 0.5 Gy, 0.75 Gy, and 1.0 Gy.The dose calculation for the single radiation field shown here took 2.3 seconds on an Intel i7-10700 CPU with 2.9 GHz.The effect of the correction factor for bone is clearly visible in the dose distribution.This can be seen when comparing the dose deposition within the jawbone or the vertebra to the deposited dose to soft tissue in the vicinity of these bony structures.As expected, the dose deposited within bone structures is significantly reduced.Furthermore, the effect of the air cavity in the center of the CT slice is clearly visible.Here, an increase of dose deposition downstream from the cavity is visible.

Discussion
The generated PBK for the mixed neutron-gamma field allow the time-efficient reproduction of DDC from Monte Carlo simulations within tens of seconds for the neutron and gamma components separately.For the first time, dose calculation on CT data is available for MEDAPP.The neutron dose component can be weighted by an individual factor to account for the relative biological effectiveness of neutrons.
For the neutron dose deposition in water, agreement within ±10% for the first 100 mm depth in water was demonstrated.Depending on the field size, this depth comprises the 70%-85% dose fall-off region.
As comparison, DDC calculated using the commercial treatment planning software for photons applied in clinical FNT at the CNTS in Seattle show deviations from measured data of less than 3% (Kalet et al 2013).It should be kept in mind that the mean neutron energy at the CNTS is significantly higher than the one at MEDAPP and therefore the depth dose characteristics also differ significantly.
As presented in this work, dose calculations for heterogeneous slab phantoms could be highly improved by the application of tissue-dependent corrections that account for the tissue composition.The performance in the two material slab phantom of soft tissue and bone shows promising improvements in comparison to the uncorrected calculations with respect to MC calculations.Nevertheless, remaining relative deviations of ±10%-20% are reported.The dose to lung is highly overestimated by the pencil beam approach.The main reason for the high deviations in the lung slab is most likely an overestimation of the scatter component by the pencil beam approach (see Söderberg et al 2003).Here, the large mean free path in lung can only be accounted for by the correction for the radiological path length on the central beam axis.The scatter contribution remains uncorrected for in the PBK algorithm so that it is expected to be overestimated.In contrast, MC methods can directly account for scattering effects towards the central beam axis.The depth-dependence of the spectral neutron fluence also influences the result of the correction factors for the PBK calculations.
Generally speaking, the results from the PBK algorithm with the heterogeneity correction applied is in good agreement with the results presented by Söderberg et al (2003).There, the agreement between MC simulations and the outcome of the photon pencil kernel dose calculation from the applied treatment planning software is also highly improved by the application of heterogeneity corrections.Similar to the findings presented here, an overestimation of dose to the lung slab by the treatment planning software is reported.Söderberg et al (2003) report relative deviation of their photon treatment planning software with respect to MC calculations for a Bragg-Gray cavity filled with water positioned in the center of the tissue slabs.For muscle, adipose tissue, and bone relative deviations between +1% and −11% and for lung tissue an overestimation of +8% of the treatment planning software with respect to MC calculations are reported.The deviations increase with increasing depth and an underestimation of more than −3% is reported for muscle and adipose tissue only for depths larger than 17 cm.

Conclusion
In conclusion, the presented time-efficient 3D dose calculation along with the separation of neutron and gamma dose deposition for our fission source provides high potential for future research.While the presented dose calculation approach is not implemented in a clinical treatment planning software, it allows dose calculation on CT data for the mixed neutron-gamma radiation field for the first time.
Further challenges that need to be addressed in the future can clearly be identified.Deviations shown in the comparison of measured data and MC simulations need to be investigated further.Here, the gamma energy spectrum is of special interest.Currently, runtime intensive MC calculations are the only way for a verification of the discussed correction method for tissue heterogeneities applied in the PBK approach.Dose measurements using radiosensitive films in heterogeneous tissue-equivalent slab phantoms are currently under investigation and need to be conducted in the verification process of the PBK calculations.

Figure 1 .
Figure 1.Pencil beam kernels for neutron component: (a) analytical fits to simulated data for direct and scattered contribution for neutrons.(b) Fitted weights for direct and scattered components to simulated dose for different field sizes with the same color coding as in (a).

Figure 2 .
Figure 2. Comparison of the measured depth dose curve (blue) and the calculated ones using MCNP (green) and the PBK algorithm (purple) for the neutron component (a) for field size 82 mm and (b) 116 mm.Relative errors of the PBK calculations with respect to MC simulations are given as black solid line and relative deviations between MC data and measurements are given as black dotted line.

Figure 3 .
Figure 3.Comparison of the measured beam profiles (blue) and the calculated ones using MCNP (green) and the PBK algorithm (purple) for the neutron component at 50 mm depth for field size (a) 82 mm and (b) 116 mm.Relative errors of the PBK calculations with respect to MC simulations are given as black solid line and relative deviations between MC data and measurements are given as black dotted line.

Figure 4 .
Figure 4. Comparison of the measured depth dose curve (blue) and the calculated ones using MCNP (green) and pencil beam algorithm (purple) for the gamma component (a) for field size 82 × 82 mm 2 and (b) for field size 116 × 116 mm 2 .Relative errors of the PBK calculations with respect to MC simulations are given as black solid line and relative deviations between MC data and measurements are given as black dotted line.

Figure 5 .
Figure 5. Dose calculation using the PBK algorithm and MC simulations for the neutron component in a slab phantom (a) of alternation soft tissue and bone in light and dark gray, respectively, and (b) with an additional layer of lung tissue in intermediate gray.The field size was set to 82 × 82 mm 2 and relative errors are indicated in black.

Figure 6 .
Figure 6.Dose distribution of neutron radiation field overlaid on CT data.

Table 1 .
Fitting parameters β for equation (2) for the neutron and gamma components of the beam.
μ γ,fit = 0.005 mm −1 in water was obtained and fitted β-values for the gamma component are again given in table 1.