Acoustic Localization of Bragg Peak Proton Beams for Hadrontherapy Monitoring †

Hadrontherapy makes it possible to deliver high doses of energy to cancerous tumors by using the large energy deposition in the Bragg-peak. However, uncertainties in the patient positioning and/or in the anatomical parameters can cause distortions in the calculation of the dose distribution. In order to maximize the effectiveness of heavy particle treatments, an accurate monitoring system of the deposited dose depending on the energy, beam time, and spot size is necessary. The localized deposition of this energy leads to the generation of a thermoacoustic pulse that can be detected using acoustic technologies. This article presents different experimental and simulation studies of the acoustic localization of thermoacoustic pulses captured with a set of sensors around the sample. In addition, numerical simulations have been done where thermo-acoustic pulses are emitted for the specific case of a proton beam of 100 MeV.


Introduction
Localization of a source is a technique in which a source is located by detecting propagated signals received in several sensors and the analysis of them [1]. There are many localization techniques proposed for wireless sensor networks [2,3]. However, in this article, a three-dimensional localization to solve the estimation of an acoustic source in a homogeneous medium is introduced. The use of acoustic sensors to locate sound sources in such practical systems is of great interest but needs further development and improved performance systems. This research has significant potential for many applications in medicine, physics, engineering, and underwater acoustics. The method to locate the tumor tissue is based on a computed tomography scan to find the area that will then be radiated by heavy particles in the Bragg peak region [4]. However, uncertainties in the patient positioning and/or in the anatomical parameters can increase the uncertainty during the radiotherapy. In these cases, acoustic source localization in medical applications has gained a lot of interest in recent years, which ought to be the necessity for improving the monitoring of tumor tissue in hadrontherapy treatments. Linear sensors can be employed for acoustic source localization in a noise environment using a time delay estimation. The method presented in this paper is based on the TDOA (time difference of arrival) [5] technique that performs very well in the localization of an acoustic event in both two-dimensional and three-dimensional spaces decreasing the error while increasing the number of sensors. The acoustic sensors in known positions and using a DAQ system to record the signal. Differences in the signal propagation path from the source lead to different phases in the detected signal. Therefore, crosscorrelation analysis is used to estimate the delays of arrivals accurately [6], even in conditions with low signal-to-noise ratio.
The pressure source localization of the Bragg peak in hadrontherapy can also be used to identify the regions of local heat released due to energy deposition. In this paper, we focus on the objective of monitoring the position for hadrontherapy through the Bragg proton beam acoustic localization. This pressure is related to the beam energy, the temporal pulse width, the size of the beam, and the number of protons by pulse, so, to some extent, it might be used in the future for dose sensitivity as well, but this aspect is out of the scope of the paper. For this reason, as a first approach, the source assessed in this article presents a pressure above the threshold of detection [7] for beams of a few million protons per spill of energies from 20 up to 200 MeV. Both in simulation and experiment, homogenous and isotropic medium is used and the wall effect is neglected since the direct signal arrives earlier to the Omni-directional receiver than the reflected signal.

Overview of Approach
Techniques based on cross-correlation and generalized correlation (GCC) [8] have been employed to determine the time difference of arrival of the signals (TDOA) given its computational cost and accuracy of the results. To obtain a better estimate of the TDOA, it is convenient to filter the signal before its integration, as shown in Figure 1 [6]. The cross-correlation between the signal and filtered by the filters and , is expressed as a function of the power spectral density , as shown below.
where {*} indicates a conjugated complex and ( ) is a frequency-dependent weight function.
Due to finite observations, we can only obtain an estimation of ( ) [9]. Therefore, to obtain the TDOA, the following expression will be used [6].
where ( ) is the obtained estimation of ( ). For each pair of sensors, the TDOA is taken as the time delay that maximizes the cross-correlation between the filtered signals of both sensors, that is: = ( ) .
A general model for three-dimensional (3-D) estimation of a source using receivers is developed. To obtain the location of the source, we start by knowing the spatial position ( , , ) of a certain number of sensors. Let ( , , , ), the position of the source to be located, the distance between the source and the i-th sensor will be: The range difference in distance between the i-th receiver and the first receiver, , is given by: The cross-correlation R x i x j between the signal x i and x j filtered by the filters H i and H j , is expressed as a function of the power spectral density G x i x j , as shown below.
where {*} indicates a conjugated complex and ϕ GCC ( f ) is a frequency-dependent weight function. Due to finite observations, we can only obtain an estimation of G x i x j ( f ) [9]. Therefore, to obtain the TDOA, the following expression will be used [6].
For each pair of sensors, the TDOA is taken as the time delay that maximizes the cross-correlation between the filtered signals of both sensors, that is: x i x j (t ) . A general model for three-dimensional (3-D) estimation of a source using M receivers is developed. To obtain the location of the source, we start by knowing the spatial position (x i , y i , z i ) of a certain number of sensors. Let (x s , y s , z s ,), the position of the source to be located, the distance between the source and the i-th sensor will be: Sensors 2019, 19,1971 3 of 13 The range difference in distance between the i-th receiver and the first receiver, d i1 , is given by: where c is the sound velocity in the medium, d 1 is the distance between the first receiver and the source, and τ i1 is the estimated TDOA between the first receiver and the i-th receiver [10]. Equation (4) considered for all the sensors form a nonlinear equation system whose solution can be found by several ways. After studying different resolution methods, it was decided to use the Newton-Raphson method since it offers very good results and computation time.

Newton-Raphson Method
To get the position of the thermoacoustic source inside the medium, we have solved the nonlinear system using the Newton-Raphson method [11] by means of partial derivatives. Consider a system of m equations and n unknowns.
This system can be written in vector form as f (x) = 0, where f is a vector of m dimensions and x is a vector of n dimensions. To solve this system of equations, we have to find a vector x such that the function f (x) equals the null vector. If we call η to the solution of the system and x r to an approximation of it, we can develop f in Taylor series around this approximation as: Since f (η) = 0 then we get, as an approximation, the following.
Now, we can define the vector x r+1 as this approximation, which is closer to the root than x r . We can continue with the iterative method to obtain approximations closer and closer to the solution. To write the iterative method, the term ∇ f (x r ) is replaced by the Jacobian of the function f , that is: which is a n × m matrix. Then, it is possible to obtain a new value of x r+1 by solving the following relationship.
Then, iteratively, we can approximate more and more the x r+1 to η until a solution error x r − x r+1 previously fixed is reached. This method has been compared with other algorithms for solving systems of nonlinear equations and has shown some advantages over the rest like good accuracy and low computing cost. However, if the initial solution value of the system differs greatly from the real solution, then the method does not converge conveniently. Figure 2 shows the convergence of the location algorithm as a function of the distances between the initial point of the method, the result of the reconstruction of the position and the real position for a total of 10,000 simulations.

Bragg Peak
The main physical advantage of heavy particles as compared to photons is their characteristic depth-dose profile, the known Bragg curve in honor of Sir William Henry Bragg who investigated the energy deposition of alpha particles, which form a Radium source in the air at the beginning of the last century [12]. While the photon dose decreases exponentially with penetration depth according to the absorption law for electromagnetic radiation, the depth-dose of heavy charged particles exhibits a flat plateau region with a low range of the particles [13]. This paper uses an analytical approach presented by T. Bortfiel in 1996 [4], and the numerical representation valid for protons with energies between 10 MeV and 200 MeV. Thus, the energy of single protons along a Zaxis in a homogeneous medium (water) is considered. The total energy released in the medium per unit mass in the Z-axis is shown below.
where ( ) is the proton flow, that is, the number of protons per , ( ) is the energy deposited on the Z-axis, and represents the mass density of the medium. The method makes use of a midpoint where a certain fraction γ of the energy released in nuclear interactions is absorbed locally while the rest is ignored. Then, the total absorbed energy ( ), will be given by: A relationship between the initial energy ( = 0) = and the range = in the medium can be approximated as = for = 1.5. This relation is valid for protons with energy close to 250 MeV. The factor is proportional to the square root of the effective atomic mass of the medium [4]. Using the inverse for ≤ 0.5 cm and assuming to be given in units of MeV, the best fit parameters for ( ) are = 1.77, = 2.2 × 10 for the proton in water [4]. The remaining energy ( ) at an arbitrary depth ≤ fails to travel the distance − according to the rangeenergy relationship shown below.
For energies above 20 MeV, there is non-negligible probability that protons may be lost from the beam due to nuclear interactions. This non-elasticity was studied and tabulated by Janni [14] as a

Bragg Peak
The main physical advantage of heavy particles as compared to photons is their characteristic depth-dose profile, the known Bragg curve in honor of Sir William Henry Bragg who investigated the energy deposition of alpha particles, which form a Radium source in the air at the beginning of the last century [12]. While the photon dose decreases exponentially with penetration depth according to the absorption law for electromagnetic radiation, the depth-dose of heavy charged particles exhibits a flat plateau region with a low range of the particles [13]. This paper uses an analytical approach presented by T. Bortfiel in 1996 [4], and the numerical representation valid for protons with energies between 10 MeV and 200 MeV. Thus, the energy of single protons along a Z-axis in a homogeneous medium (water) is considered. The total energy released in the medium per unit mass in the Z-axis is shown below.
where φ(z) is the proton flow, that is, the number of protons per cm 2 , E(z) is the energy deposited on the Z-axis, and ρ represents the mass density of the medium. The method makes use of a midpoint where a certain fraction γ of the energy released in nuclear interactions is absorbed locally while the rest is ignored. Then, the total absorbed energyD(z), will be given by: A relationship between the initial energy E(z = 0) = E 0 and the range z = R 0 in the medium can be approximated as R 0 = αE p 0 for p = 1.5. This relation is valid for protons with energy close to 250 MeV. The factor α is proportional to the square root of the effective atomic mass of the medium [4]. Using the inverse for R o ≤ 0.5 cm and assuming E o to be given in units of MeV, the best fit parameters for E o (R o ) are p = 1.77, α = 2.2 × 10 −3 for the proton in water [4]. The remaining energy E(z) at an arbitrary depth z ≤ R o fails to travel the distance R o − z according to the range-energy relationship shown below.
Sensors 2019, 19,1971 5 of 13 For energies above 20 MeV, there is non-negligible probability that protons may be lost from the beam due to nuclear interactions. This non-elasticity was studied and tabulated by Janni [14] as a function of the residual range (R 0 − z). The proton flow φ(z) can be written by using the equation below.
where φ 0 is the primary fluence and the slope parameter β was determined to be β = 0.012 cm −1 [4]. Thus, the distribution of the deposition of energy along the depth range can be expressed as the equation below.
where Γ represents the gamma function, ζ = (R 0 − z)/σ, with a σ value of 0.012R 0.935 0 , and ε represents a relatively small fraction of the fluence Φ 0 in the peak. Figure 3 shows the distribution of the dose as a function of the range for a different proton energy.
where is the primary fluence and the slope parameter was determined to be = 0.012 [4]. Thus, the distribution of the deposition of energy along the depth range can be expressed as the equation below.
where represents the gamma function, = ( 0 − )⁄ , with a value of 0.012 . , and represents a relatively small fraction of the fluence in the peak. Figure 3 shows the distribution of the dose as a function of the range for a different proton energy.

Thermoacoustic Model
In the thermoacoustic case, an excited point source emits a pressure wave proportional to the first time derivate of the excitation pulse [15]. Hence, a Gaussian excitation pulse leads to a bipolar acoustic emission consisting of a positive compression, which results in an increase in pressure. This is followed by a negative rarefaction, which is a decrease of pressure. The positive and negative pressure peaks are not only due to the heating and cooling of the medium, but the variation of the heating rate also plays a role. The medium expands or contracts according to its coefficient of thermal volumetric expansion α'. As a result, a pressure wave is observed. The pressure wave from an energy deposition in a region can be understood as the sum of the individual responses that would be observed from decomposing the spatial deposition into point sources. The resulting pressure signal depends on the time derivative of the excitation pulse. The amplitude of the wave depends on the energy deposited, the number of protons per pulse of the beam, and the temporal shape of the excitation pulse. A dose of 1 Gy generates a ~240 μK temperature increase in water [15]. Ignoring heat diffusion and cinematic viscosity, the wave equation that describes the pressure at a time and position ⃑, is shown below [16][17][18][19].

Thermoacoustic Model
In the thermoacoustic case, an excited point source emits a pressure wave proportional to the first time derivate of the excitation pulse [15]. Hence, a Gaussian excitation pulse leads to a bipolar acoustic emission consisting of a positive compression, which results in an increase in pressure. This is followed by a negative rarefaction, which is a decrease of pressure. The positive and negative pressure peaks are not only due to the heating and cooling of the medium, but the variation of the heating rate also plays a role. The medium expands or contracts according to its coefficient of thermal volumetric expansion α'. As a result, a pressure wave is observed. The pressure wave from an energy deposition in a region can be understood as the sum of the individual responses that would be observed from decomposing the spatial deposition into point sources. The resulting pressure signal depends on the time derivative of the excitation pulse. The amplitude of the wave depends on the energy deposited, the number of protons per pulse of the beam, and the temporal shape of the excitation pulse. A dose of 1 Gy generates a ∼ 240 µK temperature increase in water [15]. Ignoring heat diffusion and cinematic viscosity, the wave equation that describes the pressure p at a time t and position r , is shown below [16][17][18][19]. where c s (ms −1 ) represents the speed of sound in the middle, C p (J kg −1 K −1 ) is the specific heat capacity, and → r , t (J s −1 m −3 ) is the energy density deposited in the medium. Equation (15) can be solved using the Kirchhoff integral as shown below.
where p → r , t denotes the hydrodynamic pressure at a given place and time. The values for the thermoacoustic model were an energy of 100 MeV, a temporal profile of 1 µs, 3.4 × 10 6 protons per pulse, the beam with a size of 1 mm, and a sensor located 40 mm from the Bragg peak. The characteristics of the simulation are given by simulation results from different studies, as well as their application in clinical cases [7,15,[20][21][22][23][24][25]. The values for this case are shown in Figure 4. As a result, the pressure obtained at the reception point will be the signal that will be emitted by the piezoelectric [26] transducer to simulate a bipolar source that will be located by the sensor array. where (ms −1 ) represents the speed of sound in the middle, (J kg K ) is the specific heat capacity, and ( ⃗, ) (J s m ) is the energy density deposited in the medium. Equation (15) can be solved using the Kirchhoff integral as shown below.
where ( ⃗, ) denotes the hydrodynamic pressure at a given place and time. The values for the thermoacoustic model were an energy of 100 MeV, a temporal profile of 1 μs, 3.4 × 10 protons per pulse, the beam with a size of 1 mm, and a sensor located 40 mm from the Bragg peak. The characteristics of the simulation are given by simulation results from different studies, as well as their application in clinical cases [7,15,[20][21][22][23][24][25]. The values for this case are shown in Figure 4. As a result, the pressure obtained at the reception point will be the signal that will be emitted by the piezoelectric [26] transducer to simulate a bipolar source that will be located by the sensor array.

Experimental Setup
The experimental data measurements were made in the laboratories of the physics department at the Universitat Politècnica de València (Spain). There is a water tank with a volume of 0.64 m with a programmable 3D axis system MOCO PI MICOS arm that was programmed to move the Reson TC4014 receiver hydrophone in the tank. The hydrophone has a receiving sensitivity of −186 ± 3 dB @ 1 V μPa ⁄ and a frequency response from 15 kHz to 480 kHz. The emitter hydrophone is a Reson TC4038 with a transmitting response of 110 dB @ 1 μPa V ⁄ @ 1m and a frequency response from 50 kHz to 800 kHz. Figure 5 shows the experimental setup with the transmitter and receiver inside the tank. A National Instruments data acquisition system was used with PXI type cards to generate the signal used as input of the linear E&I A150 amplifier that feeds the transmitter. Both the receiving and the feeding signal were captured. The latter was captured with an ×100 probe to avoid overloads in the system. All signals were stored at 10 Ms/s with a duration of 500 µs.

Experimental Setup
The experimental data measurements were made in the laboratories of the physics department at the Universitat Politècnica de València (Spain). There is a water tank with a volume of 0.64 m 3 with a programmable 3D axis system MOCO PI MICOS arm that was programmed to move the Reson TC4014 receiver hydrophone in the tank. The hydrophone has a receiving sensitivity of −186 ± 3 dB @1 V/µPa and a frequency response from 15 kHz to 480 kHz. The emitter hydrophone is a Reson TC4038 with a transmitting response of 110 dB @1 µPa/V @ 1m and a frequency response from 50 kHz to 800 kHz. Figure 5 shows the experimental setup with the transmitter and receiver inside the tank. A National Instruments data acquisition system was used with PXI type cards to generate the signal used as input of the linear E&I A150 amplifier that feeds the transmitter. Both the receiving and the feeding signal were captured. The latter was captured with an ×100 probe to avoid overloads in the system. All signals were stored at 10 Ms/s with a duration of 500 µs. Two experiments have been performed. First, a calibration is done to reduce the uncertainty due to the time of arrival from which the speed of the sound has been measured. For this, 12 different reception positions were assessed in a straight line on the emitter axis. By having the 12 measurements along the line, a time-distance linear adjustment is made whose slope value corresponds to the speed of sound. In addition, the calibration allows a time correction that results in a decrease of the arrival time of the signal, according to the linear adjustment obtained in the fit.
In the second part of the experiment, 12 reception points have been set that will correspond to the 12 sensor positions of the array. The mechanical axis moves the hydrophone to each of the points and then the signal is emitted and recorded 10 times. Figure 5b shows one of the measuring points where the Reson TC4014 sensor is fixed to the mechanical arm.

Numerical Simulation
To evaluate the localization method described, the reconstruction of the location of a Gaussian pulse source of 50 μs is simulated from the reception of four sensors located on the lateral surface of different coordinates. Figure 6 shows the position of the sensors and the source in the space for the simulated model [6]. In this simulation, the sources are always into the volume covered by the coordinates of the sensors. Despite it being possible to locate the source with three sensors, the use of at least four sensors is incorporated to both improve the reliability and quality of the results. Two experiments have been performed. First, a calibration is done to reduce the uncertainty due to the time of arrival from which the speed of the sound has been measured. For this, 12 different reception positions were assessed in a straight line on the emitter axis. By having the 12 measurements along the line, a time-distance linear adjustment is made whose slope value corresponds to the speed of sound. In addition, the calibration allows a time correction that results in a decrease of the arrival time of the signal, according to the linear adjustment obtained in the fit.
In the second part of the experiment, 12 reception points have been set that will correspond to the 12 sensor positions of the array. The mechanical axis moves the hydrophone to each of the points and then the signal is emitted and recorded 10 times. Figure 5b shows one of the measuring points where the Reson TC4014 sensor is fixed to the mechanical arm.

Numerical Simulation
To evaluate the localization method described, the reconstruction of the location of a Gaussian pulse source of 50 µs is simulated from the reception of four sensors located on the lateral surface of different coordinates. Figure 6 shows the position of the sensors and the source in the space for the simulated model [6]. In this simulation, the sources are always into the volume covered by the coordinates of the sensors. Despite it being possible to locate the source with three sensors, the use of at least four sensors is incorporated to both improve the reliability and quality of the results.  To evaluate the algorithm, the volume of the cube has been modified between 27.0 × 10 m and 512.0 × 10 m . Positions of the sensors are shown in Table 1, where represents the length of the cube, whose values were 200, 300, 400, 500, and 600 mm [6].  Table 2 shows the deviation of the position of the simulated source with respect to the real position of the source. The results can be expressed as a function of the distance between the reconstructed position and the real position of the source. Figure 7 shows the results, where the To evaluate the algorithm, the volume of the cube has been modified between 27.0 × 10 −3 m 3 and 512.0 × 10 −3 m 3 . Positions of the sensors are shown in Table 1, where H represents the length of the cube, whose values were 200, 300, 400, 500, and 600 mm [6]. Table 2 shows the deviation of the position of the simulated source with respect to the real position of the source. The results can be expressed as a function of the distance between the reconstructed position and the real position of the source. Figure 7 shows the results, where the abscissa axis shows the volume in m 3 , while the ordinates axis shows the distance in mm between the prediction of the algorithm and the real position. In addition, a fitting line to the results is shown.  These reconstructed positions do not exceed 5 mm of the real position for all the studied cases. In addition, it is important to calculate the algorithm for a future application in real time. For this reason, Figure 8 shows the calculation times for 4, 6, 8, 10, and 12 sensors for a source position, depending on changes in volume. Once satisfactory results of the localization algorithm have been obtained by simulating different known source points in known sensor positions, the localization method has been evaluated experimentally in the next section.

Experimental Localization with Thermoacoustic Signals
To validate the localization method, 12 reception positions (sensors) have been configured for a single source position. The signal emitted is the pulse obtained by the simulation of the thermoacoustic model, given by Equation (15). Table 3 shows the positions of the measuring points, which have all been referenced to the lower corner of the tank. Table 3. Position of the source and sensor reception points.

Axis
Sensor Positions [cm]  Once satisfactory results of the localization algorithm have been obtained by simulating different known source points in known sensor positions, the localization method has been evaluated experimentally in the next section.

Experimental Localization with Thermoacoustic Signals
To validate the localization method, 12 reception positions (sensors) have been configured for a single source position. The signal emitted is the pulse obtained by the simulation of the thermoacoustic model, given by Equation (15). Table 3 shows the positions of the measuring points, which have all been referenced to the lower corner of the tank. Table 3. Position of the source and sensor reception points.

Axis
Sensor Positions [cm] Once satisfactory results of the localization algorithm have been obtained by simulating different known source points in known sensor positions, the localization method has been evaluated experimentally in the next section.

Experimental Localization with Thermoacoustic Signals
To validate the localization method, 12 reception positions (sensors) have been configured for a single source position. The signal emitted is the pulse obtained by the simulation of the thermoacoustic model, given by Equation (15). Table 3 shows the positions of the measuring points, which have all been referenced to the lower corner of the tank. Two inputs channels were used. Channel 1 recorded the emitted signal passed to the amplification system while channel 2 recorded the received signal of the emitter Reson TC4014. Both signals were used afterward for correlation. The receiver is fixed to the MOCO programmed axis that moves the hydrophone to each of the measurement positions from position 1 to position 12 in the direction shown in Figure 9. By doing this, we have the same information as when we using a 12-element sensor array.  Figure 9. By doing this, we have the same information as when we using a 12element sensor array. Figure 9. This figure shows the configuration of the receiver positions. The source is represented by a sphere (white) that is inside the volume generated by the sensors (cylinders of grey color with black tip) whose route is shown with the black lines and the conical red marks. In this figure, a smaller volume is represented inside the tank for its correct visualization.
The analysis of the measurements has taken into account the calibration of the measurement of the speed of sound 1492 m/s. This correction allows for a better detection precision in the TOA. In addition, the speed of the localization algorithm and the accuracy of the results will be evaluated, according to the number of reception points. Figure 10 shows an example of the signal emitted and the signal received as well as the correlation between them in order to obtain the TOA as the time of maximum amplitude of the correlation [9].
(a) (b) Figure 9. This figure shows the configuration of the receiver positions. The source is represented by a sphere (white) that is inside the volume generated by the sensors (cylinders of grey color with black tip) whose route is shown with the black lines and the conical red marks. In this figure, a smaller volume is represented inside the tank for its correct visualization.
The analysis of the measurements has taken into account the calibration of the measurement of the speed of sound 1492 m/s. This correction allows for a better detection precision in the TOA. In addition, the speed of the localization algorithm and the accuracy of the results will be evaluated, according to the number of reception points. Figure 10 shows an example of the signal emitted and the signal received as well as the correlation between them in order to obtain the TOA as the time of maximum amplitude of the correlation [9].
The analysis of the measurements has taken into account the calibration of the measurement of the speed of sound 1492 m/s. This correction allows for a better detection precision in the TOA. In addition, the speed of the localization algorithm and the accuracy of the results will be evaluated, according to the number of reception points. Figure 10 shows an example of the signal emitted and the signal received as well as the correlation between them in order to obtain the TOA as the time of maximum amplitude of the correlation [9].  Table 4 shows the estimated location of the source (mean and standard deviation) for different groups of sensors labelled, according to Figure 9 and Table 3. From Table 4, we observe that the results obtained with this configuration are satisfactory and below 1 mm. Furthermore, similar accurate results are also possible with fewer sensors such as four sensors surrounding the source and covering different directions isotropically. In cases that sensors are only in a specific region, larger offsets (several mm) may appear in specific coordinates.  10 1, 3, 4, 5, 7, 8, 9, 10, 11, 12 54.00 ± 0.17 53.00 ± 0.14 38.01 ± 0.24 10 12 All of them 54.01 ± 0.17 53.00 ± 0.12 37.99 ± 0.38

Conclusions
In this paper, the proton range of the analytical medium has been used to calculate the distribution of the dose in space. Specifically, we have addressed values of time between 1 µs and 10 µs, energies between 20 MeV and 200 MeV, and protons per pulse from 3.1 × 10 8 to 8 × 10 8 . In the same way, an analytical model to calculate the distribution of the initial pressure at a single point has been used with a computing program to find the numerical solution of the general wave equation and calculate the acoustic waves resulting from the initial pressure. Some simulation parameters include pulse widths, beam energy, and spatial and temporal configuration. Thus, it is possible to determine the parameters and frequency spectrum to select the frequency responses of the transmitter and the receiver for the experiment.
This study has described and tested a procedure for monitoring the location of a hadrontherapy acoustic source based on the detection of the signal through piezoelectric sensors and on a model for calculating the position of the energy deposition. The localization algorithm has been applied for different configurations of sensors. Thus, the results show a significant improvement when a greater number of sensors is used. For a minimal set of four sensors, the results are better if the sensors cover different directions of the space. The accuracy of the results improves as the number of sensors increases, as shown in Table 4. Although the calculation time increases with the number of sensors, the difference is not significant for any of the proposed cases. Thus, the results indicate that it would be possible to monitor in real time the hadrontherapy treatment acoustically. At first look, the case studied may be considered too simplistic since the human body is neither homogenous nor isotropic and sensors are not omnidirectional. Therefore, for a practical case in hadrontherapy, all these aspects should be considered when taking care of the geometry and tissues involved, the real response of sensors, etc. However, for most of the cases, the results in small differences in the technique proposed, and the situation considered in the paper is good, and most probably the best, for a first general assessment of the technique proposed.