diff options
Diffstat (limited to 'src/quad.c')
| -rw-r--r-- | src/quad.c | 26 |
1 files changed, 25 insertions, 1 deletions
@@ -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); |