# Deep Learning reconstruction with uncertainty estimation for $\gamma$ photon
interaction in fast scintillator detectors

G. Daniel<EMAIL_ADDRESS>M.B. Yahiaoui  C. Comtat  S. Jan  O.
Kochebina  J-M. Martinez  V. Sergeyeva  V. Sharyy  C-H. Sung  D. Yvon
Université Paris-Saclay, CEA, Service de Génie Logiciel pour la Simulation,
91191, Gif-sur-Yvette, France  Université Paris-Saclay, Inserm, CNRS, CEA,
Laboratoire d’Imagerie Biomédicale Multimodale (BioMaps), 91401, Orsay, France
Université Paris-Saclay, CEA, IRFU, Département de Physique des Particules,
91191, Gif-sur-Yvette, France

###### Abstract

This article presents a physics-informed deep learning method for the
quantitative estimation of the spatial coordinates of gamma interactions
within a monolithic scintillator, with a focus on Positron Emission Tomography
(PET) imaging. A Density Neural Network approach is designed to estimate the
2-dimensional gamma photon interaction coordinates in a fast lead tungstate
(PbWO4) monolithic scintillator detector. We introduce a custom loss function
to estimate the inherent uncertainties associated with the reconstruction
process and to incorporate the physical constraints of the detector.

This unique combination allows for more robust and reliable position
estimations and the obtained results demonstrate the effectiveness of the
proposed approach and highlights the significant benefits of the uncertainties
estimation. We discuss its potential impact on improving PET imaging quality
and show how the results can be used to improve the exploitation of the model,
to bring benefits to the application and how to evaluate the validity of the
given prediction and the associated uncertainties. Importantly, our proposed
methodology extends beyond this specific use case, as it can be generalized to
other applications beyond PET imaging.

###### keywords:

Deep Learning, Neural Networks, Uncertainty quantification, Event
reconstruction algorithms, Gamma detector, PET Imaging

††journal: Engineering Application of Artificial Intelligence

##  1 Introduction

Gamma photon detection is used in numerous industrial, medical and security
applications. It is often based on scintillating crystals coupled with a light
collection and readout system. The scintillator is either pixelated,
consisting on an array of small individual crystals, or continuous, made of a
large monolithic block. The choice between both technologies is often based on
a trade-off between sensitivity and resolution performances. For pixelated
detectors, the spatial localization of the detection is provided by the
crystal impacted by the gamma photon, whereas for continuous detectors,
designated algorithms shall be used to derive spatial information. Several
algorithms have been proposed, based either on prior knowledge of the physics
of detection or on a machine learning approach [18]. The main objective of
this work is the development of a physics informed deep learning method for a
quantitative estimation of the spatial coordinates of the gamma interaction
within a monolithic scintillator, including uncertainties. The application
framework is nuclear medicine imaging, more specifically the detection of 511
keV gamma photons in Positron Emission Tomography (PET).
PET imaging is a powerful in vivo functional imaging modality mainly used in
oncology, neurology and cardiology. It is based on the administration to the
patient of a biomarker labeled with a radionuclide that decays through
positron emission, followed by the detection in coincidence, outside the body,
of pairs of 511 keV gamma photons resulting from the annihilation of the
emitted positrons with electrons of surrounded media. The data acquisition
process is followed by the tomographic reconstruction of a three-dimensional
image of the biomarker distribution within the body. The quality of the
reconstructed image highly depends on the performances of the gamma photon
detectors in terms of sensitivity, spatial resolution, and temporal
resolution. The detector spatial resolution has a direct impact on the
contrast recovery of small structures in the image. A higher detection
efficiency translates into a higher number of detected coincidences, resulting
in a better signal-to-noise ratio (SNR) in the PET image. The improvement of
the coincidence resolving time (CRT), characterizing the capability of a pair
of detectors to resolve the difference between the times of interaction of the
two 511 keV gamma photons detected in coincidence, also helps for increasing
the SNR in the image [23]. This principle is referred to as time-of-flight
(ToF) PET.
State of the art clinical PET systems are based on pixelated detectors made of
LSO (lutetium oxyorthosilicate) or LYSO (lutetium-yttrium oxyorthosilicate)
crystals, with a pixel pitch between 3.5 and 5 mm. These PET systems have, at
best, a CRT of 210 ps [26], an intrinsic spatial resolution in the
reconstructed image of 3.5 mm [26] (about 3 mm at the detector level), and an
absolute sensitivity of $\sim$ 20 counts/sec/kBq (i.e. $\sim$ 2%) [15]. Major
efforts are being made worldwide in nuclear instrumentation research groups to
improve these parameters, in particular to decrease the CRT below 100 ps,
ideally down to 10 ps [20]. There is a trend toward the development of PET
detectors based on monolithic crystals for their higher sensitivity, since
there are no intercrystal gaps [14]. Based on the scintillation light
distribution readout, dedicated neural networks have recently been proposed to
provide a single 3-dimensional gamma photon interaction position within the
monolithic crystal (see, for example [5, 12, 8, 16, 18]). All these studies
demonstrate the benefit of using a machine learning approach for an accurate
position estimation.
In this study, we propose to address the question of the reconstruction of the
2-dimensional gamma photon interaction coordinates for a fast lead tungstate
(PbWO4) monolithic scintillator detector [27], that is not straightforward due
to the specificity of our acquisition system which is not a conventional
pixelated photodetector. Our system consists indeed in the use of a micro-
channel plate photomultiplier tube (MCP-PMT) which necessitates a dedicated
process to reconstruct gamma photon interaction parameters from the acquired
signal, as we describe in this paper. Moreover, we aim not only to perform the
reconstruction but to associate uncertainties on the reconstruction. We will
exploit and validate these uncertainties to assess the reliability of the
network prediction. This approach is new in this field and the methodology we
propose can be extended to other applications on signal processing or sensor
data analysis. In the future, these spatial coordinates uncertainties could be
used during the tomographic reconstruction and potentially improve the quality
of the PET image.
The paper is organized as follows. In section 2, we describe the gamma photon
detector based on PbWO4, the Monte Carlo simulation used to generate the
training and testing datasets, and the preprocessing of the detector raw data.
Section 3 focuses on the methodology of the deep learning approach used for
this study, including the description of a baseline method for performance
comparisons and the definition and the architecture of the selected Density
Neural Network. Results are presented in section 4. These results and the
methodology are discussed in section 5.

Figure 1:  Left: Schematic diagram of the ClearMind detection module. A 511
keV gamma-ray interaction in the crystal produces scintillation and Cherenkov
photons that are converted by the photocathode to photoelectrons. These
photoelectrons are then multiplied by the MCP-PMT and induce signals on the
transmission lines (TLs). Signals from the left and right ends of each TL are
amplified by 40 dB amplifiers and digitized by a SAMPIC module. Right:
Transmission lines Printed Circuit Board (PCB). The axis $x$ and $y$
corresponds to the coordinate system that we use to locate the interaction
position.

##  2 Materials

###  2.1 Detector description

The ClearMind gamma detector (Figure 1) is composed of a MCP-PMT sealed by a
monolithic PbWO4 crystal, acting both as the gamma conversion crystal and as
the optical window of the MCP-PMT. A high quantum efficiency photoelectric
layer is deposited on its inner face. The direct deposition of a photocathode
with a refraction index superior to the refraction index of the PbWO4 crystal
allows us to avoid total reflection at the crystal/photocathode interface,
thus maximizing the photon collection efficiency of the module [27]. The use
of this ”scintronic” crystal as an entrance window of a MCP-PMT makes it
possible to optimize the time resolution thanks to the excellent electron
transit time spread ($\sim$60 ps FWHM) to the detection anodes provided by
this type of photodetector.
The PbWO4 crystal, homogeneously doped, has a surface of
$59\leavevmode\nobreak\ \mathrm{mm}\times 59\leavevmode\nobreak\ \mathrm{mm}$
a thickness of 5 mm, and is provided by CRYTUR [9]. The photocathode deposit
and the integration of the device into a MCP-PMT structure is handled by the
PHOTEK company, based on its MAPMT-253 design [21]. We developed a signal
readout system for this device using 32 transmission lines as shown in Fig. 1,
[10]. The signals are read out at both ends of the transmission lines,
amplified and then sampled by a SAMPIC WaveShape Recorder [7].
Typically a 511 keV energy deposit in the crystal produces 185 optical photons
mostly isotropically. Out of these, $\sim$ 20% are collected by the
$53\leavevmode\nobreak\ \mathrm{mm}\times 53\leavevmode\nobreak\ \mathrm{mm}$
photocathode and generate detected photoelectrons that are collected and
amplified by the MCP-PMT photodetector.
Many processes in the signal formation involves random features. For example,
the photon direction and time of production, the photoelectron production
probability, the gain of the Micro Channel Plate used for electron
multiplication, the time transit of the electron propagation through the MCP-
PMT, the noise of readout amplifiers are best described using parametrized
random variables. These are necessary to describe the pulses shapes produced
by a single photoelectron (SPE).
Each SPE induces a signal on typically three readout lines. Thus the typical
30 SPE signals pileup at the output of the transmission lines (Figure 2), and
build the event signal registered by the SAMPIC module.

Figure 2: Set of pulses as registered by the SAMPIC waveshape recorder for a
511 keV energy deposit. For clarity purpose, only the pulses registered on one
side of the transmission lines are shown (half of the set).

###  2.2 Simulation

In order to create a sufficiently large and unbiased database to train and
test our reconstruction algorithms, we developed a detailed simulation of the
ClearMind detector [24]. The knowledge of the hidden physical processes
available in Monte Carlo simulations is necessary to provide the ground truth
(target) for the training of the Machine Learning models. It is also very
convenient to assess the intrinsic performances of our algorithms.

This simulation is based on the Gate v9.0 [17, 22] / Geant4 v7.0 [2, 3, 4]
software allowing to simulate in full details the interaction of the particle
with matter and optical photons generation and tracking. Furthermore, we have
developed specialized software to simulate the photodetector, analog, and
digital electronic components. Necessary parameters have been extracted from
dedicated measurements [11, 10]. This simulation includes the following main
parts of the detector response.

  1. 1.

The gamma interaction in the crystal accounts for three processes:
photoelectric conversion, Compton scattering and Rayleigh diffusion. The two
first processes produce relativistic electron that emits visible photons
through two mechanisms: Cherenkov radiation and scintillation ($\sim$20 and
$\sim$165 photons for 511 keV $\gamma$-quanta respectively).

  2. 2.

Each optical photon is propagated individually by the simulation program.
During the propagation all main physical effects are taken into account:
photon absorption inside the PbWO4 crystal, reflection or absorption on the
crystal borders for the different types of the crystal surface (polished,
ground, absorbing), escape of photons from the crystal into the air.

  3. 3.

Photocathode simulation includes the Fresnel reflection of visible photons at
photocathode boundaries, absorption of photons by the photocathode and
extraction of generated photoelectrons as a function of the photon wavelength.
As a result we compute, assuming a photocathode of nominal efficiency, that we
produce in average 30 photoelectrons for a 511 keV $\gamma$-ray photoelectric
conversion in the crystal and 75% of events contain at least one Cherenkov
photon converted into a photoelectron.

  4. 4.

We then simulate the propagation and the multiplication of individual
photoelectrons generated by the photocathode in the MCP-PMT and parametrize
the main PMT response features: time response, PMT gain and gain fluctuation,
signal sharing between different output anodes.

  5. 5.

Finally, we simulate the signal readout through the transmission lines with
realistic signal shapes, taking into account the possible overlay of several
photoelectrons, electronics noise and digitization sampling.

Most of the simulation parameters are adjusted to the results obtained by the
characterization of the first prototype using pulsed laser in the single-
photon regime. More details about the simulation could be found elsewhere [25,
24].

###  2.3 Waveform preprocessing and input data shaping

Depending on the energy deposited in the crystal, the SAMPIC module records
from 2 to 64 pulse shapes, typically about 30 for a deposited energy of 511
keV. This corresponds to about 10 kBytes of raw data per event. This volume of
data is to be compared with the volume of data to be reconstructed, the
properties of the gamma interaction in the crystal: 3D position, time,
deposited energy, interaction multiplicity, typically 25 bytes per event. The
acquired data are therefore redundant, but also intricate. The parameters to
be reconstructed are encoded in a complex way and are mixed in the pulse
shapes.

Thus reconstructing the properties of the gamma interaction in the crystal is
a complex task.

We first use our knowledge of the physics of the detector to calculate on the
raw data a set of statistical variables, so called "parameter observables"
highly correlated to the parameters of the gamma-ray interaction to be
reconstructed, as well as a second set of "bias observables" which monitor the
known undesirable instrumental effects (saturation of the acquisition
electronics, edge effects, etc…). The data volume of the observables is 100
bytes per event. These observables are then used as inputs of shallow, fully
connected neural networks, whose training takes only a few minutes. The
development cycle is fast, at the cost of a loss of information that is
difficult to anticipate.

In the following paragraph, we explain the parameter observables developed to
correlate to the features targeted by the neural networks presented in this
paper: the gamma interaction position $x_{\mathrm{pos}}$ (_position along the
transmission lines_) and $y_{\mathrm{pos}}$ (_perpendicular to the
transmission lines_).

First, for each transmission line $l$, that triggered the SAMPIC acquisition,
the digitized pulse shapes are acquired: $F_{l,\mathrm{Left}}(t_{j})$ and
$F_{l,\mathrm{Right}}(t_{j})$ at times $t_{j}$ for the Left and Right
transmission line, where the index $j$ corresponds to the sampling time.
Examples of acquired pulse shapes are shown in Figure 3.

We first calculate the electric charges collected (integral of the pulse shape
current over time) at both the end of the 32 transmission lines $C_{l}$.

We compute two interaction observables expected to be correlated to
$y_{\mathrm{pos}}$ :

  * •

We select the line with the largest collected charge and its two neighbors. We
fit a parabola on these three values. The position of the maximum of this
parabola is the first observable.

  * •

The second observable is the median of the distribution of line numbers
weighted by the charge collected on these lines :
$\mathrm{Med}_{L}=\mathrm{median}(l,C_{l})$. It is common that one line
carries a large fraction of the total charge collected over the 32 readout
lines registered in an event. In order to extract as much information as
possible from the surrounding line charge values, we developed an "upgraded"
median algorithm documented in section Appendix: Median algorithm on weighted
distributions. This is the "median" algorithm we will use for all the
observables calculations.

Figure 3: Waveforms registered on one trigered line $l$. Red and green lines
are the left $F_{l,\mathrm{Left}}(t_{j})$ and right
$F_{l,\mathrm{Right}}(t_{j})$, time shifted, registered pulses shapes. Black
line shows the time difference waveform
$F_{l,\mathrm{Left}}(t_{j})-F_{l,\mathrm{Right}}(t_{j})$. We identify on this
line three pulses clusters at 4.6 ns, 6 ns and 8.5 ns. For each of them, the
time difference curve shows a bipolar shape, correlated to the position of
each photoelectron charge induction along the readout line.

Quantifying observables expected to correlate to $x_{\mathrm{pos}}$ is more
complex. It turns out that on some lines, we can identify few pulses separated
in time in the recorded pulse shapes. We denote each identified pulse by the
index $p$. When this happens, we quantify the pulse detection times
($T_{l,p,\mathrm{Left}}$ and $T_{l,p,\mathrm{Right}}$) and the total charge
$C_{l,p}$ of each of these pulses. The time difference $\Delta
t_{l,p}=T_{l,p,\mathrm{Left}}-T_{l,p,\mathrm{Right}}$ is correlated to the
position of the pulse induced by this photoelectron along this line. We then
calculate the three following observables:

  * •

the _median_ of the $\mathrm{Med}_{\Delta t}=\mathrm{median}(\Delta
t_{l,p},C_{l,p})$ distribution (all lines and pulses), weighted by the
integrated charge of the each pulses. The details of the algorithm for the
weighted median are given in Section Appendix: Median algorithm on weighted
distributions;

  * •

the _mean_ of the same distribution $\mathrm{Mean}_{\Delta
t}=\mathrm{mean}(\Delta t_{l,p},C_{l,p})$;

  * •

$\Delta t_{\mathrm{lmax}}$, the time difference $\Delta t$ of the largest
pulse of the line (greatest charge collected).

We also used the shape of the digitized pulses $F_{l,p,\mathrm{Left}}(t_{j})$
and $F_{l,p,\mathrm{Right}}(t_{j})$. The
$\mathrm{Diff}_{l,p}(t_{j})=F_{l,p,\mathrm{Left}}(t_{j})-F_{l,p,\mathrm{Right}}(t_{j})$
show bipolar shapes depending on the position of the pulse injection along the
line. We integrated the $\mathbf{F}$irst and $\mathbf{S}$econd components of
the bipolar shape $\mathrm{Int}_{l,p,F}$ and $\mathrm{Int}_{l,p,S}$. We then
compute the
$\mathrm{Bipol}_{l,p}=(\mathrm{Int}_{l,p,F}-\mathrm{Int}_{l,p,S})/C_{l,p}$,
for each pulse. The saved observables are then:

  * •

the _weighted median_ of the distribution
$\mathrm{Med}_{\mathrm{Bipol}}=\mathrm{median}({\mathrm{Bipol}}_{l,p})$ (see
section Appendix: Median algorithm on weighted distributions for the details)

  * •

the _mean_ of the distribution
$\mathrm{Mean}_{\mathrm{Bipol}}=\mathrm{mean}({\mathrm{Bipol}}_{l,p})$.

Along with these observables, we also compute parameters and bias observables
relevant for reconstruction the other properties gamma interactions in the
crystal (see section 7). Some of them are also expected to be relevant for the
uncertainty estimation. These 23 observables are the inputs of the following
neural networks.

An alternative approach would be to use the full waveforms’ data as inputs for
convolutional deep neural networks. However, preliminary tests showed that the
complexity of the data required to implement deeper network structures,
involving training times of days on our GPU card, thus with slower development
cycles. We are considering this option for future works.

##  3 Methodology

###  3.1 Baseline method

The simplest approach to gamma-interaction reconstruction is using the fact
that the highest density of the detected photons corresponds to the
coordinates of the interaction point. This suggests to use the transmission
line with the maximum detected charge as the reference line. To reconstruct
the coordinate across lines, $y_{R}$, we calculate the weighted average of
center-of-line for this and two neighboring lines:

$y_{R}=\frac{\sum_{l=i-1}^{i+1}y_{l}C_{l}}{\sum_{l=i-1}^{i+1}C_{l}}\;,$ (1)

where $y_{l}$ is the y-coordinate of the line center, $C_{l}$ is the charge of
line $l$ (only the negative signal part is used for the charge calculation),
$i$ is the reference line number.

The coordinate along the lines, $x_{R}$ is reconstructed as

$x_{R}=\frac{(t_{\mathrm{Right}}-t_{\mathrm{Left}})}{2}\times
v_{\mathrm{signal}}\;,$ (2)

where $t_{\mathrm{Right}}$ and $t_{\mathrm{Left}}$ are the time measured at
the right and left ends of line $i$ respectively and $v_{\mathrm{signal}}$ is
the signal propagation speed, assumed to be 35% of speed of light in the
simulation.

This approach, based only on expert knwoledge of the physical processes is
used as a reference for comparison with the performance of our approach based
on neural networks.

###  3.2 ML approach

This section describes our approach based on a supervised Machine Learning
algorithm, that takes as input the 23 preprocessed variables described in
sections 2.3 and 7 in order to predict the interaction position. We use a
Neural Network model that is able to provide both the prediction and an
estimation of an uncertainty associated to this prediction thanks to the
paradigm of Density Neural Networks.On the contrary to the baseline method,
the use of Neural Network is mainly based on the exploitation of simulated
data to build a model that link the observable quantities to the
reconstruction of the gamma interaction position. Moreover, the baseline
method is not able to provide uncertainties associated to the reconstruction,
which is a key difference with our approach.

####  3.2.1 Density Neural Network

The conventional approach to perform a regression through a Neural Network is
to define a loss function that measures the discrepancy between the
predictions of the neural network and the expected values. With classical
notations, let $\mathcal{D}=\\{(x_{i},y_{i}),i=1,2,\ldots,N\\}$ be the
database and $\mu_{\theta}(x)$ the prediction of the expected value of the
neural network parameterized by $\theta$. By noting $||\bullet||$ the
Euclidean norm, a usual loss function is the Mean Square Error :

$\textrm{MSE}(\theta)=\frac{1}{N}\sum_{i=1}^{N}||y_{i}-\mu_{\theta}(x_{i})||^{2}.$
(3)

The learning procedure consists in finding parameters $\theta^{*}$, among the
space of possible parameters $\Theta$, that minimizes the Mean Squared Error.
It corresponds to solving the following problem:

$\theta^{*}=\underset{\theta\in\Theta}{\textrm{argmin}}\;\textrm{MSE}(\theta).$
(4)

The Density Neural Network approach proposes to make the assumption that the
expected outputs $y_{i}$ are a realization of a probability distribution such
as the normal distribution
$y_{i}\sim\mathcal{N}(\mu_{\theta}(x_{i}),\sigma_{\theta}^{2}(x_{i}))$ or a
mixture of normal distributions [6]. However, this method can be adapted to
other distributions and we propose to use a distribution that considers some
specific physical constraints of our application.

The main physical constraint on the interaction position is that this
interaction must be located in the limits of the detector. Therefore, we
decide to use a truncated normal distribution in order to insure this
constraint.

Beforehand in order to avoid confusion with $(x,y)$ coordinates of the
position of the interaction in the detector, we will note $s$ the inputs data
from the preprocessing describes in Section 2.3. So formally, our neural
network is now expected to output four values for each inputs $s_{i}$, the two
means $\mu_{\theta,x}(\mathrm{s_{i}}),\mu_{\theta,y}(\mathrm{s_{i}})$ and the
two scale parameters
$\sigma_{\theta,x}^{2}(\mathrm{s_{i}}),\sigma_{\theta,y}^{2}(\mathrm{s_{i}})$
of the assumed truncated Gaussian law. We have assumed independent normal laws
on $x_{i}$ and $y_{i}$. In this case, the hypothesis that each output
$x_{i},y_{i}$ follows a truncated normal distribution in $[a,b]$ leads to the
two following distributions relative to the two coordinates with $z=x$ or
$z=y$:

$\displaystyle
p(z_{i}|\mu_{\theta,z}(\boldsymbol{\mathrm{s_{i}}}),\sigma_{\theta,z}^{2}(\boldsymbol{\mathrm{s_{i}}}))$
$\displaystyle=$
$\displaystyle\underbrace{\frac{1}{\sqrt{2\pi\sigma_{\theta,z}^{2}(\boldsymbol{\mathrm{s_{i}}})}}\exp\left(-\frac{(z_{i}-\mu_{\theta,z}(\boldsymbol{\mathrm{s_{i}}}))^{2}}{2\sigma_{\theta,z}^{2}(\boldsymbol{\mathrm{s_{i}}})}\right)}_{\textrm{Classical
Gaussian Likelihood }}$ (5) $\displaystyle\times$
$\displaystyle\underbrace{\frac{1}{\Phi\left(\frac{b-\mu_{\theta,z}(\boldsymbol{\mathrm{s_{i}}})}{\sigma_{\theta,z}(\boldsymbol{\mathrm{s_{i}}})}\right)-\Phi\left(\frac{a-\mu_{\theta,z}(\boldsymbol{\mathrm{s_{i}}})}{\sigma_{\theta,z}(\boldsymbol{\mathrm{s_{i}}})}\right)}}_{\textrm{With
truncation in }[a,b]}$

where $\Phi$ is the Cumulative Distribution Function of the standard normal
distribution, $a$ and $b$ are the truncation boundary of the truncated normal
distribution. In our application, $a$ and $b$ represent the limits of the
detector: $a=-30\;\textrm{mm}$ and $b=30\;\textrm{mm}$. In our case, the
condition $a\leq(x_{i},y_{i})\leq b$ is assumed to be always verified : an
interaction occurs inside the detector; otherwise it is a simulation error and
we should not consider it. So, by using equation 5, we derive the likelihood
of the parameters on the whole dataset:

$L(\theta)=\prod_{i}\prod_{z\in\\{x,y\\}}p(z_{i}|\mu_{\theta,z}(\boldsymbol{\mathrm{s_{i}}}),\sigma_{\theta,z}^{2}(\boldsymbol{\mathrm{s_{i}}}))$
(6)

Maximizing this likelihood is equivalent to minimizing the negative log-
likelihood, that is easier to handle, which leads to the following loss
function:

$\displaystyle l(\theta)$ $\displaystyle=$ $\displaystyle-\log(L(\theta))$ (7)
$\displaystyle=$
$\displaystyle\sum_{i}\sum_{z\in\\{x,y\\}}\log\left(\Phi\left(\frac{b-\mu_{\theta,z}(\boldsymbol{\mathrm{s_{i}}})}{\sigma_{\theta,z}(\boldsymbol{\mathrm{s_{i}}})}\right)-\Phi\left(\frac{a-\mu_{\theta,z}(\boldsymbol{\mathrm{s_{i}}})}{\sigma_{\theta,z}(\boldsymbol{\mathrm{s_{i}}})}\right)\right)$
$\displaystyle+$
$\displaystyle\frac{1}{2}\log(2\pi\sigma_{\theta}^{2}(\boldsymbol{\mathrm{s_{i}}}))+\frac{(z_{i}-\mu_{\theta,z}(\boldsymbol{\mathrm{s_{i}}}))^{2}}{\sigma_{\theta,z}^{2}(\boldsymbol{\mathrm{s_{i}}})}$

We want to highlight that the conventional problem with MSE defined in
equation 3 is a special case of equation 7 and $a=-\infty$, $b=+\infty$ (no
truncature) and where the scale parameters
$\sigma_{\theta,x}^{2}(\boldsymbol{\mathrm{s_{i}}})$ and
$\sigma_{\theta,y}^{2}(\boldsymbol{\mathrm{s_{i}}})$ are assumed to be
constant (homoscedasticity hypothesis). Finally, the optimization problem
consists now in finding the parameters $\theta^{*}$ that minimize this new
loss function:

$\theta^{*}=\underset{\theta\in\Theta}{\textrm{argmin}}\;l(\theta)$ (8)

We can solve this problem by using a conventional gradient descent
optimization in order to train the neural network and find optimal parameters.

This approach brings two main components:

  * •

the use of a Density Neural Network allows the estimation of an uncertainty
associated to the network output through the prediction of the terms
$\sigma^{2}_{\theta,x}(\boldsymbol{\mathrm{s_{i}}})$ and
$\sigma^{2}_{\theta,y}(\boldsymbol{\mathrm{s_{i}}})$. These variances are
specific to each example and represent a higher or lower uncertainty on the
predicted positions $\mu_{\theta,x}(\boldsymbol{\mathrm{s_{i}}})$ or
$\mu_{\theta,y}(\boldsymbol{\mathrm{s_{i}}})$ according to the inputs
$\boldsymbol{\mathrm{s_{i}}}$ of the neural network. For instance, in our
application, we can expect an interaction closer to the edge of the detector
is less precisely located than an interaction at the center of the detector.
We highlight that this uncertainty is self-estimated by the neural network and
consequently must be empirically validated, we bring elements for this
validation in section 4.3. This uncertainty can be seen as aleatoric
uncertainty [19] due to the randomness of the phenomenons involved during the
detection and measurements (Cherenkov and scintillation production, loss of
optical photons, randomness in the pulse shape creation…).

  * •

the use of a truncated normal hypothesis allows us to consider the physical
constraints. In section 4, we show that this property is essential to obtain
better performances for interaction close to the edges of the detector,
compared to the use of the classical MSE approach.

####  3.2.2 Architecture and training

We develop two neural networks: one neural network associated to the
conventional Mean Squared Error loss function and one neural network
associated to our custom loss function with the truncated normal hypothesis.
For both architectures, all of the hidden layers are similar: we use six
hidden fully-connected layers with 256 neurons each and a tanh activation
function. The difference is the ouptut layer:

  * •

the output layer for the neural network associated to the MSE loss has two
neurons (one for each coordinate), with a scaled tanh function:
$\frac{b-a}{2}\tanh$ in order to insure that the prediction is located inside
the detector.

  * •

the output layer for the neural network associated to the truncated normal
loss function has four output neurons. The first ones output
$\mu_{\theta,x}(\boldsymbol{\mathrm{s_{i}}})$ and
$\mu_{\theta,y}(\boldsymbol{\mathrm{s_{i}}})$ and the activation function is
the previous scaled tanh function. The second ones output
$\sigma_{\theta,x}(\boldsymbol{\mathrm{s_{i}}})$ and
$\sigma_{\theta,x}(\boldsymbol{\mathrm{s_{i}}})$. The activation function is
$\textrm{softplus}(\alpha)+\epsilon$ where $\alpha$ is the output before
applying the activation function. $\epsilon$ is a small constant that ensures
that $\sigma_{\theta,x}(\boldsymbol{\mathrm{s_{i}}})$ and
$\sigma_{\theta,y}(\boldsymbol{\mathrm{s_{i}}})$ are not too close to 0 and
avoids divergence issues. In our case, we choose
$\epsilon=10^{-6}\;\textrm{mm}$.

The neural networks and their trainings are computed by using the Tensorflow
[1] library. The optimizer is the Adam algorithm. The training is stopped by
an early stopping procedure that monitors the validation loss; we use randomly
20 % of our data as validation data. The code has been executed on a NVIDIA
RTX Quadro 5000. The training time is a dozen of minutes.

####  3.2.3 Generated datasets

We used three datasets to train and test our supervised Machine learning
algorithms. All of them use the detector modeling explained in section 2,
corresponding to a PbWO4 scintallation crystal of 59 mm $\times$ 59 mm. The
inputs are the 23 observables described in sections 2.3 and 7 from the 32
transmission lines and the outputs are both coordinate X and Y of the gamma
interaction, which are saved thanks to the Geant4 simulation.

  1. 1.

For the training set, we simulate a Gamma photon source shaped in a 6 cm cube.
The energy spectrum has been adjusted to increase the probability of high
energy deposits in the PbWO4 and thus generate an approximately flat deposited
energy spectrum. Thus we generate 7 times more 1.2 MeV than 300 keV photons.
The gamma rays impinge the PbWO4 crystal perpendicularly and uniformly over
the entire surface of the optical window. This training batch contains 450 000
events.

  2. 2.

The first test dataset simulates an grid of $9\times 9$ 511-keV gamma ray
point sources, regularly spaced by 7 mm. Again the gamma rays impinge the
PbWO4 crystal perpendicularly to the surface of the optical window. This test
dataset contains 300 000 events.

  3. 3.

The second test dataset simulates a 6 cm cube-shaped Gamma photon source,
mono-energetic at 511 keV. The gamma rays impinge the PbWO4 crystal
perpendicularly and uniformly on the whole surface of the optical window. This
test dataset contains 600 000 events. We design this test dataset with a high
number of examples in order to conduct an accurate analysis of the
performances of the neural network considering the uncertainty estimation.

Figure 4: Simulated sources - Expected positions

##  4 Results

###  4.1 Evaluation on the grid of sources - First test dataset

The positions of the grid of sources are represented on figure 4. This
configuration helps visualising directly some properties of the different
reconstruction algorithms. The results are presented on figure 5.

Figure 5: Grid reconstruction by the different methods

They show that the algorithms are all able to reconstruct the sources in the
center area of the detector, between -20 mm and +20 mm in both directions.
However, the baseline method and the conventional Neural Network associated to
a MSE loss function are not able to give a prediction on the edge of the
detector. These predictions are even incorrectly attributed to other position
at -25 mm and +25 mm for the conventional MSE approach. This "folding" of the
predictions will imply a lower confidence of the reconstruction in this area,
as we show in section 4.2. On the contrary, the use of the truncated Gaussian
likelihood is able to provide positions close to the edge and avoid the
"folding" effect. This property shows a first advantage of the use of this
custom loss function.

The bottom right figure introduces the benefits of the uncertainty evaluation.
To produce this figure, we use a weighting on each event $i$ according to the
predicted uncertainties terms $\sigma_{\theta,x}(\boldsymbol{\mathrm{s_{i}}})$
for the X position, and $\sigma_{\theta,y}(\boldsymbol{\mathrm{s_{i}}})$ for
the $Y$ position. The weights are computed using the following equation:

$w_{i}=\frac{(\sigma_{\theta,x}(\boldsymbol{\mathrm{s_{i}}})^{2}+\sigma_{\theta,y}(\boldsymbol{\mathrm{s_{i}}})^{2})^{-1}}{F},$
(9)

where $F$ is a normalization factor such that $\sum_{i}(w_{i})$ is equal to
the number of detected events on the considered bin of the histogram.

$F=\sum_{i\in\mathrm{bin}}{(\sigma_{\theta,x}(\boldsymbol{\mathrm{s_{i}}})^{2}+\sigma_{\theta,y}(\boldsymbol{\mathrm{s_{i}}})^{2})^{-1}}$
(10)

This weighting penalizes the events with a high predicted uncertainty. The
result shows a reduction of the spreading effect, making the reconstruction
more accurate on the detector. We quantify this gain in accuracy in section
4.2.

###  4.2 Evaluation on the uniform simulation - Second test dataset

Figure 6: Reconstruction of the uniform simulation

For the global evaluation of the performances, we use the simulation of a
uniform distribution of sources (cube source) in front of the detector. The
direct reconstruction provided by the algorithms is shown on figure 6. The
results show the same properties as the grid simulation, especially the
inability of the baseline method and the conventional MSE to reconstruct the
edges of the detector. We can see a partial "folding" effect on the
reconstruction with the truncated Gaussian likelihood. This effect is
attenuated by applying the weighting as defined in equation 9.

Figure 7 shows for each true position the average 2D distance to the
reconstructed position. A high average 2D distance in this case corresponds to
a high spread of the reconstructions. At the center of the detector, the
baseline method provides less spread reconstructions, between 2 and 3 mm, than
the conventional MSE and our approach with the truncated Gaussian likelihood,
between 3 and 4 mm. However, the area of high spread, over 6 mm, is larger for
the baseline method than the machine learning method. We show later that this
default highly affects the precision of the reconstruction, even at the center
of the detector.

The use of the weighting helps the truncated Gaussian likelihood to give more
importance to the most certain events, leading to a spread between 1 and 2 mm
at the center of the detector, between 2 and 4 mm in the intermediate area.
The area with the higher spread, over 6 mm, is thiner than the other
algorithms.

Figure 7: 2D error reconstruction according to the true position,
corresponding to the spread of the reconstruction - The color scale has been
threshold at 6 mm for visualisation purpose

Figure 8 shows for each predicted position the average 2D distance to the true
position. A high average 2D distance in this case corresponds to a low
precision of the reconstructions. Black areas correspond to area with no
predicted position due to the inability of the algorithms to reconstruct some
positions, as shown on figure 6. The baseline method gives very low accurate
results, mostly over 6 mm, even at the center of the detector. By combining
with the previous information, we can give the following interpretation:

  * •

if the algorithm provides a localisation at the center of the detector, we
cannot trust the algorithm since the reconstruction error is over 5 mm, as
shown on figure 8 ;

  * •

if the true position is the center of the detector, we can trust the algorithm
because the reconstruction error is between 2 and 3 mm, as shown on figure 7.
However, for a real application, we cannot access to the true position, then
we cannot exploit this information.

Figure 8: 2D error reconstruction according to the predicted position,
corresponding to the precision of the reconstruction - The color scale has
been threshold at 6 mm for visualisation purpose

The conventional MSE shows a better precision than the baseline method, with a
reconstruction error between 3 mm at the center of the detector to 5 mm close
to the edges. The areas with very high reconstruction error, over 6 mm, are
limited. The truncated Gaussian likelihood shows similar performances as the
conventional MSE for the precision, with the ability to reconstruct the edges,
with less precision. The application of the weighting improves the results,
providing a precision between 1 and 2 mm at the center of the detector,
between 3 and 5 mm in the intermediate area and over 6 mm in a thin band close
to the edge. This band corresponds to the part of the detector without optical
layer and it is expected to have a degradation of the performances in this
area.

Finally, figures 9 and 10 show the overall results, considering the whole test
dataset. The histograms represent the global reconstruction errors, without
selecting any position. We also compute the Root Mean Squared Error (RMSE) of
the predictions.

Figure 9: Global histogram of reconstruction errors on the X position Figure
10: Global histogram of reconstruction errors on the Y position

For both coordinates, the baseline method shows the worse performances. This
result shows the advantage to use a Machine Learning approach to address this
reconstruction problem, due to the complex link between the observable
variables and the positions.

We apply a new 1D-weighting dedicated to each coordinate:

$w_{X,i}=\frac{(\sigma_{\theta,x}(\boldsymbol{\mathrm{s_{i}}}))^{-2}}{\sum_{j}{(\sigma_{\theta,x}(\boldsymbol{\mathrm{s_{j}}}))^{-2}}}$
(11)
$w_{Y,i}=\frac{(\sigma_{\theta,y}(\boldsymbol{\mathrm{s_{i}}}))^{-2}}{\sum_{j}{(\sigma_{\theta,y}(\boldsymbol{\mathrm{s_{j}}}))^{-2}}}$
(12)

On average, the conventional MSE gives slightly better results than the
truncated Gaussian likelihood in terms of RMSE. We can expect this behaviour
because the MSE approach is specifically designed to minimize the MSE, thus
its squared root the RMSE.

Table 1 summarizes several performances such as the mean error (the bias) on
the histogram, RMSE and the Standard Deviation (spread of the distribution
without bias). A small bias is observed for the reconstructions with the
truncated Gaussian likelihood methods for the X coordinate, under 1 mm. On the
contrary, the conventional MSE presents a bias for the Y coordinate
reconstruction. On the other metrics, the weighting applied to the truncated
Gaussian likelihood outperforms the other methods: the application of this
weighting helps improving the performances by giving less importance to
uncertain reconstructions. Especially, the RMSE reaches 2.61 mm for the X
coordinate and 1.91 mm for the Y coordinate.

Quantity |  Baseline method |  Conventional MSE |  Truncated Gaussian lik. |  Truncated Gaussian lik. weighting
---|---|---|---|---
Mean error X (mm) |  -0.01 |  0.01 |  -0.32 |  -0.39
RMSE X (mm) |  8.05 |  4.8 |  5.11 |  2.61
Std dev X (mm) |  8.05 |  4.8 |  5.1 |  2.58
Mean error Y (mm) |  0.00 |  0.19 |  0.04 |  0.02
RMSE Y (mm) |  5.5 |  4.72 |  5.08 |  1.91
Std dev Y (mm) |  5.5 |  4.72 |  5.08 |  1.91
Table 1: Table of performances. Best values are in bold.

###  4.3 Calibration of the uncertainties

The evaluation of the performances shows the advantage of using the
uncertainty prediction as a weighting of the uncertain reconstruction.
However, these results do not provide accurate elements to assess the
calibration of the uncertainties. In this section, we provide an evaluation of
the quality of the uncertainty estimation by using a coverage plot.

To produce this plot, we consider a probability level $\alpha$. For each event
$i$, our neural network provides a prediction of the parameters
$\mu^{(i)}_{\theta}$ and $\sigma^{(i)}_{\theta}$ for a truncated Gaussian
distribution, with bounds $a$ and $b$. Under the hypothesis of this
distribution, we compute the Prediction Intervals (PI) at level $\alpha$ for
each event $I^{(i)}(\alpha)$ such that:

$p(y_{\mathrm{true}}\in I^{(i)}(\alpha))=\alpha$ (13)

where $y_{\mathrm{true}}$ is the expected value. There is an infinite number
of possible intervals that correspond to this condition, we choose the
interval that is centered on the expected median. Then, we compute the number
of events whose expected value belongs indeed to the Prediction Interval,
which leads to the Prediction Interval Coverage Probability (PICP) at level
$\alpha$:

$\mathrm{PICP}(\alpha)=\frac{|\\{y_{\mathrm{true}}\in
I^{(i)}(\alpha)\\}_{i}|}{N}$ (14)

where $|E|$ is the cardinal of the ensemble $E$ and $N$ is the total number of
events. The PICP can be seen as the empirical frequency of the true values
that indeed belong to the Prediction Interval. Finally, we compare the
$\mathrm{PICP}(\alpha)$ to the probability level $\alpha$:

  * •

if $\mathrm{PICP}(\alpha)=\alpha$, we have a perfect calibration;

  * •

if $\mathrm{PICP}(\alpha)>\alpha$, the model is under-confident, or in other
word, conservative;

  * •

if $\mathrm{PICP}(\alpha)<\alpha$, the model is over-confident.

Figure 11 shows the results for different probability levels $\alpha$. For
both coordinates, the PICP is close to the value $\alpha$ which means the
Prediction Interval can be trusted in average. The model is slightly over-
confident for probability levels under 75 % and conservative for higher
probability levels. For instance, Prediction Intervals at 95 % are
conservative.

Figure 11: Coverage plots - The models are slightly too confident for
probability levels lower than 75 % and slightly too conservative for
probability levels higher than 75 %.

##  5 Discussion

###  5.1 Benefits of the study for the detector reconstruction and the PET
imaging

The use of Machine Learning approaches to reconstruct gamma photon position of
interaction in monolithic crystals is a hot topic nowadays. For PET
applications, most of the time it is a LYSO scintillator of several tenth of
millimeters size and about 1 cm thickness coupled with pixelated silicon
photomultiplier (SiPM). For example, in [5] the spatial resolution obtained
with a feed-forward neural network is 0.74 $\pm$ 0.01 mm for the XY plane and
an average 1.01 $\pm$ 0.01 mm for the Z coordinate. Similar results were
obtained in [12] with one multilayer perceptron with five hidden layers and
100 nodes. In [8], the sub-millimeter precision in XY plane is obtained with
an algorithm based on neural networks integrated into a second neural network
for simultaneous estimation of the event position and timestamp. Even a
preclinical PET system with an annular monolithic LYSO scintillator was
proposed by [16] with an inner diameter of 6 cm with minimum transaxial
scintillator thickness of about 1 cm covered by SiPMs. Authors of these study
obtained the spatial resolution on reconstructed interaction position of about
half of millimeter with a ten-layer deep residual-convolutional neural
networks. Sub-millimeter precise photon interaction position determination in
CeBr3 crystal and LaBr33:Ce [18] for Compton Camera imaging system was
obtained with a newly designed convolutional neural network of five layers.

As shown in these recent studies and in section 4, the Machine Learning
approaches provide promising and more accurate reconstructions than the
baseline methods, only based on the knowledge of the detector. In contrast to
prior research with LYSO or other scintillation only crystals, this current
study aims a precise reconstruction of the gamma photon interaction in a PbWO4
crystal with faster (Cherenkov effect) but fewer (lower light yield) optical
photons. Moreover, in contrary to pixelated SiPM with one readout channel per
pixel, a MCP-PMT with fewer readout channels ($2\times 32$ instead of
$32\times 32$) is used. These differences introduce additional complexity on
the precise position reconstruction. The coordinate along the transmission
lines (X-coordinate) turned out to be complex to reconstruct, but our machine
learning models were able to learn to address this complexity.

For monolithic crystals, the so-called edge-effect worsens the spatial
localization performance toward the detector borders [12]. In our study, the
use of a dedicated loss function, the truncated Gaussian likelihood, helps
recovering some dead area at the edges that are not reconstructed by the
baseline methods or the conventional MSE. This effect can be explained by the
fact that no interaction can be observed outside of the detector, so the MSE
avoids to provide predictions close to the area outside of the detector and
consequently, close to the edges. Because we include this prior knowledge in
the truncated Gaussian likelihood, the model trained with this loss function
can predict positions close to the edges.

Different from previous research, an uncertainty on the gamma photon
interaction reconstruction is provided in this study by the Density Neural
Network. This additional information is included by a weighting in the
performance evaluation. We have shown that it significantly improves the
performances. We can object that this improvement is "artificial" because we
do not improve the prediction, we give less importance in the metrics to the
uncertain prediction. However, it is very promising for the application in PET
imaging because the tomographic image reconstruction relies on the processing
of a large number of events. If some events carry less information, they can
degrade the SNR in the reconstructed image. As a consequence, accounting for
the uncertainty information in the tomographic reconstruction shall increase
the SNR in the image. The uncertainty can be included in a spatial resolution
model adapted to each individual event, enabling an event-based resolution
modeling image reconstruction. The proposed approach would be innovative for
PET imaging, with the best of our knowledge, and its impact on the
reconstructed image quality must be assessed. The combination of the proposed
gamma photon interaction reconstruction in the detector with the full PET
image reconstruction is planned for future works.

Finally, the calibration of the uncertainties show that their prediction is
satisfying, although not perfect. The disagreement between the probability
level and the empirical frequency can come from two factors:

  * •

the prediction of the neural networks can be improved by optimizing the model
and the training phase;

  * •

the truncated Gaussian hypothesis could not be the best representation of the
uncertainty. Further works are planned to use other probability distribution
functions, under the constraint that these probability distributions must not
diverge and be differentiable so that they can be used as loss functions.

###  5.2 Generalization of the methodology

Through this paper, we want to highlight the use of a specific methodology,
which can be beneficial for many applications of Deep Neural Network in
regression tasks. This approach relies on the paradigm of uncertainty
estimation, which is a growing topic in the scientific literature on
Artificial Intelligence. Several approaches are usually proposed, such as
Bayesian Neural Network or Deep Ensembles [13]. The use of the Density Neural
Network is an interesting approach when the uncertainties are dominated by
fluctuations and random phenomenon in the data, such as the source of
variability during the scintillation process in our detectors. This approach
brings additional information thanks to the uncertainties that can be used to
gain information on the prediction of neural networks and can be exploited in
numerous applications. A dedicated methodology of exploitation and validation
has been presented in this paper.

More generally, we can describe the methodology as follows:

  * •

The design of a dedicated loss function, based on a negative log-likelihood of
an assumed distribution. This distribution is often assumed as a Gaussian
distribution, but it is not restricted to this type of distribution and can be
adapted to the problem to address. In our case, we decided to add physical
constraints through a truncation of a Gaussian distribution according to the
detector edges, which requires a specific implementation of the loss function.
We show that this approach brings improvements for reconstructing specific
events, especially close to the detector.

  * •

The estimated uncertainties can be exploited to complete the prediction of the
neural network. The uncertainty can be evaluated event by event to assess the
confidence on the prediction of the model. In our case, we show that we can
statistically improve the reconstruction resolution and this feature will be
helpful in the exploitation of the detector. Uncertain events are consequently
considered as noisy and can be discriminated through this method. Applications
applying Neural Networks to analyse signals or sensors data would draw
benefits from this approach.

  * •

It must be kept in mind that the uncertainties are self-estimated by the
Neural Network, thus they are prone to the same biases as the classical
prediction. Consequently, the validation of the estimated uncertainty is
required and the use of coverage plot is a possible approach to empirically
verify the relevance and reliability of these uncertainties.

##  6 Conclusion

Our objective to predict the 2D position of gamma interaction in the ClearMind
detector is successfully achieved thanks to Machine Learning methods. Our
neural network outperforms the baseline approaches based on physical knowledge
only and used as inputs of the model. The introduction of Density Neural
Networks provides a reliable estimation of the uncertainty in the prediction
of the neural network that can be used to discriminate reconstructions that
are likely less accurate. We aim to draw benefits from this method in further
works for PET image reconstruction. Moreover, we exploit the flexibility of
the Density Neural Network approach to design a truncated gaussian loss
function based on physical constraints. This property helps for the
reconstruction of events in areas far from the center of the detector,
increasing the exploitable surface of the crystal. We want to highlight that
this methodology is generic and can be adapted to other use cases, as far as
we are able to introduce such constraints in the loss function.

The preprocessing step is also important to reduce the number of input
variables, thus to get a compact neural network that can be embedded as close
as possible to the detector. The computation time on test data on a
conventional CPU is 14 s for the 600 000 events, corresponding to 23 µs per
event.

In outlooks, we want to conduct a sensitivity analysis on the input variables.
The objective is to study the importance of the preprocessed observables for
the prediction of the positions. We will verify if we find the expected
correlation between the prediction of the positions and the identified
observables given by the physical expertise. Moreover, the analysis will be
performed on the uncertainty estimation in order to understand which
quantities have the highest influence on the uncertainty in the
reconstruction. We will compare the performances in reconstruction with the
use of the raw data, in order to assess the possible loss of information due
to this preprocessing. In further works, we will also apply this methodology
for the Depth of Interaction, energy and time reconstructions in order to
obtain complete information on the gamma ray reconstructions. Moreover, these
results will be exploited to complete PET image reconstruction in dedicated
studies, in order to improve the SNR on the image by the weighting or the
discrimination of uncertain reconstruction.

## Acknowledgments

This project has received funding from the European Union’s Horizon 2020
research and innovation program under the Marie Skłodowska-Curie grant
agreement No 800 945. We are grateful for the support and seed funding from
the CEA, Programme Exploratoire Bottom-Up, under grant No. 17P103-CLEAR-MIND,
the Cross-Disciplinary Program on Numerical Simulation of CEA for the
financial support of project AAIMME, and the French National Research Agency
under grant No. ANR-19-CE19-0009-01.

##  7 Appendix: Statistical Processing of raw Waveforms

The following paragraphs describes the preprocessed observables we developed
for the ClearMind detector in addition to the description in section 2.3.
First, for each line that triggered the SAMPIC acquisition, let’s define
$F_{l,\mathrm{Left}}(t_{j})$ and $F_{l,\mathrm{Right}}(t_{j})$, the digitized
pulse shapes (Left and Right of transmission line) at time $t_{j}$.
Then, before computing event observables, we calculate for each trigered line:

  * •

The times of the first samples of the pulse shapes $T_{0,\mathrm{Left}}$,
$T_{0,\mathrm{Right}}$

  * •

The times calculated on the rising edges of the pulses, extrapolated to the
half height between the sampled values, $T_{l,\mathrm{Left}}$ and
$T_{l,\mathrm{Right}}$.

  * •

The time lapse when the pulses have exceeded the WaveCatcher trigger
threshold.

  * •

The time lapse when the pulses have saturated the WaveCatcher.

  * •

The sum of the charges collected at both ends of the transmission lines
$C_{l}$.

The selected observables can be classified into several categories.
First, a general parameter stores the number of transmission lines that have
acquired a signal.

Then we store the observables expected to be correlated to the _interaction
time_ of the gamma photon in the crystal. These are :

  * •

The time of the first SAMPIC sample recorded on the first transmission line
that acquired a signal. This time $T_{\mathrm{First}}$ is then subtracted from
all the time variables associated with the event.

  * •

The time of the first photoelectron detected on the event. We store the lowest
photoelectron time value (defined as $0.5\times(T_{l,\mathrm{Left}}$ \+
$T_{l,\mathrm{Right}}$)) of all the lines that triggered the acquisition.

We add four bias observables, intended to help the neural network to
decorrelate possible stacking effects on the pulse shapes. These are

  * •

the time lapse beyond 50% of the amplitude on the pulse shapes
$F_{l,\mathrm{Left}}(t_{j})$) and $F_{l,\mathrm{Right}}(t_{j})$ (one
observable for Left pulse shape and one for Right)

  * •

the computed pulses’ rise time (one observable per side)

Then we store the observables expected to be correlated to the _energy
deposited_ by the interaction in the detector which is a quantity that is
planned to be reconstructed for future works. These observables are also
expected to bring information to estimate the uncertainties on the interaction
position. We store :

  * •

The sum of the charges collected over all lines $C_{T}$.

  * •

The sum of the time lapse when pulses have exceeded the triggering threshold
of the WaveCatcher $T_{\mathrm{Thres}}$ over all lines.

  * •

The sum over of the time lapse when pulses have saturated the WaveCatcher
(bias observable), $T_{\mathrm{Sat}}$.

Then we store observables correlated to the position on the axis perpendicular
to the transmission lines $x_{\mathrm{pos}}$ and along transmission line
$y_{\mathrm{pos}}$, as explained in section 2.3.

Finally, we compute quantities expected to be correlated to the _depth of
interaction_ (DOI) of the gamma in the crystal, used for future works to
achieve 3D reconstruction. We use the fact that the farther the gamma
interaction occurs from the photoelectric layer, the more photoelectrons are
produced over a large area. The depth of interaction is thus correlated to the
dispersion of the produced photoelectrons. We can also expect that this
dispersion is relevant for our neural network to estimate the uncertainty on
$x_{\mathrm{pos}}$ and $y_{\mathrm{pos}}$. To quantify this dispersion we
first calculate the distribution: $\lvert l-\mathrm{Med}_{L}\rvert$ , weighted
by the charge measured on the line l, $C_{l}$.
$\mathrm{Med}_{l}=\mathrm{median}(\lvert l-Med_{L}\rvert,C_{l})$. The two
first observables are then :

  * •

the median($\lvert l-\mathrm{Med}_{L}\rvert,C_{l}$) and

  * •

the mean($\lvert l-\mathrm{Mean}_{L}\rvert,C_{l})$

To quantify the dispersion along the transmission lines, we use the same
algorithms, but on the previously computed $\Delta t_{l,p}$ and
${\mathrm{Bipol}}_{l,p}$ distributions (section 2.3), weighted by the charge
collected in the pulses. The observables are then :

  * •

the median($\lvert\Delta t_{l,p}-\mathrm{Med}_{\Delta t}\rvert,C_{l,p}$),

  * •

the mean($\lvert\Delta t_{l,p}-\mathrm{Mean}_{\Delta t}\rvert,C_{l,p}$),

  * •

the
median($\lvert\mathrm{Bipol}_{l,p}-\mathrm{Med}_{\mathrm{Bipol}}\rvert,C_{l,p}$)

  * •

the
mean($\lvert\mathrm{Bipol}_{l,p}-\mathrm{Mean}_{\mathrm{Bipol}}\rvert,C_{l,p}$)

We did not devise until now any bias observables on the depth of interaction
parameter.

## Appendix: Median algorithm on weighted distributions

Following is our algorithm for computing the median of a binned distribution,
written in C++. It uses the STL library to sort the input values. Instead of
returning the central value of the bin where the mid-point value of the
distribution happen, we decided to make use of the values of the neighbouring
bins to qualify the result.

    bool CompPairValues(const PairValues& Pair1, const PairValues& Pair2)
    { return (Pair1.Val<Pair2.Val); }

    // Compute the mid point of a set of values VValue of weigths Vweight
    double Calc_Median_Weighted( vector<double> Vvalues, vector<double> Vweigth, bool PrintInside) {
      unsigned long NbPair = Vvalues.size();
      if(NbPair==0) { cerr<<"Error in DY_Calc_Median_Weighted: Vvalues empty" <<endl;
        exit(-1);
    }

    // We store values in pairs
      vector<PairValues> VPair(NbPair);
      double Half_SomWeight = 0.;
      for( unsigned long ind= 0; ind < NbPair; ind++) {
        VPair[ind].Val = Vvalues[ind];
        VPair[ind].weigth = Vweigth[ind];
        Half_SomWeight += VPair[ind].weigth;
      }
      Half_SomWeight *= 0.5;

      // We sort values, stable_sort ou sort
      std::stable_sort( VPair.begin(), VPair.end(), CompPairValues);

     // We find mid-weight point
      double ValPlus = 0.;  double ValMinus = 0.;
      double SomWeigPlus = 0.; double SomWeigMinus = 0.;
      unsigned long IndexPlus =0;
      for( unsigned long Index =0; Index<NbPair; Index++)
      { SomWeigPlus+=  VPair[Index].weigth;
        ValPlus= VPair[Index].Val;
        IndexPlus= Index;
        if(SomWeigPlus>Half_SomWeight) break;
        SomWeigMinus = SomWeigPlus;
        ValMinus = ValPlus;
      }

      // Compute the Median
      double ValueMedian, ValueMedian2, Median;
      double SomWeigMinDes ;
      if( (IndexPlus == 0) || (IndexPlus == (NbPair-1)) )
        Median = ValPlus;
      else {
        double Un_WeightIndPlus= 1./VPair[IndexPlus].weigth;
        // Rising Option
        double DiffSom = Half_SomWeight-SomWeigMinus;
        double DeltaVal = ValPlus-ValMinus;
        ValueMedian = ValMinus + DeltaVal*DiffSom*Un_WeightIndPlus;
        // Descending option
        double ValMinDes = VPair[IndexPlus+1].Val;
        SomWeigMinDes = Half_SomWeight*2. - SomWeigPlus;
        double DiffSom2 = Half_SomWeight - SomWeigMinDes;
        double DeltaVal2 = ValPlus-ValMinDes;
        ValueMedian2 = ValMinDes + DeltaVal2*DiffSom2*Un_WeightIndPlus;
        Median = 0.5*(ValueMedian+ValueMedian2);
      }
      if( PrintInside) // Print everything for debugging { }
      return Median;
    }

## References

  * [1] Martín Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Greg S. Corrado, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Ian Goodfellow, Andrew Harp, Geoffrey Irving, Michael Isard, Yangqing Jia, Rafal Jozefowicz, Lukasz Kaiser, Manjunath Kudlur, Josh Levenberg, Dandelion Mané, Rajat Monga, Sherry Moore, Derek Murray, Chris Olah, Mike Schuster, Jonathon Shlens, Benoit Steiner, Ilya Sutskever, Kunal Talwar, Paul Tucker, Vincent Vanhoucke, Vijay Vasudevan, Fernanda Viégas, Oriol Vinyals, Pete Warden, Martin Wattenberg, Martin Wicke, Yuan Yu, and Xiaoqiang Zheng.  TensorFlow: Large-scale machine learning on heterogeneous systems, 2015\.  Software available from tensorflow.org.
  * [2] S. Agostinelli, J. Allison, K. Amako, J. Apostolakis, H. Araujo, P. Arce, M. Asai, D. Axen, S. Banerjee, G. Barrand, F. Behner, L. Bellagamba, J. Boudreau, L. Broglia, A. Brunengo, H. Burkhardt, S. Chauvie, J. Chuma, R. Chytracek, G. Cooperman, G. Cosmo, P. Degtyarenko, A. Dell’Acqua, G. Depaola, D. Dietrich, R. Enami, A. Feliciello, C. Ferguson, H. Fesefeldt, G. Folger, F. Foppiano, A. Forti, S. Garelli, S. Giani, R. Giannitrapani, D. Gibin, J.J. Gómez Cadenas, I. González, G. Gracia Abril, G. Greeniaus, W. Greiner, V. Grichine, A. Grossheim, S. Guatelli, P. Gumplinger, R. Hamatsu, K. Hashimoto, H. Hasui, A. Heikkinen, A. Howard, V. Ivanchenko, A. Johnson, F.W. Jones, J. Kallenbach, N. Kanaya, M. Kawabata, Y. Kawabata, M. Kawaguti, S. Kelner, P. Kent, A. Kimura, T. Kodama, R. Kokoulin, M. Kossov, H. Kurashige, E. Lamanna, T. Lampén, V. Lara, V. Lefebure, F. Lei, M. Liendl, W. Lockman, F. Longo, S. Magni, M. Maire, E. Medernach, K. Minamimoto, P. Mora de Freitas, Y. Morita, K. Murakami, M. Nagamatu, R. Nartallo, P. Nieminen, T. Nishimura, K. Ohtsubo, M. Okamura, S. O’Neale, Y. Oohata, K. Paech, J. Perl, A. Pfeiffer, M.G. Pia, F. Ranjard, A. Rybin, S. Sadilov, E. Di Salvo, G. Santin, T. Sasaki, N. Savvas, Y. Sawada, S. Scherer, S. Sei, V. Sirotenko, D. Smith, N. Starkov, H. Stoecker, J. Sulkimo, M. Takahata, S. Tanaka, E. Tcherniaev, E. Safai Tehrani, M. Tropeano, P. Truscott, H. Uno, L. Urban, P. Urban, M. Verderi, A. Walkden, W. Wander, H. Weber, J.P. Wellisch, T. Wenaus, D.C. Williams, D. Wright, T. Yamada, H. Yoshida, and D. Zschiesche.  Geant4 - a simulation toolkit.  Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 506(3):250–303, July 2003.
  * [3] J. Allison, K. Amako, J. Apostolakis, H. Araujo, P. Arce Dubois, M. Asai, G. Barrand, R. Capra, S. Chauvie, R. Chytracek, G. A. P. Cirrone, G. Cooperman, G. Cosmo, G. Cuttone, G. G. Daquino, M. Donszelmann, M. Dressel, G. Folger, F. Foppiano, J. Generowicz, V. Grichine, S. Guatelli, P. Gumplinger, A. Heikkinen, I. Hrivnacova, A. Howard, S. Incerti, V. Ivanchenko, T. Johnson, F. Jones, T. Koi, R. Kokoulin, M. Kossov, H. Kurashige, V. Lara, S. Larsson, F. Lei, O. Link, F. Longo, M. Maire, A. Mantero, B. Mascialino, I. McLaren, P. Mendez Lorenzo, K. Minamimoto, K. Murakami, P. Nieminen, L. Pandola, S. Parlati, L. Peralta, J. Perl, A. Pfeiffer, M. G. Pia, A. Ribon, P. Rodrigues, G. Russo, S. Sadilov, G. Santin, T. Sasaki, D. Smith, N. Starkov, S. Tanaka, E. Tcherniaev, B. Tome, A. Trindade, P. Truscott, L. Urban, M. Verderi, A. Walkden, J. P. Wellisch, D. C. Williams, D. Wright, and H. Yoshida.  Geant4 developments and applications.  IEEE Transactions on Nuclear Science, 53(1):270–278, February 2006\.
  * [4] J. Allison, K. Amako, J. Apostolakis, P. Arce, M. Asai, T. Aso, E. Bagli, A. Bagulya, S. Banerjee, G. Barrand, B. R. Beck, A. G. Bogdanov, D. Brandt, J. M. C. Brown, H. Burkhardt, Ph. Canal, D. Cano-Ott, S. Chauvie, K. Cho, G. A. P. Cirrone, G. Cooperman, M. A. Cortés-Giraldo, G. Cosmo, G. Cuttone, G. Depaola, L. Desorgher, X. Dong, A. Dotti, V. D. Elvira, G. Folger, Z. Francis, A. Galoyan, L. Garnier, M. Gayer, K. L. Genser, V. M. Grichine, S. Guatelli, P. Guèye, P. Gumplinger, A. S. Howard, I. Hřivnáčová, S. Hwang, S. Incerti, A. Ivanchenko, V. N. Ivanchenko, F. W. Jones, S. Y. Jun, P. Kaitaniemi, N. Karakatsanis, M. Karamitros, M. Kelsey, A. Kimura, T. Koi, H. Kurashige, A. Lechner, S. B. Lee, F. Longo, M. Maire, D. Mancusi, A. Mantero, E. Mendoza, B. Morgan, K. Murakami, T. Nikitina, L. Pandola, P. Paprocki, J. Perl, I. Petrović, M. G. Pia, W. Pokorski, J. M. Quesada, M. Raine, M. A. Reis, A. Ribon, A. Ristić Fira, F. Romano, G. Russo, G. Santin, T. Sasaki, D. Sawkey, J. I. Shin, I. I. Strakovsky, A. Taborda, S. Tanaka, B. Tomé, T. Toshito, H. N. Tran, P. R. Truscott, L. Urban, V. Uzhinsky, J. M. Verbeke, M. Verderi, B. L. Wendt, H. Wenzel, D. H. Wright, D. M. Wright, T. Yamashita, J. Yarba, and H. Yoshida.  Recent developments in Geant4.  Nucl. Instrum. Meth. A, 835:186–225, 2016.
  * [5] M.V. Belov, V.A. Kozlov, V.S. Tskhay, and M.V. Zavertyaev.  Resolution estimation in different monolithic PET detectors using neural networks.  Physica Medica, 106:102527, February 2023.
  * [6] Christopher M. Bishop.  Mixture density networks.  Neural Computing Research Group Report, NCRG/94/004, February 1994\.
  * [7] D. Breton, C. Cheikali, E. Delagnes, J. Maalmi, P. Rusquart, and P. Vallerand.  Fast electronics for particle time-of-flight measurement, with focus on the sampic asic.  Nuovo Cimento C, 43:7, 2020.
  * [8] Pietro Carra, Maria Giuseppina Bisogni, Esther Ciarrocchi, Matteo Morrocchi, Giancarlo Sportelli, Valeria Rosso, and Nicola Belcari.  A neural network-based algorithm for simultaneous event positioning and timestamping in monolithic scintillators.  Physics in Medicine & Biology, 67(13):135001, July 2022.
  * [9] CRYTUR, spol. s r.o., Na Lukách 2283 511 01 Turnov Czech Republic.  https://www.crytur.com, 2020.
  * [10] M. Follin, R. Chyzh, C.-H. Sung, D. Breton, J. Maalmi, T. Chaminade, E. Delagnes, K. Schäfers, C. Weinheimer, D. Yvon, and V. Sharyy.  High resolution mcp-pmt readout using transmission lines.  Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 1027:166092, 2022.
  * [11] M. Follin, V. Sharyy, J.-P. Bard, M. Korzhik, and D. Yvon.  Scintillating properties of today available lead tungstate crystals.  Journal of Instrumentation, 16(08):P08040, aug 2021.
  * [12] Marta Freire, Sara Echegoyen, Andrea Gonzalez-Montoro, Filomeno Sanchez, and Antonio J. Gonzalez.  Performance evaluation of side-by-side optically coupled monolithic LYSO crystals.  Medical Physics, 49(8):5616–5626, August 2022.
  * [13] Jakob Gawlikowski, Cedrique Rovile Njieutcheu Tassi, Mohsin Ali, Jongseok Lee, Matthias Humt, Jianxiang Feng, Anna M. Kruspe, Rudolph Triebel, Peter Jung, Ribana Roscher, Muhammad Shahzad, Wen Yang, Richard Bamler, and Xiao Xiang Zhu.  A survey of uncertainty in deep neural networks.  CoRR, abs/2107.03342, 2021.
  * [14] Andrea Gonzalez-Montoro, Antonio J. Gonzalez, Shirin Pourashraf, Robert S. Miyaoka, Peter Bruyndonckx, Garry Chinn, Larry A. Pierce, and Craig S. Levin.  Evolution of PET Detectors and Event Positioning Algorithms Using Monolithic Scintillation Crystals.  IEEE Transactions on Radiation and Plasma Medical Sciences, 5(3):282–305, May 2021.
  * [15] Alexander M. Grant, Timothy W. Deller, Mohammad Mehdi Khalighi, Sri Harsha Maramraju, Gaspar Delso, and Craig S. Levin.  Nema nu 2-2012 performance studies for the sipm-based tof-pet component of the ge signa pet/mr system.  Medical Physics, 43(5):2334–2343, 2016.
  * [16] Gangadhar Jaliparthi, Peter F Martone, Alexander V Stolin, and Raymond R Raylman.  Deep residual-convolutional neural networks for event positioning in a monolithic annular PET scanner.  Physics in Medicine & Biology, 66(14):145008, July 2021.
  * [17] S. Jan, G. Santin, D. Strul, S. Staelens, K. Assié, D. Autret, S. Avner, R. Barbier, M. Bardiès, P. M. Bloomfield, D. Brasse, V. Breton, P. Bruyndonckx, I. Buvat, A. F. Chatziioannou, Y. Choi, Y. H. Chung, C. Comtat, D. Donnarieix, L. Ferrer, S. J. Glick, C. J. Groiselle, D. Guez, P.-F. Honore, S. Kerhoas-Cavata, A. S. Kirov, V. Kohli, M. Koole, M. Krieguer, D. J. van der Laan, F. Lamare, G. Largeron, C. Lartizien, D. Lazaro, M. C. Maas, L. Maigne, F. Mayet, F. Melot, C. Merheb, E. Pennacchio, J. Perez, U. Pietrzyk, F. R. Rannou, M. Rey, D. R. Schaart, C. R. Schmidtlein, L. Simon, T. Y. Song, J.-M. Vieira, D. Visvikis, R. Van de Walle, E. Wieërs, and C. Morel.  GATE: a simulation toolkit for PET and SPECT.  Physics in Medicine & Biology, 49(19):4543, September 2004.
  * [18] M Kawula, T M Binder, S Liprandi, R Viegas, K Parodi, and P G Thirolf.  Sub-millimeter precise photon interaction position determination in large monolithic scintillators via convolutional neural network algorithms.  Physics in Medicine & Biology, 66(13):135017, July 2021.
  * [19] Alex Kendall and Yarin Gal.  What uncertainties do we need in bayesian deep learning for computer vision?  CoRR, abs/1703.04977, 2017.
  * [20] Paul Lecoq, Christian Morel, John O Prior, Dimitris Visvikis, Stefan Gundacker, Etiennette Auffray, Peter Križan, Rosana Martinez Turtos, Dominique Thers, Edoardo Charbon, Joao Varela, Christophe de La Taille, Angelo Rivetti, Dominique Breton, Jean-François Pratte, Johan Nuyts, Suleman Surti, Stefaan Vandenberghe, Paul Marsden, Katia Parodi, Jose Maria Benlloch, and Mathieu Benoit.  Roadmap toward the 10 ps time-of-flight PET challenge.  Physics in Medicine & Biology, 65(21):21RM01, October 2020.
  * [21] MAPMT-253 Multi-Anode MCP-PMT Datasheet, from Photek Inc.  https://www.photek.com/pdf/datasheets/detectors/ds034-auratek-mapmt253-detector-datasheet.pdf, 2020\.
  * [22] David Sarrut, Mateusz Bała, Manuel Bardiès, Julien Bert, Maxime Chauvin, Konstantinos Chatzipapas, Mathieu Dupont, Ane Etxebeste, Louise M. Fanchon, Sébastien Jan, Gunjan Kayal, Assen S. Kirov, Paweł Kowalski, Wojciech Krzemien, Joey Labour, Mirjam Lenz, George Loudos, Brahim Mehadji, Laurent Ménard, Christian Morel, Panagiotis Papadimitroulas, Magdalena Rafecas, Julien Salvadori, Daniel Seiter, Mariele Stockhoff, Etienne Testa, Carlotta Trigila, Uwe Pietrzyk, Stefaan Vandenberghe, Marc-Antoine Verdier, Dimitris Visvikis, Karl Ziemons, Milan Zvolský, and Emilie Roncali.  Advanced Monte Carlo simulations of emission tomography imaging systems with GATE.  Physics in Medicine & Biology, 66(10):10TR03, 2021.
  * [23] Dennis R Schaart.  Physics and technology of time-of-flight PET detectors.  Physics in Medicine & Biology, 66(9):09TR01, apr 2021.
  * [24] C. H. Sung et al.  Detailed simulation for the ClearMind prototype detection module and event reconstruction using artificial intelligence.  ArXiv e-prints 2209.11587 [physics.ins-det], Submitted to Nucl. Instrum. Meth. A, 9 2022.
  * [25] Chi-Hsun Sung.  Simulation and artificial intelligence for a gamma detector for high resolution PET imaging.  PhD thesis, Université Paris-Saclay, 2022.
  * [26] Joyce van Sluis, Johan de Jong, Jenny Schaar, Walter Noordzij, Paul van Snick, Rudi Dierckx, Ronald Borra, Antoon Willemsen, and Ronald Boellaard.  Performance Characteristics of the Digital Biograph Vision PET/CT System.  Journal of Nuclear Medicine, 60(7):1031–1036, 2019.
  * [27] D. Yvon, V. Sharyy, M. Follin, J.-P. Bard, D. Breton, J. Maalmi, C. Morel, and E. Delagnes.  Design study of a “scintronic” crystal targeting tens of picoseconds time resolution for gamma ray imaging: the ClearMind detector.  Journal of Instrumentation, 15(07):P07029–P07029, jul 2020.