Multi-channel, multi-template event reconstruction for SuperCDMS data using machine learning

TL;DR

This paper introduces a novel machine learning-based algorithm for event reconstruction in SuperCDMS detectors, effectively accounting for correlated noise and pulse shape variations to enhance energy resolution and position estimation.

Contribution

The paper presents the N×M filter algorithm and PCA-GBDT analysis pipeline, improving energy resolution and enabling position reconstruction in SuperCDMS data analysis.

Findings

Reduced energy resolution for mono-energetic events by over three times.

Achieved energy resolution below 50 eV for 71Ge K-shell electron-capture peaks.

Demonstrated improved noise handling and pulse shape analysis using machine learning.

Abstract

SuperCDMS SNOLAB uses kilogram-scale germanium and silicon detectors to search for dark matter. Each detector has Transition Edge Sensors (TESs) patterned on the top and bottom faces of a large crystal substrate, with the TESs electrically grouped into six phonon readout channels per face. Noise correlations are expected among a detector's readout channels, in part because the channels and their readout electronics are located in close proximity to one another. Moreover, owing to the large size of the detectors, energy deposits can produce vastly different phonon propagation patterns depending on their location in the substrate, resulting in a strong position dependence in the readout-channel pulse shapes. Both of these effects can degrade the energy resolution and consequently diminish the dark matter search sensitivity of the experiment if not accounted for properly. We present a new algorithm for pulse reconstruction, mathematically formulated to take into account correlated noise and pulse shape variations. This new algorithm fits N readout channels with a superposition of M pulse templates simultaneously - hence termed the N$\times$M filter. We describe a method to derive the pulse templates using principal component analysis (PCA) and to extract energy and position information using a gradient boosted decision tree (GBDT). We show that these new N$\times$M and GBDT analysis tools can reduce the impact from correlated noise sources while improving the reconstructed energy resolution for simulated mono-energetic events by more than a factor of three and for the 71Ge K-shell electron-capture peak recoils measured in a previous version of SuperCDMS called CDMSlite to $<$ 50 eV from the previously published value of $\sim$100 eV. These results lay the groundwork for position reconstruction in SuperCDMS with the N$\times$M outputs.

Tables

Table 1. Table 1 : Comparison of resolution between N × \times M with and without position correction and 1 × \times 1 OF.

Estimator	Peak resolution in %
1$\times$1 OF	2.56% $\pm$ 0.14%
N$\times$M OF without position correction (Elin-est)	2.39% $\pm$ 0.13%
N$\times$M OF with $G_{\text{BDT}}$ position correction (EN×M)	0.73% $\pm$ 0.04%

Table 2. Table 2 : Comparison of energy resolution of the N × \times M filter with other SuperCDMS energy estimators in SuperCDMS along with the 1 σ \sigma statistical uncertainties from the Gaussian fits in Figure  8 . NF stands for the non-stationary filter. The distribution of the resolutions from 500 bootstrapped samples are used to compare the K-, L- and M-shell peak resolutions.

Estimator	K shell $\sigma$ in eVee	L shell $\sigma$ in eVee	M shell $\sigma$ in eVee
E${}_{1\times 2}{}^{}{}_{}{}^{\mathrm{corrected}}$	100.42 $\pm$ 0.45	37.67 $\pm$ 0.16	16.33 $\pm$ 0.14
E${}_{\mathrm{NF}}{}^{}{}_{}{}^{\mathrm{corrected}}$	98.44 $\pm$ 0.46	35.73 $\pm$ 0.15	12.17 $\pm$ 0.09
EN×M	45.73 $\pm$ 0.15	33.43 $\pm$ 0.15	12.11 $\pm$ 0.06

Equations

$${\chi }^2(a)=\sum _{j=1}^S{\left(\frac{d_j-a{T}_j}{\sigma }\right)}^2.$$

$${\chi }^2(a)=\sum _{j=1}^S\sum _{k=1}^S(d_j-a{T}_j)(d_k-a{T}_k)(V^{-1}{)}_{jk}.$$

$${\chi }^2(a)=\sum _{f\ge 0}\frac{|{\tilde{D}}_f-a{\tilde{T}}_f{|}^2}{{\sigma }_f^2}.$$

$$\hat{a}=\frac{{\sum }_f{\tilde{T}}_f^{\ast }{\tilde{D}}_f/{\sigma }_f^2}{{\sum }_f{\tilde{T}}_f^{\ast }{\tilde{T}}_f/{\sigma }_f^2}.$$

$${\tilde{T}}_f(\delta )=e^{-2\pi if\delta }{\tilde{T}}_f(0).$$

$$\mathbf{V}={\mathbf{V}}_{\mathrm{noise}}+{\mathbf{V}}_{\mathrm{shape}}.$$

$$d_{\mathrm{total},j}=\sum _{c=1}^N{w}_c{d}_{c,j},$$

$${\chi }^2(E)=\sum _{c=1}^N{\chi }_c^2=\sum _{c=1}^N\sum _{f\ge 0}\frac{|{\tilde{D}}_{c,f}-a_c(E){\tilde{T}}_{c,f}{|}^2}{{\sigma }_{c,f}^2}.$$

$$\left[\begin{array}{cccc}\mathrm{cov}(\Re (z_1),\Re (z_1))& \mathrm{cov}(\Re (z_1),\Im (z_1))& \mathrm{cov}(\Re (z_1),\Re (z_2))& \mathrm{cov}(\Re (z_1),\Im (z_2))\\ \mathrm{cov}(\Im (z_1),\Re (z_1))& \mathrm{cov}(\Im (z_1),\Im (z_1))& \mathrm{cov}(\Im (z_1),\Re (z_2))& \mathrm{cov}(\Im (z_1),\Im (z_2))\\ \mathrm{cov}(\Re (z_2),\Re (z_1))& \mathrm{cov}(\Re (z_2),\Im (z_1))& \mathrm{cov}(\Re (z_2),\Re (z_2))& \mathrm{cov}(\Re (z_2),\Im (z_2))\\ \mathrm{cov}(\Im (z_2),\Re (z_1))& \mathrm{cov}(\Im (z_2),\Im (z_1))& \mathrm{cov}(\Im (z_2),\Re (z_2))& \mathrm{cov}(\Im (z_2),\Im (z_2))\end{array}\right].$$

$${\chi }^2(a)=\sum _f\sum _{c_1}\sum _{c_2}({\tilde{D}}_{c_1,f}-a{\tilde{T}}_{c_1,f})({\tilde{D}}_{c_2,f}-a{\tilde{T}}_{c_2,f}{)}^{\ast }{V}_{c_1{c}_2,f}^{-1}.$$

$$p_{c,s}=\sum _{i=1}^M{a}_{ci}{T}_{ci,s}({\delta }_{ci})=\sum _{i=1}^M{a}_{ci}{T}_{ci,s}({\delta }_c).$$

$${\vec{p}}_c={\vec{a}}_c\cdot {\mathbf{T}}_c({\delta }_c),$$

$${\overrightarrow{\stackrel{~}{P}}}_c={\vec{a}}_c\cdot {\tilde{\mathbf{T}}}_c\cdot \mathbf{R}({\delta }_c).$$

$${\Delta }_{c,f}={\tilde{D}}_{c,f}-{\tilde{P}}_{c,f}={\tilde{D}}_{c,f}-({\vec{a}}_c\cdot {\tilde{\mathbf{T}}}_c\cdot \mathbf{R}({\delta }_c){)}_f.$$

$${\chi }_c^2=\sum _f{\chi }_{c,f}^2=\sum _f\frac{|{\Delta }_{c,f}{|}^2}{{\sigma }_{c,f}^2}=\sum _f\frac{[\Re {({\Delta }_{c,f})]}^2+[\Im ({\Delta }_{c,f}){]}^2}{{\sigma }_{c,f}^2}$$

$${\chi }^2=\sum _f{\chi }_f^2=\sum _f\sum _{j=1}^N\sum _{k=1}^N\sum _{\alpha =\mathrm{Re},\mathrm{Im}}\sum _{\beta =\mathrm{Re},\mathrm{Im}}{\Delta }_{j\alpha ,f}{\Delta }_{k\beta ,f}({\mathbf{V}}^{-1}{)}_{j\alpha ,k\beta ,f}.$$

$${\tilde{D}}_f^T=\left(\begin{array}{c}\Re (chan1\mathrm{FFT})\\ \Im (chan1\mathrm{FFT})\\ \Re (chan2\mathrm{FFT})\\ \Im (chan2\mathrm{FFT})\\ ⋮\\ \Re (\mathrm{chanN}\mathrm{FFT})\\ \Im (\mathrm{chanN}\mathrm{FFT})\end{array}\right).$$

$${\tilde{P}}_f^T=\left(\begin{array}{c}\Re (chan1\mathrm{prediction}\mathrm{FFT})\\ \Im (chan1\mathrm{prediction}\mathrm{FFT})\\ \Re (chan2\mathrm{prediction}\mathrm{FFT})\\ \Im (chan2\mathrm{prediction}\mathrm{FFT})\\ ⋮\\ \Re (\mathrm{chanN}\mathrm{prediction}\mathrm{FFT})\\ \Im (\mathrm{chanN}\mathrm{prediction}\mathrm{FFT})\end{array}\right).$$

$$\vec{a}=(a_{11},a_{12},...a_{1M},a_{21},a_{22},...a_{2M},...a_{N1},a_{N2},...a_{NM})$$

$${\tilde{\mathbf{T}}}_f=\left(\begin{array}{c}\begin{array}{cccccccc}\Re ({\tilde{T}}_{11})& \Im ({\tilde{T}}_{11})& & & & & & \\ \Re ({\tilde{T}}_{12})& \Im ({\tilde{T}}_{12})& & & & \text{0}& & \\ ⋮& ⋮& & & & & & \\ \Re ({\tilde{T}}_{1M})& \Im ({\tilde{T}}_{1M})& & & & & & \\ & & \Re ({\tilde{T}}_{21})& \Im ({\tilde{T}}_{21})& & & & \\ & & \Re ({\tilde{T}}_{22})& \Im ({\tilde{T}}_{22})& & & & \\ & & ⋮& ⋮& & & & \\ & & \Re ({\tilde{T}}_{2M})& \Im ({\tilde{T}}_{2M})& & & & \\ & & & & \ddots & & & \\ & & & & & \Re ({\tilde{T}}_{N1})& \Im ({\tilde{T}}_{N1})& \\ & & & & & \Re ({\tilde{T}}_{N2})& \Im ({\tilde{T}}_{N2})& \\ & \text{0}& & & & ⋮& ⋮& \\ & & & & & \Re ({\tilde{T}}_{NM})& \Im ({\tilde{T}}_{NM})& \end{array}\end{array}\right).$$

$$\begin{array}{rl}{\mathbf{R}}_f(\vec{\delta })& ={\mathbf{R}}_f({\delta }_1,{\delta }_2,\cdots ,{\delta }_N)\\ & =\left(\begin{array}{cccc}\begin{array}{cc}\cos 2\pi {\delta }_1& \sin 2\pi {\delta }_1\\ -\sin 2\pi {\delta }_1& \cos 2\pi {\delta }_1\end{array}& & & \text{0}\\ & \begin{array}{cc}\cos 2\pi {\delta }_2& \sin 2\pi {\delta }_2\\ -\sin 2\pi {\delta }_2& \cos 2\pi {\delta }_2\end{array}& & \\ & & \begin{array}{cc}\ddots & \ddots \\ \ddots & \ddots \end{array}& \\ \text{0}& & & \begin{array}{cc}\cos 2\pi {\delta }_N& \sin 2\pi {\delta }_N\\ -\sin 2\pi {\delta }_N& \cos 2\pi {\delta }_N\end{array}\end{array}\right).\end{array}$$

$$\begin{array}{rl}{\chi }^2& =\sum _f{\chi }_f^2\\ & =\sum _f{\left({\overrightarrow{\stackrel{~}{D}}}_f-{\overrightarrow{\stackrel{~}{P}}}_f\right)}^T{\mathbf{V}}_f^{-1}\left({\overrightarrow{\stackrel{~}{D}}}_f-{\overrightarrow{\stackrel{~}{P}}}_f\right)\\ & =\sum _f{\left({\overrightarrow{\stackrel{~}{D}}}_f-\vec{a}\cdot {\tilde{\mathbf{T}}}_f\cdot {\mathbf{R}}_f({\delta }_c)\right)}^T{\mathbf{V}}_f^{-1}\left({\overrightarrow{\stackrel{~}{D}}}_f-\vec{a}\cdot {\tilde{\mathbf{T}}}_f\cdot {\mathbf{R}}_f({\delta }_c)\right).\end{array}$$

$${\vec{a}}^T={\left(\sum _f\left({\tilde{\mathbf{T}}}_f\cdot {\mathbf{R}}_f({\delta }_c)\right)\cdot {\mathbf{V}}_f^{-1}\cdot {\left({\tilde{\mathbf{T}}}_f\cdot {\mathbf{R}}_f({\delta }_c)\right)}^T\right)}^{-1}\sum _f\left({\tilde{\mathbf{T}}}_f\cdot {\mathbf{R}}_f({\delta }_c)\right)\cdot {\mathbf{V}}_f^{-1}\cdot {\overrightarrow{\stackrel{~}{D}}}_f^T.$$

$$E_{\mathrm{lin}-\mathrm{est}}=\sum _{i=0}^{NM-1}{P}_i\cdot {\alpha }_i,$$

$${\chi }^2=\sum _j{\left[E_{\mathrm{true}}^j-\sum _i\left(P_i^j\cdot {\alpha }_i\right)\right]}^2,$$

$${\alpha }_i=\frac{{\sum }_j{P}_i^j\cdot E_{\mathrm{true}}^j}{{\sum }_j{\left(P_i^j\right)}^2}.$$

$$G_{\text{BDT}}=\frac{E_{\mathrm{true}}-E_{\mathrm{lin}-\mathrm{est}}}{E_{\mathrm{lin}-\mathrm{est}}}.$$

$$E_{\mathrm{N}\times \mathrm{M}}=E_{\mathrm{lin}-\mathrm{est}}\cdot \left(1+G_{\text{BDT}}\right).$$

$$E_t=E_r+n_{eh}eV.$$

$$E_{\text{NF}}^{\mathrm{corrected}}=E_{\text{NF}}^{\mathrm{uncorrected}}\cdot \left(1+C\cdot [A_{\mathrm{fast}}-A_{\mathrm{fast}}^{\mathrm{mean}}]\right),$$

Full text

1 Introduction

Experimental data often take the form of time series. For example, the SuperCDMS experiment [1, 2] digitizes phonon signals in Transition Edge Sensors (TESs) that are deposited on the surfaces of cryogenic semiconductor detectors [3]. Some SuperCDMS detectors also digitize charge signals produced by ionization inside the detectors. These digitized time series contain pulses whose amplitudes are proportional to the amount of energy deposited in the detector. The pulse shapes are determined by the location of the energy deposit inside the detector, the physics of phonon propagation, TES behaviour, and the response of the electronics. We seek to fit templates/shapes to the time series data to determine underlying parameters of the initiating event. Traditionally, the SuperCDMS experiment estimates the energy deposited by particles scattering in its detectors by fitting a fixed pulse shape to the digitized time series and then converting the fitted amplitude to an equivalent energy – the so-called optimal filter [4]. In this paper, we generalize this approach to the case of pulse shapes that vary with the position of the energy deposits which are measured in multiple sensors containing both independent and common-mode (correlated) noise components. The general approach can be adapted to any experiment that analyzes time series data with varying pulse shapes and/or multiple sensors.

Section 1 of the paper introduces optimal filtering and describes the complications that arise when dealing with multiple channels, pulses with a priori unknown start times, and shape variations. Section 2 describes a computationally tractable method for addressing these complications. The remaining sections of the paper describe how this method is applied in detail to SuperCDMS, including validation with simulation and data.

1.1 Fitting a time series with an optimal filter

Let $d_{j}$ denote a discrete time series with $S$ equally spaced time samples at times $t_{j}=j\Delta t$, where $\Delta t$ is the sampling interval. Let $T(t)$ denote a template: a pulse shape of normalized amplitude that we wish to fit to the time series to determine the amplitude of the measured pulse. In the case of Gaussian white noise, the measurement uncertainty $\sigma_{j}$ on each measured point is constant and statistically independent of the uncertainty on any other point. In this restricted circumstance, one can write a simple time-domain $\chi^{2}$ for the fit as

[TABLE]

Here $T_{j}\equiv T(t_{j})$, $a$ is an unknown amplitude we wish to fit for, and the sum is over all $S$ time samples. The best-fit amplitude $\hat{a}$ is the value of $a$ that minimizes this $\chi^{2}$. In general, different points in the time series may have correlated uncertainties, i.e. their covariance matrix $V_{jk}\equiv{\rm cov}(d_{j},d_{k})$ is not necessarily diagonal. In this case, one can generalize the $\chi^{2}$ to account for a non-diagonal covariance matrix as follows:

[TABLE]

If the covariance matrix $V$ between measured time samples is known and if it can be efficiently inverted, this formulation may suffice. In practice, for a time series with $S$ samples, the covariance matrix, of size $S\times S$, may be too large to practically invert. Instead, it is usually computationally more efficient to formulate the problem in the frequency domain. Let $\tilde{D}_{f}$ and $\tilde{T}_{f}$ denote the discrete Fourier transforms of the time series $d_{j}$ and $T_{j}$. It is common for the statistical properties of the measurement noise to be invariant as a function of time (“stationary”). If the measured noise is stationary, then the Fourier amplitudes at different frequencies will be statistically independent: ${\rm cov}(\tilde{D}_{f_{1}},\tilde{D}_{f_{2}})=0$ for $f_{1}\neq f_{2}$. In this case a $\chi^{2}$ expression in the frequency domain takes a simple form:

[TABLE]

Here the sum is over frequencies $f=(1/\Delta t)k$, where $k=0,1,...~(S/2)$, and the term in the denominator, $\sigma_{f}^{2}$, is the power spectral density (PSD) of the measurement noise. Note as well that although a Fourier transform yields both positive- and negative-frequency components, for a real time series $d_{j}$ the negative-frequency components of its Fourier transform are redundant with their positive frequency counterparts, related through complex conjugation $(\tilde{D}_{-f}=\tilde{D}_{f}^{*})$. Hence, the sum in the $\chi^{2}$ can be restricted to non-negative frequencies with no loss of information.

Working in the frequency domain inherently requires the use of complex numbers, yet their incorporation into the $\chi^{2}$ to yield a real value must be done with care. The form of the numerator in Equation 3 can be better understood by recalling that, for a complex number $z$, $|z|^{2}=[\Re{(z)}]^{2}+[\Im{(z)}]^{2}$. For $z=\tilde{D}_{f}-a\tilde{T}_{f}$ being the difference between the data and scaled template in the frequency domain, Equation 3 can be interpreted as saying that the real and the imaginary components of this difference can be treated as two statistically independent data points that each contribute to the $\chi^{2}$. The statistical independence of the real and imaginary components follows from the fact that the noise is assumed to be stationary. One can also understand that nothing is lost in restricting the sum in the $\chi^{2}$ to positive frequencies: the number of independent real and imaginary degrees of freedom in $\tilde{D}_{f}$ for $f\geq 0$ is exactly equal to the number of time samples $S$ in the time domain111One might worry that because there are $S/2+1$ frequencies in the sum, each of which contributes a real and an imaginary part, Equation 3 contains $S+2$ random variables. However for a real time series the first and last terms in the sum will always have imaginary components of zero, i.e. $\Im{(D_{0})}=\Im{(D_{f_{\mathrm{max}}})}=0$, thus reducing the number of free components to the expected $S$..

Because Equation 3 is quadratic in $a$, the value of $a$ that minimizes the $\chi^{2}$ can be found analytically by taking its derivative with respect to $a$ and setting it equal to zero, yielding

[TABLE]

This expression for $\hat{a}$ is referred to as the optimal filter (OF) estimate for the amplitude $a$. It is optimal since, by the Gauss-Markov theorem, it is an unbiased estimator for $a$ with the least possible statistical variance of all linear functions of the data.

1.2 Dealing with time offsets between the data and the template

The above discussion assumed that the fitted template $T(t)$ is aligned in time with the pulse in the time series data. If however the template is offset by some time shift $\delta$ relative to the pulse, one can account for the shift by making the template a function of $\delta$, writing $T_{j}(\delta)=T(t_{j}-\delta)$ in the time domain, or $\tilde{T}_{f}(\delta)={\rm DFT}[T_{j}(\delta)]$ for the frequency domain version, where DFT is the Discrete Fourier Transform. This shifted template $\tilde{T}_{f}(\delta)$ can be substituted into Equation 3, with the $\chi^{2}$ then becoming a function of both $a$ and $\delta$. Unfortunately, the dependence on $\delta$ is non-linear, and the minimum cannot be found analytically. It can however be found numerically by scanning over a range of $\delta$ and minimizing the $\chi^{2}$ at any particular value of $\delta$ with respect to $a$. Equation 4 still holds provided that one substitutes $\tilde{T}_{f}(\delta)$ for $\tilde{T}_{f}$, with $\hat{a}$ now becoming a function of $\delta$. The best-fit $\hat{a}$ is then the one corresponding to the global minimum of $\chi^{2}$ with respect to $a$ and $\delta$.

For the case of $\delta$ being an integer multiple of the sampling period $\Delta t$, $\tilde{T}_{f}$ has a particularly simple form:

[TABLE]

In this case, the Fourier transform of the time-shifted template can be readily obtained from the Fourier-transformed unshifted template simply by rotating each component by a complex phase $-2\pi f\delta$. It is important to emphasize that this relation only holds for the case that $\delta$ is an integer multiple of $\Delta t$. Moreover, Equation 5, which can be derived from the expression of the discrete Fourier transform, implicitly assumes that $T(t)$ is periodic with period $S\Delta t$. It has the effect of “wrapping around” the template — any part of the template shifted beyond the end of the sampling period appears back at the start of the template, as shown in the bottom part of Figure 1. This does not cause problems if the templates go to zero at both ends (top of Figure 1), but can result in template discontinuities otherwise (bottom of Figure 1). In cases where the time shift is not an integer multiple of $\Delta t$ or where Equation 5 introduces wraparound effects for time shifts of interest, one must instead directly calculate $T(t_{j}-\delta)$ and explicitly take its Fourier transform.

1.3 Dealing with pulse shape variation

Section 1.1 assumed a fixed pulse shape with no variation, for which the overall amplitude was the only free parameter. If the template shape $T(t)$ assumed in the fit is not an accurate representation of the average pulse shape, a simple one-parameter fit for the amplitude will typically be biased. If the pulse shape varies from event to event, then fitting every event with a single fixed pulse shape will tend to broaden the amplitude resolution due to the unmodeled variation in shape. In the context of the SuperCDMS experiment, dependence of the phonon pulse shape on the position of the energy deposit inside the detector produces a noticeable smearing of the energy resolution if no effort is made to account for the position variation.

The approach in Section 1.1 can be readily generalized to the case of template shapes that depend on additional parameters. In this case the $\chi^{2}$ can be minimized with respect to all free parameters. Provided that the parameterized templates accurately model the variation of pulse shapes in the data, the additional parameters provide additional degrees of freedom that help reduce the deviation between the fitted and true pulse shape.

If the dependence of the pulse shape model on the fitted parameters is not linear, the best-fit parameters that minimize the $\chi^{2}$ cannot generally be solved for analytically, although numerical minimization is usually possible. The addition of a relative time shift between the data and fitted template, described in Section 1.2, is an example in which the model is linear in one parameter (amplitude $a$) but non-linear in a second (time shift $\delta$). In some cases, a simple reparameterization can make a problem linear. For example, fitting a sine wave with unknown phase and amplitude can either be done as a non-linear fit $T(t)=a\sin(\omega t+\phi)$ for $a$ and $\phi$, or equivalently as a linear fit $T(t)=A\sin\omega t+B\cos\omega t$ for $A$ and $B$.

An alternative to attempting to fit for a varying pulse shape is to consider any variation in pulse shape to be a systematic uncertainty in the fit and to include that systematic uncertainty as part of the covariance matrix. For example, the covariance matrix $V_{jk}$ for time samples in Equation 2 can be expressed as a noise contribution and a contribution from pulse shape variation:

[TABLE]

Here, ${\bf V}_{\mathrm{shape}}$ contains the covariances between time bins due to any pulse shape variation. Equation 6 can be substituted into Equation 2 and the resulting $\chi^{2}$ is minimized. This approach does not attempt to fit for any parameters that describe the shape of the pulse, but effectively “deweights” those features of the time series that are uncertain due to the underlying variation in pulse shapes. The additional contribution to the covariance matrix can be thought of as an additional uncertainty that will tend to broaden the fitted distribution of $a$ relative to the case in which there is no variability in pulse shape. However, the broadening from including $\bf{V}_{\mathrm{shape}}$ is typically smaller than the broadening in resolution caused by fitting a fixed template shape to pulses that do not all follow that shape.

While stationary noise results in a diagonal noise matrix expressed in the frequency domain, the additional covariance contribution from ${\bf V}_{\mathrm{shape}}$ can have any form, and in general the combined covariance matrix will no longer be diagonal when the fit is reformulated in the frequency domain. That is, the pulse shape uncertainty can be thought of as a source of non-stationary noise. Equation 3, which assumes that there is no covariance between the measurement uncertainties of the real and imaginary parts of the Fourier amplitudes at different frequencies, thus does not hold, and there is little benefit from working in the frequency domain rather than the time domain in this case.

The SuperCDMS experiment has traditionally used a combination of Equations 2 and 6 as its most reliable simple energy estimator, and refers to this as a “non-stationary optimal filter” [5, 6]. This combination is found to give better energy resolution than a simple stationary optimal filter which assumes a fixed pulse shape.

Computationally, approaches in which the pulse shape variation can be captured by a linear function of a small number of fit parameters are much more efficient than either a non-linear fit or attempting to absorb the effects of shape variation as an additional uncertainty that results in a non-diagonal covariance matrix.

1.4 Dealing with multiple channels

Data may consist of multiple time series, either from separate detectors (e.g. jointly fitting gravitational wave data from physically separate detectors) or from multiple channels within a detector (e.g. the segmented phonon sensors on a SuperCDMS detector, each covering a different portion of the detector’s surface). If the time series from different channels are aligned in time (or can be straightforwardly aligned by shifting to correct for “time of flight” effects), a simplistic approach may be to sum the time series and then fit the summed pulse:

[TABLE]

where $d_{c,j}$ is the time series data for channel $c$, $w_{c}$ is a weight factor for that channel, and $N$ is the number of channels. This summed time series $d_{\mathrm{total},j}$ can then be used in any of the fits given in previous sections. This approach, while easy to implement, is not guaranteed to give optimal results, especially if the weights $w_{c}$ are not chosen carefully.

SuperCDMS has traditionally calculated its energy estimators using the summed time series from all of a detector’s phonon channels. Not only is this approach easy, it also tends to mitigate (although not eliminate) energy resolution broadening due to position dependence across the detector, because adding the signals from all channels tends to average away this position dependence. This is especially true in the non-stationary optimal filter, where a summed pulse is created first to mitigate position dependence to first order and then the rest of the position dependence is treated as correlated (non-stationary) noise [5].

A better approach is to form a $\chi^{2}$ as in Equation 3 for each channel, and then to sum these $N$ $\chi^{2}$ terms to form a combined $\chi^{2}$ that can be minimized as a function of any free parameters. Imagine that each channel has a separate template $T_{c}$ and amplitude $a_{c}(E)$ that depends on a free parameter $E$, which we can think of as being an energy. We can then generalize Equation 3 to the case of multiple channels as

[TABLE]

This expression accounts for the possibility that the template shape, amplitude, and noise PSD $\sigma_{f}^{2}$ may differ from channel to channel.

The approach of simply adding the $\chi^{2}$ expressions from separate channels implicitly assumes that there are no correlated uncertainties between different channels. However, it may not be appropriate when combining data from different channels inside one detector. For example, noise induced by mechanical vibrations or electronic interference is likely to appear as correlated noise across all channels, and in this case Equation 8 will not correctly account for these correlations. Simply summing the raw time series as per Equation 7 has the same issue: while uncorrelated noise will tend to average out in the sum, noise that is correlated between channels will not cancel in the sum and does not benefit from any averaging effect.

To properly account for correlated noise, one must include the covariances between measured quantities. However, in the context of Equation 3, in which the measured quantities are expressed as the complex Fourier amplitudes $\tilde{D}_{f}$, particular care must be taken. Consider two complex random variables $z_{1}$ and $z_{2}$. Fully characterizing the covariances of these variables requires ten different real numbers:

[TABLE]

In contrast, problems arise if one tries to directly calculate covariances between complex numbers, i.e. ${\rm cov}(z_{i},z_{j})$. The proper definition of variance of a complex random variable is ${\rm E}[|z-{\rm E}[z]|^{2}]$ rather than ${\rm E}[(z-{\rm E}[z])^{2}]$, and is thus a real number [7]. The covariance of two random variables is defined as ${\rm E}[(z_{1}-{\rm E}[z_{1}])(z_{2}-{\rm E}[z_{2}])^{*}]={\rm E}(z_{1}z_{2}^{*})-{\rm E}(z_{1}){\rm E}(z_{2}^{*})$. A $2\times 2$ covariance matrix $V_{ij}={\rm cov}(z_{1},z_{2})$ formed using these definitions from two complex numbers can be parameterized by only four real numbers, illustrating that this covariance matrix does not fully capture all of the possible correlations. To fully characterize all of the covariances between complex random variables, one must also include pseudo-covariances, defined as ${\rm pseudocov}(z_{1},z_{2})\equiv{\rm cov}(z_{1},z_{2}^{*})={\rm E}[(z_{1}-{\rm E}[z_{1}])(z_{2}-{\rm E}[z_{2}])]$. Together, the $2\times 2$ covariance and pseudo-covariance matrices fully capture all possible covariances among components of the complex variables. Unfortunately, the form of a likelihood or $\chi^{2}$ in this general case is quite complicated [7].

In the special case of complex variables that have circular symmetry (i.e. are uniformly distributed in their complex phase), the pseudo-covariance matrix will be identically zero, and there is a simple form for the $\chi^{2}$ including correlations among channels as was derived in [6]:

[TABLE]

Here, the summation is over frequency bins and over channels indexed by $c_{1}$ and $c_{2}$. Note that the covariance matrix $V_{ij,f}={\rm cov}(\tilde{D}_{c_{1},f},\tilde{D}_{c_{2},f})$ is an $N\times N$ covariance matrix that is a function of frequency, and that $V_{ij,f}=V_{ji,f}^{*}$ (ensuring that the above expression is purely real). Equation 10 assumes that the noise is uncorrelated between frequencies, i.e. at any given frequency the noise may be correlated between channels, but the noise at any given frequency is assumed to be uncorrelated with the noise at any other frequency.

In order to avoid relying on the assumption of circularly symmetric noise, and to avoid the complexities of dealing with both covariance and pseudo-covariance matrices, in this paper we will formulate our $\chi^{2}$ expressions using only real numbers, treating each complex random variable as two real random variables. Therefore, if we have a vector of $N$ complex amplitudes, we rewrite this as a real vector with $2N$ components and express their covariances with a $2N\times 2N$ matrix along the lines of Equation 9.

1.5 The approach in this paper

A general solution to the problems of pulse shape variation and correlated noise is to formulate a generic $\chi^{2}$ that encodes all possible correlations, and to allow the fitted signal templates to be functions (possibly non-linear) of the fit parameters. Unfortunately, both aspects of this approach present considerable computational difficulties. With $S$ samples and $N$ channels of data, the covariance matrix capturing any possible correlations between data points would have a size of $SN\times SN$. For anything but very short time series the inversion of this matrix will likely be computationally intractable.

In this paper, we make the simplifying assumption that noise is uncorrelated between frequencies but allow for correlations between channels at any given frequency. This simplification effectively means that rather than having to invert a single matrix of size $SN\times SN$, we instead have to invert $\mathcal{O}$($S/2$) matrices of size $2N\times 2N$ (one for each non-negative frequency in the discrete Fourier transform), where each covariance matrix is a purely real matrix along the lines of Equation 9, encoding the correlations among the real and imaginary components of $\tilde{D}_{f}$. Due to the much smaller matrix size, this approach is computationally tractable, even if it must be done for each frequency.

Accounting for pulse shape variation requires a different approach. An accurate parameterization of the pulse shape as a function of the underlying physics parameters would be ideal, but developing this parameterization may be non-trivial, especially if the detector physics is complicated and not analytically calculable. Furthermore, if the pulse shape varies non-linearly with these parameters, then the $\chi^{2}$ cannot generally be minimized analytically, increasing the computational load. In this paper, we take the approach of empirically modelling the varying pulse shapes as linear combinations of different templates, derived from data, which capture the most significant variations. The end result is that we fit time series from $N$ channels with linear combinations of $M$ different templates, resulting in a total of $NM$ fitted amplitudes222Throughout this paper we use $N\times M$ to refer to the fitter, while $NM$ is the number of fitted amplitudes.. Section 2 describes the mathematical formalism to perform this fit while taking into account the measured covariances between channels at each frequency.

The $N\times M$ amplitudes encode information about the characteristics of the underlying pulse, such as its energy and position. However, the parameterization employed is empirical, and we need to relate the fitted amplitudes to parameters of physical interest, such as the energy and position of the event inside a SuperCDMS detector. Section 3 describes how machine learning techniques can be used to extract meaningful pulse templates and to map the fitted amplitudes to the energy and position of the energy deposit inside the detector.

2 The N$\times$M optimal filter

Suppose that we have $N$ channels to which we fit $M$ templates with different shapes. We can write the prediction $p_{c,s}$ in the time domain of the data of channel $c$ for the $s^{\textrm{th}}$ time bin sample as:

[TABLE]

Here, the sum over $i$ denotes a sum over templates $T_{ci,s}$, and $a_{ci}$ and $\delta_{ci}$ are the amplitude and time delay for the $i^{\textrm{th}}$ template. From a physics standpoint, we expect that the templates of the same channel share the same starting point, which is dependent on the event interaction time and signal propagation speed. Consequently, in the last equality we impose the constraint that templates of the same channel have identical time delays, thus $\delta_{ci}=\delta_{c}$.

Similar to the single-channel, single-template OF case, the $\chi^{2}$ has a nonlinear dependence on $\delta_{c}$. As the number of channels increases (beyond four), the $\chi^{2}$ minimization involving $\delta_{c}$ becomes computationally challenging. For instance, if each of twelve channels allows five possible time-delay values, a brute-force grid search would require evaluating $5^{12}$ combinations. We address this issue by first extracting the time delay for each channel using standard single-channel, multiple-template optimal filter. This process can be solved analytically using an inverse Fourier transform, as described in [8]. Ideally, the template would be shifted by the amount of the time delay. However, we opted to shift the pulse data instead to reduce memory consumption. For time delay shifts smaller than 40 $\mu$s, the two operations are effectively equivalent for frequencies below 25 kHz, where the majority of the pulse information resides.

Equation 11 implies an equal number of templates ($M$) across channels, albeit this is primarily a mathematical formality. In practice, one can first identify the maximum number of templates required for any channel and set $M$ to this maximum number. Channels with fewer templates can be accommodated by including additional templates whose amplitudes are fixed to be zero.

The final output of the N$\times$M optimal filter is therefore $NM$ amplitude coefficients $a_{ci}$. The correlation among these coefficients can be very complicated, and the ultimate energy or position estimators will be expressed in a non-linear mapping of $a_{ci}$, as demonstrated in Section 3.

Equation 11 can be written in matrix form as

[TABLE]

where $\vec{p}_{c}$ is a $1\times S$ vector containing the predicted time series of channel $c$, $\vec{a_{c}}$ is an amplitude vector of size $1\times M$, and ${\bf{T}}_{c}$ is a template matrix of the size $M$ $\times$ $S$ containing the $M$ templates for channel $c$, each containing $S$ samples. Taking a Fourier transform of Equation 12 and expressing the time shift of the templates explicitly using Equation 5 (for a time shift adequately approximated by an integer number of sampling periods) yields

[TABLE]

Here, $\vec{\tilde{P}}_{c}$ is the Fourier transform of the predicted time series for channel $c$, with a length given by the number $F$ of independent frequencies in the Fourier transform333For $S$ even, $F=S/2+1$, as discussed in Section 1.1.. The $M\times F$ matrix $\bf{\tilde{T}}_{c}$ contains the Fourier transforms of the $M$ templates evaluated at a time shift of $\delta_{c}=0$. $\bf R(\delta_{c})$ is an $F\times F$ diagonal matrix that applies the time shift to the templates by rotating each frequency component by $e^{-2\pi if\delta_{c}}$, i.e. $R_{f_{1}f_{2}}=\delta_{f_{1}f_{2}}e^{-2\pi if_{1}\delta_{c}}$. All of the matrices in Equation 13 are complex matrices.

Let us define the difference between measurement and prediction as $\Delta$, so that at some frequency $f$

[TABLE]

The $\chi^{2}$ for channel $c$ can be written as

[TABLE]

where $\sigma^{2}_{c,f}$, a scalar, is the (variance) noise of channel $c$ at frequency $f$. Expanding the $\chi^{2}$ to include multiple channels means that the correlations between different channels, as well as between the real and imaginary components of the same channel, need to be taken into account.

For a given frequency, a real covariance matrix of size $2N$ $\times$ $2N$ is needed to encode the covariances among the real and imaginary components of the data from the $N$ channels. Note that some off-diagonal elements of this matrix may be zero, e.g. between the real and imaginary terms of the same channel if the noise is stationary444The general observation from data is that these terms are close to but not exactly zero – either because the stationary noise assumption is not strictly met or because of statistical variation.. Let us denote this covariance matrix as ${\bf V}_{j\alpha,k\beta}$, where $\alpha$ and $\beta$ can represent either the real or imaginary parts of the data for the channels $j$ and $k$, respectively. Likewise, we denote the real and imaginary parts of $\Delta_{c}$ with a Greek subscript. The total $\chi^{2}$ can then be expressed as

[TABLE]

Ignoring the sum over frequencies, the remaining summations are over channels (indices $j$ and $k$) and the real/imaginary components of the Fourier amplitudes (indices $\alpha$ and $\beta$). Note that this expression corresponds to the inverse covariance matrix ${\bf V}_{f}^{-1}$ (dimension $2N\times 2N$, an entirely real matrix) multiplied on either side by a vector of length $2N$ representing the real and complex components of the residuals $\Delta_{f}$. Equation 16 accounts for all correlated noise between channels at any frequency, but assumes that noise at different frequencies is uncorrelated.

The exact format of $\Delta$ and ${\bf V_{\it f}}$ in Equation 16 can be arbitrarily defined in principle. Here we illustrate one possible option which writes everything in terms of purely real quantities. Let $\tilde{D}_{f}$ be a row vector of length $2N$, which contains data from all channels at frequency $f$ in sequential order, alternating between real and imaginary components. Its transpose $\tilde{D}_{f}^{T}$ looks like

[TABLE]

The transpose of the corresponding $\tilde{P}_{f}$ row vector is

[TABLE]

Similarly, we can redefine terms in Equation 13 to accommodate multiple channels. First we define an amplitude vector $\vec{a}$ that contains all amplitudes ($NM$ amplitudes in total):

[TABLE]

We define ${\bf\tilde{T}}_{f}$ to be a matrix of size $(NM\times 2N)$ that holds the Fourier transform (real and imaginary parts) of the templates, so that $\vec{a}\cdot{\bf\tilde{T}}_{f}$ is a $1\times 2N$ vector that contains the predicted values for the real and imaginary parts of the Fourier amplitudes for each channel, in the same format as Equation 18. Each row of ${\bf\tilde{T}}_{f}$ corresponds to one of the $NM$ amplitudes and contains the real and imaginary Fourier amplitudes at frequency $f$ for each channel resulting from that amplitude. Given how the amplitudes are defined in terms of channels and templates, only one channel (two columns) of each row is non-zero:

[TABLE]

Note that the templates here are assumed to be those before any time shifts are applied.

Finally, we define a $2N\times 2N$ time shift matrix ${\bf R}_{f}(\vec{\delta})$, which is a function of the time delay vector $\vec{\delta}$ containing the time shifts for all channels. This is a block diagonal matrix of size $2N\times 2N$ where each block is a $2\times 2$ rotation matrix and all other entries outside these $2\times 2$ blocks are zero:

[TABLE]

Each $2\times 2$ block is determined by a $\delta_{i}$, which is the time delay for that channel. ${\bf R}_{f}$ simply applies the phase rotation to each channel.

With the above definitions, we can write the $\chi^{2}$ from Equation 16 in matrix form:

[TABLE]

This $\chi^{2}$ can be analytically minimized with respect to the amplitudes, with the minimum at

[TABLE]

3 Applications in SuperCDMS

SuperCDMS employs cryogenic detectors (cylindrical Si or Ge crystals, approximately 10 cm in diameter and $\sim$3 cm thick) for direct detection of dark matter [1, 2]. Quantized vibrations, called phonons, produced by an incoming particle scattering off an electron or nucleus in the crystal are detected as the primary source of information about the interaction energy. Phonon signals are collected bolometrically in multiple sensors on the top and bottom of the detectors, using TESs to convert small thermal excursions into measurable currents which are later amplified and then digitized.

Event reconstruction via optimal filters described in the previous sections is our standard method to extract these quantities from a pulse occurring in the waveforms. Many factors affect our reconstruction ability, several of which were described in Section 1. The rest of this paper discusses how to apply the N$\times$M filter introduced in Section 2 within the context of SuperCDMS. Our ultimate analysis goal is to obtain a position-corrected energy estimator. The discussion covers the topics of template generation for the N$\times$M filter and using these templates to fit for a total of $NM$ amplitudes as given by Equation 23. We will also present the use of machine learning techniques to remap the amplitudes to meaningful physical quantities of interest: energy and position. The overall data processing pipeline is shown in Figure 2 with more details in the following sub-sections.

3.1 Template generation for the N$\times$M filter

The N$\times$M filter allows us to fit multiple ($M$) shapes to each channel (see Equation 23). There are two aspects to template generation that require consideration. The first is how one ought to generate the templates. The templates generated must (a) be reasonably independent of detector noise, and (b) capture the variation in pulse shape across the detector. These requirements help prevent unwanted biasing of the filter outputs by ensuring that the templates can accurately represent pulse shapes.

Furthermore, the template shapes carry useful information about the location of events in the detector. Charge and phonon transport processes result in primary and secondary populations of phonons arriving at the sensors at different times555For the purposes of time alignment, all templates are constructed with a common starting time, representing the arrival time of the earliest phonons.. Prompt phonons preserve information about the location of the original energy deposit in the crystal and contribute a fast component to the pulse shape. Position reconstruction in SuperCDMS detectors is achieved by translating the resulting pulse-shape differences into quantities that correlate with position. In the past, SuperCDMS has used a type of optimal filter called the two-template optimal filter (1$\times$2 OF) which effectively separated the energy and position information in a pulse into amplitudes derived from slow and fast templates, respectively [6, 9]. The amplitude of the fast component primarily provided position information, while the amplitude of the slow template provided an energy estimate. For this fitter, the slow template was constructed from the average pulse shape, while the fast template was calculated as the average of residuals between many measured pulses and the slow template after inverting those residuals with negative peaks.

The N$\times$M filter generalizes and systematizes this approach, and so ought to be able to capture the energy and position information more effectively than other methods. The N$\times$M filter also addresses the issue of correlated noise between channels and the fact that in newer detector designs the pulse shapes cannot be so cleanly separated into a fast and a slow component.

The second aspect of template generation concerns the number of templates. In principle, one can imagine fitting a very large number of templates to a single pulse to obtain a near-perfect match to the measured shape. Such a fit would be limited by computational resources and possibly by degeneracies among the fit parameters. Furthermore, physical interpretation of the results may not be straightforward when using a large number of empirically derived templates.

With both of these aspects in mind, we present the use of principal component analysis (PCA) [10] to derive effective templates for an N$\times$M fit. The advantage of this approach is two-fold: (a) it presents an elegant mathematical structure that extracts the most important features in the data, and (b) it provides a convenient way to determine the number of useful templates. PCA enables identification of the most significant dimensions in a multi-dimensional data set (often termed “feature identification”). In addition, it provides a framework to achieve dimensionality reduction by decomposing a data set into components (“feature reduction”). We exploit PCA feature identification to obtain templates and feature reduction to restrict the number of templates used in the N$\times$M filter.

Principal component analysis starts with a data set of $n$ events, each having $p$ numbers associated with it. A basis vector in this $p$-dimensional parameter space is derived that maximizes the variance of the projection of the events onto this vector. This procedure is equivalent to finding the direction in the parameter space that explains the largest amount of variance in the data. After this first principal component is found, its projection is subtracted from each event to remove the variance due to that component, and then the procedure is repeated to find the next basis vector which maximizes the remaining variance in the data. This process is continued iteratively until a full set of $p$ basis vectors is found. The procedure is equivalent to finding the eigenvalues of the covariance matrix of the data. The eigenvectors are called the principal components (PCs), and their associated eigenvalues, arranged in descending order, identify how much of the overall variance in the data can be associated to each component. In this paper, we will refer to the first component (the one with maximum information or variance) as component 0, the second component as component 1 and so on.

Data from the new detectors built for SuperCDMS SNOLAB consist of digitized waveforms (called traces) sampled at a frequency of 625 kHz (2.5 MHz) for phonon (charge) channels. Each time bin of the waveform contributes information about the position, energy, and noise of an event. A waveform containing $p$ time bins therefore has $p$ features. To determine templates for the N$\times$M filter, we apply PCA to a sample of $n$ traces each with a trace length of $p$ samples, aligned with respect to the start point of the pulse in the trace. The start point is typically configured to occur near the midpoint of the trace. Alignment of traces is necessary to avoid capturing artifacts in the PCs due to small differences in pulse start times. To ensure that the PCs capture pulse-shape variation across the full detector, we compute a common set of templates (to be used for all detector channels) from digitized events that include traces from all detector channels. We normalize the templates by the maximum height of their respective PCs.

PCA produces $p$ PCs for a data set with $p$ features. $p$ corresponds to the trace length in our case and it can be as large as 32768 (131072) samples for phonon (charge) traces corresponding to a 52 ms trace. It is therefore important to determine the optimal balance between number of templates and amount of useful information in the templates. We produce a scree plot [11] which displays the amount of information in each component as a function of component number. The $y$ axis is the fractional contribution of the eigenvalue or variance of that component to the total variance in the data set. Examples of PCs and the corresponding scree plot for a simulated data set (which is described in Section 4.1) are shown in Figure 3. We normalize each PC in Figure 3 so its maximum height is equal to 1. We perform a qualitative assessment of the PCs and the scree plot to determine a final template selection.

We aim to combine the amplitudes across different channels to obtain a position-corrected energy estimator. We achieve this via a two-step process which first calculates a linear estimator from the N$\times$M amplitudes and then corrects for higher-order effects with supervised machine learning. This procedure is described in Sections 3.2 and 3.3.

3.2 Energy reconstruction with linear estimator

The $M$ normalized PCs obtained from the procedure described in Section 3.1 are used to generate amplitudes for each of the $N$ channels using the N$\times$M filter. In all, the N$\times$M filter produces a total of $NM$ amplitudes which represent the relative contribution of each PC. Each amplitude contains information about both the energy and position of the event, with each channel’s component 0 containing most of the energy information. A weighted sum of the channel amplitudes provides a reasonable, first-order estimate of the reconstructed energy:

[TABLE]

where $P_{i}$ is the $i^{th}$ amplitude from the N$\times$M fit and $\alpha_{i}$ is the corresponding weight for that component. (Note that this sum is over both channels and principal components.) This estimate corrects for linear dependencies, such as relative calibration666Relative calibration factors are weights applied to each channel to account for unequal TES response owing to differences in bias points, manufacturing tolerance, etc. of the channels, by appropriately weighting the contribution of each component as follows. The weights, $\alpha_{i}$, are obtained by minimizing the following $\chi^{2}$ summed over a set of training events:

[TABLE]

where $E_{\mathrm{true}}^{j}$ is the true energy for event $j$ (for example, simulated events of known energy or events selected in a calibration peak). The subscript $i$ and superscript $j$ on $P_{i}^{j}$ denote the $i^{th}$ amplitude component for the $j^{th}$ event. Equation 25 can be analytically minimized to give the weight of the $i^{th}$ amplitude component:

[TABLE]

These weights can be calculated provided we have a data set with events whose true energies are known. Calibration data sets with readily identifiable features can be used for this purpose. Generally speaking we use a mono-energetic set of events, such as from a calibration peak, to determine a set of weights.

When applied to a set of mono-energetic events, Equations 24 and 25 implicitly include a partial correction for position dependence in the estimated energy. Position dependence in the fitted amplitudes will tend to broaden the energy resolution. By minimizing the $\chi^{2}$ in Equation 25, we produce a set of weights that yields a linear estimator that minimizes the energy resolution, including an effective linear correction for any position dependence of the amplitudes. These weights will also reflect the noise levels in the detector, since amplitudes that are particularly affected by noise will be de-weighted as a result of this procedure.

It follows that the optimal weights (those which minimize Equation 25) are energy dependent: ${\alpha(E)}=\alpha_{0}+\beta(E)$, where $\alpha_{0}$ is a constant and $\beta(E)$ is some function of energy. This energy dependence results from the fact that the relative contributions of noise and position dependence vary with energy, with position-dependent effects dominating over noise at higher energies. The $\beta(E)$ term accounts for these position effects. Because the linear weights depend on both the noise levels in individual channels and on position dependence of the pulse amplitudes and shapes, which are not known a priori, the weights to be applied for any given event must be derived from data. We have devised an iterative method to approximate the weights from data using calibration peaks with known true energies (such as mono-energetic signatures resulting from decays of 71Ge). A preliminary energy estimate is calculated by setting the weights to one for every event. This is sufficient to allow us to identify peaks in the energy spectrum. Then, using events from the calibration peaks, we estimate the weights at those energies. We assume that the same weights can be applied to events within 5$\sigma$ of the calibration peak. The preliminary energy estimates for events outside the peaks are then used to linearly interpolate the weights, and those interpolated weights are then used in Equation 24 to generate an improved energy estimate. Events whose preliminary energy estimates fall above or below the 5$\sigma$ region of the first and last calibration peaks, respectively, are assigned the same weights as the nearest calibration peak. The use of only a small number of known energy points for the linear interpolation of the weights introduces systematic uncertainty, which we acknowledge as a limitation of the technique that can be studied in more detail with future simulation datasets. Note that we have introduced an additional requirement – the need for multiple calibration peaks in the data set – to obtain optimal energy resolution from the N$\times$M filter.

3.3 Machine learning reconstruction of energy and position from $N\times M$ amplitudes

It may be possible to further improve the resolution performance relative to Equation 24. The linear estimate is corrected for correlated noise implicitly with the non-diagonal covariance matrix in Equation 22, and the resolution broadening induced by position dependence is corrected to first order as described above. However, the reconstructed energies can manifest non-linear dependencies on the event location inside the crystal. The nature of these higher-order terms is nearly impossible to determine analytically even with the availability of accurate simulations whose data are indistinguishable from experimental data. In cases such as this, where the model that maps a set of variables to an output is unknown, one can use machine learning to evaluate an empirical model. We use Gradient Boosted Decision Trees (GBDTs) [12] from the LightGBM package [13] to estimate non-linear corrections to energy from position dependence. The use of PCA templates constructed in the manner described in Section 3.1 allows us to capture position information in the amplitudes. If these amplitudes are used as input variables, the GBDT combines the amplitudes in an optimal fashion which can improve the resolution of a chosen calibration peak.

Although decision trees are well known in the context of classification problems, here we use them to perform regression tasks. Let us consider the case of a one-dimensional mapping of a variable $X$ to another variable $Y$. Assume that the exact nature of the relationship between $X$ and $Y$ is unknown and therefore needs to be calculated empirically. GBDTs work by finding the ideal way to split the given data in $X$ into categories such that each category minimizes the prediction error for $Y$. Mean Squared Error is the commonly used cost function to evaluate the efficacy of each split. In multi-dimensional data, i.e. a mapping of a set, $\{X\}$, to an output, $Y$, the GBDT algorithm attempts to divide the multi-dimensional space so as to minimize the cost function summed over events in the training set. An individual tree that performs such a regression does not necessarily produce the best results and is often termed a weak predictor. Therefore, the method is expanded to an ensemble of decision trees, each of which on its own is a weak predictor, but jointly act as a strong predictor [14]. Each successive tree learns from the errors of its predecessors via gradient descent. GBDTs have the general advantage of being more interpretable compared to neural networks, even though visualizing a tree for high-dimensional data may not be straightforward.

LightGBM GBDT is effective on comparatively small data sets [15], which makes it a suitable algorithm to employ in a rare event search experiment such as SuperCDMS. We use LightGBM GBDT in the following way. The $NM$ amplitudes from the N$\times$M filter are rescaled to be of similar order of magnitude to ensure that no single amplitude is unfairly valued by the algorithm. The rescaled amplitudes serve as the features that the GBDT trains on. The goal is to map the $NM$ amplitudes to an energy correction factor. We call this energy correction factor $G_{\textrm{BDT}}$, which is meant to supplement the typical linear position gain corrections employed in SuperCDMS [16]. $G_{\textrm{BDT}}$ is defined as the fractional deviation of true energy from $E_{\mathrm{lin-est}}$:

[TABLE]

We assume that this correction factor is independent of $E_{\mathrm{true}}$, corresponding to the case that the correction results from non-uniformities in phonon collection efficiency that are independent of the true energy. This assumption allows us to generalize a sample trained on a single calibration peak to other energies. Use of training samples with a larger range of energies, which are currently not available, would allow us to relax this assumption. Among the hyperparameters that define the structure of the GBDT are the number of trees in the ensemble, the maximum number of leaves in the trees, the depth of the tree, and the minimum number of samples in a child node for the next split. A grid search on a validation sample was employed to find the optimal number of trees, leaves and depth.

Calibration peaks provide events with known energies which can be used to derive $G_{\textrm{BDT}}$-trained models. Obtaining position labels is not possible for real data sets, as we do not know the true positions of events in the detector. However, SuperCDMS simulations based on GEANT4 can produce both position- and energy-labeled data sets. Although it may not be fully representative, we can use the simulated data set (described in Section 4.1) to demonstrate a proof of concept. In real data, the training is performed on events from one calibration peak, validated on a statistically independent sample of events from the same calibration peak, and then tested on a completely different calibration peak to evaluate generalizability across the full energy spectrum. The final corrected energy for a new data point after training is given by

[TABLE]

The fundamental assumption in training any machine learning model is that the training and test data sets are independently derived from the same underlying distribution. In our stable detector operation modes, such as in a dark matter data set, we fix important physics parameters such as the TES bias points and the detector voltage. Although the noise environment could potentially change with time, the covariances and weights used in the N$\times$M filter can readily be updated to optimally handle the noise. However, the operational bias points and pulse shapes of each detector vary based on the particulars of its design, material, and fabrication variances, which in turn manifest themselves in our input $NM$ amplitudes. Therefore, a separate GBDT needs to be trained for each detector; however, once trained, the same model can be applied to all data sets acquired under the same operating conditions for the same detector. Moreover, the hyperparameters that define the structure of each GBDT need to be tuned separately for each detector.

For most real data sets, our inputs are the $NM$ amplitudes and the output to be trained is the value of $G_{\textrm{BDT}}$. In the case of simulated data, we know the true position of the events in addition to their energies, which allows us to map the amplitudes to position as well. This is not the case if training on real data whose position labels are unknown. For simulated data, one can slightly alter the model to map the $NM$ features to a $G_{\textrm{BDT}}$ value and three position parameters corresponding to the three geometric axes of the detector. However, a model trained on simulated data is not guaranteed to produce best results on real data unless the simulation parameters are tuned for that detector.

The pipeline in Figure 2 was tested on data drawn from simple analytic models and on simulated events and finally applied to real data from CDMSlite [17], a study done with a previous generation of SuperCDMS detectors operated in the Soudan mine. In the following sections, we describe results from applying the data processing pipeline to GEANT4 simulations and CDMSlite data.

4 Validation of the $N\times M$ analysis framework

In Section 4.1 we present results from applying the N$\times$M data processing pipeline to simulated events, followed by results from applying the pipeline to CDMSlite data in Section 4.2.

4.1 Validation with GEANT4 simulations

We used the SuperCDMS Detector Monte Carlo (DMC) framework to simulate 10,000 mono-energetic (0.5 keV) events distributed uniformly throughout the volume of a SuperCDMS SNOLAB-style High Voltage (HV) detector operated at 100 V ($\pm$50 V on either face of the detector) [18]. HV detectors produce an amplified phonon signal proportional to the applied voltage ($V$) via the Neganov-Trofimov-Luke (NTL) effect [19, 20, 21],

[TABLE]

Here, $E_{r}$ is the original recoil energy of an electron for a measured energy $E_{t}$ and $n_{eh}=E_{r}\cdot Y(E_{r})/\epsilon$ is the number of electron-hole pairs produced in the interaction777$\epsilon$ is the ionization energy. $Y(E_{r})\equiv 1$ for electron interactions and is an energy-dependent factor $<1$ for nuclear recoils.. HV detectors are projected to achieve excellent sensitivity to low mass WIMPs. They contain 12 phonon channels and do not have any direct charge readout. The simulation uses GEANT4 [22] to generate particle hits and G4CMP [23] to simulate phonon and charge transport in the crystal, and combines these with the detector response to produce simulated events. The noise PSD from a real data set is used to generate noise traces to add to the simulated pulses. The noise described by a PSD does not account properly for correlations between channels and hence the noise data used in this example is not correlated. The DMC runs an extensive pipeline which produces raw traces in the data format used by the experiment.

4.1.1 Energy Reconstruction

The raw traces produced from the simulation are run through the full processing pipeline shown in Figure 2. Following the pipeline, we first determine the N$\times$M templates. Figure 3 shows that only the first three PCA components contain significant information about the pulse shape. Therefore, we perform a 12 channel $\times$ 3 template fit to obtain 36 amplitudes. The N$\times$M amplitudes produced are first used to construct $E_{\mathrm{lin-est}}$, and then the data set is split 90/10 to perform supervised learning. A training sample consisting of 90% of the simulated events is used to train the GBDT, with the remaining 10% used to assess the performance of the model in terms of the width of the reconstructed peak. Note that unlike a real data set, here we are restricted to a single mono-energetic peak for training and testing. The resolution of the peak on the test sample is estimated by fitting a Gaussian to a histogram of the reconstructed energies. A comparison to the width of the peak obtained from the simple 1$\times$1 OF described in Section 1 is shown in Table 1. We find that the N$\times$M filter without GBDT corrections is statistically equivalent to the simplest OF applied to the sum of phonon channels. This observation is attributed to the fact that the events used in the simulation do not have any added correlated noise. Hence the statistical treatment of correlated noise in Equation 23 does not provide significant enhancements in energy resolution compared to the simple 1$\times$1 OF. While E${}_{\textrm{lin-est}}$ corrects for linear effects of position dependence, this appears to be a sub-dominant contribution in this particular data set. Subsequent application of $G_{\textrm{BDT}}$ using Equation 28 provides a factor of $\sim$3.3 improvement in resolution of the test data showing that position effects in these data are predominantly non-linear in nature. A comparison of the spectrum in each case is shown in Figure 4.

4.1.2 Position Reconstruction

It is straightforward to verify that the N$\times$M amplitudes contain position information. Figure 5 shows the dependence of some of the N$\times$M amplitudes on position. A channel map of the simulated detector is shown in Figure 6. In general, an event located closer to a certain channel deposits more energy in that channel compared to others. This effect is clearly visible in Figure 5. For example, events near the outer ring channel (PAS1) deposit the least energy in the central channel (PFS1) and vice versa (see Figure 6). The high level of non-linear dependence of amplitudes on position is also evident in these plots. Furthermore, the correlation of each component amplitude with true event position is non-linear and different for different channels, which justifies our use of the GBDT to perform a non-linear correction to the N$\times$M energy estimator.

An additional advantage of working with simulated data is that, unlike real data, we can train on position labels to reconstruct true position from the N$\times$M amplitudes. We use the same GBDT but supply it with 3 additional labels to reconstruct position: X, Y and Z. We find that the GBDT is able to reconstruct X and Y position with excellent accuracy (within 0.52 and 0.55 mm, respectively) as shown in Figure 7. We also find that Z position reconstruction is poor (prediction error is about 5.49 mm), which is not surprising given the relatively low level of correlation in Figure 5 between the amplitudes and true Z positions. This lack of correlation is expected in HV detectors because the NTL effect rapidly produces additional phonons along a column in the Z direction, which smears out Z-position sensitivity. Inclusion of phonon arrival times from an appropriately defined time delay quantity has been shown to help reduce this effect. Furthermore, the effectiveness of applying a machine learning model trained on simulated data to a real data set remains to be verified. We emphasize once again that at this point in time, the most realistic simulations available for SuperCDMS detectors are not fully representative of real pulses. Additional work is required to extend the application of this position reconstruction approach to SuperCDMS SNOLAB.

4.2 Validation with CDMSlite data

CDMSlite was an experiment conducted as part of SuperCDMS Soudan. It involved running the conventional SuperCDMS detector at the time, the iZIP, which had both phonon and charge readout, in “CDMSlite mode” [16, 17, 24, 25]. The detector was grounded on one side, and a high voltage of up to $V=-75$ V was applied on the other side to produce NTL amplification.

Three successful runs with iZIP detectors in this mode inspired the design and development of HV detectors for SuperCDMS SNOLAB [16, 24, 25]. We have applied the N$\times$M pipeline in Figure 2 to data from CDMSlite Run 3 [25]. CDMSlite Run 3 operated a Ge detector for a period of about 5 months. 70Ge, when exposed to a neutron source, produces 71Ge, which decays via electron capture with a half life of 11.43 days [24, 26]. The subsequent x-ray and Auger-electron emissions generate three distinct peaks corresponding to the shell from which the electron was captured: the K-shell peak at 10.37 keV, the L-shell-peak at 1.30 keV and the M-shell peak at 160 eV. These peaks are used for the energy calibration of the detector.

CDMSlite used the Non-Stationary Filter (NF) described in Section 1.3 as the final energy estimator for its dark matter analysis, after correcting for various environmental effects (e.g. temperature) with gain factors and a simple correction for position dependence using the amplitude of the fast amplitude from the 1$\times$2 OF (described briefly in Section 3.1) [16]. The main part of the analysis which differs for the N$\times$M energy estimator is the position correction. In the CDMSlite dark matter analysis, the NF estimator, $E_{\textrm{NF}}$, was found to correlate with the fast amplitude and therefore a linear correction was devised to decorrelate the two, thus effectively correcting for position effects. This correction took the form

[TABLE]

where $A_{\mathrm{fast}}$ denotes the fast amplitude with the superscript “mean” denoting the mean fast amplitude obtained from the K-shell events, $C$ is a constant, and the superscript “uncorrected” denotes the estimator prior to position correction. The final optimized estimator, $E^{\mathrm{corrected}}_{\textrm{NF}}$, was used in the dark matter analysis. Similarly, for comparison in this study, we also construct the 1$\times$2 OF estimator, E${}^{\mathrm{corrected}}_{1\times 2}$, after similar corrections are applied.

Here we attempt to compare our best position-corrected estimator from the N$\times$M, i.e., $G_{\textrm{BDT}}$-corrected $E_{\mathrm{lin-est}}$ (EN×M), with two other energy estimators: NF ($E^{\mathrm{corrected}}_{\textrm{NF}}$) and 1$\times$2 OF (E${}^{\mathrm{corrected}}_{1\times 2}$), used in CDMSlite. The sensors of the Soudan iZIP detectors were sectioned into 8 phonon readout channels, compared to 12 in SuperCDMS SNOLAB detectors (see Figure 6), of which 4 were not read out in the CDMSlite measurements to accommodate application of a large voltage bias on the detector. PCA analysis showed that most of the information was contained in the first 4 PCs. Therefore, we performed a 4$\times$4 fit to the CDMSlite data. The 16 amplitudes produced from this processing were only corrected for environmental effects that showed a correlation with them. They are then combined using Equation 24 with weights estimated from Equation 26.

The final step is to correct for position dependence. Since we do not have position labels to train the GBDT, we do not attempt to reconstruct the events’ absolute positions. Energy corrections with $G_{\textrm{BDT}}$ were calculated with the GBDT. Unlike the simulated data set in Section 4.1, in this case we have multiple independent energy peaks to train, validate, and test our machine learning model. The smaller number of features in this data set (16 amplitudes) allows us to obtain reasonable results even with only thousands of training data points. We select events in the Ge K-shell activation peak in the spectrum to train the GBDT. The events in the peak are selected following the same procedure used in the published analysis [25], resulting in 1984 events. Of the selected events, we reserve about 100 events ($\sim$5% of the data) to perform an independent assessment of the K-shell peak resolution after training.

We can increase the size of the training sample by rotating the amplitudes for each event among the channels. We make use of the rotational symmetry of our detector (see Figure 6) by rotating the data of the three inner channels by $0^{\circ}$, $120^{\circ}$, or $240^{\circ}$, and effectively tripling the training data size. The augmented data set thus has rotational symmetry built into it which allows the GBDT to learn this important property of our data. For SuperCDMS SNOLAB, one can use rotational symmetries combined with symmetry about the vertical axis to flip sides of the detector, thereby increase the training data size by up to a factor of 6 for HV detectors. The trained model is also applied to the other two calibration peaks which served as test data sets.

A caveat in this procedure is the limited sample size in the L- and M-shell peaks. To truly assess the efficacy of this technique, the validation must be performed on independent test samples. Test samples must be independent to validate the machine learning and, separately, to validate the weighted linear-sum estimation of the N$\times$M pipeline. While sufficient training data points ($\sim 5100)$ are available for the K-shell peak, significantly fewer events are available for the L- and M-shell peaks for weight calculation (a factor of ten lower than the K-shell events). Therefore, to produce a fair comparison of the techniques, we use the “bootstrap” method to produce 500 data sets by sampling the event with replacement. In each bootstrapped data set, 100 events are reserved as an independent test set and the remaining are used to calculate the weights $\alpha_{i}$ according to Equation 26. Each bootstrapped sample is then run through the full N$\times$M pipeline and the final calibration peak from the test sample is fit with a Gaussian each time. The resulting distributions of the widths of the calibration peaks are shown in Figure 8. One of the estimators used for comparison is the non-stationary filter, NF, which was explained earlier. The other estimator is the slow component from the 1$\times$2 OF which predominantly contains energy information (described in Section 3.1). Both of these estimators are corrected for linear position dependence using the residual amplitude estimated from the 1$\times$2 filter (see Equation 30).

In Table 2, we compare the mean energy resolution of the calibration peaks obtained from the fitted Gaussian mean of the bootstrapped distribution of energy resolutions for the three mentioned estimators. We find that the BDT-corrected N$\times$M filter shows the best resolution: about a factor of 2.6 better than the 1$\times$2 estimator for the K-shell peak, up to $\sim$8% better for the L-shell peak and about 3% better for the M-shell peak. The small amount of training data available to calculate the weights (a factor ten lower than the K-shell training sample) potentially results in underfitting of the weights for the L and M shells. With larger training samples, resulting in better determination of the weight factors, the resolution performance of the N$\times$M OF for the L- and M-shell peaks could likely be improved. The position uncertainty in CDMS detectors scales as a function of energy and therefore the amount of broadening from position effects decreases with energy. Consequently, the amount of improvement from the BDT correction is expected to decrease with energy.

4.3 Comment on training samples

Most supervised learning algorithms are data hungry: deep learning typically calls for very large data sets to produce a usable model [27]. Rare event search experiments such as SuperCDMS unfortunately are often limited in the size of the calibration and search data sets available for training. While we are able to demonstrate significant improvements in energy resolution by our choice of a data-efficient algorithm, we are still limited by small training samples at low energies. Moreover, the training procedure we described in this paper results in data being discarded from the training peak, which is not ideal. Although simulated data sets could in principle be an option to provide independent training data, they are not guaranteed to be fully representative of real-world dark matter data because (a) differences arising from built-in assumptions can skew the simulated training data away from real detector data, and (b) tuning simulation input parameters to match real data sets is an involved process that itself may require large validation data sets. Model training in machine learning requires that the test data set be similar in characteristics and features with the training data for reliable performance. Current research to integrate machine learning into the SuperCDMS data analysis pipeline follows a two-pronged approach: (a) develop more realistic simulated data that closely reproduce detector physics, and (b) explore applicability of machine learning models trained on currently available simulated data sets to existing real data sets.

5 Conclusion

We have developed a method for accounting for correlated noise between channels and position-dependent variations in pulse shape in an optimal filter energy estimator. These effects can be captured by simultaneously fitting data from $N$ channels with $M$ templates derived from principal component analysis, yielding $NM$ amplitudes in total. A combination of a linear estimator based on these fitted amplitudes and a gradient boosted decision tree can be used to generate energy or position estimates for events in a SuperCDMS-style detector. We have demonstrated the validity of this method with simulated data and with data from an older SuperCDMS detector, showing significant improvements in energy resolution.