Characterization of the X-ray spectrum of a linear accelerator

The Nuclear Measurement Laboratory (LMN) at CEA Cadarache in France uses a high-energy electron linear accelerator, LINAC (9-21 MeV), to characterize nuclear waste drums. It enables to explore new examination modalities, such as active photon interrogation or dualenergy CT to scan large concrete objects with diameters up to 140 cm. These techniques require precise awareness of the photon spectrum emitted by the LINAC. However, direct measure of this photon energy spectrum cannot be achieved because of the accelerator pulses causing detector saturation. During the last few years, a large number of indirect methods has been developed. From an experimental point of view, the simplest indirect method for spectrum estimation method is ransmission measurements. Because it can be set up easily and accurately using an ionization chamber as well as an appropriate screen. The obtained transmission curve depends on the photon energy spectrum, which can be estimated using inverse models. In this paper, we present the development of a numerical model to determine the energy spectrum from an attenuation curve via transmission measurements which combines two types of inverse models: a continue model and a discrete model. We validate this tool using a test spectrum and its transmission curve obtained via Monte-Carlo simulation. This qualification allowed us to determine its sensitivity (signal-to-noise ratio, SNR) in order to have a good convergence. We show that if the SNR is less than 4%, we have a good estimation of the photon energy spectrum. Then, it was experimentally tested with a transmission curve obtained at the laboratory. Keywords—Linear electron accelerator, High energy computed tomography, Spectrum characterization, LINAC


I. INTRODUCTION
OR the last few years, the characterization of energy spectra from linear accelerators becomes a major issue for medical sector (radiotherapy) or industry (to develop new technics of non-destructive testing). The Nuclear Measurement Laboratory (LMN) at CEA Cadarache has an imaging system currently implemented in the CHICADE experimental platform allowing radiographic and tomographic analyses on large objects (Fig.  1). A major upgrade of this system was launched involving the provision of a new high-energy accelerator called SATURNE, with photon energy up to about 15 MeV [1]. The latter is made of new elements and elements from an older accelerator (such as accelerating waveguide), consequently its characteristics, like the electron spectrum or the photon maximal energy, are not accurately know.
This high-energy accelerator allows exploring new examination modalities, such as active photon interrogation (IPA) [2,3] or dual-energy CT [4]. These techniques require precise awareness of the photon spectrum emitted by the LINAC. However, direct measurement of high-energy linear accelerators photon energy spectrum (for example, with a hyper-pure germanium spectrometer) is difficult to implement. The SATURNE accelerator used by the laboratory is a pulsed accelerator with 5 µs pulses during which the photon flux produced is very high. Spectrometers placed directly in the output beam saturate during the pulse, making impossible the direct measurement.
Consequently, to characterize the photon spectrum of highenergy accelerators, a large number of indirect methods have been developed.
One indirect method employs Monte-Carlo simulation to determine without experimental protocol the accelerator photon energy spectrum [5]. Nevertheless, this method requires a precise knowledge of the accelerator characteristics including electron beam energy spectrum, which is difficult to assess accurately experimentally or numerically. A second indirect method is based on an experimental protocol using energy spectrum from the spectroscopy of Compton-scattered photons [6]. Photons emerging from the accelerator interact by Compton scattering in a screen. The Compton spectrum measurement at a specific scattering angle can be connected directly to the photon spectrum by a linear correlation between the incident photon energy and the recorded photon position. The main drawback of this method is that it requires a complex experimental set-up. Indeed, a massive shielding is necessary in order to avoid detector saturation and the consequent loss of resolution. Another experimental method, based on photoactivation, is based on the photonuclear reaction cross-section of different materials according to energy and radioisotopes produced by these reactions [7]. Screens of different materials, with well-known photonuclear cross-section and nuclear reaction, are irradiated. The type and quantity of radioactivity produced by this irradiation for the different screens depends on the photon energy spectrum. This method requires a large number of screens of known composition, not always available at the laboratory.
The simplest indirect method from an experimental point of view is the spectrum estimation method from transmission measurements. It consists on dose rate measurements through more or less attenuating known materials, leading to the construction of a transmission curve [8,9]. With this curve and attenuation coefficients depending on the energy for the material used, the energy spectrum can be estimated using inverse models. Attenuation measurements can be set up easily and accurately using an ionization chamber of small volume as well as an appropriate screen. Moreover, spectrum estimation from transmission measurements requires only the use of an adapted mathematical inverse model. It is for these reasons that the spectrum estimation from transmission measurements was chosen in our study.

II. NUMERICAL MODELS
Relative transmission T(xi), at attenuator thickness xi can be expressed by (1). (1) With S(x0) the detector signal without an attenuating screen. If the energy is partitioned into m discrete energy channels Ej, then it may be possible to trace the energy spectrum of the photons, , for example via an inverse process, from the transmission measurements, with the equation (2).
With the value of energy spectrum for energy Ej, the detector response and the linear attenuation coefficient of the attenuator screen at energy Ej [10].
Two methods for energy spectrum estimation from the transmission measurements have been selected: -discrete numerical model, -continuous numerical model.
Discrete models are the most widely used in the literature [8]. The searched spectrum is initialized via a given function (often triangular shape) which is then corrected by an iterative process in order to minimize the cost function summing the quadratic deviations between the experimental transmission curve and the estimated spectrum transmission. The final spectrum will then be estimated by solving the linear transmission equations from (2).
Continuous models [5] are based on mathematical equations and the fact that different photon spectra obtained from linear electron accelerators always have a similar shape (Bremsstrahlung spectrum), where only the energy of the peak flux and the width of the spectrum vary. The photon flux initially increases monotonically with the energy to a maximum (peak flux), and then decreases monotonically to the maximum photon energy (Fig. 2) equal to the maximum electron energy. Some work has proposed the use of a continuous model where these changes can be represented by a lognormal function [5]. Then, the energy spectrum can be estimated from the lognormal function equation in discrete form for the different energy channels via (3).
With the expectation and the standard deviation of the lognormal function which will be updated at each iteration. The estimated energy spectrum is used to calculate the corresponding transmission curve through the attenuating screen thicknesses and compare it with the experimental transmission curve.
For both models, the estimated spectrum is obtained by minimizing the cost function defined as the sum of the squared deviations between the reference and the estimated spectrum transmission curve for each attenuator screen thickness.

A. Advantages and drawbacks of the two models
The discrete model requires a hypothesis on the initial shape of the spectrum, based on the knowledge of the peak flux and the maximum energy of the photons, which can strongly influence the result. Indeed, the accelerators used in the literature (especially in the medical field) are standard and very well characterized. The maximum energy of the photons as well as the energy of the flux peak are known. For example, in Orman Zaid Assatel's thesis [8], the initialization is done with a triangular distribution which requires to fix the peak energy of this function corresponding to the peak flux energy of the searched energy spectrum. Since in this case the energy of the flux peak is known, the energy of the peak of the triangular function can be fixed and a monotonicity condition can be added to constrain the spectrum to be monotonically increasing before this peak and monotonically decreasing after this peak. In this way, the peak flux energy and the maximum photon energy are fixed when the spectrum is initialized and will not change during the estimation. The model used without the monotonicity condition then gives too much freedom in the estimation of the spectrum and therefore provides spectra whose shape does not make physical sense. However, when initialized with a function close to the spectrum to be estimated, this model gives very good results thanks to a channel by channel adjustment.
With the continuous model, initializing the spectrum with a lognormal function avoids hypotheses about the peak flux energy and allows realistic spectra to be obtained that are close to the simulated reference spectra. However, it imposes a lognormal distribution shape. The continuous model therefore gives an initial idea of the general shape of the spectrum without any prior hypothesis. It gives good results at low energies and provides the energy of the peak flux of the photon spectrum.

B. Combination of both models
The challenge of our study is to characterize the energy spectrum of an accelerator without first knowing the shape of the spectrum to be estimated. In view of the disadvantages and advantages of the two models, a combination of the two would allow this characterization to be carried out. The continuous model is run first in order to have a first spectrum shape and the energy of its flux peak. The discrete model is then initialized with the spectrum estimated via the continuous model (instead of the triangular function). Then, the condition on monotonicity before and after the flux peak can be added.
The thickness of the Bremsstrahlung target being very well characterized (5 mm), it allows a first estimation of the maximum energy of the photons as a function of the measured dose from Monte-Carlo calibration simulations where the electrons are monoenergetic. As this energy is not accuratly known, it cannot be fixed but will vary randomly during the different iterations within a limited range compared to the estimate (approximately ± 2 MeV, i.e., 8 channels of 0.43 MeV). Furthermore, the peak energy found with the continuous model will also vary within a small range around the initial value (about ± 0.5 MeV or 2 channels of 0.43 MeV).

C. Selection of the attenuating screen
Spectrum evaluation method based on transmission curve requires the use of different thicknesses of attenuating screen. The choice of the attenuator screen can be made on the basis of its total linear attenuation coefficients [8]. The variation of the linear attenuation coefficient of the chosen material must be monotonous over the energy range of the spectrum of interest in order to avoid two energies having the same linear attenuation coefficient value. In our case, the attenuating screen chosen is aluminum 5754 alloy. Indeed, its linear attenuation coefficient is monotonous over the energy range of the accelerator (up to 17 MeV) and it is inexpensive so easily procurable in quantity. Moreover, its density of 2.68 g.cm -3 limits the total thickness required. It takes about 60 cm to obtain 1% transmission at 15 MeV instead of 150 cm of PMMA (Polymethyl methacrylate) for example.

A. Model parameters
In order to validate the method with conditions close to the experimental ones, a modelling of an operating point of the SATURNE accelerator is carried out with the MCNP6.1 Monte-Carlo calculation code [11] allowing to obtain a known photon spectrum. The electrons are monoenergetic with an energy of 12 MeV corresponding to the maximum energy of the photon spectrum. From this simulated spectrum, a transmission curve is generated via (2) through 30 plates of 2 cm thick aluminum 5754 alloy, i.e., 60 cm. As the resulting curve has 30 measurement points, the spectrum must have a maximum of 30 energy channels.
The spectrum is estimated over a range from 0.2 MeV to 13.1 MeV with channels of 0.43 MeV (30 channels). Indeed, the lower limit of 0.2 MeV is justified by the fact that below 0.4 MeV, the photons have a contribution to the spectrum of less than 4.10 -3 %, which is 40 times smaller than the average uncertainty of our experimental dose measurements of 0.26 % (cf paragraph 4). The final channel will vary randomly between 9.66 MeV and 13.1 MeV (i.e. 8 channels of 0.43 MeV).
One of the parameters in equation (2) of the transmission is the ionization chamber response, . Experimentally, the detector used for the dose measurements is a small volume (0.125 cm 3 ) ion chamber of the PTW brand (type 31010) [12]. The response of this ionization chamber is determined by Monte Carlo simulations (MCNP6.1) where, for each incident photon energy, the average deposited energy is calculated. This response of the ionization chamber, shown in Fig. 3, is used in our algorithms. However, we have checked that the convergence of the method depends only slightly on the response used, but the response must be consistent with the transmission measurements. The convergence is the same within 10% even we use a monotonic increasing or flat response.

B. Results and validation by simulation
The method described in section 2.2 is launched 5000 times to obtain 5000 estimated spectra. For each of these spectra, 2000 iterations are performed during the channel by channel modification. The number of iterations has been fixed so that convergence is always obtained. Indeed, the convergence of the method always happens before 2000 iterations as can be seen on the example in Fig. 4. The energy of the final channel is set randomly at the beginning of each spectrum estimation. Among all spectra obtained, 10 % (i.e. 500 spectra) with the lowest cost function, and therefore the best convergence, are kept (the choice to keep only 10 % of the spectra is detailed in paragraph 5). The energy of the end channel of the estimated final spectrum, , is determined via (4) with N equal to 10 % of the 5000 estimated spectra, i.e. 500 here, and the energy of the end channel of spectrum i. The end channel of the estimated spectra is defined as the mean last non-zero channel.
. ∑ The average energy of the final channel is here 12.65 ± 0.02 MeV. The energy of the end channel of the reference spectrum is 12 MeV so the average energy found with our method is correct within 2 channels (0.65 MeV).
The final spectrum is estimated by averaging the spectra whose final channel corresponds to the channel of or the channel closest to this energy. In the end, 213 spectra corresponding to this condition are averaged to obtain the estimated final spectrum, shown in blue in Fig. 5   For each channel the uncertainty is represented by the green error bars in Fig. 5. The relative average uncertainty per channel with respect to the average spectrum is 1.3.10 -4 . This is defined by (5), where m is the number of channels in the energy of the spectra, is the flux value of channel j for spectrum i, ̅ is the average flux of channel j for the n number of retained spectra, here n=213.

(in linear
The spectra are therefore relatively close to each other. Moreover, the estimation has converged well since graphically the reference transmission curve from the models and that of the estimated average spectrum (Fig. 6. (a)) are very close: the relative deviations between these two curves (defined by (6), with T the transmission) are between -1.6.10 -4 and 2.5.10 -4 and the average relative deviation per point is -2.1.10 -5 (Fig. 6. (b)).

C. Adding Gaussian noise
In order to finalize the validation of our method, a Gaussian noise is added to the transmission curve resulting from the modelling. The standard deviation of this normal distribution is fixed at σ = 2 % in order to be higher than the maximum uncertainty on the experimental measurements which is 1.35 % with a mean uncertainty of 0.26 % (cf. paragraph 4).
The method described in section 2.2 is performed 5000 times. Among all spectra obtained, 10 % with the lowest cost function are chosen (i.e. 500 spectra). The average energy of the end channel is 11.75 ± 0.04 MeV (defined by (4)). The estimated average final spectrum is shown in Fig. 7 (in linear Fig. 7. (a) and logarithmic Fig. 7. (b) scales). The red spectrum is the reference spectrum from Monte-Carlo simulations. For each channel the uncertainty is represented by the green error bars in Fig. 7. The average uncertainty per channel with respect to the average spectrum is 3.5.10 -4 (defined by (5)), which is twice as large as for the noiseless estimate of the transmission curve. The noisy reference transmission curve and the estimated average spectrum curve are plotted in Fig. 8. (a). The relative deviations of the two transmission curves (defined by (6)) are between -9.6.10 -2 and 9.1.10 -2 , with an average relative deviation per point of 3.4.10 -2 ( Fig. 8. (b)). This deviation, which is much greater than for the previous part without noise (which was 6.1.10 -5 ), is explained by the noise of the reference transmission curve, which is therefore no longer monotonous.
The Gaussian noise added to the reference transmission curve with a standard deviation of 2 % encompasses the uncertainty in the experimental measurements (which is at most 1.35 %). Nevertheless, the estimated spectrum has a shape close to the expected spectrum. The estimation method developed is therefore validated for a spectrum derived from Monte-Carlo simulations and can be used to estimate the spectrum of the SATURNE accelerator.

IV. EXPERIMENTAL SET-UP
The transmission measurements are performed with the SATURNE accelerator through 30 aluminum 5754 alloy plates of 2 cm thickness and 10 cm side. The experimental setup is shown in Fig. 9. Dose measurements are performed with a small volume (0.125 cm 3 , type 31010) PTW ionization chamber [12] surrounded by its 3 mm PMMA cap, located after the attenuator screen plates. In order to correct for any variation in dose rate during the measurements, a CdTe probe is integrated next to the target. The measured signal is used to normalize the signal obtained with the ionization chamber (beam monitoring). At the output of the accelerator, two tungsten blocks spaced 4 mm apart allow the beam collimation. The ionization chamber is also collimated with lead, including two 20 cm blocks separated by 1 cm in the direction of beam incidence.
For each aluminum thickness, 120 dose measurements are made and averaged. The average uncertainty of the 3600 measurements made is 0.26 % with a maximum of 1.35 %. In order to make the measurement faster by reducing the number of trips to remove each aluminum plate, a step, made up of 5 plates staggered by one centimeter, is placed on a translation plate that can be remotely controlled. Thus, moving the step by one centimeter simulates the removal of a plate. Every measurements, 5 aluminum plates are removed. The experimentally obtained transmission curve is shown in Fig. 10 (with the measurement uncertainty for each thickness in light blue). The method described in section 2.2 is executed and 5000 spectra are obtained. They are sorted according to their increasing cost function in order to select only the 10 % (i.e. 500 spectra) with the best convergence. The end channel energy of the estimated final spectrum is taken to be equal to the average end channel energy of these 500 spectra.
The average energy of the final channel is 11.61 ± 0.4 MeV (defined by (4)). The spectrum is estimated by averaging the spectra whose final channel corresponds to the channel of the mean energy or to the channel closest to this energy. In the end, 77 spectra match this criterion and are averaged to obtain the estimated final spectrum shown in Fig. 11. For each channel, the uncertainty is represented by the green error bars in Fig. 11. The average uncertainty per channel with respect to the average spectrum is 7.1.10 -4 (defined by (5)).
This value is close to the average uncertainty found when estimating the spectrum from the noisy reference transmission curve (which was 3.5.10 -4 ). The experimental transmission curve and the estimated average spectrum curve are plotted in Fig. 12. (a). The relative deviations of the two transmission curves (defined by (6)) are between -5.4.10 -2 and 1.4.10 -2 , with an average relative deviation per point of -4.9.10 -3 ( Fig. 12. (b)). This deviation is also close to that found with the noisy reference transmission curve (which was 3.4.10 -2 ).
V. DISCUSSION Fig. 13 provides a visual comparison of the estimated spectrum from the experimental measurements and the spectrum from Monte-Carlo simulations. The estimated spectrum looks different from the simulated spectrum with a much slower decay after the flux peak. However, as the SATURNE accelerator is composed of an assembly of elements of different origins, its characteristics are not precisely known, so the models probably do not reflect reality. Moreover, it is unlikely that the electrons are monoenergetic. For this study, 10% of the spectra are selected for the estimation of the energy spectrum shape. This value is the result of an optimization between the execution time of the method to have a good convergence and the accuracy of the results. For example, for the spectrum estimated from the transmission curve from the simulations: -By keeping 5 % of the spectra instead of 10 %, there are 102 spectra kept for the estimation of the spectrum shape instead of 213, i.e. 70 % less spectrum for the same execution time. But the average value of the cost function of the estimated spectrum decreases by 50 % (better accuracy). -By retaining 15 % of the spectra instead of 10%, there are 321 spectra retained for spectrum shape estimation instead of 213, i.e. 40 % more spectrum for the same execution time. But the average value of the cost function of the estimated spectrum increases by 34 % (worse accuracy).
Finally, the discrete model used in our method has also been implemented with the simulated annealing method [13] which gives equivalent results with a 5 times higher convergence time. For this reason, this method has not been retained in our study.

VI. CONCLUSIONS
In this paper, the results of the development of an experimentally simple and robust tool for characterizing the photon spectrum of a high-energy accelerator via transmission measurements are presented. The combination of a discrete and a continuous numerical model enables the estimation of the energy spectrum as well as the maximum photon energy (to the nearest 0.65 MeV, compared to the reference spectrum derived from Monte Carlo simulations) of an accelerator whose characteristics are not precisely known. The method has been developed and tested using Monte-Carlo simulations. The shape of the estimated spectrum is physically realistic and very close to the reference spectrum used for comparison. The addition of noise, corresponding to the uncertainty of the experimental measurements, on the transmission curve resulting from the simulations allows to validate our method under conditions close to the experimental ones. This method was then tested on a transmission curve obtained with the SATURNE accelerator. The energy spectrum of the accelerator was estimated as well as the maximum energy of the photons, which is 11.6 ± 0.4 MeV.