# Neutron and Gamma Pulse Shape Discrimination by Robust Determination of the Decay Shape
## 2 Materials and Methods
### Discrimination Methods
A cooling factor \(\beta\) with a value close to one (such as \(\beta=0.999\)) has to be defined as well, taking into account that the closer to one the slower the optimization process (or, in the analogy, the cooling process of the solid hopefully becoming a crystalline network). Finally, a small probability value \(p\) (such as \(p=0.01\)) has to be selected to allow for the occasional selection of worse solutions (configurations of higher energy in the analogy). With the initial solution \(x_{0}\), the residuals are computed by means of Equation (9) and the corresponding value for the objective function \(\rho(r)\) is obtained by means of Equation (11). These values are stored as \(x_{previous}\) and \(\rho_{previous}\), respectively, to start the following iterative process. They are also stored as \(x_{best}\) and \(\rho_{best}\). The initial movement amplitude is stored as \(\sigma_{previous}\). 1. A new solution \(x_{current}\) is obtained by taking the previous one as a base for a movement of amplitude \(\sigma_{current}\), where \[\sigma_{current}=\beta\sigma_{previous}\] (12) \[\Delta x\text{-}N(0,\sigma_{current})\] (13) \[x_{current}=x_{previous}+\Delta x\] (14) and where the increment \(\Delta x\) has been taken from the multivariate normal distribution with zero mean and \(\sigma_{current}\) standard deviation. It has to be checked that \(x_{current}\) belongs to the search domain, otherwise a new solution in the search domain needs to be obtained again by Equations (13) and (14). 2. The residuals are computed by Equation (9) for the current solution and then, by Equation (11), the value for the objective function, \(\rho_{current}\), is computed. If this value is smaller than the previous one \(\rho_{previous}\), the current solution is selected as base solution for the next iteration (that is \(x_{current}\) becomes \(x_{previous}\) and \(\rho_{current}\) becomes \(\rho_{previous}\)). Furthermore, if \(\rho_{current}\) is smaller than \(\rho_{best}\) then \(x_{current}\) is also stored as \(x_{best}\) and \(\rho_{current}\) as \(\rho_{best}\). Otherwise, if \(\rho_{current}\) is not smaller than the previous one \(\rho_{previous}\), then with the small probability \(p\), the new solution is accepted as base solution or, otherwise, it is rejected so that the solution of the previous iteration is kept for the next iteration. The process is stopped if the current movement amplitude \(\sigma_{current}\) is smaller than the desired final movement amplitude \(\sigma_{final}\). Otherwise, the algorithm turns back to step 1 above. 3. If the cooling process has been slow enough, at the end of the algorithm execution, the current solution must be coincident with the best solution (up to negligible differences of the order of \(\sigma_{final}\)). Furthermore, it is a necessary (though not sufficient) condition to attain the global minimum that repeated executions of the algorithm yield the same results (again, up to negligible differences of the order of \(\sigma_{final}\)). More details on the application of the simulated annealing to the computation of this robust estimator can be found in [18]. Furthermore, the computer code given in the Supplementary Materials includes this algorithm.