VARIAN CLINAC 6 MeV Photon Spectra Unfolding using a Monte Carlo Meshed Model

Energy spectrum is the best descriptive function to determine photon beam quality of a Medical Linear Accelerator (LinAc). The use of realistic photon spectra in Monte Carlo simulations has a great importance to obtain precise dose calculations in Radiotherapy Treatment Planning (RTP). Reconstruction of photon spectra emitted by medical accelerators from measured depth dose distributions in a water cube is an important tool for commissioning a Monte Carlo treatment planning system. Regarding this, the reconstruction problem is an inverse radiation transport function which is ill conditioned and its solution may become unstable due to small perturbations in the input data. This paper presents a more stable spectral reconstruction method which can be used to provide an independent confirmation of source models for a given machine without any prior knowledge of the spectral distribution. Monte Carlo models used in this work are built with unstructured meshes to simulate with realism the linear accelerator head geometry..


Introduction
Detailed knowledge of the spectrum emitted by the linear accelerator is essential, due to the fact that dosimetric factors depend directly on the beam energy.The principal objective of accurate dose calculation is delivering a prescribed dose to a target volume while minimizing radiation damage to the surrounding healthy tissues.
However, the energy spectrum of a Bremsstrahlung photon beam used in cancer treatments is difficult to be characterized by direct measurements because the emitted photon flux is very high and it produces the detector`s saturation [1].
For this reason, a methodology which guarantees a realistic shape of the spectrum emitted by LinAc is necessary.This work presents a robust and practical procedure to obtain reconstructed spectra by using depth dose measurements in a water phantom and Monte Carlo simulations [2].
A first order Fredholm integral equation allows to relate the depth dose curve measurements in the water phantom`s central axis with the Bremsstrahlung spectrum that forms it, but the solution to this equation is not simple due to the ill-conditioned nature of the mathematical problem.
The work is focused on solving the inverse problem based on reconstructing a spectrum from its experimental depth dose curve.The reconstructed spectrum will be used as source spectrum in a Monte Carlo simulation to compare depth dose curves between experimental measures and simulated results, validating the performed procedure.
It is here where it becomes important to have a detailed, accurate and realistic Monte Carlo model of the linear accelerator head to run a simulation which takes into account each material and geometric measures according to manufacturer data.
To that, this study has used the last version of MCNP6.1.1Monte Carlo code to simulate the LinAc head and water phantom irradiation.MCNP v.6 is the first Monte Carlo code which allows introducing meshed geometries into simulations [3].
The use of meshed geometries to generate the threedimensional models presents the advantage of a higher accuracy in the geometry modelling, which is reflected in a precision improvement of the absorbed dose estimation.
According to this, the developed procedure presents an alternative method to calculate a linear accelerator photon beam spectrum based on indirect measurements.Specifically, the study has been focused on the reconstruction of a 6 MeV primary photon beam issued by the medical linear accelerator Varian Clinac available at the Hospital Universitari i Politècnic La Fe, Valencia, Spain.

Spectrum reconstruction basic theory
The unfolding methodology used in this work for spectra reconstruction was developed and applied in previous works obtaining accurate results in other LinAc machines and for different energy beams [4].

Linear Problem and Discretization
The photon attenuation from the treatment beam is governed by an exponential law whose attenuation coefficients vary with the energy of these photons, therefore the first order equation of Fredholm can relate the energy deposition to different depths of photons in the water phantom.This equation can be expressed as: where d(x), is the depth dose, x is the water phantom depth, ‫)ܧ(߯‬ the incident energy spectrum and ‫ݎ‬ ா ‫)ݔ(‬ a response matrix that contains simulated depth dose curves for different consecutive mono-energetic beams.
If the preceding equation is inverted, it makes it possible to calculate the spectrum which produces the depth dose curve.Nevertheless, the inverse matrix problem is ill-conditioned due to the nature of the problem.For this reason, to solve the equation it is necessary to discretize the problem into different intervals of depth positions.
The energy is discretized as well in different energy bins of the spectrum ‫,)ܧ(߯‬ where ߯ is the photon spectrum with energy between ‫ܧ‬ ିଵ and ‫ܧ‬ .
On the other hand, ‫ݎ‬ ா ‫)ݔ(‬ can be approximated by a M x N matrix, A: where ‫ܣ‬ is the deposited dose in a position i, for a given energy spectrum j.
Thus the equation becomes as follows: where ߯ ⃗ = (߯ ଵ , … , ߯ ) is the original unknown beam and ݀ ⃗ = (݀ ଵ , … , ݀ ) ் is the depth dose curve which has been obtained experimentally at the Hospital Universitari i Politècnic La Fe .
To obtain the ‫ܣ‬ matrix (M × N) sequential Monte Carlo (MCNP6.1.1)simulations were executed acquiring depth dose curves for different mono-energetic incident beams with different energies.These doses have been calculated as a sum of the dose due to electrons and the dose due to photons.
Figure 1 shows the simulated results which configure the response matrix.This response matrix (dim.300x350) with energy range between 0-7 MeV (energy bin 20 KeV) and the depth stepped by millimeter were calculated by simulation taking the contribution of photons and electrons.
A small line of noise can be seen in this figure, it is probably due to the fact that these values have been obtained from a simulation using low number of simulated particles and this produces statistic errors higher than rest.This fact had no relevance or significant influence in the calculation of the reconstructed spectrum.

Fig. 1. Response matrix with energy range between 0 -7 MeV, A є RM×N
For the 6 MeV energy range we are considering, depth dose curves from mono-energetic spectra were calculated using the MCNP6 Monte Carlo code.The Source Surface Distance (SSD) from a point source on the central axis to a homogeneous water phantom has been 100 cm.The used field size has been 10 cm × 10 cm.Each Monte Carlo simulation was run until the uncertainty in all evaluated points of computed depth dose curves was less than 0.3%.To register the contribution of photons and electrons, the FMESH tally has been used with its respective flux-to dose conversion factors, implemented as the DE/DF target in MCNP.
Once the response matrix is obtained, the mathematical procedure is applied to reconstruct spectrum.The discretization of the problem offers a linear system equation matrix: The original photon spectrum has been defined here (Eq.6) as the solution vector ߯.

Singular Value Decomposition Method (SVD)
In this case, the matrix equation ‫ܣ‬ • ߯ = ݀ allows to obtain the original spectrum if the response matrix is known but several error affect the spectrum´s determination as explained above and it can´t be obtained directly, hence we use the singular value decomposition method (SVD): (7) where ‫ݑ‬ and ߥ are the orthonormal singular vectors and ߪ the singular values.The same equation can be expressed by matrix form: where ߑ is a diagonal matrix with its singular values ߪ in the diagonal and the rest zeros, while ܷ and ܸ are the matrix of the singular vectors.
The SVD algorithm works very well for noise-free data.For data affected by noise, direct SVD reconstruction is extremely inaccurate but there are several regularization schemes that improve its performance.
Tikhonov methodology has been used to solve this matrix system [5], [6].With this numeric algorithm, an optimum solution ߯ can be obtained by generating a new response matrix, removing the parts of the solution corresponding to the smallest singular values.Tikhonov regularization is the most common and well-known form of numerical stabilization.This methodology corresponds to making the normal system of the least-squares problem ‫ܣ‬ ் ‫ݔܣ‬ = ‫ܣ‬ ் ݀ better conditioned.
The idea is to define the general regularization solution xλ as the minimum of the following weighted combination of residual norms and constraint: where L is the discrete derivative operator and x* is the residual solution.L is a general matrix with column rank complete, and which is usually related to some derivation operator.The classic standard form of Tikhonov regularization considers that L is the identity matrix, so the solution is: Hence, the equation is equivalent to solving the linear system: where the regularization parameter ߣ controls the weight granted to the minimization of the restriction of the part relative to the minimization of the residual norm.Clearly, a high λ (equivalent to a large amount of regularization) favors a small semi-norm solution and a large residual norm, while a small λ (i.e. a small amount of regularization) has the opposite effect.Thus, the regularization parameter is an important parameter that controls the properties of the regularized solution and therefore should be carefully selected.
One way to determine the optimal regularization parameter, λ, is to plot the L curve, which in an ideal environment has a relatively vertical segment and a relatively horizontal line.For this curve, the norm-2 of the solution vector is drawn, ‫‖ݔ‖‬ ଶ (ordinate), against the norm-2 of the residual vector, ‫ܣ‖‬ • ‫ݔ‬ − ݀‖ ଶ (abscissa) for different values λ.As it has been explained, it is recommended to place the truncation parameter λ, at the value closest to the point of maximum curvature of this curve, that is, by selecting the value corresponding to the corner that conforms to L.

Monte Carlo Simulation
The Varian CLINAC head was modeled in detail according to the manufacturer's specification and each element is assumed to have the typical natural abundance of different nuclides.
The LinAc of the present study is the Varian CLINAC 6 MeV photon beam configuration installed at the Hospital Universitari i Politècnic La Fe which is equipped with an adjustable jaw photon collimator as well as with an multi leaf collimator (MLC) with 120 leafs.The accelerator geometry has been accurately modeled as shown in figure 2.

Fig. 2. Three-dimensional model of the Varian CLINAC
The model has been generated using 3D Modeling Software for Engineering ANSYS SpaceClaim creating quickly and easily the solid model CAD.
Later, the 3D solid model was meshed with Abaqus/CAE [7], which is a software used for finite element modeling and analysis of mechanical components and assemblies.It has been used to define the meshed geometry that will be introduced in the MCNP6.1.1 simulation.
Meshed geometry modeling tool is a novelty of MCNP version 6.1.1,(older MCNP versions and other Monte Carlo codes do not present this capability) and it is an important advantage when modelling complex geometries with precision.

Fig. 3. Varian CLINAC meshed with Abaqus/CAE
The use of meshed geometries, presents the advantage of a higher accuracy in the geometry modelling which is reflected in an increased absorbed dose estimation precision.Unstructured meshes enable the use of cells of different types and/or different size, allowing complex geometries and optimizing the number of cells used according to the needs.That is why, for the present work, we have used unstructured grids imported from Abaqus/CAE.Finally, the Varian CLINAC simulation has been performed using the reconstructed spectrum as a source.
The simulation was divided in two parts in order to speed up the calculations.The first simulation corresponds to the Varian CLINAC head part and the particles are stored in a surface source (phase space) located behind the MLC.This phase space was generated by MCNP6.1.1 simulation.The second simulation corresponds to the water phantom and starts from the previous surface source generated with more than 10 7 collected particles.[8] To speed up the calculations without compromising the accuracy of the results, the MCNP6 (version 6.1.1)code has been parallelized in SENUBIO ISIRYM research group´s cluster Quasar using the MPI protocol with 26 processors.Finally, the computational real time real time was about 48 hours.
Concerning radiation transport, a detailed physics treatment for a coupled photon and electron mode simulation has been considered, taking into account the following physical processes: photoelectric effect with fluorescence production, Compton and Thomson scattering, and pair production in the energy range between 0.001 and 7 MeV.
Figure 4 shows the photon flux obtained by simulation passing through the whole Varian CLinac head.
An important mesh advantage is the easy results visualization.The results supplied in ASCII by means of huge matrix can be converted easily using some code tools or FORTRAN routines to obtain a compatible threedimensional visualization format that can be opened with most of 3D visualization software such as Paraview or Visit [9].By this way, an interactive and complete threedimensional view is achieved and dose or flux results can be shown quickly and easily in any point of the simulated geometries.
Figure 5 shows 3D dose distribution in a water phantom and isodose curves obtained by MCNP6 simulation of a 6 MeV treatment beam transport.

MCNP6 water phantom dose and its isodose curves
visualized with Paraview.

Experimental procedure
In order to validate the Monte Carlo simulation and the unfolded spectrum, experimental measures were obtained using a water phantom to determine the real depth dose curve and compare it with the simulated depth dose curve.These experimental measures were taken at the Hospital Universitari i Politècnic La Fe in Valencia and were performed with a high resolution detector (registering the dose contribution of both electrons and photons) placed in a motorized guide inside the 50 cm side phantom ("PTW 31010 Semiflex").
The semiflex chambers are designed for therapy dosimetry, mainly for dose distribution measurements in motorized water phantoms.They have a short stem for mounting and a flexible connection cable.The nominal useful energy range is from 140 kV to 50 MV photons and 10 ... 45 MeV electrons.
The wall material is graphite with a protective acrylic cover.The guard rings are designed up to the measuring volume.Air density correction is required for each measurement, and a radioactive check device is available as an option.
Chambers are shaped cylindrically with an inner diameter of 5.5 mm.The 0.125 cm 3 chamber is ideal for 3D dosimetry in a water phantom, since the measuring volume is approximately spherical resulting in a flat angular response over an angle of ± 160° and a uniform spatial resolution along all three axes of a water phantom.The water cube was irradiated with a 10 cm x 10 cm field size using a 6 MeV photon beam maintaining the source to surface distance (SSD) equal to 100 cm, source to collimator distance (SCD) equal to 37.9 cm, collimator angle was zero degrees, temperature was 20 degrees Celsius and 1013.25 mbar of pressure.
Depth dose curves and profile curves can be measured accurately in this phantom because precise positioning of high resolution detectors can be easily accomplished using a guide driven by reinforced toothed belts.
Only have been used depth dose curves to validate the spectrum and profile curves have not been used because the curves obtained by mono-energetic beams used in the unfolding process are depth dose curves too.

Results
The reconstructed photon spectra obtained using unfolding Tikhonov method for 6 MeV photon beams is shown in Figure 7.The figure above represents the reconstructed 6 MeV spectrum and its respective polynomial fit.It has been calculated according deconvolution procedure explained.The method is based on a previously simulated response matrix, as well as an experimental depth dose curve.
The reconstructed spectrum has been validated comparing depth dose curves obtained by Monte Carlo simulations including the Tikhonov reconstructed spectrum as a source with experimental depth dose curves measured at the hospital.This comparison is shown in Figure 8; it shows that simulated dose curve fits accurately the measured depth dose curve, being the average relative error is 1.89 %.Experimental values have been compared with the results of the MCNP simulation using as a source the Skeikh-Bagheri spectrum, spectrum commonly publicized and used in several works [10].It can be appreciated in figure 9, that results are worst, the peak is not coincident, and the dose curve using Sheikh-Bagheri Spectrum is below the experimental spectrum.The average relative error in this case is 7.1, higher than the one obtained using the spectrum reconstructed with the presented unfolding technique.
The importance of obtaining a reconstructed spectrum and not using a theoretical spectrum extracted from the literature is because each radiotherapy unit will have slight variations in its spectrum due to minor modifications in the installation of its components and the life of the machine.Thus, if we want a real calculation of the dose given to the patient, we will need a precise calculation of the real spectrum used in the treatment.

Conclusions
The procedure of this study is an important approach that allows characterizing the primary photon spectrum emitted by a LinAc.A computer model of a Varian CLINAC head has been generated to be used by MCNP6.1.1.The model was used to calculate the depth dose curves generated by reconstructed spectra and validate them by comparison with experimental depth curves.
This paper also presents the innovation of using unstructured mesh geometries in Monte Carlo calculations which means a direct impact on calculations because it allows to define the geometry with more accuracy.
Obtained results prove that the methodology described estimates accurately the beam energy spectrum emitted by LinAc.The validation was performed comparing depth dose curve in a water tank between simulated and measurements.

Fig. 8 .
Fig. 8. Depth dose curve obtained with simulation of 6 MeV photon spectrum by Tikhonov method

Fig. 9 .
Fig. 9. Depth dose curve obtained with simulation of 6 MeV photon spectrum from Skeikh-Bagheri by Tikhonov method