1. Introduction
The R3 B setup (Reactions with Radioactive Relativistic Beams setup)
is a multi-purpose experimental setup used to study nuclear structure properties of short-lived isotopes [1]. The setup will be located
at the high-energy branch of the Super FRagment Separator (SuperFRS) [2] of the Facility for Antiproton and Ion Research (FAIR) [3]
near Darmstadt in Germany.
The heart of the R3 B setup consists of a fixed target, at which a
secondary beam generated in the Super-FRS is directed. The general
goal of the R3 B experiment is to provide a kinematically complete
reconstruction of all particles participating in the reaction [1], so that
the nuclear structure of the beam isotope (which could be very shortlived) can be studied. In order to accomplish this, the fixed target is
surrounded by many different detector systems (see Fig. 1).
For most R3 B experiments, the target is surrounded by silicon strip
vertex detectors [4] to measure the charged-particle tracks close to the
target. The CALIFA (CALorimeter for In-Flight detection of gamma-rays
and high energy charged pArticles) [5] system encloses the vacuum
chamber housing the silicon detectors. CALIFA measures gamma-rays
from the decay of excited nuclei and light-charged particles that are
produced at large angles. The heavier charged-particles produced at
the reaction are usually (strongly) boosted in the forward direction
and travel through the superconducting dipole magnet GLAD (Gsi
LArge Dipole magnet) [6] downstream of the target. This magnet
allows for spectrometric analysis of these charged particles. A complex
tracking system of many detectors located both upstream and downstream of GLAD [7] is used to measure the tracks of these heavier
charged-particles. Several of these detectors are located inside the
vacuum chamber downstream of GLAD, but some are also located
at the end of the tube downstream of the vacuum chamber. Finally,
neutrons produced at the target, which are generally also very forwardboosted, travel right through GLAD and are detected by NeuLAND (Neu
Large-Area Neutron Detector) [8].
NeuLAND [8] is a Time-of-Flight spectrometer for the detection of
fast neutrons in the range 200 MeV − 1000 MeV [10]. The detector has a
modular design and consists almost entirely of active RP408 scintillator
material [11], which is a type of BC408 [12]. The neutrons that should
be detected can undergo hadronic scattering within the scintillators
of NeuLAND, which can lead to the production of secondary charged
particles. These secondary charged particles are then detected through
their scintillation light [10] and their energy and time of arrival are
registered by the corresponding electronics.
Each NeuLAND module, called a double-plane (dp), consists of 50
horizontally oriented scintillators of 5 cm × 5 cm × 250 cm, followed by
50 such scintillators in vertical orientation. Each scintillator is read out
by two Photo-Multiplier Tubes (PMTs) (one at each endpoint) and is
held in place with two 5 cm wide and 0.5 mm thick aluminum strips
(see Fig. 2). From these specifications, it follows that each module
(dp) has an active area of 2.5 m × 2.5 m and a thickness of 10 cm.
The NeuLAND design goal is to have 30 dp. However, in the present
situation (at the publication date of this work), funding is only secured
for 16 dp and only 8 dp are currently mounted in the R3 B setup (see
Fig. 3).
The problem that we wish to address in this paper, is that a single
high-energy neutron can create a complex shower of secondary particles (both charged and uncharged) in NeuLAND, which can produce
signals in any number of scintillators between one and a hundred. This
raises the issue which of these signals corresponds to the head of the
shower; the first signal generated by the particle shower. This shower
head is the information we need from NeuLAND [13], because together
with the time and position of the reaction at the target, the time and
position of the shower head allows us to reconstruct the 4-momentum
vector of the neutron [8].
The problem of finding the shower head(s) among all the scintillator
signals in NeuLAND is challenging. Especially in the situation where
multiple neutrons have to be detected in coincidence, solutions are
far from trivial because of two reasons: (1) it is not (always) known
a priori how many neutrons have impinged on the detector and (2)
showers from distinct neutrons tend to overlap quite often. These
problems were already addressed in Refs. [8,13]. However, the solutions proposed there have typical shortcomings in neutron detection,
mainly in terms of efficiency (how often one succeeds in correctly
identifying the shower heads). A good efficiency for identifying the
shower heads is crucial, because this determines the required beamtime for an experiment. For this reason, we propose a newly-developed
Machine Learning (ML) algorithm [14] to find the shower heads in
NeuLAND.
The methodology of the new ML algorithm is discussed in Section 2.
The performance of this algorithm (mainly in terms of efficiency)
and its systematic uncertainties are discussed in Sections 3 and 4. In
Section 5, future possibilities for further improvements on top of the
new algorithm are discussed. Finally, conclusions and a summary are
presented in Section 6.
2. Methodology
2.1. Data generation and preparation
In order to explain the implementation of the ML algorithm in the
NeuLAND data-analysis procedure, a brief overview of this procedure
is given. This procedure is explained in Fig. 4. Central in this scheme
are the NeuLAND ‘‘Hits’’. These are tuples of 5 numbers: {𝐸, 𝑡, 𝑥, 𝑦, 𝑧}.
Within one event, each scintillator that produced a signal is assigned
exactly one such a 5-tuple. 𝑡 and 𝐸 stand for the energy deposition and
the time of the scintillator hit, as reconstructed from the PMT pulses.
The position coordinate (either 𝑥 or 𝑦) along the scintillator is also
obtained from these PMT pulses. The other two position coordinates
are obtained from the geometric position of the scintillator. Together,
these 5-tuples (NeuLAND ‘‘Hits’’) contain all relevant information from
NeuLAND.
The ‘‘Hits’’ can either be obtained from the experiment, or from
simulation. In both cases, the information from all PMT pulses generated in NeuLAND is first converted to QDC ((Q)charge-to-Digital
Convert) values and TDC (Time-to-Digital Convert) values. Next, a
‘‘Hit’’ is produced for each scintillator in NeuLAND where a PMT pulse
was generated at both ends. In such a ‘‘Hit’’, 𝐸 is obtained as the
geometrical mean of the two QDC values of those two PMT pulses [8].
𝑡 is obtained as the standard mean of the two corresponding TDC
values. The position coordinate along the scintillator is obtained from
the difference of those same TDC values.
To obtain the ‘‘Hits’’ from simulation, the PMT pulses were phenomenologically obtained from the outcome of the Monte Carlo transport. The corresponding TDC values were calculated with respect to the
beginning of the event, which is exactly known in a simulation. This
process is called digitization and it accounts for several experimental
effects in the PMT pulses, such as time resolution (𝜎 = 150 ps), energy
resolution (𝜎 = 50 keV), detection threshold (1 MeV of deposited energy
for virtually all experiments), saturation effects, Birk’s law, and light
attenuation. For more details on how the PMT pulses were calculated,
the interested reader is referred to Sect. 7.2 in Ref. [15]. The resolution
effects result in uncertainties in the ‘‘Hit’’-values of 𝜎 = 212 ps for the
𝑡-value and 𝜎 = 3 cm for the position coordinate along the scintillator.
In this work, the simulations were carried out by R3BRoot [16,17],
which is an integral framework for simulation and data analysis of
the R3 B experiment. For our simulations, the NeuLAND geometry, as
discussed in Section 1, was used (with a distance of 14 m between
NeuLAND and the target position). Individual scintillators were modeled in agreement with Refs. [10,13]. For the physics list, one of the
following two reference physics lists was used: QGSP_INCLXX_HP or
QGSP_BERT_HP. Both of these physics lists have been benchmarked
against experimental data [13] and reality was found to be about
halfway between them. As event generator, the traditional NeuLAND
simulation files of the 132 Sn breakup reaction with a relative energy of
500 keV were used [13]. The relative energy between the neutrons of a
single event is defined as the invariant mass of those neutrons together
minus the sum of the individual neutron masses:
The invariant mass is obtained by first adding the 4-momenta of all
neutrons within the event, and then taking the length of the resulting
4-momentum vector. As such, the relative energy is a measure for the
excitation energy of the recoil nucleus (although not exactly the same).
When the ‘‘Hits’’ are obtained from the experiment, the TDC values
are measured with respect to a common-stop signal. However, as not all
data channels have exactly the same processing speed and the commonstop signal usually does not correspond to the beginning of the event,
the TDC values have to be synchronized before they can be used to
obtain ‘‘Hits’’ that can be compared to the simulation. Since for all
NeuLAND experiments, rates are expected to be low enough (< 1 MHz)
so that different events do not overlap, this synchronization can be
handled event-by-event. Likewise, QDC values have to be calibrated
(event-by-event) before they can be used to obtain such ‘‘Hits’’. This
calibration and synchronization procedure is also used to correct for
experimental effects that were not considered in the digitization process, such as biases and walk effects. This procedure to first apply such
corrections to the experimental data and then obtain ‘‘Hits’’ for further
data processing, is what is typically used in NeuLAND experiments such
as Ref. [18].
A full discussion of the two considered physics lists
QGSP_INCLXX_HP (INCLXX for short) and QGSP_BERT_HP (Bertini
for short) used is beyond the scope of this work, but we shall briefly
discuss the most relevant differences between them. NeuLAND is a
neutron detector of BC408 [8], a CH2 -organic scintillator. Hence, the
most important parameters in the physics list are the neutron-tohydrogen and neutron-to-carbon cross sections. These cross sections
have been tabulated in Table 1 for various neutron energies. These cross
sections were obtained by directly extracting them from simulations.
They are not dominated by one or two specific reaction types, but are
a sum of many different reactions.
The Bertini physics list (QGSP_BERT_HP) is based on the Bertini
cascade model developed to simulate high-energy particle physics [19].
On the other hand, the INCLXX physics list (QGSP_INCLXX_HP) is an
experimental/phenomenological physics list developed for simulating
nuclear physics experiments in the intermediate energy range [20].
As INCLXX is dedicated to nuclear physics, it explicitly takes the
production of secondary light nuclei into account, such as deuterium,
tritium and helium-3. Bertini does not do this, which is the main reason
why the neutron-to-carbon cross sections are significantly lower for
Bertini. Both Bertini and INCLXX do take the production of secondary
alphas into account. Energy spectra of the produced secondary protons
do not differ very much. However, as our benchmark in Ref. [13]
is about halfway between the two physics lists, we conclude that
INCLXX overestimates the production of light nuclei, while Bertini
underestimates this.
2.2. Multiplicity determination
After the ‘‘Hits’’ are identified, they are first clustered together
according to the procedure outlined in Refs. [13,15]. This procedure
assigns ‘‘Hits’’ that have 𝛥𝑥 ≤ 7.5 cm, 𝛥𝑦 ≤ 7.5 cm, 𝛥𝑧 ≤ 7.5 cm and
𝛥𝑡 ≤ 1.0 ns to the same cluster. All these clusters (even those of a single
‘‘Hit’’) are then considered in the subsequent analysis. The next step
(see Fig. 4) is to determine the multiplicity of the event (the number of
neutrons that came from the target). For this step, we propose the use of
ML. We have chosen to use of a Deep Neural Network for this, because
this seemed to be the fastest option in terms of CPU processing time.
A Convolutional Neural Network (CNN) has also been investigated for
the multiplicity determination [21], but while this method was found to
have about the same performance (in terms of efficiency) as our work
(see Appendix), its CPU processing time was much slower.
The DNN was equipped with three input neurons per scintillator.
These input neurons contain the ‘‘Hit’’ energy deposition 𝐸, the time of
the scintillator ‘‘Hit’’ 𝑡 and the position coordinate of the ‘‘Hit’’ along
the scintillator (either 𝑥 or 𝑦). For scintillators without a hit, three
zeros were given as inputs. This way, all ‘‘Hit’’ information (which
is all relevant information NeuLAND has to offer [8]) is fed into the
DNN. The information of the other two position coordinates from the
‘‘Hits’’ (which are determined from the geometrical position of the
scintillator) is contained in the positions of the input neurons in the
network, because the same input neuron always correspond to the same
scintillator. Since it was known from Ref. [21] that adding the total
energy deposition in the full detector per event and the total number
of clusters per event as extra inputs to the DNN improved efficiency,
these two quantities were also given to the DNN as inputs. This resulted
in 9002 input neurons for the NeuLAND design case (30 dp of 100
scintillators each) [8]. All inputs were linearly scaled to [0, 1] before
they were fed into the DNN. 𝑡-values (‘‘Hit’’-times) above a certain
threshold were artificially put to zero. This was done, because a small
fraction of the ‘‘Hits’’ (< 0.1%) has very large 𝑡-values of the order of
104 ns. Such rare and very large 𝑡-values were found to be extremely
disruptive to the DNN training procedure and, therefore, had to be
removed. The proper time threshold depends on the simulated scenario.
For the design-case of 30 dp, a neutron energy of 600 MeV and a
distance between NeuLAND’s front-end and the target of 14 m, 100 ns
was found to be a good threshold (in this situation the neutrons need
72 ns to travel to the rear end of NeuLAND).
As this DNN is supposed to classify a given event according to its
multiplicity (the number of neutrons that came from the target), a
separate output neuron for each possible multiplicity was implemented.
Typical NeuLAND simulations consider multiplicities up to five [8,13],
which would result in five output neurons. These output neurons were
equipped with the SoftMax activation function [22], as this transforms
the network output into a probability per multiplicity. As the best
method for multiplicity determination is still unknown at the time
of this work, a densely connected network was considered with the
ReLU activation function on each hidden neuron, as this a-priori allows
any combination of the input data to be considered as a possibility
to determine the output. The ReLU activation function is defined as
𝑓 (𝑥) = 𝑥 for 𝑥 ≥ 0 and 𝑓 (𝑥) = 0 for 𝑥 < 0 [23].
The network was implemented with the Keras [23] user-interface
to the TensorFlow [24] framework. Training was done with Supervised
Learning (SL) [25] using the ADGRAD [26] minimization algorithm
(learning rate 0.001, 𝜖 = 0 and zero decay) and the Categorical Cross
Entropy [27] minimization function. The number of hidden layers was
2, with 9000 (first layer) and 1200 (second layer) neurons each. These
parameters (minimization algorithm, minimization function, number
of hidden layers and number of neurons per hidden layer) were determined by optimization. This procedure is discussed in detail in the
Appendix.
One million simulated events were generated for the training
(200.000 per multiplicity). During the training, only events were used
where the number of neutrons detected by NeuLAND was identical to
the number of neutrons that came from the target. This was done,
because it would be incorrect to train the DNN on events where it
will never be able to determine the correct answer, which is the case
when not all information is available: the situation where not all
neutrons are detected by NeuLAND. The remaining events (that passed
this condition) were subdivided into batches of 1000 events and the
network weights were updated after each of these batches. The number
of one million events and the batch size of 1000 were chosen as such,
because it was observed that these numbers were sufficient to reach
a saturated network accuracy under almost all circumstances within
2 epochs. Since the multiplicity condition implies that events with
higher multiplicities are more often discarded, the undiscarded highermultiplicity events were given larger weights in the minimization
function to prevent the DNN from becoming biased. This bias could
result in a loss from roughly 25% (30 dp) to several hundred percent
(few dp) in the proper reconstruction of multiplicity-five events.
2.3. ‘‘Hit’’ selection
In the last step in Fig. 4, the ‘‘Hit’’ selection, the actual shower heads
have to be identified. For this, we also propose the use of ML. The network design for the ‘‘Hit’’ selection is based on the successes obtained in
the multiplicity determination. The network for the ‘‘Hit’’ selection was
chosen as a densely connected DNN with 14 input neurons, two output
neurons with the SoftMax activation function and 12 hidden layers of
200 neurons with the ReLU activation function each. An event is then
sent though this network cluster-by-cluster and a different copy of this
network was trained for each multiplicity under consideration. The 14
input neurons contain different cluster properties like Time-Of-Flight,
relativistic beta, energy deposition, number of ‘‘Hits’’ in the cluster,
etc. (all inspired by Ref. [13]). The two output neurons give the two
probabilities of the cluster containing at least one shower head, and of
not containing any shower heads at all, respectively.
Using the network output, a score is computed for each cluster as the
difference between the two output neurons. Subsequently, all clusters
in the event are sorted according to this score. Then, all clusters are
discarded, except for the ones with the highest score. The number
of clusters kept is decided by the multiplicity DNN. Subsequently,
the ‘‘Hits’’ in the clusters which have the shortest TOF define the
shower heads. The Supervised Learning for the ‘‘Hit’’ selection was done
similarly to the multiplicity determination.
The full simulation and network training as described in all sections
above required a computational time of roughly 6 h when the configuration of 30 dp with a neutron beam energy of 600 MeV was used. The
generation of simulated data was done using an i7-8750H CPU and 16
GB of RAM memory. The network training was done using a 4 GB GPU
(NVIDIA GeForce GTX 1050 Ti 4 GB). The disk-storage requirements
for this computation were about 250 GB. The computer code of the
new ML algorithm (both the multiplicity determination and the ‘‘Hit’’
section) is available as public-domain software [28].
2.4. Reference algorithms
To evaluate the performance of our new ML algorithm (denoted
DNN algorithm in the following), its efficiency is compared to two
‘‘Perfect Tracking’’) [13] rely on the use of simulation data to optimize
certain decision parameters. Since simulations generally do not perfectly agree with experimental conditions, this dependence introduces
systematic uncertainties to the application of all algorithms. These
systematic uncertainties come from three different sources (the three
main components in a GEANT4 simulation): the geometry, the event
generator and the physics list. The systematic uncertainties from the
first two sources can generally be suppressed by closely matching
the geometry and the event generator to the experimental conditions,
although one should keep in mind that this requires a re-optimization
(or re-training) for each new experiment.
However, finding a physics list that closely matches the reality
with the desired precision is still an open issue for NeuLAND experiments [10,13]. The alternative to using simulation data would
be to perform measurement of reactions such as the break-up of the
deuteron which would produce neutrons with known energies and
angles once the proton is tagged with a high precision. However, these
measurements have to be done for a large number of energies to be
useful and as such are not a very viable alternative for the large
range of energies desired. Because of the fact that the geometry and
the event generator contributions to the systematic uncertainties can
generally be made much smaller than the physics list contribution,
we have only concentrated on the quantitative study of the physics
list contribution to the systematic uncertainties. For this reason, we
repeated all simulations for this work four times and changed the
physics list between QGSP_INCLXX_HP and QGSP_BERT_HP at both
the optimization (or training) and the validation level. The mean of the
four outcomes was then considered as the algorithm performance and
the largest difference between the mean and each of the four individual
outcomes was used as our estimation of the physics list uncertainty.
As reality was determined to be halfway between two physics
lists [13], one could, in principle, also compose a realistic training set
by mixing events of the two physics lists. However, as the benchmarking of the two physics lists in Ref. [13] was only done for neutron
energies of 110 MeV and 250 MeV and only for the NeuLAND detection
efficiency, this procedure is unreliable. The appropriate weight factors
of the two physics lists are indeed roughly 50% in that situation, but
there is no experimental evidence that these weight factors remain the
same at higher neutron energies, or for other observables. For this reason, we have decided not to pursue this idea and estimate the physics
list uncertainties as described above. However, once appropriate weight
factors have been measured over the full range of neutron energies,
event mixing may prove to be a powerful tool in the future to reduce
the physics-list uncertainties and/or to improve the physics-list models.
other shower-head identification algorithms in the upcoming sections:
the TDR algorithm and the ‘‘Perfect Tracking’’ algorithm. Note that
from now on, the notion ‘DNN algorithm’ stands for the combination
of multiplicity determination and ‘‘Hit’’ selection, both handled by ML.
The ‘‘Perfect Tracking’’ algorithm relies on information from the Geant4
simulation. It traces the particle showers created by each individual
neutron to find out whether this shower produced any NeuLAND ‘‘Hits’’
(see Fig. 4). If ‘‘Hits’’ are produced, the neutron is assigned a shower
head: the first produced ‘‘Hit’’. Otherwise, the neutron is assigned no
shower head at all.
Since the ‘‘Perfect Tracking’’ algorithm relies on the ‘true’ particle
showers from the Geant4 simulation, it will always provide the best
possible shower-head identification. For this reason, its outcome was
used as the ‘correct’ output for the Supervised Learning of the DNNs.
However, despite its name the ‘‘Perfect Tracking’’ algorithm does not
always result in perfect efficiency. This is because sometimes, a neutron
may just not interact (sufficiently) with NeuLAND to produce any
‘‘Hits’’ [13]. Hence, a comparison between the DNN algorithm and
the ‘‘Perfect Tracking’’ algorithm can reveal which part of the DNN
algorithm efficiency is caused by imperfections in the algorithm and
which part is caused by external factors (the lack of proper neutron
interactions leading to ‘‘Hits’’).
The TDR algorithm refers to the algorithm originally proposed in
the NeuLAND Technical Design Report (TDR) [8]. This procedure is
illustrated in Fig. 5 for a situation of 30 dp and 600 MeV neutrons.
Here, the multiplicity is determined by plotting the number of clusters
per event versus the total energy deposition in that event and imposing
linear ‘‘decision’’ cuts (the envelopes defined by the diagonal black lines
in Fig. 5). The slope is the same for all cuts, but the distance between
the cuts is different for different multiplicities. The slope and distances
are optimized using simulation data. Subsequently, an event is assigned
a multiplicity depending on where it is located between the ‘‘decision’’
cuts.
We would like to note that the situation where all neutrons within
a single event have (almost) the same energy is actually quite common [8]. In most R3 B experiments., the beam consists of a heavier
nucleus of interest with an energy of several hundred MeV per nucleon.
This nucleus is then studied in inverse kinematics by impinging it on a
very light fixed R3 B-target. As such, the knocked-out neutrons typically
have energies of several hundred MeV as well, but with variations of
the order of only a few MeV (determined by the shell from which
they are removed). Hence, for such experiments, our algorithms can
be trained/optimized with mono-energetic neutrons.
The ‘‘Hit’’ selection of the TDR method is exactly as discussed for
the DNN algorithm in Section 2.3, only clusters are not sorted according
to some DNN-computed score, but according to their so-called 𝑅-value.
