Addendum: Measurement of charged particle yields from PMMA irradiated by a 220 MeV/u ¹²C beam

The characterization of the secondary charged particles production induced by ion beams of interest for particle therapy (PT) was the subject of different experimental campaigns in the past (Agodi et al 2012, Henriquet et al 2012, Gwosch et al 2013), due to its fundamental importance in assessing the feasibility of beam range monitoring in PT exploiting the charged secondary emission.

The analysis of the secondary charged particles produced in the interaction of a 220 MeV/u fully stripped $\newcommand{\car}{^{12}{{\rm C}}} \car$ ion beam impinging on a $20 \times 5 \times 5$ cm³ PMMA has been the subject of a dedicated publication in 2014 (Piersanti et al 2014). The published data were collected in 2012 at the GSI facility in Darmstadt (Germany).

During an experimental campaign performed in February 2014 at the Heidelberg Ion-Beam Therapy Center (HIT, Heidelberg, Germany), an heterogeneous light response of the LYSO crystals used as secondary charged particles detector was observed. This evidence triggered a re-analysis of the published data to evaluate the impact on the measurements of the proton (p), deuteron (d) and triton (t) populations in the reconstructed track sample.

In this paper, the inhomogeneity of the crystal light response is presented in section 2, as well as its impact on the Particle IDentification (PID) results. The new emission velocity ($\beta = v/c$) spectra and particle yields obtained with the improved analysis are presented in sections 3 and 4, respectively. While re-analyzing the data collected at GSI, we have also refined the calculation of the detector and geometrical efficiencies applied to the reconstructed track sample, explicitly taking into account the energy dependence of the charged fragments. The new efficiencies calculation is discussed in section 4.1.

As the total flux of charged particles ($p + d + t$) produced in the interaction of a PT ion beam with a PMMA target is of the same order of magnitude than the one published in Piersanti et al (2014), the implications for PT applications discussed in the previous paper still hold and are not discussed nor presented again hereafter.

Figure 1 left shows a sketch of the experimental setup (not to scale) of the 2012 data taking (GSI data in the following). A 1.1 mm thick plastic scintillator (start counter, SC), placed before the PMMA target, was used to count the primary ions. A drift chamber (DC) and a LYSO detector were placed both at $\theta = 60^\circ$ or $90^\circ$ with respect to the beam impinging direction, at a distance of $41~{\rm cm}$ and $65.8~{\rm cm}$, respectively, from the PMMA center. The DC, composed by 12 wire layers (alternated X and Y views) with 3 cells per layer and filled with a Ar/CO2 (80/20) gas mixture, was used to select and track secondary charged particles, exploiting the distribution of the number of fired cells (N_hits). The LYSO detector was a $2 \times 2$ matrix of LYSO crystals, $1.5 \times 1.5 \times 12.0$ cm³ each. Each crystal was wrapped in Tyvek^®16 by 5 sides and the four crystals (shown in figure 1(right)) were held together by black tape and read by a single photomultiplier tube (PMT). The response signal from the LYSO detector was used in time coincidence with the SC to form the experiment trigger and for the PID.

Zoom In Zoom Out Reset image size

A $10.0 \times 10.0 \times 0.2$ ${\rm cm}^{3}$ plastic scintillator (VETO in figure 1(left)) was placed between the target and the DC to stop electrons with energies below 0.55 MeV. The 220 MeV ${\rm u}^{-1}$ fully stripped $\newcommand{\car}{^{12}{{\rm C}}} \car$ beam was Gaussian-shaped, with a standard deviation in the transverse plane of  ∼1 cm, measured with  ∼$1\%$ relative uncertainty by the GSI beam monitor chambers. The beam rate was ranging from hundreds of kHz to  ∼1 MHz, with a temporal spectrum of the beam typical of a synchrotron machine. The beam was pulsed, with a triangular shape spill $5{\rm{-}}10$ s long and a peak rate occurring $2{\rm{-}}3$ s after the spill start. The beam spill was interleaved with a no beam time window of $2{\rm{-}}3$ s. More details on the experimental setup and on the Data Acquisition (DAQ) can be found in Piersanti et al (2014).

In 2014 the LYSO crystals used for the GSI experiment were calibrated at the HIT center using protons of several energies (50 MeV $\div$ 210 MeV, HIT data in the following). The calibration data allowed to measure the quenching effect in the LYSO crystals and the saturation of the LYSO PMT taking place at energies of interest for PT applications. The response of the different crystals appeared to be significantly different, as shown in figure 2 for protons of 70 MeV (a) and 91 MeV (b), and was characterized by the presence of three clearly distinguishable peaks (two crystals happened to have nearly the same light response).

Zoom In Zoom Out Reset image size

Therefore, a re-analysis of the GSI data has been performed in order to include the effect of the crystal inhomogeneity when looking at the energy deposited in the LYSO detector (E_LYSO) as a function of the Time of Flight (ToF) of charged secondary particles, as defined in Piersanti et al (2014). The charged particles signal selection exploiting the DC detector, asking for $N_{\rm hits} \geqslant 8$, has not been changed.

Figure 3 left shows the E_LYSO versus ToF distribution published in Piersanti et al (2014), while figure 3 right highlights the contribution of each crystal using different colors. The crystal regions (see figure 1) were identified for each detected fragment tracked using the DC.

Zoom In Zoom Out Reset image size

The limited knowledge about the actual crystal (and hence its energy response) hit by each fragment, due to the DC tracking resolution when extrapolating the DC reconstructed tracks forward to the LYSO detector entrance face, is taken into account when assigning the systematic uncertainty associated to the PID efficiency and the definition of the selection bands (see section 4.3).

In the GSI data taking, the observed LYSO4 energy response is significantly different with respect to the other ones, possibly due to the optical coupling with the PMT. Instead, the (E_LYSO, ToF) regions of the other three crystals are approximately overlapping. The handling of the crystals, their transportation and re-assembly between the two different experiments can very well account for the different results (the three distinct peaks observed in the HIT data appear as two main populations in the GSI data).

Protons and deuterons impinging on the LYSO4 crystal where thus identified, using the PID selection bands defined in the published work, as deuterons and tritons.

A new PID procedure has been defined to evaluate the charged particle yield and velocity distribution for each particle species taking into account the specific response of each crystal.

Fragments collected by the LYSO4 crystal, the one that showed a significantly different response from the other ones, were excluded from the data sample as their contribution to the decrease of the total statistical uncertainty would have been balanced by the increased systematic uncertainty related to the Monte Carlo (MC) simulation description of the peculiar crystal response (see section 4).

Electrons, that can be easily identified as the cluster of points centered around  ∼2.0 ns in the (E_LYSO, ToF) plane, were rejected by requesting ToF $>$3.5 ns and ToF $>$3.0 ns for the $90^\circ$ and $60^\circ$ angular configuration data samples, respectively. The correct geometrical acceptance, computed excluding the LYSO4 data, and the systematic error coming from the electrons rejection cut were taken into account when evaluating the yield, as described in section 4.

Figure 4 shows that only few tritons are present. Therefore, we decided to neglect such contribution and to perform only the study of protons and deuterons fluxes and velocity spectra. To implement the PID selection, a full MC simulation of 220 MeV/u $\newcommand{\car}{^{12}{{\rm C}}} \car$ ion pencil beam impinging on the PMMA target was generated using the FLUKA code (Ferrari et al 2005, Battistoni et al 2007, 2015), taking into account a detailed description of the primary beam characteristics, of the experimental geometry and the detector resolution, for both $90^\circ$ and $60^\circ$ angular configurations, as it was done in the previous analysis.

Zoom In Zoom Out Reset image size

The MC LYSO energy response was tuned over the data distribution in order to achieve a good data-MC agreement. The simulation was used to identify the charged particle species and compute the selection efficiencies (see section 4.1). Figure 5 left shows the E_LYSO versus ToF measured distribution (Data, black dots), superimposed to the MC distribution of Protons (red dots), Deuterons (blue dots) and Tritons (green dots) for the $90^\circ$ angular configuration.

Zoom In Zoom Out Reset image size

To draw the PID selection functions, shown as solid lines in figure 5(right), the MC E_LYSO versus ToF distribution has been divided in ToF slices of 0.5 ns width. Figure 6 shows the E_LYSO distribution in a ToF bin ($8.0{\rm{-}}8.5$ ns) for the $90^\circ$ case, fitted with a double Gaussian function: the first peak (mean value $\mu 1$) is related to the proton fragments and the second peak (mean value $\mu 2$) to the deuteron fragments.

Zoom In Zoom Out Reset image size

The functions used to separate the proton, deuteron and triton populations (shown in figure 5(right)) were obtained from the fitted central values of the proton and deuteron peaks as a function of ToF and by moving 5 and 4 standard deviations, respectively, from the central value ($\mu1 + 5\sigma1$, $\mu2 + 4\sigma2$). 'Hard' and 'soft' selections were defined in order to compute a systematic uncertainty on the PID selection, by setting the separation 4 (3) and 6 (5) σ from the p (d) central value respectively (see section 4.3).

With the new PID selection, the β spectra have been re-obtained following the same procedure described in Piersanti et al (2014). Figure 7 shows the new β spectra (red triangles) superimposed to the old published ones (black circles) for p (left) and d (right) fragments detected at $90^\circ$ (top) and $60^\circ$ (bottom), normalized to the number of primary carbon ions.

Zoom In Zoom Out Reset image size

As expected the measured spectra have the same shape, but different normalization with a depletion of the deuteron population and an increase of the proton one.

The measurement of the charged secondary particle yield depends on the fragment detected kinetic energy ($E_{\rm kin}^{\rm ToF}$) and reconstructed emission position (x_PMMA) through the calculation of the detection efficiency. $E_{\rm kin}^{\rm ToF}$ was evaluated from the fragment ToF and x_PMMA was evaluated from the backward extrapolation of the secondary track reconstructed by the DC to the (x,y) plane containing the beam line inside the PMMA. While performing a reanalysis of the data, the explicit dependence of the efficiency on the detected kinetic energy $E_{\rm kin}^{\rm ToF}$ and x_PMMA was accounted for, unlike what was done in the previous publication where an efficiency averaged over the full energy production range was used. The current strategy implies a correction which significantly affects the determination of the production yields, especially at low energies where the probability that a fragment is absorbed in the PMMA is higher. The yield of charged secondary particles was computed using the following equation:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Phi_{\rm ch} = \frac{1}{4\pi \times N_{\rm C} \times {\rm corr}_{\rm GT} \times \epsilon_{\tau} \times \epsilon_{\rm sel}^{\rm ch} } \sum_{i,j} \left(\frac{N_{\rm ch}}{\epsilon(x_{\rm PMMA}, E_{\rm kin}^{\rm ToF})} \right)_{ij} \label{eq:flux} \nonumber \end{align} \tag{ 1 }$$

where N_C is the number of primary carbon ions as defined in Piersanti et al (2014) and $\newcommand{\e}{{\rm e}} \epsilon_{\tau}$ is the DAQ dead time efficiency. N_C was corrected by the factor corr_GT that takes into account the underestimation of the incoming ion total number: multiple carbon ions interacting with the SC detector within the discrimination time of the SC signal (80 ns) cannot be detected. The corresponding value of the correction at 90 and 60 degrees are corr$_{\rm GT}(90^\circ) = 1.12 \pm 0.07$ and corr$_{\rm GT}(60^\circ) = 1.37 \pm 0.19$. More details on the corr_GT evaluation can be found in Mattei et al (2015). The symbol $(N_{\rm ch})_{ij}$ represents the number of charged secondary particles (identified as p or d) in each $E_{\rm kin}^{\rm ToF}$–x_PMMA (i–j) bin. Figure 8 shows the measured $E_{\rm kin}^{\rm ToF}$ versus x_PMMA distribution of proton selected particles detected in the $90^\circ$ (left) and $60^\circ$ (right) angular configuration. The $\newcommand{\e}{{\rm e}} \epsilon(x_{\rm PMMA}, E_{\rm kin}^{\rm ToF})_{ij}$ is accounting for the DC signal selection and the geometrical and detector efficiencies and $\newcommand{\e}{{\rm e}} \epsilon_{\rm sel}^{\rm ch}$ is the PID selection efficiency. The following (section 4.1) describes the evaluation of each efficiency term.

Zoom In Zoom Out Reset image size

4.1. Efficiencies

To evaluate the DAQ dead time efficiency $\newcommand{\e}{{\rm e}} \epsilon_{\tau}$, the DAQ dead time τ was computed as a function of the instantaneous event by event rate, taking into account the event size measured during the data acquisition, as detailed in Mattei et al (2015). In the data taking in the $90^\circ$ configuration, $\newcommand{\e}{{\rm e}} \epsilon_{\tau}$ is  ∼$80\%$ while at $60^\circ$ is  ∼$70\%$, depending on the DAQ event rate. The detection and geometrical efficiency $\newcommand{\e}{{\rm e}} \epsilon(x_{\rm PMMA}, E_{\rm kin}^{\rm ToF})$ was computed using a dedicated MC simulation. Since $\newcommand{\e}{{\rm e}} \epsilon(x_{\rm PMMA}, E_{\rm kin}^{\rm ToF})$ relies only on the capability of the MC software to simulate the secondary fragments transport (and not on the fragments production nuclear models), a flat simulation was used in which secondary protons and deuterons were generated inside the PMMA target with uniform energy and production position spectra. The spatial distribution of secondary fragments in the z–y plane was modeled with a Gaussian function centered along the PMMA central axis with a 1 cm RMS width equal to the primary ion beam width. The primary ion beam range is $8.90 \pm 0.03$ cm, determined from a MC simulation (GEANT) (see Piersanti et al (2014)) and therefore the BP is expected to be at  ∼$2.45$ cm. The production along the PMMA axis (x axis) was continuous in the range $x_{\rm PMMA}= [-10, 15]$ cm, larger than the PMMA dimensions (going from $-6.45$ cm to $13.55$ cm) in order to fully account for the whole spatial distribution of the reconstructed tracks, while the kinetic energy distribution was continuous in the 10–250 (20–250) MeV range for the p (d) samples.

The efficiency matrix $\newcommand{\e}{{\rm e}} \epsilon(x_{\rm PMMA}, E_{\rm kin}^{\rm ToF})$ for protons detected in the $90^\circ$ ($60^\circ$) angular configuration is shown in figure 9left (right) and is computed as the MC detected $E_{\rm kin}^{\rm ToF}$ versus x_PMMA distribution of the simulated protons normalized to the number of secondary protons produced in each bin.

Zoom In Zoom Out Reset image size

Given the energy loss of secondary fragments inside the PMMA, in order to take into account and correct for range boundary effects in the $\Phi_{\rm ch}$ evaluation, the efficiency values in the edge bins in energy (10–20 MeV and 220–230 MeV) were extrapolated to the lower and higher energy values. Few charged secondary tracks are reconstructed outside the PMMA region due to the multiple scattering effect and reconstruction algorithm errors. Therefore, to correct for x_PMMA range boundary effects on $\Phi_{\rm ch}$, the efficiency values in the edge bins ($-7.0{\rm{-}}-6.5$ cm and 13.5–14.0 cm) were extrapolated to the lower and higher x_PMMA values.

The selection efficiency $\newcommand{\e}{{\rm e}} \epsilon_{\rm sel}^{\rm ch}$ was computed from a full MC simulation by measuring the probability that a fragment of type u could be measured in the region v ($u, v = p, d$ for protons and deuterons respectively) and taking the diagonal values of the resulting mixing matrix: $\newcommand{\e}{{\rm e}} \epsilon_{\rm mix}^{uv} = N_{uv}/N_u$. $N_{uv}$ is the number of fragments u measured in the region v normalized to the total number of fragments N_u in the (E_LYSO–ToF) plane.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \varepsilon_{\rm mix} = \left(\begin{array}{@{}ll@{}} \varepsilon^{\,pp} & \varepsilon^{\,pd} \nonumber \\ \varepsilon^{dp} & \varepsilon^{dd} \end{array} \right). \label{mixmat} \nonumber \end{align} \tag{ 2 }$$

The u and v regions are defined as previously described in section 2. The mixing matrix defined in equation (2) is used to estimate the selection efficiency (diagonal elements) and to correct the yield estimate for the wrong PID assignment (cross feed of different PID hypothesis) whenever the off diagonal elements are statistically different from zero. The matrix has been computed both for the $90^\circ$ and $60^\circ$ angular configurations:

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \varepsilon_{\rm mix}^{90} = \left(\begin{array}{@{}ll@{}} 99.65 \pm 1.37 \% & 0.33 \pm 0.06 \% \nonumber \\ 5.73 \pm 0.59 \% & 93.31 \pm 3.20 \% \end{array} \right)\nonumber \\ \varepsilon_{\rm mix}^{60} = \left(\begin{array}{@{}ll@{}} 99.72 \pm 0.68 \% & 0.28 \pm 0.03 \% \nonumber \\ 7.90 \pm 0.46 \% & 90.00 \pm 2.04 \% \end{array} \right). \nonumber \end{align*}$$

The only off diagonal element that give a significant contribution to the yield evaluation is $\varepsilon^{dp}$, that takes into account the d population that is mis-reconstructed as p as a consequence of the partial fragment containment inside the LYSO detector.

4.2. Proton and deuteron production yields

The measured yields, computed using equation (1), for the p and d fragments populations shown in figure 8 and detected in the $90^\circ$ and $60^\circ$ angular configurations are listed hereafter:

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{r@{}lll} \Phi_p^{90^\circ} &= (4.90 \pm 0.06_{\rm stat} \pm 0.57_{\rm sys}) \times 10^{-3}~{\rm sr}^{-1} \nonumber \\ \Phi_d^{90^\circ} &= (0.70 \pm 0.02_{\rm stat} \pm 0.12_{\rm sys}) \times 10^{-3}~{\rm sr}^{-1} \nonumber \\ \Phi_p^{60^\circ} &= (11.28 \pm 0.05_{\rm stat} \pm 2.30_{\rm sys}) \times 10^{-3}~{\rm sr}^{-1} \nonumber \\ \Phi_d^{60^\circ} &= (2.15 \pm 0.02_{\rm stat} \pm 0.44_{\rm sys}) \times 10^{-3}~{\rm sr}^{-1}\end{array}. \nonumber \end{align*}$$

The calculation of the systematic uncertainties is detailed in the following (section 4.3).

4.3. Systematic uncertainties

Several sources of systematic uncertainty were considered, related to the different terms that enter equation (1), used to compute the particle yields.

In order to account for the limited precision on the measurement of the number of incoming carbon ions, the term $\sigma_{\rm N_{\rm C}}$ has been computed. Such contribution to the systematic uncertainty is related to the uncertainty on the correction factor corr_GT. It was evaluated by means of a dedicated MC simulation of the beam rate and spill shape, varying the MC beam shape models and mean rate conditions. More details can be found in Mattei et al (2015).

The systematic uncertainty $\newcommand{\e}{{\rm e}} \sigma_{\epsilon_{\tau}}$ is related to the $\newcommand{\e}{{\rm e}} \epsilon_{\tau}$ evaluation. Such term takes into account the different conditions in which the DAQ was operated, during the data taking and during the laboratory sessions dedicated to the τ (DAQ dead time) measurement. $\newcommand{\e}{{\rm e}} \sigma_{\epsilon_{\tau}}$ was computed considering the maximal and minimal τ values compatible with the extrapolation of the τ laboratory measurements at the event by event rate and event size measured during the data taking. $\newcommand{\e}{{\rm e}} \epsilon_{\tau_{\rm min}}$ and $\newcommand{\e}{{\rm e}} \epsilon_{\tau_{\rm max}}$ were calculated, respectively, and the corresponding yields ($\Phi_{\rm max}$ and $\Phi_{\rm min}$) were obtained. The systematic uncertainty $\newcommand{\e}{{\rm e}} \sigma_{\epsilon_{\tau}}$ was defined as the semi-dispersion of the corresponding yield results (Mattei et al 2015).

The systematic error $\sigma_{\rm PID}$ accounts for the uncertainty on the PID efficiency and DC LYSO crystals mis-identification. It was evaluated by varying the particle identification functions used to identify protons and deuterons in the (E_LYSO–ToF) plane. Figure 5 right shows the different PID functions determined for the hard and soft selections of the proton and deuteron fragments with respect to the default PID functions (see section 2). The semi-dispersion of $\Phi_{\rm ch}$ was computed applying the hard and soft selections and it was taken as an estimate of the PID related systematic uncertainty.

As $\varepsilon_{\rm mix}$ was evaluated from the FLUKA MC simulation with a sample of limited statistics, the statistical uncertainty of the matrix elements was propagated as systematic uncertainty when performing the yield calculation. This contribution to the total systematic error is negligible (absolute value of the order of $10^{-7}{\rm{-}}10^{-8}$ for both p and d yields measured at $90^\circ$ and $60^\circ$).

The yield relative systematic error from each systematic source is given in table 1.

Table 1. Relative systematic error for each source of systematic uncertainty contributing to the evaluation of the proton (p) and deuteron (d) yields measured at $90^\circ$ and $60^\circ$ with respect to the beam incoming direction.

	$\boldsymbol{90^{\circ}}$		$\boldsymbol{60^{\circ}}$
	$\boldsymbol{p}$	$\boldsymbol{d}$	$\boldsymbol{p}$	$\boldsymbol{d}$
$\boldsymbol{\sigma_{\rm PID}}$	$1\%$	$12\%$	$2\%$	$4\%$
$\boldsymbol{\sigma_{\rm N_{\rm C}}}$	$6\%$		$14\%$
$\newcommand{\e}{{\rm e}} \boldsymbol{\sigma_{\epsilon_{\tau}}}$	$10\%$		$15\%$

The discovery of an heterogeneous light response of the different scintillating LYSO crystals, used for the charge secondary particles identification procedure in the analysis reported in Piersanti et al (2014), required a full data sample re-analysis to correct the published fragment yields measurements.

The production yields and emission velocity distributions of the different particle species produced in the interaction of a 220 MeV/u $\newcommand{\car}{^{12}{{\rm C}}} \car$ ion beam impinging on a polymethyl methacrylate (PMMA) target were recomputed. In the new analysis, the description of the detector and geometrical efficiencies for the yield evaluation was improved, taking into account the energy dependence of the emission spectrum of secondary fragments, as well as the DAQ dead time efficiency evaluation.

While the observed triton production yield is now compatible with zero, within the experimental uncertainties, the estimated total charged secondary yield ($p+d$) is however close to the one published in Piersanti et al (2014), with an overall increase due to the improved efficiency calculation that takes into account the energy dependence of the fragments absorption probability inside the PMMA target.

The implications for beam range monitoring applications in Particle Therapy that exploit charged fragments, discussed in Piersanti et al (2014), are hence holding true and are supported by the updated results.