# Neutron and Gamma Pulse Shape Discrimination by Robust Determination of the Decay Shape
## 2 Materials and Methods
### Discrimination Methods
The steepest slope indicates a gamma-ray pulse. for the unknown vector \(x\)  
\[\mathbf{x}= \begin{pmatrix}d\mathbf{\tau}\\ \delta\mathbf{A}\end{pmatrix}  \tag{6}\]  
and the vector of observation residuals \(r\) to be determined after the adjustment. Then the characteristic decay time and residual amplitude of the pulse will be  
\[\mathbf{\tau}=\mathbf{\tau}_{approx}+\mathbf{d}\mathbf{\tau} \tag{7}\]  
and \(\delta A\), respectively. In Equations (4) and (5), the amplitudes and times for data points 1, 2, \(\ldots n\) in a given pulse are used, that is, \(A_{1}\), \(A_{2}\), \(\ldots A_{n}\), and \(\Delta t_{1}\), \(\Delta t_{2}\), \(\ldots\Delta t_{n}\), respectively. In this linearization, an approximate value \(\tau_{approx}\) has been used. This approximate value can be set either by experience, or by inspection of the pulse shape in a first determination, which can also be subsequently improved by iteration: one simple strategy is to use a certain time \(\Delta t\) and the corresponding amplitude \(A\) (for example, one like the lower threshold in Figure 2) and then obtain \(\tau\) from Equation (2) neglecting \(\delta A\). Other values can be obtained by evaluating a few points just to obtain a suitable value to start the algorithm. Thus, for the current dataset, an approximate value of 12 ns is obtained. The selection of a poor approximate value will have an adverse effect on the results of the least squares adjustment but not so much on the results after robust estimation, due to its high insensitivity to systematic errors in the model. The solution by _least squares_ is simple and well-known:  
\[\mathbf{x}= \begin{pmatrix}\mathbf{A}^{T}\mathbf{A}\end{pmatrix} ^{-1}\mathbf{A}^{T}\mathbf{k} \tag{8}\]  
\[\mathbf{r}=\mathbf{A}\mathbf{x}-\mathbf{k} \tag{9}\]  
\[\mathbf{C}_{\mathbf{x}}=\frac{\mathbf{r}^{T}\mathbf{r}}{\mathbf{n}-\mathbf{2}} \begin{pmatrix}\mathbf{A}^{ T}\mathbf{A}\end{pmatrix} ^{-1} \tag{10}\]  
where \(\mathbf{C}_{\mathbf{x}}\) is the covariance matrix of the unknowns, which has in the diagonal the squared values of the estimated standard errors of the unknowns, that is, \(\sigma_{\tau}\) and \(\sigma_{\delta A}\), respectively. For solution by robust estimation, it is proposed to use the minimum L\({}_{1}\)-norm computed by a global optimization method: the simulated annealing method, which is based on the analogy with the self-construction process of crystalline networks. Instead of minimizing the sum of squared residuals, now, the objective is to find a solution that minimizes the L\({}_{1}\)-norm of residuals, that is, the function \(\rho\) defined as  
\[\rho(\mathbf{r})=\sum_{i=1}^{n}\lvert\mathbf{r}_{i}\rvert \tag{11}\]  
To start the algorithm, an initial solution \(x_{0}\) is selected (sensible choices are the approximate value for tau \(\tau_{approx}\) mentioned before and zero for \(\delta A\), or their corresponding values obtained after the least squares solution). Then a search domain centered on this initial solution is defined: a certain number of times, 3 or 5, the standard deviation of the least squares solution could be a possible choice to define the boundaries of the search domain. It has to be taken into account that a too large search domain will require a slower, more thorough, search, while an excessively small search domain might leave outside the global optimum that is sought for. An initial movement amplitude \(\mathbf{\sigma_{0}}\) (a vector with the standard deviation of the least squares solution for each unknown or similar values that could be used to explore all regions of the search domain in a few steps) is also set up. A final movement amplitude \(\sigma_{final}\) is also selected, with a size that can be considered negligible for the problem to solve (\(\sigma_{final}=\sigma_{0}/100\) may be a reasonable choice).