Theoretical analysis of angular distribution of scattering in nozzle components using a response-function method for proton spot-scanning therapy

In spot-scanning proton therapy, highly precise beam control is required in the treatment nozzle such that the proton beam does not spread out during transportation by restraining the divergence of the beam angle and spot size, simultaneously. In order to evaluate the beam-broadening behaviour induced by passing through the various nozzle components, we have developed a new method to calculate the angular divergence profile of a proton beam in the nozzle. The angular divergence of the proton beam for each nozzle component is calculated by the Monte Carlo simulation code, Geant4, assuming that the initial beam has no divergence. The angular divergence profiles generated in the various nozzle components are then fitted by the analytic function formula with triple Gaussian distributions. The fitted profiles can be treated like analytic response functions and the angular divergence profile in the nozzle can be easily and systematically calculated by using a convolution theorem. The beam-broadening behaviour during transportation in the nozzle is carefully evaluated. The beam profiles are well-characterized by the proposed angular divergence analysis, i.e. triple Gaussian profile analysis. The primary Gaussian part of the beam profile is mainly generated by air and dose monitors with plate electrode components. The secondary and tertiary Gaussian parts are so-called wide-angle scattering and generated mainly by spot-position and profile monitors with metal window and wire components. The scattering of the nozzle component can be analysed using the proposed response function method for the angular distribution. Multiple convolved angular scattering can be determined from the response function of the individual nozzle components. The angular distribution from small to large angle regions can then be quantitatively evaluated by the proposed method. The method is quite simple and generalized, and is a straightforward way to understand the nozzle and component characteristics related to the beam-broadening behaviour.


Introduction
In recent years, the number of facilities adopting spot-scanning proton beam therapy has risen (PTCOG 2016). The spot-scanning method was developed to deliver highly conformal doses to the target tumours while sparing surrounding healthy tissues by using a well-controlled narrow proton beam (Chu et al 1993, Pedroni et al 1995, Gillin et al 2009. With regards to creating a high-quality therapeutic beam, the treatment nozzle design is one of the most important factors in manufacturing the beam transport system. Typically, several basic components are incorporated in the nozzle, such as X and Y-direction scanning magnets, a beam profile monitor (PRM), a dose monitor (DM) and a spot position monitor, among others.
Highly precise beam control is required in the treatment nozzle, so that the proton beam can maintain an accurate spot position and its reproducibility. Furthermore, it is essential that the beam should be prevented from spreading out during transportation by restraining the divergence of the beam angle and spot size, simultaneously. However, when the beams pass through the various nozzle components, a large angle dispersion caused by Theoretical analysis of angular distribution of scattering in nozzle components using a response-function method for proton spot-scanning therapy high Z materials may occur, as well as the dispersion of the small angle. It is important to identify how each nozzle component affects the beam broadening along the beam line by understanding the performance of the nozzle.
In previous studies, several attempts have been made to understand the beam-broadening behaviour in the nozzles and to evaluate their performance. Sawakuchi et al conducted a Monte Carlo (MC) simulation of a singlespot proton beam passing through a nozzle, and reported that vacuum windows and PRM are the main causes of large angular scattering (Sawakuchi et al 2010a). Gottschalk et al proposed categorizing beam-broadening into four parts, i.e. the core, halo, aura and spray, moving outwards from the beam's centre. They reported that the spray is mainly caused by the tungsten wires of the PRM (Gottschalk et al 2015). The beam profiles, after passing through nozzle components, were calculated by using model distribution functions, such as a summation of three Gaussian functions, also known as triple-Gaussian model function (TGF). It was possible to accurately describe the beam profile (Lin et al 2014). These studies are crucial to nozzle design and performance validation; however, they are not sufficient to evaluate the systematic contribution of each nozzle component on the beambroadening process.
We have developed a new analysis method to evaluate the beam-broadening process in angular distribution space. We define an angular distribution response function using model TGFs that represent the angular distribution response function calculated by the MC method for each component. In this approach, the angular divergence profile of any combination of the components in the nozzle can be calculated by the convolution of each component response function. Because the angular distribution response function of a component is independent from its spatial position in this method, it is fairly easy to calculate the overall response function of the nozzle, even when a nozzle component is added or modified; recalculations of the response function specific to the component and the overall convolution are straightforward. Without calculating the spatial beam profile, we can evaluate and compare the basic beam-broadening characteristics of the nozzle components and its combinations. The method is quite user-friendly, yet powerful and systematic, and it is straightforward to understand the beam-broadening mechanism, both qualitatively and quantitatively, through the analysis of component response functions and its convolutions.
In this paper, we introduce the response function method for evaluating angular scattering in nozzles and its application to a typical scanning nozzle. We then analyse the characteristics of each nozzle component related to beam broadening and its propagation to the following components using this method. Finally, their effects on the total nozzle performance are discussed.

Materials and method
2.1. Geometry of the scanning nozzle at Hokkaido University In this paper, the nozzle geometry for the proton beam therapy centre at Hokkaido University is used for the evaluation, which is a compact nozzle dedicated for the spot-scanning method manufactured by Hitachi, Ltd (Tokyo, Japan). The nozzle design is similar to the one at University of Texas MD Anderson Cancer Center (Smith et al 2009). This is composed of many components, including a vacuum window at the end of the proton-transport beamline, a PRM, a helium chamber, DM, a spot-position monitor (SPM) and a protection film, as illustrated in figure 1. It also has a pair of scanning magnets for the X-and Y-directions, none of which introduce materials in the beam path other than air or helium gas. Each component consists of subcomponents as follows. The vacuum window is made of thin titanium. The PRM have Kapton thin-film windows and copper and aluminium thin plates. They also have thin tungsten-wire electrodes, which are arranged in a three-layer structure with wire spacers. The helium chamber and the helium window consist of Kapton film, copper and aluminium to suppress proton-beam scattering caused by the air. The two main and sub DMs have plate electrodes composed of Kapton film and copper. The configuration of the SPM is almost the same as that of the PRM with five layers of tungsten wires at equal intervals. At the exit of the nozzle, there is a protection film consisting of Kapton, copper and aluminium.

MC simulation tool for proton transfer evaluation in the nozzle
The MC simulation toolkit Geant4.10.01.p01 (Agostinelli et al 2003, Collaboration 2015 was used for the proton transfer evaluation in the nozzle. The typical nozzle geometry of the proton beam therapy facility dedicated to the spot scanning was used for geometrical modelling. For the physics model, the standard parameters for particle therapy are used, as described in table 1. The parameter choices affect the calculation results, as discussed by various authors (Goudsmit and Saunderson 1940, Lewis 1950, Urban 2006. The result of MC simulation was validated by the measured data obtained at our facility; the range uncertainty was within 0.05 mm in percent depth dose and that of the beam size was within 0.2 mm in lateral profile. The relative dose errors for the in-water percentage-depth dose and lateral profiles were mostly within ±1%.

Evaluation of the angular distribution of the proton beam by MC simulation
At first, the angular distributions are calculated by the MC simulation code, Geant4, for all nozzle components in place, with a broad incident beam that is homogeneous with no angular divergence and no energy broadening. The calculations are performed from the beam-transport window or the entrance of the nozzle, all the way down to the isocentre, where events are recorded. The ten initial energies, Es, of the proton beam are selected from the energies used for spot-scanning between 70.2 MeV and 220 MeV. Using the Cartesian coordinate system, the beam-travelling direction is defined as the Z-axis, with the X and Y-axes running perpendicular to it. The angle θ = θ x , θ y , with respect to the Z-axis, is defined as (1), using the proton momentum p = p x , p y , p z , The events are recorded as a function of θ x , θ y and the angular distributions are obtained as o MC θ x , θ y . This function is normalized by the initial proton number, N 0 , and the bin widths, ∆θ x and ∆θ y . The final number of protons recorded by the detector, N, is obtained by the summation: Note that the angular distribution is scored with a detector with a sufficiently large solid angle, Ω. ,or N to satisfy N ≈ N 0 ; additionally, o MC θ x , θ y can be regarded as a probability-density distribution. It is also assumed that the energy lost when passing through the nozzle component is negligible.
For further procedures, it is better to use 1D histograms; therefore, the angular distribution is checked for isotropy and all 2D data are reduced to one dimension. The 2D histogram is summed in the θ x and θ y directions, which are 2D data distributions projected into the X and Y directions, respectively. After confirming that the agreement is reasonable, the 1D angular distribution for each component i

Fitting procedure of the angular distribution of the proton beam by the TGF
The 1D angular-distribution functions obtained by the MC simulations are fitted by the TGF described above for both the individual components and the complete nozzle. The TGF as a function of θ is expressed as (3) for the component i.: The first term is called primary Gaussian, second one secondary and third one tertiary, and denoted by j in the followings such as g i,j , A i,j and σ i,j . The fitting is performed in the following manner; (i) Fit with only the first term, g fit i,1 (θ) and obtain A i,1 , σ i,1 ; (ii) Fit with only the first two terms, g fit i,1 (θ) + g fit i,2 (θ) , with the initial conditions set to A i,1 , σ i,1 for the primary Gaussian, but allow the parameters to vary within ±10% of the initial parameters and obtain A i,1 , σ i,1 , A i,2 and σ i,2 . Check whether A i,1 and σ i,1 stay well inside the boundaries; (iii) Fit with all the three Gaussians, g fit i,1 (θ) + g fit i,2 (θ) + g fit i,3 (θ) , with the initial conditions set to A i,1 , σ i,1 , A i,2 and σ i,2 and, once again, set the boundaries to ±10% of the initial values and obtain all six parameters. This procedure is repeated for all nozzle components, including the air contribution. Eventually, g fit i (θ) functions for the individual components and the g fit All (θ) function for the complete nozzle are obtained. One of our fitting results for the PRM is shown in figure 2 as an example.

Energy dependence of the angular distribution of the proton beam
It was perceived that there might be a universal behaviour to the energy dependence of the beam profiles. Therefore, it was verified whether or not this was the case. If there is a universal behaviour, we can predict the beam profiles at any energy simply by calculating the profile at any given energy and scaling it to others.
The angular distribution, o MC i (θ), obtained by an MC simulation, as well as the fitted TGF, g fit i (θ), and the obtained parameters, σ i,j , have energy dependences. It is determined whether the σ i,j parameters are proportional to pv 70.2 MeV /pv, where p is the momentum of protons, v is their velocity and v 70.2 MeV is the minimum constant velocity in our calculation. It is also verified whether the amplitudes, A i,j parameters, obtained by the fitting have no energy dependences that are proportional to the number of protons belonging to one of the three Gaussians.

Convolution of the angular-distribution-response functions
Having obtained all the response functions, g fit i (θ), for the nozzle components, it is fairly straightforward to calculate the final angular distribution in the analytical formula by the convolution method for the beam passing through all the components of the nozzle, as shown in (4): To calculate the convolution, the Fourier transform of each g fit i (θ) is first calculated to obtain G fit i (k). The convolution, g conv All (θ), is expressed as the inverse Fourier transform of the product of all G fit i (k), as shown in (6): where F −1 denotes the inverse Fourier transform of the function in the parentheses. We determined whether the obtained g conv All (θ) is consistent with the corresponding MC simulation result, o MC All (θ), with all components.

Analysis of the primary Gaussian part of the angular distribution of the proton beam
The contributions from the primary-Gaussian part of several selected components (SC) to the overall TGF profile are analysed by convolution calculation. The convolution of the selected primary-Gaussian parts is calculated as We check whether g conv SC, 1 (θ) can reproduce main part of the final angular distribution.

Analysis of the secondary-and tertiary-Gaussian parts of the angular distribution of the proton
The contributions from the secondary-Gaussian part of several SC to the overall TGF profile can be analysed by the following method: first, the convolution of the SC with the secondary-Gaussian part missing is calculated; it is then convolved with the full TGF for the remaining components to obtain g conv SC,2 (θ), which is shown in (8): Finally, the convolved function, g conv SC,2 (θ), is subtracted from g conv All (θ) to obtain the difference, g conv SC,2 (θ) : Before conducting the above calculation, in order to choose the SC, the secondary-Gaussian part for each component is evaluated by using the same method described above, i.e. only one component is chosen as the SC and g conv SC,2 (θ) is calculated for all the component. A few components that have large secondary parts are then selected.
Likewise, for the tertiary-Gaussian part, we use (10): and the difference:

Angular distribution of the proton beam calculated using MC simulation
As described in section 2.3, the 2D angular distribution, o MC All θ x , θ y , at the isocentre for all nozzle components, including the air contribution, was calculated by an MC simulation. Using the same framework, o MC i θ x , θ y was also calculated for each component. One such calculated result for the PRM is shown in figure 3(a).
The anisotropy in the 2D angular distribution was checked by comparing the 1D data, as a function of x, integrated over the Y-direction, and as a function of y, integrated over the X-direction. As shown in figure 3(b), the difference between the integrated intensities for the θ x and θ y directions for the PRM was very small. This indicates that there was almost no anisotropy down to the order of 10 −4 . It will be discussed in more details in the discussion session. Therefore, in the following analysis, only the 1D data obtained in such a way were used.
1D angular distributions, with all nozzle components in place were calculated for each component using the MC simulation code. Selected results for 70.2 MeV, 99.9 MeV, 140.8 MeV and 220.0 MeV are shown in figure 4. As shown in the figures, when the energy increases, the angular distributions become narrower as expected. In addition, the major components that contributed to the broadening of small-angle regions near the centre were air scattering and the DM. For the wide-angle-tail parts, at a scale of approximately 10 −2 , the main components that contributed to them were the PRM, air scattering, DM and SPM. 3.2. The results of fitting of the angular distribution of the proton beam by the TGF As described in section 2.4, the results have been fitted with TGFs using a least-squares method. The results of the fitting for an energy of 140.8 MeV are shown in figure 5 for (a) all nozzle components in place, (b) air, (c) the PRM and (d) the DM. Each Gaussian contribution is indicated by coloured curves, together with the MC simulation result indicated by triangular markers. The fitting precisions for all three Gaussian parts are quite good, as can be seen from the figures. We assessed the precision of the fitting using the maximum distance, D, of the Kolmogorov-Smirnov test (Bellinzona et al 2015). It is evident that, even in the worst case, D, is less than 0.005, which is satisfactory.
In the case of the air, the contributions of the primary and secondary Gaussian parts are relatively large; in the case of the PRM, it shows a small broadening effect in the primary part, but the contributions to the secondary and tertiary Gaussians are relatively large; in the case of the DM, it shows modest broadening in the primary part, but the contribution to the secondary part is as large as that for the PRM. To check the fitting uncertainty for PRM furthermore, we calculated the ratio g fit PRM (θ) /o MC PRM (θ). It was within 1.0 ± 0.2 in the angular range of θ < 10 mrad. Even in the tertiary-Gaussian region, which is θ > 10 mrad, the averaged ratio was also around 1.0 ± 0.2, but because of the poor statistics of the result of MC simulation in the region, it was as high as 1.0 ± 0.4.
All the parameters thus obtained (

Result of energy-dependence analysis of the angular distribution of the proton beam
As described in section 2.5, the energy dependences of the A i,j parameters have been checked first. The graphs of A i,j parameters for the primary, secondary and tertiary parts are shown in figure 6, and it is apparent that they have low energy dependences. Secondly, the σ i,j parameters, normalized by pv 70.2 MeV /pv are shown in figure 7, also show little energy dependences, implying that our above assumption of the energy dependence is a reasonable one. The reason why we have chosen pv 70.2 MeV /pv as the normalization will be discussed in the following session. There are discrepancies in the normalized σ i,3 values for all the nozzle component and air cases. These are understandable, because the tertiary-Gaussian part for air is relatively small compared with the secondary part and when compared with the cases with other components, resulting in some errors in the lower energy cases. The case with all the nozzle components was also affected by the air contribution. Since there is very little energy dependence in the A i,j values and the normalized σ i,j values, or the normalized overall profiles, the results for only one energy are shown in the following figures.

Results obtained by the convolution of the angular-distribution-response functions
Following the procedure described in section 2.6, we checked whether an analytically convolved TGF, g conv All (θ) , can well describe the MC simulation results when all the nozzle components are in place. The former cases are indicated by the red curves in figure 8, whereas the latter results are denoted by markers for the four energies: 70.2, 99.9, 140.8 and 220.0 MeV. The convolved functions overlap with the corresponding results quite well at all energies.

Results of analysis of the primary Gaussian part of the angular distributions of proton beam
Following the method of section 2.7, the angular distributions of the proton beam have been analysed in terms of the primary, secondary and tertiary parts for all nozzle components in place and for each component.
For the primary part, the peak normalized g fit i,1 (θ) functions for each individual component have been obtained for an incident-proton-beam energy of 140.8 MeV, which is shown in figure 9(a), together with the simulation results with all the nozzle components in place, o MC All (θ). The functions having a hat, '^', in the following section are normalized by the peak. The air scattering shows the largest contribution to the primary part and slightly smaller contributions are obtained from the two DMs. The curves for the two DMs clearly coincide.
Based on the above results, air scattering and these two DMs were chosen as the SC in (7). Convolution of the primary-Gaussian functions for the SC that have large primary Gaussians parts, were calculated and normalized by the peak. The convolved function for the SC, g conv (Air,DM1,DM2),1 (θ), thus obtained is shown by a dashed line in figure 9(b), which roughly reproduces, albeit slightly smaller, the major part of the primary Gaussian contrib utions. The latter, g conv All,1 (θ), is the result of the convolution obtained when all the nozzle components are in place. It is a little bit smaller than that obtained by the MC simulation, o MC All (θ), which can be attributed to the contrib utions of the secondary and tertiary parts.

Results of the analysis of the secondary-and tertiary-Gaussian parts of the angular distribution of the proton beam
Following the method of section 2.8, the wide-angle parts of the angular distribution of the proton beam passing through all nozzle components were analysed. For the secondary parts, the SC were the PRM, two DMs and the air contribution, which have large secondary parts, in equation (8), while the PRM, DM, SPM and air were chosen as tertiary parts in equation (10).
The effects of the secondary and tertiary parts could be obtained by subtracting g conv SC,2 (θ) or g conv SC,3 (θ) from g conv all (θ). Note that g conv SC,2 (θ) and g conv SC,3 (θ) are the products of the response functions for the SC with missing secondary or tertiary contributions and the response functions of all other components with Gaussian contributions.

Secondary part
For the secondary parts, g conv SC,2 (θ) was calculated for the SC for an incident-proton-beam energy of 140.8 MeV, as shown in figure 10(a), together with the simulation results with all nozzle components in place, o MC All (θ).  Table 2. The parameters obtained by the triple-Gaussian fitting for 70.2 MeV: the integrated intensity for the jth Gaussian part ( A i,j ), the standard deviation of the corresponding Gaussian part (σ i,j ) and the Kolmogorov-Smirnov-test parameter (D). The PRM shows the largest contribution to the secondary Gaussian part, with a peak intensity of around 10 −2 . Slightly smaller contributions from the two DMs and air have been obtained. The overall secondary Gaussian contributions from the above four components are shown as dashed lines in figure 10(b), while that obtained with all the nozzle components in place are shown by a solid line. As can be seen from the figures, the secondary parts, g conv All,2 (θ) , mostly consist of the above four components, g conv (PRM,Air,DM1,DM2),2 (θ), in the intensity region between 10 −2 and 10 −3 . Table 3. The parameters obtained by the triple-Gaussian fitting for 99.9 MeV: the integrated intensity for the jth Gaussian part ( A i,j ), the standard deviation of the corresponding Gaussian part (σ i,j ) and the Kolmogorov-Smirnov-test parameter (D).

Tertiary part
For the tertiary parts, g conv SC,3 (θ), was calculated for the SC for an incident-proton-beam energy of 140.8 MeV, and is shown in figure 11(a) together with the simulation results with all nozzle components in place, o MC All (θ). The PRM showed the largest contribution to the tertiary Gaussian, with a peak intensity at around or below 10 −3 . The contributions of the SPM, DM1 and DM2 are approximately two to three times smaller than that of the    PRM, while that of air is about 3 to 4 times smaller. The overall contribution from all the five components are shown as a dashed line in figure 11(b) and the contribution for all the nozzle components in place is indicated by the red curve. As can be seen from the figures, the tertiary Gaussians, g conv All,3 (θ), mostly consist of the above five components, g conv (PRM,DM1,DM2,SPM,Air),3 (θ), in the intensity region below about 10 −3 . The red curve shows the convolution result for all the nozzle components in place, g conv All,1 (θ), the dashed line shows the one for major components, two DMs and air scattering contributions, g conv (Air,DM1,DM2),1 (θ), and the triangular markers show the MC simulation result for all the components in place, o MC All (θ). All the curves are normalized by their peak.

Discussion
Using the proposed method, the angular distribution of a proton beam passing through the nozzle can be calculated by multiple convolution of each response function of the component, as expressed by equation (6). The response function for the nozzle under consideration is defined as an individual combination of triple-Gaussian functions. Once we use the response function method, it is fairly straightforward to estimate the effect of changing or modifying one of the components in the nozzle. We can simply calculate the response function of the component and convolve it with those of the remaining components. Furthermore, the angular distribution passing through the specific part of the nozzle can be evaluated by convolving the response functions for the SC. The proposed method, therefore, can provide quite a powerful yet intuitive method for understanding the broadening behaviour and to better design or modify treatment nozzles. The response functions of the components are well represented by the triple-Gaussian functions. We can assume that the primary part of the function can be mainly attributed to multiple-Coulomb scattering (MCS) effects. The secondary and tertiary parts of the function can be attributed to wide-angle scattering, including wide-angle Rutherford scattering, although there is not a clear one-to-one correspondence between the distribution and the scattering effects. We can recognize the characteristics of each component's effect on the beambroadening behaviour from its response function and evaluate how it happens. If a component consists of many light nuclei, such as air or polymer windows, the main contribution is reflected in the primary Gaussian part. On the other hand, if it consists of heavy nuclei, such as metal windows, electrodes, or wires, the main contributions give rise to the secondary and tertiary Gaussian parts. Therefore, from response-function analysis, we can obtain much information regarding the components and how they affect the behaviour of the nozzle under consideration.
There is another reason for the wide-angle-scattering Gaussian part, which is caused by wire-scattering in PRM. A proton beam that goes through windows of the PRM produces a sharp Gaussian distribution caused by the MCS and wide-angle scattering part. In addition to this, a small number of protons that actually hit the wires of the PRM produce a broader Gaussian distribution because of the wire material, which is tungsten. Since the rate of hitting the wires is of the order of 10 −2 , it overlaps with the sharp Gaussian distribution caused by the window materials, and the wide-angle component, at the level of 10 −2 , caused by the wire. The main parts of both components should be caused by the MCS rather than the wide-angle Rutherford scattering. The latter contributions should contribute to much wider angled scattering parts.
Moreover, we need to consider that the wide-angle Rutherford scattering, which occurs scarcely, should be convolved with the MCS in the other materials in the component, but also with the same material that caused the wide-angle scattering. Therefore, the wide-angle parts in the vicinity of the primary Gaussian part should be expressed by a convolution of the wide-angle Rutherford scattering function and the Gaussian functions caused by the MCS.
The energy dependence analysis conducted in section 3.3 shows reasonable results that can be explained by the above discussions. The proton beam broadening induced by MCS is approximated by the Highland approximation (Highland 1975  We would also like to discuss about anisotropy in the wide-angle or in very low intensity region. In contrast to our simulation results, Lin et al (2014) reported anisotropic beam divergence in the 2D intensity distribution in the 0.01% intensity region. Note, however, that MC simulation modelling is quite different for our case and for Lin's case. In the former case, simulation starts from isotropic beam divergence at the nozzle entrance and in the latter case, starting particle transport from the phantom surface particle transport from the phantom surface allowing anisotropy in beam divergence. Also, in our case the profiles are normalized to the area, which is equivalent to the number of protons, whereas in Lin's case, intensity is normalize to its peak, results in 0.01% in Lin's results (Lin et al 2014) corresponds to below 10 −4 in our case. Therefore, it is reasonable that we observed no anisotropy in our simulation. There may be inaccuracy of multiple Coulomb scattering modelling in the simulation tools in the vicinity of and below 10 −4 region in our case, as have been reported by Grevillot et al (2010) and Sawakuchi et al (2010b), which should not affect our calculation.
(12) Here, θ 0 is the width of the angular distribution assuming Gaussian distribution, x/X is the thickness of the medium in a radiation length unit. As can be seen from the equation, it has the energy dependence of 1/pv. Therefore, it is reasonable to assume that the energy dependence of the primary Gaussian, which is a result of MCS, can be expressed by 1/pv. Our analysis clearly shows that this assumption holds quite well for the given energy range. It is, however, quite interesting to see that it also holds for the secondary and tertiary Gaussian parts. It is reasonable that PRMs, which have wires as components, show a 1/pv dependence because the main part of their angular profile is caused by MCS. For the other components, the wide-angle parts are also well characterized by the 1/pv dependence. This is a kind of empirical result produced by the MCS. Therefore, we do not know of a clear reason for this, but as wide-angle Rutherford scattering should always be convolved with MCS that has the energy dependence of 1/pv, there is a possibility of the overall wide-angle scattering having nearly a 1/pv dependence, at least for the vicinity of the primary Gaussian part.

Conclusion
We have developed a new analysis method that can calculate angular divergence profiles in the treatment nozzle of a spot-scanning proton beam therapy system using an angular distribution response function defined for each nozzle component. The angular divergence profile of any combination of the components in the nozzle can be calculated by the convolution of each component response function. As the response angular profile of a component is independent of its spatial position in this method, we can easily calculate the overall response function of the nozzle, even when a nozzle component is added or modified. Following the analysis, we can identify the components that broadened the beam most or contributed the most to wide-angle scattering. We can effectively use such information while designing a nozzle, adding a new component or when performing quality   assurance. The method is quite simple and generalized, and is a straightforward way to understand nozzle and component characteristics related to the beam-broadening behaviour.
The angular divergence of the proton beam for each nozzle component is calculated by the MC simulation code, Geant4, assuming that the initial beam has no divergence. Subsequently, the angular divergence profiles generated in the various nozzle components are fitted by the analytic function formula with TGF distributions. The fitted profiles can be treated like analytic response functions and the angular divergence profile in the nozzle can be easily and systematically calculated by using a Fourier transformation and its convolution theorem. The beam-broadening behaviour during transportation in the nozzle is carefully evaluated.
We have analysed the angular divergence profile of the proton beam in a typical scanning nozzle. The divergence profiles through the nozzle are calculated by the convolution of those response functions and are wellcharacterized by the proposed angular divergence analysis. The primary Gaussian part of the beam profile is mainly generated by air and DM with plate electrode components. The secondary and tertiary Gaussian parts correspond to so-called wide-angle scattering and generated mainly by SPMs and PRMs with wire electrode components.
We conclude that the proton beam scattered by the nozzle component can be analysed using the proposed response function method for the angular distribution. Multiple convolved angular scattering can be determined from the response function of the individual nozzle components. Subsequently, the angular distribution from small to large angle regions can be quantitatively evaluated by the proposed method.