Elena E Romero

John J Martinez

Christophe Bérenguer
email: christophe.berenguer@grenoble-inp.fr

Degradation of a wind-turbine drive-train under turbulent conditions : effect of the control law

This work compares the simulated degradation produced at the drive-train of a given wind turbine when it is functioning at optimal and sub-optimal operating points. The simulation considers different classes of wind conditions and supposes the use of the same Maximum Power Point Tracking (MPPT) algorithm but calibrated with different feedback gains. The dissipation at the drive-train is modeled using contact mechanics principles, and it is intended for modeling the mechanical fatigue due to changes in wind speed and turbulence. The paper presents the proposed model for the drive-train degradation and the obtained power curves for a simulated variable-speed 2M W (100m rotor diameter) wind turbine, with fixed gear box and horizontal-axis.

I. INTRODUCTION

Nowadays, the growing interest in the renewable energy market has driven research in wind turbines. This sector grew up near 53% in 2020, installing more than 93 GW of wind power, and some experts expect that over 469 GW of new wind capacity will be installed in the next five years. Nevertheless, some disbelief is still associated with the numerous random factors that can affect the profitability of this technology [START_REF] Lee | GWEC -global wind report 2021[END_REF]. The nature of the wind can be the origin of unexpected fatigue damage by the aerodynamic loads that create deformations of the drive-train, reducing the useful life of the turbine and raising the cost of the energy [START_REF] Bianchi | Wind Turbine Control Systems[END_REF].

The wind turbine drive-train is a complex mechanical system where many load modes can induce different problems. Among all the faults of wind turbines, the transmission faults occur more frequently and usually lead to the most extended maintenance works and excessive maintenance cost [START_REF] Lu | Wind turbine planetary gearbox fault diagnosis based on self-powered wireless sensor and deep learning approach[END_REF]. In this regard, the analysis of these parts of the system has been the high interest for the industry and the academia, providing numerous methods to know the torsional loads and fatigue behaviors [START_REF] Li | Fatigue reliability analysis of wind turbine drivetrain considering strength degradation and load sharing using survival signature and FTA[END_REF], [START_REF] Pan | Performance degradation assessment of wind turbine gearbox based on maximum mean discrepancy and multi-sensor transfer learning[END_REF]. Nevertheless, many designed models may result in unstable mechanical modes, because they ignore the dynamics of the flexible shaft [START_REF] Rahimi | Drive train dynamics assessment and speed controller design in variable speed wind turbines[END_REF], resulting in inaccurate estimation. The dissipated energy has been suggested as an indicator of damage produced by tower top displacements [START_REF] Berglind | Wind turbine control with active damage reduction through energy dissipation[END_REF]. A similar degradation principle could be used to model degradation on other wind turbine components.

On the other hand, steady wind speed can cause less vibration on a wind turbine shaft compared to the case with the same average but high turbulence intensity wind speed. The latter may cause sudden overloads to the shaft leading to its degradation. Therefore, to obtain a complete image of the effects of wind conditions on the shaft degradation, the automatic control system must be taken into account [START_REF] Ma | Flexible wind speed generation model: Markov chain with an embedded diffusion process[END_REF].

This paper presents a novel shaft model, for Variable-Speed Fixed-Pitch wind turbines, to simulate the drive-train degradation produced by the system's operation at points below and above optimal feedback control. Different wind speeds and turbulence intensities are taken into account to illustrate the obtained degradation due to the random effect of the wind conditions. The proposed simulation relies on a wind speed model estimated from real data measurements of a laminar wind speed and two classes of turbulent wind speed data obtained from stochastic differential equations. Furthermore, a mechanical contact principle is used for modeling the transmission damping, which is leveraged to estimate the power and the energy dissipated at the shaft and consequently to estimate an image of the fatigue in the drive-train.

II. WIND TURBINE DYNAMICAL MODEL

A. Mechanical model

For modeling the degradation phenomena that have a significant influence on the lifetime of the drive-train in a wind turbine, we consider a simplified drive-train representation shown in Fig. 1. The nomenclature is presented in Table I. It is possible to represent the drive-train system as two rigid bodies linked by a flexible shaft that is being deformed with an angle θ s when the rotor speed ω r is slightly different to the generator speed ω g . The rigid bodies are an image of all mechanical devices located at each side of the effective shaft. The considered system corresponds to a Variable-Speed Fixed-Pitch turbine (VS-FP); in this type of turbine, it is common to use a reduced model of the wind energy conversion system taking into account the resonance mode of the drive-train, and the simplified system can thus be described by the following equations, as presented in [START_REF] Bianchi | Wind Turbine Control Systems[END_REF] :

θs ωr ωg = 0 1 -1 -Ks Ir -Bs Ir Bs Ir Ks Ig Bs Ig -Bs Ig θs ωr ωg + 0 0 1 Ir 0 0 -1 Ig ( τr τc ) (1)
where τ r is the aerodynamic torque applied to the rotor by the wind, and can be calculated as:

τ r = 1 2 ρπR 3 C p (λ, β) λ V 2 (2)
with V the wind speed, ρ the air density and R the rotor radius, C p (λ, β) is the power coefficient, β is the pitch angle and λ the tip speed ratio, defined as:

λ = ω r R V (3)
The power coefficient C p is an indicator of the useful power in the wind flow and it is a function of the pitch angle β and λ. The power coefficient can be defined using different methods, including the theory of aerodynamics, blade element momentum (BEM) theory, computational fluid dynamics (CFD), fuzzy logic, or generalized dynamic wake (GDW) models [START_REF] Asghar | Estimation of wind turbine power coefficient by adaptive neuro-fuzzy methodology[END_REF], [START_REF] Wang | The effect of reduced power operation of faulty wind turbines on the total power generation for different wind speeds[END_REF], [START_REF] Siavash | Prediction of power generation and rotor angular speed of a small wind turbine equipped to a controllable duct using artificial neural network and multiple linear regression[END_REF]. Nevertheless, a numerical approximation of the aerodynamic power coefficient is often accepted for achieving enough simplicity and accuracy [START_REF] Dai | Research on power coefficient of wind turbines based on SCADA data[END_REF], [START_REF] Garcia | Wind Energy Systems[END_REF]. This type of model has the form:

C p (λ, β) = c 1 c 2 λ i - c 3 β -c 4 λ i β -c 5 β x -c 6 e - c 7
λ i +c 8 λ (4) In [START_REF] Saint-Drenan | A parametric model for wind turbine power curves incorporating environmental conditions[END_REF], Eq.( 4) was optimized, obtaining the parameter λ i and the different constants as: The torque and power coefficients are of special interest for control purposes. The power coefficient C p has its maximum at (λ 0 , β 0 ), with β 0 being a very small angle, ideally zero, and maximum conversion efficiency is accomplished at λ 0 . To realise the potential benefits of the variable-speed operation, the rotational speed must be adjusted initially in proportion to the wind speed to maintain an optimum tipspeed-ratio, [START_REF] Bianchi | Wind Turbine Control Systems[END_REF].

λ i -1 = (λ + c 9 β) -1 -c 10 (β 3 + 1) -1 (5)
An optimal control law, of VS-FP turbines, considers the generator torque control τ c as a function of the rotor speed, see [START_REF] Johnson | Control of variablespeed wind turbines: Standard and adaptive techniques for maximizing energy capture[END_REF], as :

τ c = K c (ω r ) 2 (6)
where K c is an optimal feedback control gain given by:

K c = 1 2 AR 3 C pmax λ 0 3 ( 7
)
where A is the rotor swept area, C pmax is the maximum power coefficient and λ 0 is the corresponding tip-speed ratio when C pmax occurs. Thus, the generator power will be:

P g = τ c ω r (8)
and, by consequence the generated energy will be:

E g = t 0 P g dt (9)

B. Proposed Degradation Model

This section aims to present a novel drive-train model that considers contact mechanics principles to simulate the degradation in the system, using dissipated energy at the turbine shaft, denoted E d , as an indicator of deterioration.

A two-mass model usually represents the drive train of a wind turbine. Fig. 1 shows a system of this type, where the low-speed mass of the turbine is connected to the high-speed mass of the generator through a flexible shaft modeled as a spring and damper. Here, we consider that the drive train system is subject to various phenomena producing fatigue, as fatigue by vibration, by friction, by impacts, and by other aerodynamics factors that have to be considered in the deterioration model. In addition, cyclic torque fluctuations potentially reduce the useful life of the drive-train [START_REF] Bianchi | Wind Turbine Control Systems[END_REF], [START_REF] Rahimi | Drive train dynamics assessment and speed controller design in variable speed wind turbines[END_REF], [START_REF] Cetrini | Online fatigue alleviation for wind turbines by a robust control approach[END_REF]. For this reason, we consider a flexible shaft modeled as a spring and damper based on contact mechanics, adopting the contact model proposed in [START_REF] Hunt | Coefficient of restitution interpreted as damping in vibroimpact[END_REF]. Thus, we consider here a damping coefficient that can be modeled as a nonlinear function of the angular deformations.

In this work, the damping coefficient B s is considered to be a function of the torsion angle θ s , of the constant stiffness of the transmission K s , and of a constant parameter α that depends of the material. That is,

B s (θ s ) = 3 2 θ s αK s (10)
Therefore, the damping torque can be obtained as:

τ d = B s (θ s )(ω g -ω r ) (11)
hence providing the following dissipated power:

P d = τ d (ω g -ω r ) (12)
Consequently, the amount of energy that is dissipated by the drive-train will be:

E d = t 0 P d dt (13)
which is considered here as an image of the drive-train deterioration. Remark that this dissipated energy is a function of both the amplitude of the angular shaft torsion and the square of the relative velocity (ω g -ω r ).

III. SIMULATION PROCESS

This section illustrates the behavior of a wind turbine for different wind conditions. In particular, we present the obtained levels of dissipated energy at the turbine shaft that are compared to that actually generated by optimal and suboptimal control feedback gains.

A. Wind Speed Generation Model

To obtain a complete analysis of a wind turbine, it is necessary to consider the variations of wind speeds, along the time, that can affect the efficiency and durability of the mechanical parts of the turbine.

Different types of models have been developed in the literature to reproduce or predict wind speed in short-term periods, see for instance [START_REF] Manwell | Wind Energy Explained[END_REF], [START_REF] Liu | A new hybrid ensemble deep reinforcement learning model for wind speed short term forecasting[END_REF]. By considering the description of relative motion on fluids proposed by Reynolds [START_REF] Reynolds | On the dynamical theory of incompressible viscous fluids and the determination of the criterion[END_REF], equation ( 14) allows modeling the wind speed V (t), at any instant t, by taking into account its mean value V (t) and its fluctuation v(t):

V (t) = V (t) + v(t) (14)
The term V (t) is often considered as an output of a simple low-pass filter corresponding to the daily, monthly, season or annual mean behavior, and the fluctuation v(t) can be considered as an output of a high-pass filter, see for instance [START_REF] Calif | PDF models and synthetic model for the wind speed fluctuations based on the resolution of Langevin equation[END_REF].

In [START_REF] Ma | Flexible wind speed generation model: Markov chain with an embedded diffusion process[END_REF], a modelling approach is proposed, which considers that the wind speed dynamics can be modelled as a diffusion process following a stochastic differential equation, defining a so-called Ornstein-Uhlenbeck (OU) process:

dV (t) = a(V (t), t)dt + b(V (t), t)dW (15)
where a(V (t), t) and b(V (t), t) are the drift and diffusion terms, while dW is the standard Wiener process (or standard Brownian motion, a continuous process whose increments are normally distributed).

Turbulence can be classified into three different classes, where the first two include 99% of the wind speed sequences [START_REF] Calif | PDF models and synthetic model for the wind speed fluctuations based on the resolution of Langevin equation[END_REF], [START_REF] Calif | Characterization and Stochastic Modeling of Wind Speed Sequences[END_REF]. The method presented in [START_REF] Ma | Flexible wind speed generation model: Markov chain with an embedded diffusion process[END_REF] allows generating different classes of wind speed, by using the stochastic equation ( 16) as a particular case of model [START_REF] Johnson | Control of variablespeed wind turbines: Standard and adaptive techniques for maximizing energy capture[END_REF], by suitably choosing the model parameters â, b, and û.

dV (t) = -â(V (t) -û)dt + b dW (t) (16)
In our experiment, we used the parameters â, b, and û estimated from real wind speed records on a specific period and presented in [START_REF] Ma | Flexible wind speed generation model: Markov chain with an embedded diffusion process[END_REF].

The wind speed simulation has thus been performed using parameters presented in [START_REF] Ma | Flexible wind speed generation model: Markov chain with an embedded diffusion process[END_REF] providing the following models:

dV (t) = -0.0314(V (t)-10.0245)dt+0.2517 dW (t) (17) dV (t) = -(V (t) -10.0245)dt + 0.6459 dW (t) (18)
During the simulation process, it was considered different wind speed scenarios with an air density of 1.1 kg/m 3 :

• The first scenario corresponds to laminar wind flow where real measurements of wind speed were used to feed the model, see Fig. 3a). In addition, turbulent regimen was also simulated, as it was mentioned before, with two classes of turbulent wind:

• Equation ( 17) was used to generate wind flow with a low variance, see Fig. 3b) • Equation ( 18) was used to generate wind flow with a high variance, see Fig. 3c). Using these three types of wind conditions makes it possible to complete the different situations that a turbine may be subjected to in a more realistic environment.

B. Simulation Setting

As was mentioned before, it is common to use a simplified representation for the VS-FP turbines, as shown in Fig. 1, where whole transmission system is considered as a system of two rigid bodies connected by a flexible shaft. In this simulation, we consider a VS-FP turbine of 2 MW with a 100 m rotor diameter, with fixed gear box and horizontalaxis. The simulation is performed to obtain a degradation process in the transmission shaft in a short period of 30 minutes (t = 1800s), where K s has a value of 1e8 and α was arbitrarily chosen considering the suggested range presented in [START_REF] Hunt | Coefficient of restitution interpreted as damping in vibroimpact[END_REF] with a value of 0.5.

For this turbine, the power coefficient curve, C p versus λ, is presented in Fig. 4, where the value of C pmax that can be obtained is equal to 0.4615 at λ 0 equal to 6.4.

Thus, the optimal feedback control gain [START_REF] Berglind | Wind turbine control with active damage reduction through energy dissipation[END_REF], will be K c = 9.5065e5. Nevertheless, one part of this work aims to illustrate the performance of the wind turbine sub-optimal control gains. For this reason, two additional scenarios were proposed to complete the analysis:

• System controlled at 10% below the optimal K c • System controlled at 10% above the optimal K c .

C. Results and discussions

The dynamical system (1) was used to simulate the variations in the torsion angle in the transmission shaft by considering an optimal control for three wind speed cases, see Fig. 5. The obtained torsion angles have a small magnitudes but depending of the wind conditions, the speed of the fluctuations are more or less important.

Figure 6a shows the dissipated energy for laminar flow with optimal control gain and it can be compared to the generated one, see Fig. 6b.

Figure 2 shows the generated power versus rotational speed curves for different control gains. As expected, the generated power is a cubic function of the rotor speed. In Fig. 7: Generated power with a sub-optimal control gain for turbulent wind conditions. the case that the control gain is 1.1K c , the produced power is greater than that produced by the optimal one and for 0.9K c , at same rotor speeds.

Figure 7 illustrates the effect of the turbulence in the produced energy. When the turbine is submitted to a lower wind turbulence it can produce more power than in the case of wind turbulence with greater variance.

A complete analysis can be made by considering both the dissipated and generated energy, and comparing with respect to the optimal control scenario:

• Concerning the dissipated energy, note that if the control gain is below the optimal point, the magnitude is smaller than the other two cases. Nevertheless, when the wind flow has a more significant variance (turbulent), the dissipated energy increases notably for a higher gain and decreases considerably for a smaller gain, see Fig.

8.

• Regarding the generated energy, as expected, the amount of generated energy is lower than the optimal one in all cases of wind conditions. Nevertheless, a control gain higher than the optimal one always leads to an energy production higher than a control gain smaller than the optimal one. Note that the difference in the energy production between the tuning 0.9K c and 1.1K c is reduced when the flow is more turbulent, see Fig. 9.

Note that the amount of energy that can be dissipated in a short period when the wind conditions are in high turbulence is significant. As a result, the degradation in the system can be accelerated, increasing the probability of a failure in the transmission.

IV. CONCLUSIONS

This paper proposes a novel model of shaft degradation based on dissipated energy. The proposed model is based on contact mechanics and allows us to estimate through simulation the dissipated energy at the shaft for different wind conditions and different control gains (optimal and suboptimal). The proposed model was tested using real data of wind speed measurements concerning laminar flow, but also using simulated data of turbulent wind conditions obtained from stochastic differential wind models. The presented simulation provides a complete panoramic result about the possible situations affecting the turbine degradation.

The simulation results illustrate the impact of persistent variations in the shaft angle when the system is submitted to wind speed with high variances. Additionally, if the system works in sub-optimal control conditions, the results shown two possible situations: more generated energy but paying the price of more degradation, or less dissipation energy (in by consequence less degradation), but providing less generated energy. In all the cases, the optimal control gain always provide the maximum generation power, with a "nominal" degradation. This work is a first step towards a degradationaware control approach that would allow to find dynamically the optimal trade-off between the generated energy and the turbine degradation (dissipated energy) taking into account the actual wind conditions.

Fig. 1 :

 1 Fig. 1: Drive-train with a flexible shaft representation.

Fig. 2 :

 2 Fig.2: Generated power with respect to rotor speed produced by optimal and sub-optimal torque controllers.

Fig. 3 :

 3 Fig. 3: Considered wind speed conditions: (a) laminar and (b), (c) turbulent ones.

Fig. 4 :

 4 Fig. 4: Considered power coefficient C p versus tip-speed ratio λ.

  Turbulent flow with high variance.

Fig. 5 :

 5 Fig. 5: Obtained torsion shaft angle for different wind conditions.

Fig. 6 :

 6 Fig.6: Generated and dissipated shaft energies for optimal feedback control.

  Turbulent with low variance.

  Turbulent with high variance.

  Turbulent flow with high variance.

Fig. 8 :

 8 Fig. 8: Dissipated energy at the flexible shaft versus time for different wind conditions.

  Turbulent flow with low variance.

  Turbulent flow with high variance.

Fig. 9 :

 9 Fig. 9: Generated energy versus time for different wind conditions.

TABLE I :

 I Nomenclature

	Symbol	Units	Physical meaning
	Ks	[kg/s 2 ]	Stiffness of the transmission
	Bs	[kg/s]	Damping of the transmission
	Ir	[kgm 2 ]	Inertia of the rotor
	Ig	[kgm 2 ]	Inertia of the generator
	τr	[N m]	Aerodynamic torque
	τc	[N m]	Generator torque
	ωr	[rad/s]	Rotor speed
	ωg	[rad/s]	Generator speed


Adrien Laviron
email: adrien.laviron@csnsm.in2p3.fr

Valentin Gourlaouen

Clarisse Hamadache

Corentin Hiver

Jürgen Kiener

Jean Peyré

Vincent Tatischeff

Optimization of CeBr 3 position-sensitive calorimeter module

Keywords: Gamma--Ray Astronomy, Compton Imaging, Artificial Neural Network, Scintillating Calorimeter, Cerium Bromide, Position Reconstruction, Spatial Resolution

For the next generation of MeV range gamma-ray telescopes, position sensitive calorimeters based on a monolithic scintillator coupled to a pixelated photodetector could be an important building block. In this paper, we present the optimization of the position reconstruction algorithms using machine learning, for a detector based on a 51 × 51 × 10 mm 3 CeBr 3 crystal. For that purpose, we used an automated test bench and collimated radioactive sources to generate experimental data of known energy and position by irradiating the detector with gamma rays. We found in these data different gamma-ray interaction morphologies for which position reconstruction algorithms perform differently, and we developed an algorithm to automatically classify them. We also conducted an extensive optimization of the artificial neural networks that perform the 3D position reconstruction using the Keras Python library with Theano backend.

We found that at 662 keV, 90% of events have a morphology that facilitates position reconstruction. The optimised position reconstruction algorithms give for those events a rms error in the plane of the detector of 1.8 mm on each axis.

The rms error in the depth of the crystal is found to be 2 mm.

Introduction

Despite a large number of existing science topics, from nucleosynthesis and nuclear spectroscopy to multi-messenger astronomy and the extreme universe, gamma-ray astronomy in the 100 keV to 10 MeV energy range is understudied due to the lack of sensitivity of current generation instruments. The best concept for next generation instruments is thought to be a space-borne Compton imager (with pair production imaging capabilities up to 100 MeV) [START_REF] De Angelis | The e-ASTROGAM mission. Exploring the extreme Universe with gamma rays in the MeV -GeV range[END_REF]. Constellations of nano satellites are also envisioned, because that concept would allow all-sky monitoring and provide a powerful insight in high-energy transient sources (see e.g. [START_REF] Werner | CAMELOT: Cubesats Applied for MEasuring and LOcalising Transients mission overview[END_REF] and references therein).

Compton imaging relies on the kinematics of Compton scattering which is the dominant process of interaction of gamma rays with matter in this energy range. Imaging gamma-rays using this principle was first proposed for astronomy in the MeV range [START_REF] Schonfelder | A telescope for soft gamma ray astronomy[END_REF] and nuclear medicine [START_REF] Todd | A proposed γ camera[END_REF] in the 1970's. In an instrument as proposed in [START_REF] De Angelis | The e-ASTROGAM mission. Exploring the extreme Universe with gamma rays in the MeV -GeV range[END_REF], the incoming photon undergoes a scattering in one or several layers of position-sensitive silicon strip detector before being absorbed in a position-sensitive calorimeter based on inorganic scintillators. The measurement of both positions and energy deposits provides information on the photon's source direction [START_REF] Zoglauer | First light for the next generation of Compton and pair telescopes : Development of new techniques for the data analysis of combined Compton and pair telescopes and their application to the MEGA prototype[END_REF].

The localization accuracy and detection sensitivity achievable with this technique are also of interest for the localization of gamma-ray emitting radioactive material on Earth. Thus our results are also relevant for the ComptonCAM project [START_REF] Tatischeff | ComptonCAM[END_REF], which aims at the development of an ultra sensitive, wide-field gamma camera for localization and characterization of low-intensity radioactive waste. This study follows [START_REF] Gostojic | Application of artificial neural network in 3D imaging with lanthanum bromide calorimeter[END_REF] and [START_REF] Gostojić | Characterization of LaBr 3 :Ce and CeBr 3 calorimeter modules for 3D imaging in gamma-ray astronomy[END_REF] on the development of a position-sensitive calorimeter module based on Cerium Bromide inorganic scintillator. In this paper we first introduce the experimental setup used to perform this study. We then present the characterization of the datasets, before explaining the methods used to reconstruct the position of interaction and discuss the results.

Experimental setup and measurements

Cerium bromide calorimeter module

The calorimeter module used in this study is based on a square 51 × 51 × 10 mm 3 Cerium Bromide (CeBr 3 ) inorganic crystal (manufactured by SCIONIX company) coupled to a multi-anode photomultiplier tube (MAPMT). We then have with this setup a pixelated gamma-ray detector. Six calorimeter modules of this type have been characterized.

Because of the hygroscopic nature of the scintillating crystal used, it is embedded in an air-tight aluminium casing. One side of the crystal is coupled to a 1.5 mm thick quartz optical window, itself coupled to the MAPMT using optical grease. The other sides of the crystal are wrapped in an optically diffusive layer of Teflon. The MAPMT used is a square, two-inches, 8 × 8 anodes Hamamatsu H12700A. Its photocathode side is 48.5 mm. It is read by a self-triggered electronic system ROSMAP (developped by IDEAS company), using four 16-channels ASICs and an FPGA to output the data in an Ethernet compatible format. Figure 1 shows a photography of the three components of the module.

In order to use this calorimeter module for Compton imaging, we need both a position and an energy measurement of each gamma-ray interacting in the crystal. We define an event as one or several interactions of a single gamma ray that deposits enough energy to trigger the read-out electronics.

For each event, the MAPMT measures the scintillation light distribution that enables the reconstruction of the gamma-ray's position of first interaction.

The integral of the scintillation light distribution provides the energy deposit.

For the interaction position measurements, the key idea of this setup is to analyze the shape of the scintillation light distribution on the 64 pixels (64 charge values collected over the 64 anodes). Figure 2 shows an example of recorded event. From the shape of the light distribution we expect to determine the 3D coordinates (x,y (surface), z (depth)) of the gamma-ray first interaction in the crystal. Position reconstruction is based on machine learning algorithms, and will be discussed extensively in this paper.

Test bench

To assess the position reconstruction ability of the calorimeter module, we developed a test bench constituted of a radioactive source and a calorimeter module. The radioactive source is collimated to generate a gamma-ray beam directed towards this module.

For this study we used a combination of two 7.6 MBq 137 Cs sources emitting 662 keV gamma rays, and a 424 kBq 241 Am source emitting 59.5 keV gamma rays. In addition, we used two uncollimated sources, one 2.4 MBq 133 Ba source emitting 81 keV and 356 keV gamma-rays, and one 37 kBq 207 Bi source emitting 570 keV and 1064 keV gamma-rays. The 137 Cs source is collimated using a 10 cm outer diameter, 10 cm long tungsten collimator for which the inner diameter can be changed. A good compromise between localization precision and detected source activity has been found at 1.3 mm inner diameter. The 241 Am source is collimated using a 10 cm outer diameter, 1 mm inner diameter, 1 cm long copper collimator. The localization error associated with these collimators width contribute marginally to the calculated position reconstruction errors. It will be included in all results presented hereafter. or by using a full tungsten cylinder instead of the collimator. In order to reduce the background activity, the test bench is shielded by 4 cm thick lead bricks below it as well as on its four sides. A 5 mm copper plate can be added on the top of the calorimeter module to further reduce background at low energy.

Background radioactivity measurements can be made by removing the source

The calorimeter module can be moved automatically to any position in a plan orthogonal to the beam direction by two worm gears driven by stepper motors controlled by an Arduino board. The Arduino is controlled through a serial port emulation by a LabVIEW program that also record the data from the calorimeter module. The module can also be rotated by 90°for the beam to irradiate either its front or its side. Figure 3 shows both irradiation modes. Figure 4 is a photography of the test bench in side irradiation configuration.

These properties of the setup aims at providing events of known coordinates and of known energy to perform the calibration of the calorimeter module. We call hereafter the coordinates of the beam the "mechanical positions". The coordinates system used in this study sets the origin in the center of the crystal.

The xy plane spans the front of the detector, and the z axis represents the depth of interaction (see Fig. 3). These coordinates are defined in the frame of the scanning bench by measuring the count rate in 0.5 mm steps to find the edges of the crystal.

Data characterization

Spectra and energy selection

Figure 5 shows a spectrum of a non collimated 137 Cs source aquired with this calorimeter module, obtained by simply summing the ADC values read by each anode. On the spectrum we can see the full-energy peak near channel 150, 000 corresponding to 662 keV, the Compton front below channel 120, 000 and the backscattering peak near channel 45, 000. Since we are designing a calorimeter, we will mostly be interested in the full-energy events for reference. However, simply selecting events in a given energy window is not optimal for two reasons The maximum of full-energy peak can be noticeably shifted from a mechanical position to another, as shown in Figure 6 A broader than necessary energy selection allows more background radioactivity events in the datasets For these reasons, a finer energy selection, depending on the mechanical position, has been performed on the datasets. For each mechanical position, a gaussian function is fitted to the full-energy peak and events of energy inside a 2σ energy window are selected.

Corrections to improve spectral resolution

The energy resolution obtained by simply summing the charge of the 64 anodes is about 8.3% at 662 keV. This unexpected high value for a CeBr 3 scintillator is due to the position-dependent full energy peak spectral shift shown in Figure 6.

That effect has already been noticed at a smaller scale on smaller scintillators [START_REF] Ulyanov | Localisation of gamma-ray interaction points in thick monolithic CeBr 3 and LaBr 3 :Ce scintillators[END_REF]. It can be explained by the loss of scintillation photons at the interface between the crystal and the reflector. Because of absorption and diffusion (studied in detail in [START_REF] Gostojić | Characterization of LaBr 3 :Ce and CeBr 3 calorimeter modules for 3D imaging in gamma-ray astronomy[END_REF]), the collected charge is higher near the center than close to the edges, and therefore the measured energy and energy resolution are lower there. This shift in measured energy can be quantified by fitting, for each mechanical position, a gaussian function to the full-energy peak and retrieving the position of the center of each fitted gaussian. These fitted energy shifts are shown on Figure 7 for two different energies.

The gains differences between the 64 anodes of the multi-anode PMT can also explain part of these shifts, since those gains can vary by a factor of 2. However, this effect is found to be less important than the former, as shown on Figure 7, on which the spectral shifts are represented before and after the gain correction. It shows that the gains have a large influence on the asymetry of the detector response, but not on its non-uniformity. This correction will not be applied in the following, unless it is explicitly stated.

Let's now consider a reference mechanical position, the center of the crystal for example, so that the spectral shift is defined as a shift with respect to this position's spectrum. The shift can now be expressed as the ratio of the full-energy peak's maximum at the reference position to the full-energy peak's maximum It should be noted that this method takes into account part of 145 the anode gain differences, even though the correction is not applied. Indeed, the scintillation light distribution can be modelled as a continuum and a peak (see Figure 2). The anode gain correction has no effect on the contribution of the continuum to the measured energy, and the peak position is largely determined by the position of inter- An example of this method's results is shown on Figure 8, where the spectral resolution improves from 8.3% to 4.7% at 662 keV. The same correction with the same coefficient gives good results at other energies, see e. g. Figure 9 for the 570 keV and 1064 keV lines of 207 Bi for which we had respectively 8.3% and 7.8% FWHM spectral resolution and attain respectively 4.8% and 4.4%. In this last case, this method is probably limited by our ability to precisely reconstruct the position of interaction at high energies (see Section 4.2).

Morphology of events

Data analysis by visual examination of the events' scintillation light distributions showed a variety of morphologies. Those morphologies can be classified in three main categories as illustrated in Fig. 10. The impact of these morphologies on position reconstruction will be discussed in Section 4.2.5.

The first category presents a single, localised energy deposit. The second category presents two, identifiable and localised energy deposits. The third category presents an energy deposit that is wide and continuous. It should be noted that this classification discretizes a continuous physical phenomenon and therefore that some events are borderline cases. A procedure to automatically recognise events of the first category has been developed. Its core idea is to assess if a significant amount of light is spread or clustered. For that purpose we choose, amongst all pixels, those that have received the most light. The number of pixels considered is a parameter of this algorithm. This selection is made after applying to the pixels the correction of the gains of the anodes. Then we define a shape (typically a square) that we consider acceptably clustered. This is also a parameter of the algorithm. Finally, we check that the considered pixels can fit in the shape.

This check relies on a two-dimensional convolution of the scintillation light distribution, masked with chosen pixels, by the shape. If such a position exists, the event is selected.

The event selection procedure uses this algorithm two times sequentially.

First, we select the events for which the 4 considered pixels are found in a 2 × 2 pixels square. The events that are not selected by this first selection go through a second round, during which are selected those for which the 5 considered pixels can be found in a 3 × 3 pixels square. This two-stage selection is necessary because the second round will select events that have been rejected by the first one but have a shifted or truncated scintillation light distribution, especially on the borders of the crystal. But it should be noted that events selected during the first round could have been rejected by the second. Indeed, choosing too many pixels includes in the analysis pixels from the light continuum that have the highest values because of photon noise, discarding events that have an acceptable morphology. This procedure selects between 89% and 94% of events (depending on the calorimeter module characterized) in datasets recorded with a 137 Cs source and selected by the fine energy selection around 662 keV described in Section 3.1, and between 85% and 93% of events selected the same way around the 356 keV line of a 133 Ba source. This result shows that most of events at these energies have a localised energy deposit in the crystal. This procedure is meaningless for events below a certain energy because not enough scintillation light is produced. At low energy, as shown on Figure 11 with an event of energy 59.5 keV, most events will have scintillation light distribution compatible with a single, localised energy deposit, because gamma-rays interact in CeBr 3 mostly via photoelectric absorption.

Position reconstruction

Artificial neural networks

An artificial neural network (ANN) is a network of simple algorithms that aims at mimicking the operation of a biological neuron. Those algorithms execute a weighted sum of several inputs, and apply to this number an R → R function called activation function. We use ANNs to reconstruct the position of the first interaction of a gamma ray from the measured scintillation light distribution. The general idea of this position reconstruction method has been proposed in [START_REF] Delorme | Use of a neural network to exploit light division in a triangular scintillating crystal[END_REF] for positron emission tomography scanners.

The machine learning used in this study is supervised learning. It means that ANNs are first trained to output the correct answer, using input data with known output, by minimizing an error function. For our case, the input data are the 64 channels of our MAPMT, the output is the real position of interaction of the gamma ray in the crystal. The mechanical positions approximate the known outputs. We choose for our error function a mean square error. Once ANNs used in this work are multi-layer perceptrons: neurons are organised into one or several layers, that are executed sequentially. Such non-cyclic neural networks have been described in 1943 in [START_REF] Mcculloch | A logical calculus of the ideas immanent in nervous activity[END_REF]. The number of layers and their size (number of neurons) are the parameters of the ANN chosen by the programmer. Layers outside of the input and output ones are called "hidden layers". The activation function is also such a parameter, as well as the algorithm used for training. These parameters define the neural network's architecture.

During training, the neurons weights are adjusted by iterations, called epochs [START_REF] Ball | Data Mining and Machine Learning in Astronomy[END_REF], to minimize the error from the reconstructed gamma-ray interaction coordinates to the mechanical position. After a certain number of epochs, this error decreases without improving the general ANN's performance. This phenomenon is called overlearning. In this case the neural network tends to reconstruct better the events from the training dataset, but looses performance for any other event. To prevent overlearning, the error function is also calculated on a second dataset, not used for the training and called validation dataset. The training is stopped whenever this second error stops decreasing. All this work has been done using the Python Keras framework [13] with Theano [START_REF] Theano Development | A Python framework[END_REF] backend.

In this part we will present the effect of various event selection algorithms and neural network architectures on the performance of position reconstruction, first in 2D in the xy plane (see Fig. 3), and then in the depth of the crystal along the z axis. Gamma rays interacting in the border of the crystal are more likely to escape after a Compton diffusion. So, for a same number of triggers in the detector, the number of events selected by the fine energy selection can be twice less. In that case, the weight evaluation procedure during the neural network training would be biased, underestimating the importance of near-border events. This would result in degraded position resolution near the borders of the detector.

Therefore, the datasets we used contain a fixed number of events per mechanical position, typically 50. This allows for a more uniform position reconstruction performance as well as contained training durations. The morphology selection described in Section 3.3 is only applied once the ANN has been trained and is used for the position reconstruction, since we found that applying it before the ANN training yielded no significant improvement to the ANN performance.

Performance indicators

This study also investigates the effects of the ANN's architecture on the position reconstruction performance. For the xy reconstruction, the performance of an ANN is measured by way of its mean uncertainty σ XY , defined as the quadratic average on all mechanical positions p of the σ XY,p ,

σ XY = 1 N p N p σ 2 XY,p (1)
where N p is the number of events at mechanical position p. The quantity σ XY,p is defined as the square root of the mean squared distance from the reconstructed x and y coordinates (in the frame defined on Fig. 3) to the known ones,

σ XY,p = 1 N p Np i=1 (x p -x rec i ) 2 + (y p -y rec i ) 2 , (2)
where N p is the number of events at the mechanical position p, x p (respectively y p ) is the known x (resp. y) coordinate of the mechanical position p and x rec i (resp. y rec i ) is the x (resp. y) coordinate of event i reconstructed using the ANN.

In order to quantify any systematic offset between the known position and the reconstructed position, an offset error σ o XY,p has been defined as

σ o XY,p = (x p -x) 2 + (y p -y) 2 , (3)
where x (resp. y) is the mean of the values of the x rec i (resp. y rec i ) of events at this mechanical position (see Section 2.2).

The σ XY,p can be expressed as the quadratic sum of σ o XY,p and a point spread error

σ s XY,p = 1 N p Np i=1 (x rec i -x) 2 + (y rec i -y) 2 , (4)

Exploration of ANN's architecture parameters space

The analysis has been performed for all combinations of the following archi- Such output has been studied also for the mean point spread error σ s XY (see Eq. 4) and the mean offset error σ o XY (see Eq. 3) as shown on Figure 13. This analysis shows that the offset error is systematically lower than the point spread error. Some architectures, using elu or softplus activation functions and more than 16 neurons per hidden layer, perform better than the others on this test.

Moreover, some architectures such as those using the softsign display a point spread error better than most, at the cost of a larger offset error. Architectures using more neurons do not seem to overlearn because of the validation dataset monitoring. However, they require more calculations, wich may be of critical importance depending on the application.

Energy dependency

This study has been performed in two different energy ranges, one around the full-energy peak of a 137 Cs source at 662 keV and one around the 59.5 keV line of a 241 Am source. In this second energy range, the same fine energy selection has been applied. No morphology selection was made as discussed in Section 3.3.

To assess the energy dependency of neural networks, the same analysis has been conducted with 50 events per mechanical position, recorded at energies near 662 keV and near 59.5 keV. This study showed that using two neural networks, each trained with dataset from a single 59.5 keV or 662 keV energy window, does not perform significantly better than using a single neural network trained with a dataset composed with both energy windows.

Finally, we conclude that, for the xy position reconstruction and the set of architectures studied here, there is no significant improvement over 2 layers and ≈ 19 neuron per hidden layer (see Fig. 12). The best training algorithm for our purpose appear to be adamax, and both elu and softplus activation functions give good results. In the following, we use one ANN for all energies and the elu Colorscale represents the performance of the ANN σ XY in mm (lower is better, see Eq. 1).

(colour online) We see that most gamma rays are reconstructed within a short distance of their associated mechanical position, and with a good consistency throughout the detector. The events at mechanical positions near the border of the scintillating material are harder to reconstruct because their scintillation light distribution is truncated or misshapen. We still obtain σ XY = 5.0 mm on that dataset.

Some events are reconstructed very far from their associated mechanical position. They are interpreted as background radioactivity. It should be noted that all calculations of the ANN performance σ XY (as well as σ s XY and σ o XY ) using Eq. 1 (respectively Eq. 4 and Eq. 3) include some events from the background radioactivity. To assess its contribution to the position resolution, we recorded a dataset with the 137 Cs source at a single mechanical position (x = -3.3 mm, y = -2.6 mm) with a 1 mm inner diameter collimator. The background radioactivity has been estimated by placing a thick tungsten absorber between the source and the detector, as described in Section 2.2, and 335 recording data for the same duration (two hours). Events were selected using a large energy window between ≈ 580 keV and ≈ 700 keV. A histogram of the reconstructed gamma-ray positions of interaction can be seen on Figure 15. The blue lines represent the reconstructed coordinates of all events, the green lines the reconstructed coordinates of morphologically selected events, and the red 340 lines the background radioactivity. The blue and green lines meet the red one at large distance from the position of interaction, which validates the interpretation that the events that are reconstructed far enough from the known position of interaction are background radioactivity.

When the background radioactivity is subtracted bin to bin to the data in the histogram, a new estimation of the performance of the ANN σ XY,p can be Morphology selection

σ X,p σ Y,p σ XY,p σ X,p σ Y,p σ XY,p
Not applied 3.6 mm 3.5 mm 5.0 mm 2.9 mm 2.6 mm 3.9 mm Applied 3.0 mm 2.9 mm 4.2 mm 1.9 mm 1.7 mm 2.6 mm

σ XY,p = 1 N n bin=1 N bin × (x known -x rec bin ) 2 + (y known -y rec bin ) 2 ( 5
)
where N is the number of events in the dataset after subtraction, n the number of bins in the histogram, N bin the number of events in a given bin, x known and y known the mechanical position of the dataset (x = -3.3 mm, y = -2.6 mm),

x rec bin and y rec bin the coordinate of the center of the bin bin. It gives, at this position p, σ XY,p = 3.9 mm for all events and σ XY,p = 2.6 mm for morphologically selected events. The results are presented in Table 1 for each axis. They point at an average localization error of 1.8 mm. These σ XY,p and σ XY,p should not be interpreted as the standard deviation of a gaussian distribution. Figure 16 shows the histograms of the difference from the reconstructed to the known coordinates (as in Figure 15 but in linear scale), in the dataset used to evaluate σ XY , together with a gaussian fit that provides a poor description of the data, and especially of the distribution tails. The standard deviations of the fitted gaussians are ≈ 1.55 mm at 662 keV and ≈ 1.7 mm at 59.5 keV for both the x and y axes.

As shown on Figure 17, the variations of the position resolution σ XY,p (for which background radioactivity has not been subtracted) is dominated by the loss of performance in the ≈ 5 mm near the borders of the crystal.

Depth of interaction reconstruction

Depth of interaction reconstruction is the reconstruction of the position of interaction along the z axis (see Fig. 3). Since the depth of interaction is unknown when irradiating the front of the crystal, acquiring data that allow neural networks training requires new procedures, which will be presented in Section 4.3.1. New performance indicators will be defined (Section 4.3.2) to be used in the other architecture's parameter space explorations that have been conducted, presented in Section 4.3.3. Finally, the results of artificial neural networks will be compared against methods based on the known absorption length of radiation in the crystal in Section 4.3.4.

Data acquisition and selection

For that part of the study, in order to know the z coordinate of the interaction, the detector is irradiated from the side (see Fig. 3 in Section 2.2).

We use 662 keV photons from the 137 Cs source to be able to explore all the crystal, since gamma rays below ≈ 300 keV do not penetrate enough in the crystal. We then obtain data of different energies by selecting events in the full-energy peak and in the Compton plateau (see Fig. 5). In the Compton plateau, events mostly result from a single interaction of a gamma ray in the crystal. Therefore, selection on event morphology is not needed. Since the gamma rays have an initial energy of 662 keV, absorption is no longer a problem. This method allows to probe all x, y and z positions in the crystal with known z position at various energies. However, this method approximates the scintillation light distribution of a lower energy photon depositing all its energy by the scintillation light distribution of a Comptonised electron of the same energy. The difference should be negligeable, given that the mean free path of the generated low energy electrons in the CeBr 3 is < 0.1 mm and that the photodetector pixels are 6 mm wide. Three energy ranges were defined in addition to the full-energy peak, around 300 keV between the Compton front and the backscattering peak, around 130 keV below the backscattering peak, and near 50 keV. In those energy ranges, no fine energy selection (as defined in Section 3.1) could be applied because there is no peak to fit. No selection on event morphology was performed because the selection on event energy already selects single interaction events.

For each of these energy ranges, the dataset used for training spans the whole yz plane of the detector with a 1 mm step. The gamma-ray absorption in CeBr 3 at these energies is such that events may interact at any x coordinate for each yz known coordinates. After the energy selection (and fine energy selection in the energy range near 662 keV), 200 events are selected for each mechanical position, out of which 30% will be used for validation.

Performance indicators

To assess the performance of different methods on the position reconstruction along the z axis, a new performance indicator σ Z has been defined as the quadratic mean on all mechanical position p of all σ Z,p , themselves defined as the one-dimensionnal square root of the mean squared distance,

σ Z,p = 1 N p Np i=1 (z p -z rec i ) 2 , (6)
where z p is the known z coordinate of the mechanical position p and z rec i is the z coordinate reconstructed using a given method.

Just as σ s XY and σ o XY were defined in Section 4.2.2 (see Eq. 4 and Eq. 3), we define a mean point spread error σ s Z and a mean offset error σ o Z as the quadratic mean on all mechanical positions p of the

σ s Z,p = 1 N p Np i=1 (z rec i -z) 2 (7)
σ o Z,p = |z p -z| (8)
where z is the mean of the z rec i values for all events recorded at given mechanical position p (see Section 2.2).

Even though the y coordinate is known for each of these events, it will not be reconstructed here. For the neural network, first investigations on the reconstruction of the depth of interaction showed that asking for the reconstruction of both y and z coordinates led to poor results in z reconstruction, because y is both easier to reconstruct and spans a larger range of values, so the training algorithm tends to prioritize that coordinate. Therefore, the y coordinate will not be used in this part.

Exploration of ANN's architecture parameters space

Given the weak dependency of the performance of neural networks on the total number of neurons, shown by preliminary tests, this part of the analysis focused on the activation function and training algorithm optimization. The analysis has been performed for the following set of architecture parameters:

4 hidden layers 80 neurons per hidden layers 10 activation functions [13]: sigmoid, hard sigmoid, tanh, elu, relu, softplus, softmax, softsign, exponential, and linear

training algorithms: adam, nadam and adamax

For our purpose, the training algorithm that performed the best was adamax.

The most performing activation functions were elu and relu. We will use elu in the following. The energies at which depth of interaction is best reconstructed are around 662 keV and 300 keV, with σ Z ≈ 2 mm. In the 130 keV energy range, ANNs are significantly worse with σ Z ≈ 2.4 mm, and in the 50 keV energy domain they perform even worse with σ Z ≈ 2.6 mm. This could be explained by the lower number of scintillation photons generated by the interaction, which gives a noisier measured scintillation light distribution, as shown in Figure 11. One should also note the presence in selected events of source X-rays from the cesium Kα line at 31 keV, disturbing the ANN training.

Depth of interaction reconstruction at low energies

At energies lower than 100 keV, the attenuation length of gamma rays in CeBr 3 is lower than 2 mm, therefore using a neural network that provides an uncertainty on position reconstruction larger than 2 mm is suboptimal. Instead, two main methods could be used: the simplest one consists in assigning to every event, independently of its scintillation light distribution, a fixed z coordinate. Below ≈ 100 keV, an appropriate coordinate is z = -5 mm, that is the coordinate of the entry face of the crystal. A more complex method is to assign to every event a depth of interaction that depends only on its measured energy, and equal to the mean of several simulated depth of interaction in the crystal. This simulation approach allows to easily account for the finite size of the crystal. Depths are calculated using the NIST XCOM photon cross-section database data [START_REF] Berger | Xcom: Photon cross section database (version 1.5)[END_REF].

Figure 22 represents the σ Z error (see Eq. 6) achieved with different depth of interaction reconstruction methods. Apart from the ANN results, these data are simulated with a Monte-Carlo method. The green and black curves show the results of assigning to every event the coordinate z = -5 mm and z = 0 mm (center of the crystal), respectively. The choice of z = 0 mm is better at high energies because then the absorption is not dominant. The red line represents the error achieved by the "mean depth of interaction" method. At high energies, since the mean goes to z = 0 mm, we find back the results of this method. At low energies, it is always better than the z = -5 mm approach because it uses a better model of the reality. We note that these methods rely on the assumption that a low energy deposit is due to a gamma ray coming from the front of the detector and being absorbed, which may or may not be the case depending on the situation. Finally, neural networks show a better performance than any other method above ≈ 150 keV.

Results and discussion

We studied the performance of a gamma-ray detection module made of a CeBr 3 monolithic scintillating crystal coupled to a 64-channels MAPMT.

We developed algorithms to optimize the spectral resolution and 3D position resolution.

We obtain a spectral resolution ehancement down to 4.7% FWHM thanks to a position-dependent energy correction, as presented in Section 3.2. This correction does not require a strict calibration of the ADC to charge relation nor a correction for the different gains of the 64 anodes, and it is fairly lightweight from a computational perspective once the gamma-ray position of interaction is reconstructed.

The selection on event morphology allows a better position reconstruction as shown on Figure 15. It also leads to a better understanding of the data. A 662 keV gamma ray in CeBr 3 will in most cases interact via Compton scattering. However, at this energy, photons are mostly scattered in the forward direction. This may explain why a majority of events shows a localised energy deposit resulting of several, closely packed gamma-ray interactions.

It can also explain the large continuous energy deposits. The third identified morphology, showing several distinct energy deposits, can then be interpreted as several Compton interactions. The fractions of selected events at 662 keV presented here varies between 89% and 94%, except for one module for which 85% of events were selected. At 356 keV, this fraction varies from 85% to 93%, and is 81% for this one module (detailed results are presented in Table 2).

The behavior of that module may be explained by a small anomaly in the optical coupling between the crystal housing and the PMT. The preferred forward scattering depends on the energy, so this variation of the fraction of selected events was expected. Future work could further investigate this variable. This event selection based on morphology allow a position resolution improvement from σ XY,p = 3.9 mm to σ XY,p = 2.6 mm at 662 keV. It may or may not be used in the context of a Compton space observatory, depending on the expected flux of the source.

Investigations on the data that were discarded for their morphology shows that neural networks probably cannot reconstruct them correctly. Indeed, when neural networks are trained using only these events, all tested architectures output σ XY > 7 mm. However, neural networks trained only with events having that unaccepted morphology and a maximum of the scintillation light distribution located near the mechanical position give great results, down to σ XY = 4 mm (including background radioactivity events). This cannot be implemented in an actual Compton imager, since the real position of interaction is not known. But from that analysis, we can assume that neural networks reconstruct these events at the position of the highest energy deposit, which may or may not be the position of first interaction.

The position reconstruction in the xy plane has been optimised down to an uncertainty of 1.8 mm for all tested energies. Following investigations showed that this method can be generalized to other calorimeter modules with the same geometry and electronics, even though the gain of individual PMT anodes may vary by a factor 2. Quantitative results for the six modules that have been characterized are reported in Table 2. We found the same ANN architecture to perform very well for all calorimeter modules: 2 hidden layers, 19 neurons per hidden layer using the elu activation function and trained using the adamax algorithm. However, that requires a new neural network training with data from that module.

Preliminary investigations showed that a single neural network could also be used with several calorimeter modules, if it has been trained with a combined set of data from each of those modules, with minimal performance loss, even though the individual anode gains of the PMT were not corrected in this study. At 662 keV, the neural network dedicated to the module 1 yielded σ XY = 4.57 mm, and the one dedicated to the module 2, σ XY = 4.50 mm, while the neural network trained with data from both modules gave for each module, respectively, σ XY = 4.56 mm and σ XY = 4.52 mm. For comparison, neural networks that were trained with a single module give consistently σ XY > 6 mm on the data from other modules. Future work will investigate further this result, as it may have a significant impact on the dataflow design of future instruments. Future work could also investigate the ability of such software to reconstruct the position of interaction in thicker crystals, since those crystals would have larger gamma-ray absorption ability, which may be useful in a low count rate application such as a space telescope.

The results obtained in position and spectral resolution fulfill the requirements of the calorimeter of next generation space-borne Compton telescopes.

Moreover, it has a position resolution noticeably finer than its pitch (1.8 mm for 6 mm square detectors), reducing the number of channels to read at the cost of some computation. The effect of this tradeoff should be discussed for embedded applications such as space observatories and nano-satellites in which both electrical power and computer resources for data management are limited.

The characteristics of the module we developed make it also an interesting detector for a portable gamma camera aimed at locating and characterising radioactive waste following the dismantling of nuclear facilities [START_REF] Tatischeff | ComptonCAM[END_REF].

Depth of interaction reconstruction is less critical than xy position reconstruction in the context of a Compton telescope, and especially since our calorimeter modules are not very thick (10 mm). However, its importance increases as the scattering detector and the calorimeter get closer. The neural network approach yielded good results at the cost of an increase in ANN complexity for the medium and high energy events (> 150 keV). The optimised z ANN architecture uses 4 hidden layers and 80 neurons per hidden layer using the elu activation function, and is trained using the adamax training algorithm. At low energy, the simple method based on the interaction length of gamma rays in CeBr 3 gives satisfactory results.

Conclusions

For the next generation of MeV range gamma-ray telescopes, position sensitive CeBr 3 scintillating detectors could be an important building block. We studied six detectors using that technology and developed a position-dependent energy correction, an event selection procedure that recognises event morphologies and we optimised the position reconstruction algorithms. We attained an energy resolution of 4.7% FWHM at 662 keV with the energy correction.

The 3D position reconstruction of a gamma-ray interaction relies on two artificial neural networks, one for the plane xy of the detector and one for its These results make that calorimeter module technology a very appealing choice for the next generation of MeV range gamma-ray space borne observatories.

Figure 1 :

 1 Figure 1: Photography (exploded view) of the calorimeter module. From left to right, we can see the scintillating crystal, the MAPMT and the read-out electronics

Figure 2 :

 2 Figure 2: Example of an event recorded by the calorimeter module. Both parts of the figure shows the same data, one in parallel perspective and the other in colorscale. One can read on the z axis or the color scale the quantity of light collected by each pixel of coordinates (x, y).

Figure 3 :Figure 4 :

 34 Figure 3: Scheme of the calorimeter module. The gamma-ray interacts in the CeBr 3 that emits scintillation light. This light is detected by a pixelated photodetector optically coupled to the scintillating crystal. The shape of the scintillation light distribution provides information on the 3D position of the gamma-ray interaction. Top: front irradiation, z in unknown. Bottom: side irradiation, x in unknown

Figure 5 :

 5 Figure 5: Spectrum of a 137 Cs source (662 keV)

Figure 6 :

 6 Figure 6: Spectra of a collimated 137 Cs source at different mechanical positions, (0, 0) being the center of the crystal, (14, 14) an intermediary position and (22, 22) a position close to a corner.

Figure 7 :

 7 Figure 7: Relative detector spectral response through its surface, as given by the position of the center of a gaussian fitted to a spectral line, for 25×25 mechanical positions with a 2 mm pitch. Top: measured for the 662 keV line of 137 Cs. Bottom: measured for the 59.5 keV line of 241 Am. Left: before correction of the variable gain of the 64 anodes. Right: after correction. (colour online)

Figure 8 :

 8 Figure 8: Left: simple spectrum, with spectral resolution of 8.3% at 662 keV. Right: corrected spectrum, with spectral resolution of 4.7% at 662 keV

Figure 9 :

 9 Figure 9: Left: raw spectrum of the 570 keV (top) and 1064 keV (bottom) lines of a 207 Bi radioactive source. Right: corrected spectrum, showing the two lines with respectively 4.8% and 4.4% FWHM spectral resolution.

Figure 10 :

 10 Figure 10: Three events with different morphologies. Data are represented twice, once in parallel perspective (left) and once in colorscale (right). The red mark in the colorscale representation indicates the mechanical position. Top: single, localised energy deposit. Center: two energy deposits. Bottom: wide, continuous energy deposit.

Figure 11 :

 11 Figure 11: Event recorded from a 241 Am radioactive source.

  the training is completed, reconstructing unknown events consists in executing the neural network program with the parameters calculated during the training process. The weights of the neurons are free parameters adjusted during the training of the ANN.

4. 2 .

 2 Optimization of 2D position reconstruction 4.2.1. Datasets The aim of this part is to investigate the performances attainable with ANN based xy position reconstruction. For that purpose, datasets were acquired with the automated test bench. A first dataset, used for training, contains mechanical positions spanning the whole detector with a 2 mm step along both x and y axes. The validation dataset contains mechanical positions every 3 mm, in order to have data at different mechanical positions. Both datasets result from the fine energy selection discussed in Section 3.1.

  Figure 12 shows an output example of such an analysis. Colorscale represents the value of σ XY , in mm, attained after training, calculated on a set of N = 14450 events evenly distributed on the xy plane of the detector.

Figure 12 :

 12 Figure 12: Output example of a systematic ANN's architecture parameter space exploration. On the x axis is represented the number of neurons per layer. On the y axis are indicated the number of hidden layers (HL), the activation function and the training algorithm used.

Figure 13 :

 13 Figure 13: Ouput of an exploration of some architectures for neural networks with 2 hidden layers and adamax training algorithm. Left: colorscale represents σ s XY Right: colorscale represents σ o XY (colour online)

Figure 14 :

 14 Figure 14: Example of a dataset's position reconstruction by our ANN. Mechanical positions change by steps of 4 mm along x and y axes, and 200 events of energy 662 keV were selected for each mechanical position. Events were selected by energy and morphology. Mechanical positions are indicated by the red markers, and gamma-rays reconstructed coordinates shown as colored dots. A single (but not unique) color is associated to a mechanical position. Dotted black lines show the borders of the pixels.

Figure 15 :

 15 Figure 15: 1D histograms of reconstructed x and y coordinates for a single mechanical position x = -3.3 mm, y = -2.6 mm. Left: projection along the x axis. Right: projection along the y axis. The blue dotted lines represent the reconstructed coordinates of events in the dataset acquired with the source. The green solid lines represent the reconstructed coordinates of events from the same dataset that have been selected for their morphology. The red dashed lines represent the reconstructed coordinates of events from the background radioactivity dataset (source behind absorber). The black dashed-dotted vertical lines represent visually the known coordinates.

Figure 16 :

 16 Figure 16: 1D histograms of the difference between the reconstructed and known coordinates. Left: x coordinates. Right: y coordinates. The blue solid lines represent the 662 keV dataset, for which the selection on event morphology has been applied. The green dashed lines represent the 59.5 keV dataset. The red dotted-dashed line show a gaussian fit of the distribution at 662 keV.

Figure 17 :

 17 Figure 17: Position resolution estimator σ XY,p uniformity through the detector x axis, calculated using morphologically selected events at 662 keV with one data point every two millimeters. Blue: at y = 0 mm. Green: at y = 14 mm

Figure 18 :

 18 Figure 18: Histogram of reconstructed depths of interaction at 662 keV. Yellow dotted: at y = 10 mm, z = -4 mm. Blue solid: at y = 0 mm, z = 0 mm. Green dashed: at y = 10 mm, z = 4 mm. Vertical lines show visually the known depth of interaction

Figure 20

 20 Figure20shows the histograms of the difference between the reconstructed and known z coordinate. As energy decreases, the reconstruction of the depth of interaction becomes less precise and more asymetric towards over-estimated depths.The combination of the two optimised ANN allow 3D position reconstruction, as shown on Figure21.

Figure 19 :

 19 Figure 19: σ Z,p (top) and σ s Z,p (bottom) represented in colorscale in mm throughout the yz plane. (colour online)

Figure 21 :

 21 Figure 21: 3D position reconstruction using two optimised ANN for 5 xy mechanical positions.

Figure 22 :

 22 Figure 22: σ Z (in mm) versus the energy of incoming gamma-ray for various methods. ANN's results are shown in blue. The black curve represents the error achieved when all events are reconstructed at z = 0 mm. The green curve represents the error achieved when all events are reconstructed at z = -5 mm. The red curve represents the error achieved when all events of a given energy are reconstructed at their mean depth of interaction.

depth z .

 . The optimised xy ANN architecture uses 2 hidden layers and 19 neurons per hidden layer using the elu activation function. The optimised z ANN architecture uses 4 hidden layers and 80 neurons per hidden layer using the elu activation function. Both were trained using the adamax training algorithm. We found the optimised ANN architectures to perform equally on different calorimeter modules for a wide energy range (from 60 keV to 662 keV). The z position can be reconstructed with a ≈ 2 mm error. The xy position resolution depends on the morphology of events, but a 1.8 mm error has been attained for ≈ 90% of events.

Table 1 :

 1 Summary of the position resolution measured at the position x = -3.3 mm, y =

	-2.6 mm
	calculated for a single mechanical position p as

Table 2 :

 2 Summary of results obtained in xy position reconstruction and morphology analysis for the six modules characterised.

Acknowledgements

This work benefited from the Investments for the future program of the French Government and operated by the French National Radioactive Waste Management Agency (ANDRA).


Tuan Anh Tran

Carine Jauberthie

Louise Travé-Massuyès

Hung Quoc

Lu

published or not. The documents may come

An Interval Kalman Filter enhanced by lowering the covariance matrix upper bound