Position Reconstruction in the DEAP-3600 Dark Matter Search Experiment

This method is based on the pattern of charges measured by the 255 PMTs around the AV. The position is calculated by maximizing a likelihood function, constructed with a probability distribution determined with a MC simulation of the detector.

For a given simulated event position, the number of photons measured by each of the PMTs is calculated with the full optical model in our MC simulation [18] which includes the refractive index, group velocity, and Rayleigh scattering length for VUV and visible photons as well as their dependence on thermodynamic properties. The simulation also takes into account the solid angle of the PMTs as viewed from the vertex, light absorption, re-emission at the TPB wavelength shifter, diffusion, and reflections at surfaces. The probability of detecting photons by a PMT is a highly non-linear function of the event’s coordinates and energy, which changes rapidly as the interaction point comes closer to the PMT and has to be determined with the simulation. To determine the event’s location $\vec{x}$, the following likelihood is maximized with a simplex algorithm:

where $P(n_{i};\lambda_{i})$ is the Poisson probability of observing $n_{i}$ PE in PMT $i$ over the full 10 $\mu$s time window. The coordinate system has its origin $\vec{x}=(0,0,0)$ at the center of the spherical AV and the $z$ axis pointing in the vertical direction towards the detector’s neck. The parameter

describes the expected number of PE hitting the PMT $i$ located at position $\vec{r}_{i}$, given the event position $\vec{x}$, an angle $\theta_{i}$ between $\vec{x}$ and $\vec{r}_{i}$, and the total number of detected photons $n_{\rm PE}$ (see Figure 1 for a visual description of the variables). The factorized total number of photons $n_{\rm PE}$ provides an absolute normalization for $\lambda_{i}$. The simulation assumes a completely filled detector and the simulated events are generated on a grid with 20 points along the positive $x$ axis, 20 along the positive $y$ axis, and 40 symmetrically on the positive and negative $z$ axis. This choice results in 20 values for the radius $r=\sqrt{x^{2}+y^{2}+z^{2}}$. Assuming the spherical symmetry of the detector, events with $x<0$ and $y<0$ are modeled as the corresponding positive coordinate values. The $z$ direction is probed in both directions to account for the asymmetry introduced by the presence of the detector neck. The points are chosen to minimize their total number while probing the volume uniformly. Since these points lay on the three axes, the influence of individual direction-dependent effects can be tested. After the simulation, we define 20 cubic spline interpolated functions $\lambda_{i}(\cos\theta;r_{j}/r_{0})$ with the index $j$ corresponding to each of the radius values. Here $r=|\vec{x}|$ and $r_{0}=85$ cm is the inner radius of the vessel. For a given angle, a second cubic spline $\lambda(r/r_{0};\cos\theta)$ is defined to interpolate in the radius variable.

The time-of-flight algorithm for position reconstruction calculates the arrival time of a photon at each PMT. The algorithm uses a simplified optical model, illustrated in Figure 2.

For a scintillation interaction event taking place at event time $t_{0}$ and position $\vec{x_{0}}$ in LAr, the arrival time of a photon at a PMT can be calculated as

in which $T_{f}(\vec{p},\vec{x_{0}})$ is the time-of-flight from the event vertex $\vec{x_{0}}$ to position $\vec{p}$, the center of the light guide front surface. The time $\tau$ is a random variable based on three physical processes: singlet decay of argon with a time constant of 2.2 ns, TPB re-emission with a time constant of 1.5 ns, and a resolution term of 0.6 ns that includes PMT response, resolution, and variability in the propagation time through the light guides. These values were determined by matching the model to the measured time response based on reconstructed variables. The vector $\vec{p}$ rather than the physical PMT position is used both for the determination of time calibration constants and for reconstruction. When a VUV photon is wavelength-shifted by TPB at a point $\vec{x^{\prime}}$ on the acrylic surface, the resulting visible photon can be emitted into a light guide (in which case $\vec{p}=\vec{x^{\prime}}$) or be emitted back into the inner detector volume and travel with a greater group velocity before reaching a light guide front point $\vec{p}$. The time-of-flight can then be calculated,

where $\rm{v_{uv}}$ is the group velocity (133 mm/ns) of the 128 nm VUV light [19] and $\rm{v_{vis}}$ is the group velocity (241 mm/ns) of 420 nm visible light in LAr [20]. A test with the values of group velocities shifted by $\pm 5\%$ resulted in a negligible effect on the position reconstruction resolution. Reflections and scattering of visible photons in TPB and Rayleigh scattering in the LAr are neglected in this algorithm.

$P_{t}(t;r,\sin\theta)$ is the probability density function (PDF) that corresponds to the flight time from a particular position, with $r$ and $\sin\theta$ defined in Figure 2. The probability of measuring a photon in a PMT-light guide assembly at position $\vec{p}$ starting from a position $\vec{x_{0}}$ between times $t_{1}$ and $t_{2}$ is described by the following equations:

Eq. 3.5 gives the overall probability $P_{t}$ and assumes that the scintillation light and the TPB fluorescence emission are isotropic. The functions $H$ and $G$ impose the conditions that the time-of-flight must be between $t_{1}$ and $t_{2}$ and that the light from the TPB fluorescence strikes the light guide of interest. $d_{\mathrm{uv}}$ and $d_{\mathrm{vis}}$ are the distances traveled by VUV and visible light, respectively. $r_{\mathrm{lg}}=9.5\,\rm{cm}$ is the radius of the face of the light guide. Eq. 3.6 is solved to determine $d_{\mathrm{uv}}$ and to impose the geometry of the TPB surface at a radius of $r=85$ cm, which marks the termination of the VUV light. Eq. 3.7 is the distance between the TPB and PMT. For the time-of-flight calculation, we approximate that the distance to the center of the light guide can always be calculated. This is an approximation because the time-of-flight in the detector depends on where the photon strikes the face of the light guide. However, the largest uncertainty arises from fluorescence occurring directly in front of the PMT of interest, with a maximum value of approximately 1/3 ns $\left(r_{\mathrm{lg}}/\mathrm{v_{vis}}\right)$. This uncertainty is accounted for in the time resolution of the PMT-light guide system.

For the calculation of the solid angle $G$, we use the distance traveled by visible light between the TPB fluorescence and the acrylic vessel at the face of the light guide, $d^{\prime}_{\mathrm{vis}}$, and check if the point of impact is within the radius $(r_{\mathrm{lg}})$ of the face of the light guide.

The likelihood $\mathcal{L}(t_{0},\vec{x_{0}})$ for a given event time $t_{0}$ and position $\vec{x_{0}}$ is

where we sum over the $N$ pulses in the first 35 ns, with the detection time $t_{i}$, distance $r_{i}$, and angle $\theta_{i}$, based on the relevant PMT and vertex position. A time-of-flight correction is applied if the path from vertex to PMT passes through the gas layer on the top of the detector. The likelihood in Eq. 3.10 is maximized using the TMinuit2 fitter in ROOT [15].

Since the calculation of the PDF ($P_{t}$) requires significant time, we employ interpolation techniques to evaluate the PDF while fitting. The detailed calculation is performed at discrete points $r=0,1,2,...,170~\rm{cm}$, and at $\sin\theta=0,0.01,0.02,...,1.0$, and in 0.25 ns time bins. The PDFs are then convolved with a time distribution that accounts for the random variable $\tau$ from Eq. 3.3. We then use a cubic spline interpolation in time, followed by a bilinear interpolation in $r$ and $\sin\theta$ during the fit.

The left panel of Figure 3 shows the convolved PDF in a two-dimensional representation where the vertical axis is the time-of-flight, the horizontal axis reports the identification numbers (ID) of spatial cells defined as ID $=100\times\sin\theta+\textrm{\tt floor}(r[\textrm{mm}]/10)\times 100$ and the colour is the value of the PDF. The right panel shows a zoomed view around $r$ = 850 mm representing a $\sin\theta$ range from 0.00 to 0.99. For large values of $\sin\theta$, visible light arrives earlier than most VUV light (for cell numbers greater than 8560) because these cells are close to the surface.

The non-linearity of the position reconstruction problem in LAr points towards a solution employing machine learning (ML). The good performance of the hit pattern algorithm-based likelihood reconstruction suggests a similar approach using a supervised training algorithm. Another strong motivation for considering a ML approach is the possibility of training it with scintillation events in the LAr volume and events in other parts of the detector. In this context, the ML training includes events originating in the neck of the detector which constitute a relevant background for dark matter searches.

We investigated several algorithms, from Naive Bayes [21] to more complex models based on boosted decision trees [22] and neural networks. Our choice converged on a fully connected feed-forward neural network (FFNN) [23] for its better performance, relative simplicity, and training speed [24].

For the position reconstruction in the DEAP-3600 detector volume, the algorithm must realize a mapping between the 255 input charge values measured by the PMTs and the $(x,y,z)$ position coordinates as output. In practice, three single-output FFNN models were trained: one for each output Cartesian coordinate. Including time-of-flight information as additional FFNN inputs was also tested, but no significant improvements in performance were observed. The final network architectures of the selected FFNNs is reported in Table 1. The number of layers and neurons was determined by iteratively enlarging the network until no increase in the performance was observed.

The training of the FFNNs was performed using MC simulations of $10^{5}$ $\beta$-decays from ³⁹Ar, $10^{5}$ nuclear recoils on ⁴⁰Ar, and $10^{5}$ shadowed $\alpha$-decays from ²¹⁰Po happening on the acrylic argon flow-guides in the neck. These samples are divided into training and validation sets with proportions of $70:30$. The training process for each Cartesian coordinate minimizes the mean squared distance (MSD) between the predicted position $y_{\rm pred,i}$ and the true simulated position $y_{\rm true,i}$ for each event $i$:

where $N$ is the number of input events presented to the neural network during the training stage. Training for reconstructing the $x$ and $y$ coordinates required 400 epochs, while the $z$ coordinate required 40 epochs.

To compare the three algorithms in their ability to reconstruct the position of scattering events from $\beta$-decays of ³⁹Ar and ⁴⁰Ar recoils in LAr, we apply an appropriate selection to data and MC generated events. We consider events that have (i) good quality waveforms recorded by the digitizers, e.g. well-formed baseline voltage and valid trigger time, (ii) physics trigger signals only, excluding calibration/cosmics triggers, (iii) no pulse multiplicity consistent with two or more coincident physics events. To avoid light leakage from a previous event, the number of pulses registered in the first 1600 ns of the waveform must be less than or equal to three and the time between the event and the previous one must be greater than 20 $\mu$s. The reconstructed number of PE is selected to be within 90–200 PE, matching the WIMP region of interest (ROI) [7].

The left panels of Figure 4 show the normalized cubed radius distribution of reconstructed ${}^{39}\rm{Ar}$ ER events (data and MC) and ${}^{40}\rm{Ar}$ nuclear recoil events (MC) with the three algorithms and the right panels show the $z$ distributions. The hit pattern algorithm produces a more uniform distribution in the bulk region and is biased near the AV surface. In contrast, the distribution obtained with the time-of-flight algorithm is less smooth in the bulk while less biased near the surface. The enhancement of the distribution is more pronounced for the hit pattern algorithm since the probability distribution used in the likelihood is strongly non-linear when an event is very close to a PMT which collects most of the light.

Figure 5 shows the contained LAr mass for Hit Pattern and Time-of-Flight based algorithms for ³⁹Ar data and MC along with the curve calculated from the geometric volume of the detector.

Figure 6 shows the uniformity of reconstructed positions from the hit pattern algorithm for ${}^{39}\rm{Ar}$ data and ${}^{40}\rm{Ar}$ nuclear recoil MC data.

Time-of-flight and hit pattern position reconstruction algorithms use the hypothesis that an event comes from a single flash of isotropic light in LAr. Thus, the two algorithms should agree when the hypothesis is true, such as for WIMPs or ${}^{39}\rm{Ar}$ ER events. However, for neck events, the light is emitted above the LAr surface and shadowed by the flow-guides: the result is a very different time-of-flight distributions and a hit pattern concentrated at the bottom of the detector.

Figure 7 shows the difference between the two algorithms’ reconstructed $z$ positions and vertex positions. ${}^{39}\rm{Ar}$ data, simulated WIMP-like events (${}^{40}\rm{Ar}$ NR events), and simulated alpha events from the neck region are selected according to the WIMP ROI. It is shown that the reconstructed vertical positions agree within 35 mm for 50% of ³⁹Ar events in data and simulated WIMPs. For simulated neck alpha decays, the reconstructed vertical positions of the two algorithms are remarkably different, with the time-of-flight algorithm systematically reconstructing these events much closer to the top of the detector than the hit pattern algorithm.

In order to compare the three models in reconstructing the position of events in LAr and in the neck, the difference between the MC “true” position $\vec{x}_{\rm{true}}$ and the reconstructed position $\vec{x}_{\rm{rec}}$ are calculated. Events in LAr are defined as the region $-850\,\rm{mm}<z_{\rm{true}}\leq 550\,\rm{mm}$ while the neck region is $z_{\rm{true}}>850\,\rm{mm}$. The results are presented in Figure 8. The new ML algorithm is comparable to the existing hit pattern and time-of-flight algorithms in reconstructing the position in LAr, while it outperforms the other two in the reconstruction of events in the neck.

The possibility to employ a radioactive $\gamma$ source offers additional opportunities for testing the position reconstruction algorithms. In particular, the source signal can be precisely tagged, effectively eliminating other background events. Moreover, the source can be placed above the LAr surface for investigating a region which is not covered by ³⁹Ar decays. The DEAP-3600 detector has a calibration polyethylene tube (known as CAL F [5]) running around the steel shell holding the AV. The source is housed in a stainless steel capsule and can slide in the pipe and be placed at 9 pre-determined positions (see Figure 9) with an ex-situ measured uncertainty of $\sim 1$ cm.

The three $\gamma$s emitted almost simultaneously by the source allow for a tagging scheme where two photons are detected by LYSO crystals (32 photons/keV light yield) read out by PMTs on both sides of the source, while the third one is directed towards the LAr. The tagging algorithm relies on the energy deposited in the LYSO crystals, particularly on the timing between the event time in the LAr volume and the time when one of the tagging PMTs triggers. More details about this source are given in [5].

The energy spectrum of the $\gamma$s from the ²²Na source, interacting with the LAr volume, shows variations under three different tagging conditions (Figure 10).

1. 1.

   No tagging by the calibration PMTs: this allows the observation of the ³⁹Ar beta decay spectrum on top of the ²²Na events;

2. 2.

   Tagging by only one of the calibration PMTs: this occurs for instance, when the 1.27 MeV $\gamma$ is detected by one of the calibration PMTs, one 511 keV $\gamma$ escapes detection, and the other 511 keV $\gamma$ reaches the LAr volume, resulting in the observed peak at 511 keV ($\sim$3200 PE in Figure 10, green spectrum);

3. 3.

   Double tagging: this represents the most stringent condition, where both 511 keV $\gamma$s are tagged by the calibration PMTs, and the 1.27 MeV interacts with the LAr ($\sim$8000 PE in Figure 10).

The full energy peak of the 1.27 MeV $\gamma$ is visible for all the tagging conditions while the 511 keV peak is clearly visible only for the single tag configuration.

To study the performance of the hit pattern and time-of-flight reconstruction algorithms, MC simulations were performed with the ²²Na source at positions 2 and 4 (Figure 9). The chosen positions are below and above the LAr surface.

A total of 1.2$\times$10⁷ and 4.5$\times$10⁶ events were generated from the ²²Na capsule in position 2 and 4, respectively. The difference in the number of simulated events arises from variations in the number of events that trigger the detector, which depend on their position due to the detector’s geometry

Reconstructed coordinates for positions 2 and 4 are shown in Figures 11 and 12, for data and MC, both reconstructed with the hit-pattern and time-of-flight algorithms (normalizing to the same number of events).

A general good agreement between data and MC is observed for all coordinates but differences between the two algorithms can be noticed. From Figures 11(d) and 12(d), it is clear that the time-of-flight algorithm tends to reconstruct events closer to the center of the detector with respect to the hit pattern-based algorithm.