update quad() to not abort if the matrix is singular

2 files changed, 30 insertions, 5 deletions

diff --git a/src/fit.c b/src/fit.c
index ef22d44..8fa3fad 100644
--- a/src/fit.c
+++ b/src/fit.c
@@ -5248,17 +5248,18 @@ int fit_event2(event *ev, double *xopt, double *fmin, int *id, size_t n, double
     size_t npeaks;
     size_t *result;
     size_t nvertices;
+    int status;
 
     /* Create the minimizer object. */
     opt = nlopt_create(NLOPT_LN_BOBYQA, 4+3*n);
     nlopt_set_min_objective(opt,nlopt_nll2,&fpars);
 
     /* Guess the position and t0 of the event using the QUAD fitter. */
-    quad(ev,pos,&t0,10000);
+    status = quad(ev,pos,&t0,10000);
 
-    if (t0 < 0 || t0 > GTVALID || NORM(pos) > PSUP_RADIUS) {
-        /* If the QUAD fitter returns something outside the PSUP or with an
-         * invalid time we just assume it's at the center. */
+    if (status || t0 < 0 || t0 > GTVALID || NORM(pos) > PSUP_RADIUS) {
+        /* If the QUAD fitter fails or returns something outside the PSUP or
+         * with an invalid time we just assume it's at the center. */
         fprintf(stderr, "quad returned pos = %.2f, %.2f, %.2f t0 = %.2f. Assuming vertex is at the center!\n",
                 pos[0], pos[1], pos[2], t0);
         pos[0] = 0.0;
diff --git a/src/quad.c b/src/quad.c
index c869aed..848a12e 100644
--- a/src/quad.c
+++ b/src/quad.c
@@ -14,9 +14,16 @@
 #include <gsl/gsl_sort.h>
 #include <gsl/gsl_cblas.h>
 #include "misc.h"
+#include <gsl/gsl_errno.h>
 
 char quad_err[256];
 
+void my_handler(const char *reason, const char *file, int line, int gsl_errno)
+{
+    fprintf(stderr, "gsl: %s:%d: %s: %s\n", file, line, "ERROR", reason);
+    return;
+}
+
 /* Returns an approximate vertex position and time using the QUAD fitter.
  *
  * The QUAD fitter was originally developed for SNO and estimates the event
@@ -83,6 +90,8 @@ int quad(event *ev, double *pos, double *t0, size_t npoints)
     static size_t p2[MAX_PMTS];
     double tmin;
     size_t max_tries;
+    gsl_error_handler_t *old_handler;
+    int status;
 
     wavelength0 = 400.0;
     n_d2o = get_index_snoman_d2o(wavelength0);
@@ -156,7 +165,22 @@ int quad(event *ev, double *pos, double *t0, size_t npoints)
         gsl_vector_view n_view = gsl_vector_view_array(n,3);
 
         gsl_linalg_LU_decomp(&m.matrix, p, &s);
-        gsl_linalg_LU_invert(&m.matrix, p, &minv.matrix);
+
+        /* Ocassionaly the matrix is singular and we can't invert it. Since we
+         * don't want our program to quit when this happens, we install a new
+         * gsl error handler, do the matrix inversion, and then restore the old
+         * handler. If the matrix inversion failed, then we just continue with
+         * the next quad point. */
+
+        /* save original handler, install new handler */
+        old_handler = gsl_set_error_handler(&my_handler);
+
+        status = gsl_linalg_LU_invert(&m.matrix, p, &minv.matrix);
+
+        /* restore original handler */
+        gsl_set_error_handler(old_handler);
+
+        if (status) continue;
 
         gsl_blas_dgemv(CblasNoTrans,1.0,&minv.matrix,&k_view.vector,0.0,&g_view.vector);
         gsl_blas_dgemv(CblasNoTrans,1.0,&minv.matrix,&n_view.vector,0.0,&h_view.vector);