A torus source and its application for non-primary radiation evaluation

Smaller (and thus cheaper) synchrotrons, with diameters down to 5 m have been developed in the last decade (Paganetti et al 2021). Compact synchrotrons are widely used as a main accelerator by manufacturers such as P-Cure (Balakin et al 2018, Pryanichnikov et al 2018), Protom (Addendum 2009, Balakin et al 2021), NewRT (Wang et al 2020), APTron (Xu et al 2013), and Hitachi (Hiramoto et al 2007, Umezawa et al 2015). This manuscript investigates the PT system provided by P-Cure as a case study, which can accelerate protons to 30–330 MeV with a 4.5 m diameter ring. However, 250 MeV is sufficient for maximum clinical requirements.

To calculate an accurate absorbed dose distribution, a dedicated geometry with a tandem injector, injection, a synchrotron ring, extraction, a BTL, and a scanning nozzle is constructed and simulated in FLUKA. The layout can be found in figure 2, and the beam loss is shown in table 1.

Zoom In Zoom Out Reset image size

The major difference compared with our previous work is beam loss model of beam pipes through bending magnets (Wang et al 2020). For NewRT system, to simulate the proton-transport process, synchrotron and BTL dipole magnets are set with a given magnet field in their central zone, then the protons loss behavior is controlled by magnet field. The disadvantage of this method is obvious.

• (1)  
  The beam loss direction will be biased to one side if the magnet field is not properly set up.
• (2)  
  Proton of momentum when $p=Gq$ (G: Magnetic rigidity, q: charge of a proton) can travel forever in the accelerator ring without loss, this will underestimate the dose level.
• (3)  
  There is no proton loss on up or down inner pipe shell.

The process of non-primary radiation evaluation is similar to radiation shielding evaluation, as both need to calculate dose generate by all the sources from system and target. However, for shielding, radiation dose outside the bunker is not that sensitive to models from sources but very sensitive to energy and weight. For nonprimary, radiation dose inside the bunker is considered which is sensitive to both source models, and energy weight. Thus, the old model would be not suitable for non-primary radiation evaluation as the IEC regions are close to magnets of synchrotron and BTL. Thus, new models of torus source beam loss model and uniform straight pipe beam loss model are adopted for current research.

There are seven non-primary radiation sources in this case:

• (1)  
  The target itself: the target is a water tank at the isocenter in place of the patient. Planned dose irradiation of patient generates large amounts of secondary patient, including neutrons, photons, protons and so on (Robert et al 2013). All these secondary particles contribute to non-primary dose around the target. The proton beam loss at the target is 31.7%. A 2 Gy spread-out Bragg peak (SOBP) was modulated to approximate clinical conditions as accurately as possible (Wang 2019a). The range depth with different weights ranges from 22 to 32 cm, as shown in figure 3. Weight adaption-based SOBP method is adopted here again same as in our previous work in reference (Wang et al 2020) figure 7. The only difference is the range extend from 8 to 10 cm.
• (2)  
  The scanning nozzle consists of two scanning magnets, a vacuum chamber, and several ionization chambers (ICs). The beam path of the scanning nozzle is in a vacuum chamber. There is approximately 2.7% proton loss in the beam path during scanning, and a dedicated Monte Carlo (MC) nozzle model is implemented in this case (Wang 2019b). The energy ranges from 70 to 250 MeV.
• (3)  
  The BTL is short and consists of several magnets. The center of the BTL is a vacuum tube, and there is approximately 1.6% proton loss evenly distributed on the beam pipe. The energy ranges from 70 MeV to 250 MeV. The uniform straight pipe beam loss model demonstrated by Xu et al (2016) is adopted in this study.
• (4)  
  Extraction: there is 34.4% proton loss during extraction as a concentration point source. The energy ranges from 70 to 250 MeV.
• (5)  
  The synchrotron ring consists of a magnetic system and a radio frequency acceleration system. The magnetic system of the synchrotron is formed by four identical quadrants separated by large free gaps connected by beam pipes. Each quadrant is formed by four C-shaped iron blocks with parallel poles. The four magnets with a homogeneous field in each quadrant are arranged in a pairwise common winding configuration. There is 18.6% proton loss in the center of each quadrant, which can be described by the torus source beam model. The torus source beam loss model will be discussed in the next section. The proton energy in the synchrotron ring ranges from 2 to 250 MeV.
• (6)  
  Injection: the injector is designed for the initial production of protons and their acceleration to an energy of approximately 2 MeV. The injector consists of a pulsed arc source of ions with a pulsed hydrogen inlet, electrostatic lenses, and a tandem high-voltage accelerator. There is 10.9% proton loss in this process, but these protons contribute little to non-primary radiation, as 2 MeV protons mostly generate x-rays.
• (7)  
  Concrete reflection and scattering: secondary particles such as neutrons, photons, and protons from (1) to (6) could be reflected or scattered by the shielding concrete wall around the PT system. This part is difficult to distinguish from the sources of (1)–(6); thus, a full MC model including a concrete wall should be established for non-primary radiation evaluation as shown in figure 4.

All sources are summarized in table 1, and there are two large categories: one is the target at the isocenter in place of the patient, and the other is the system including sources (2)–(7). These two categories are set because there is an advantage of easily distinguishing how the PT system impacts the non-primary dose.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To describe beam loss in the four quadrants as close to realistically as possible for a synchrotron, a torus source model is developed in this work. A torus can be made by revolving a small circle (radius $r$) along a larger circle (radius $R$). As shown in figure 5, a torus beam pipe is defined by two radii $r$ and $R$ and thickness $t.$ Here, $r=5\,{\rm{cm}},$ $R=15\,{\rm{cm}},$ and $t=1\,{\rm{cm}}$ for example. Inside the pipe is a vacuum, and a uniform magnetic field maintains the beam orbit in the center. The vector of loss proton $\vec{{\boldsymbol{V}}}$ (dashed blue line) has a small angle $\theta$ to the beam orbit. The pipe shell is iron, and there are air and magnet components around the pipe. When loss protons hit the pipe shell, very large secondary particles such as neutrons, gamma rays and electrons are immediately generated, which will greatly contribute to non-primary radiation around the target in the treatment area.

Zoom In Zoom Out Reset image size

The torus surface is expressed in parametric equation (1) (Shiohama and Takagi 1970, Georgiev 2019, Colley 2002)

$$\begin{eqnarray}&&\left\{\begin{array}{c}x={x}_{0}+\left(R+r\,\cos \,u\right)\,\cos \,v\\ y={y}_{0}+\left(R+r\,\cos \,u\right)\,\sin \,v\\ z={z}_{0}+r\,\sin \,u\end{array}\right.,\end{eqnarray} \tag{ 1 }$$

where $v\in \left[0,\,2\pi \right)$ is the angle about the z axis for the large circle and $u\in \left[0,\,2\pi \right)$ is the angle for the small circle. The center of the large circle is at $\left({x}_{0},\,{y}_{0},\,{z}_{0}\right).$ This model assumes that the loss protons form a cone beam (blue) and that angle $\theta$ determines the interaction area on the pipe shell, which is marked with red lines in figure 6. There are large red interaction regions when angle $\theta \geqslant 45^\circ ,$ and the whole pipe shell will be directly activated by protons. There are small red interaction regions when angle $\theta \lt 45^\circ ,$ and the activated area will be limited to the outer side of the pipe shell. Thus, the dose distribution will be concentrated on the horizontal plane for the real case of C in figure 6, as the synchrotron can always bound the proton beam in the center orbit.

Zoom In Zoom Out Reset image size

For case C, the thickness of the proton beam that goes through the pipe shell is $T:$

$$\begin{eqnarray}&&T=\displaystyle \frac{t}{\sin \,v}=\displaystyle \frac{t}{\sin \left[{\cos }^{-1}\left(\displaystyle \frac{R}{R+r}\right)\right]}.\end{eqnarray} \tag{ 2 }$$

For the P-Cure PT system, $r=4\,{\rm{cm}},$ $R=180\,{\rm{cm}}$ and $t=0.3\,{\rm{cm}}.$ Then, the thickness can be calculated: $T=1.447\,{\rm{cm}}.$ High-energy protons can easily pass through this thickness and irradiate magnet components. Thus, C-shaped iron blocks for magnets could work as a self-shielding material and should be constructed in the MC model, as shown in figure 2.

The parametric equation (1) define the position domain of any loss proton on the torus surface; correspondingly, another group of parametric equations needs to be derived to describe the travel direction of the protons. The major task here is to decompose vector $\vec{{\boldsymbol{V}}}$ into the x, y and z axes.

Let an off-orbit proton go from point M on the center orbit to point P on the inner pipe surface along vector $\vec{{\boldsymbol{V}}},$ as shown in figure 7(a)). P is on a new auxiliary plane $soz.$ The first step is to decompose $\vec{{\boldsymbol{V}}}$ into ${V}_{1}$ in the $soz$ plane and ${V}_{2}$ in the $xoy$ plane as expressed in equation (3):

$$\begin{eqnarray}&&\left\{\begin{array}{c}{V}_{1}=V\,\sin \,\theta \\ {V}_{2}=V\,\cos \,\theta \end{array}\right..\end{eqnarray} \tag{ 3 }$$

Zoom In Zoom Out Reset image size

Then, ${V}_{1}$ is decomposed into components along auxiliary axis $s$ and axis $z,$ as shown in figure 7(b)) and expressed in equation (4)

$$\begin{eqnarray}&&\left\{\begin{array}{c}{V}_{1s}=V\,\sin \,\theta \,\cos \,u\\ {V}_{1z}=V\,\sin \,\theta \,\sin \,u\end{array}\right..\end{eqnarray} \tag{ 4 }$$

${V}_{1z}$ is a final component of ${V}_{1}$ along axis $z.$ However, ${V}_{1s}$ should be further decomposed into the components along axis $x$ of ${V}_{1sx}$ and axis $y$ of ${V}_{1sy},$ as shown in figure 7(c)) and expressed in equation (5).

$$\begin{eqnarray}&&\left\{\begin{array}{c}{V}_{1sx}=V\,\sin \,\theta \,\cos \,u\,\cos \,v\\ {V}_{1sy}=V\,\sin \,\theta \,\cos \,u\,\sin \,v\\ {V}_{1z}=V\,\sin \,\theta \,\sin \,u\end{array}\right.\,.\end{eqnarray} \tag{ 5 }$$

${V}_{2}$ in the $xoy$ plane is easily decomposed into the components along axis $x$ of ${V}_{2x}$ and axis $y$ of ${V}_{2y},$ as shown in figure 7(d)) and expressed in equation (6). The $z$ axis component is 0 for ${V}_{2}.$

$$\begin{eqnarray}&&\left\{\begin{array}{c}{V}_{2x}=V\,\cos \,\theta \,\cos \,v\\ {V}_{2y}=V\,\cos \,\theta \,\sin \,v\\ {V}_{2z}=0\end{array}\right.\,.\end{eqnarray} \tag{ 6 }$$

The final step is to combine the components in the Cartesian coordinate system from equations (5) and (6), that is, equation (7):

$$\begin{eqnarray}&&\left\{\begin{array}{c}{V}_{x}=V\,\sin \,\theta \,\cos \,u\,\cos \,v-V\,\cos \,\theta \,\cos \,v\\ {V}_{y}=V\,\sin \,\theta \,\cos \,u\,\sin \,v+V\,\cos \,\theta \,\sin \,v\\ {V}_{z}=V\,\sin \,\theta \,\sin \,u\end{array}\right.\,.\end{eqnarray} \tag{ 7 }$$

Thus, equation (1) describe the position domain of any loss proton and equation (7) describe the direction of the off-orbit proton bombardment on the inner shell of the torus pipe in the Cartesian coordinate system. We call the mathematical equations (1) and (7) the torus source model.

Equations (1) and (7) work as a sampling function to generate a torus proton distribution for MC simulations. In order to intuitively investigate the influence of different parameters on the distribution of lost protons, the trajectories are shown in figure 8 for different angles: (A) $\theta =80^\circ ,$ (B) $\theta =45^\circ$ and (C) $\theta =5^\circ .$ In the center it is not a synchrotron ring, but a mathematic torus. Here, $r=5\,{\rm{cm}},$ $R=15\,{\rm{cm}}$ and ${x}_{0}={y}_{0}={z}_{0}=0$ in equation (1). On the torus surface trajectory is dense. The horizontal distribution has a radial shape, and the vertical distribution is similar to the wings of a butterfly. For the case of angle $\theta =5^\circ ,$ the wings are compressed into a flat shape. The torus source model is easy to adjust beam loss in vertical direction by changing angle $\theta ,$ this is another advantage compared to the old one.

Zoom In Zoom Out Reset image size

For a real synchrotron case, as the off-orbit travel direction of the protons is tangent to the orbit circle, angle $\theta =1^\circ$ is adopted in the following calculation for non-primary radiation (Xu et al 2016). This small angle can also guarantee a conservative estimation for non-primary radiation distribution for target as most dose will be concentrated at the isocentric level as shown in the vertical of figure 8. If this worst case non-primary dose can comply IEC standard, then a larger angle should be safe.

FLUKA, Geant4 and MCNPX are general MC toolkits (Böhlen et al 2010, Vanaudenhove et al 2014). The GEANT4 code was originally developed for applications in high energy physics. Some experimental validations of the hadronic physics integrated in the code have been reported (Dedes et al 2014, Pinto et al 2016). Conversely, the FLUKA code has been developed and successfully applied both in the high and the low energy ranges. Extensive validations of experimental data have been published (Ferrari et al 2002, Agosteo et al 2007, Böhlen et al 2014, Kozłowska et al 2019). In the framework of dose calculations in hadrontherapy, the codes were also compared with experimental data. For protons, Parodi et al (2007b) (Parodi et al 2007, Koch et al 2008) showed a good agreement between simulated and experimental depth dose profiles.

FLUKA is adopted to calculate the non-primary radiation dose distribution (Ferrari et al 2005, Böhlen et al 2014, Wang et al 2018, Kozłowska et al 2019, Ahdida et al 2022). SimpleGEO and Flair, which are advanced interfaces for FLUKA, are designed to create a complex MC model, as shown in figures 2 and 4 (Theis et al 2006, Vlachoudis 2009). Non-primary radiation is contributed mostly by secondaries during the beam transport process such as neutron, gamma/xray photon and electron. Thus, PRECISIOn is set as FLUKA suitable physics models for our specified problem which include (Fasso' 2003, Ballarini et al 2004):

• (1)  
  EMF on is used to request a detailed transport of electrons, positrons and photons. Even if the primary particles are not photons or electrons, photons are created in high-energy hadron cascades.
• (2)  
  Low energy neutron transport on down to thermal energies included. Fully analogue absorption for low-energy neutrons.
• (3)  
  Particle transport threshold set at 100 keV, except neutrons.
• (4)  
  Multiple scattering threshold at minimum allowed energy, for both primary and secondary charged particles.

The torus source model and uniform straight model work as sampling functions for the synchrotron and BTL, respectively. Different from the standard circle in section 2.3, there are four identical quadrants in the synchrotron ring. Thus, the torus source model should be modified into four parts with the same weight:

$$\begin{eqnarray*}&&{\rm{Quadrant}}\,1:u\in \left[0,\displaystyle \frac{1}{2}\pi \right),{\rm{where}}\,{x}_{0}=36.29\,{\rm{cm}},{y}_{0}=93.7\,{\rm{cm}},{z}_{0}=42.79\,{\rm{cm}};\end{eqnarray*}$$

$$\begin{eqnarray*}&&{\rm{Quadrant}}\,2:u\in \left[\displaystyle \frac{1}{2}\pi ,\pi \right),{\rm{where}}\,{x}_{0}=36.29\,{\rm{cm}},{y}_{0}=93.7\,{\rm{cm}},{z}_{0}=-42.79\,{\rm{cm}};\end{eqnarray*}$$

$$\begin{eqnarray*}&&{\rm{Quadrant}}\,3:u\in \left[\pi ,\displaystyle \frac{3}{2}\pi \right),{\rm{where}}\,{x}_{0}=-36.29\,{\rm{cm}},{y}_{0}=93.7\,{\rm{cm}},{z}_{0}=-42.79\,{\rm{cm}};\end{eqnarray*}$$

$$\begin{eqnarray*}&&{\rm{Quadrant}}\,4:u\in \left[\displaystyle \frac{3}{2}\pi ,2\pi \right),{\rm{where}}\,{x}_{0}=-36.29\,{\rm{cm}},{y}_{0}=93.7\,{\rm{cm}},{z}_{0}=42.79\,{\rm{cm}}.\end{eqnarray*}$$

USRBIN is adopted to record the absorbed dose for each simulation. USRBIN is a FLUKA command that scores the spatial distribution of absorbed dose. As recommended by FLUKA lectures, a relative dose error of 10% is a generally reliable quality for a simulation. To guarantee a non-primary dose statistical error of less than 10% around the target, the primary run number is set as $9.5\times {10}^{8}$ per simulation in this study.

Proton beam sampling process is shown in figure 9, which composed by source model, sampling random number PE (probability of proton energy) and energy. Gray source is target generated by SOBP as shown in figure 3. Yellow sources such as nozzle, extraction and injection are generated simply by concentration beam loss model. Green sources BTL and synchrotron ring are generated by uniform straight pipe beam loss model and torus source beam loss model correspondingly. Sampling random number PE is the summary of the previous energy loss weight. For example, p3 = 0.317 for target, p4 = p3 + 0.5% for 250 MeV on nozzle, p16 = p15 + 10.9% = 1 for injection. All these information should be coded into FLUKA SUBROUTINE source file, and before simulations the source file should be compiled properly.

Zoom In Zoom Out Reset image size

The typical time consumption for a ${10}^{6}$ primary run number is nearly 2 h on a standalone computer (Intel Core (TM) i7-8550U CPU @ 1.8 GHz and 16 GB memory). Approximately hundreds of hours will be required for one simulation, which is approximately one week. Thus, a parallel method should be employed to reduce the time consumption. The Amazon Web Service (AWS) High-Performance Computing (HPC) cluster is adopted for simulations (Wang et al 2022, 2023).

The non-primary dose distribution from the target is shown in figure 10. The left figure shows the $XOZ$ plane, and the right figure shows the $YOZ$ plane. The planned dose is 2 Gy at the isocenter ($x=215\,{\rm{cm}},y=155\,{\rm{cm}},z=555\,{\rm{cm}}$) in a water tank beside the nozzle. The target is a $10\,{\rm{cm}}\times 10\,{\rm{cm}}\times 10\,{\rm{cm}}$ region. The water tank also records the non-primary dose around the target. Two evaluation lines are set to investigate the 1D non-primary dose distributions in the $X$ and $Y$ directions for the target. The dose rapidly decreases as the distance from the isocenter increases.

Zoom In Zoom Out Reset image size

Figure 11 shows the 1D doses extracted from the 2D dose distributions of figure 10 along evaluation lines x and y. The flat top dose is the planned dose of 2 Gy. The dose exponentially decays outside of the 1 liter target region. The dose distribution is symmetric with respect to the isocenter.

Zoom In Zoom Out Reset image size

Table 2 shows the non-primary dose extracted from figure 11 in different IEC regions defined by figure 1. The doses in all these regions of interest are smaller than the threshold ratio of 0.5% or 0.1%.

Figure 12 shows the non-primary dose from the PT system in the XOZ plane (left) and YOZ plane (right). The Cartesian origin point is at the center of the synchrotron. From the left figure, there are 3 hot dose points for every quadrant as high-energy protons pass through the pipe shell between magnets. The magnets are made of pieces of iron, and they function as self-shielding to reduce the contribution to the non-primary dose in the target region. Hot lines emerge from the large free gaps between the quadrants. In particular, the hot line from between Q1 and Q2 threatens the patient area. Although the weights of the nozzle and BTL are much lower than that of the synchrotron, these sources are close to the patient area, and their contribution to the non-primary dose could be even higher (Wang et al 2020).

Zoom In Zoom Out Reset image size

From the right figure, extraction is a hot source due to its high weight, and the secondary particles irradiate the ground concrete. Then, the reflected and scattered particles could also contribute to the non-primary dose in the target area.

Figure 13 shows the 1D dose distributions extracted along evaluation line X (left) and line Y (right) for the PT system. The highest dose is ∼6.5E-6 Gy. This is a low level, and there is no symmetry regardless of the X or Y direction, as the sources are asymmetrically distributed. The hot line from between Q1 and Q2 contributes a bulge after X = 400 cm. Concrete reflection and scattering contribute little growth in the X < 50 cm and X > 400 cm regions.

Zoom In Zoom Out Reset image size

Table 3 shows the non-primary doses extracted from figure 13 in different IEC regions defined by figure 1. The doses in all these regions of interest are maintained at a very low level of around 1E-6 Gy and are smaller than the threshold ratio of 0.5% or 0.1%.

For direct comparison with the IEC requirements, non-primary doses from the subsources for tables 2 and 3 are summarized in tables 4 and 5. Table 4 shows the summarized non-primary doses for the region of 15–50 cm. The ratios of summarized non-primary doses to the planned dose of 2 Gy are all smaller than the IEC requirement. Moreover, the ratios are one to two orders of magnitude lower than 0.5%.

Table 5 shows the summarized non-primary doses for the region of 50–200 cm. The ratios of summarized non-primary doses to the planned dose of 2 Gy are all smaller than the IEC requirement. Moreover, the ratios are three to four orders of magnitude lower than 0.1%. Thus, the P-Cure PT system is quite safe for patients according to IEC requirements. It is a clean system, and there is no need for inner shielding concrete between the accelerator and target.

The weight of the non-primary dose from the target is more than 80% in the 15–50 cm region, whereas the PT system dominates up to 80% when the distance increases to 200 cm. This occurs because the sources of the PT system are distributed over a large scale, and the secondary particles travel all around in the shielding bunker.

Although it is a clean PT system, it can be further optimized by small efforts such as (1) adding local shielding between Q1 and Q2 to reduce the hot line and (2) painting a thin layer of neutron absorbing material on the concrete. This would be patient-friendly and helpful for reducing secondary malignancies.