update fit.c to fit multiple vertices

This commit adds a new function fit_event2() to fit multiple vertices. To seed
the fit, fit_event2() does the following:

  - use the QUAD fitter to find the position and initial time of the event
  - call find_peaks() to find possible directions for the particles
  - loop over all possible unique combinations of the particles and direction
    vectors and do a "fast" minimization

The best minimum found from the "fast" minimizations is then used to start the fit.

This commit has a few other updates:

- adds a hit_only parameter to the nll() function. This was necessary since
previously PMTs which weren't hit were always skipped for the fast
minimization, but when fitting for multiple vertices we need to include PMTs
which aren't hit since we float the energy.

- add the function guess_energy() to guess the energy of a particle given a
position and direction. This function estimates the energy by summing up the
QHS for all PMTs hit within the Cerenkov cone and dividing by 6.

- fixed a bug which caused the fit to freeze when hitting ctrl-c during the
fast minimization phase.

1 files changed, 5 insertions, 5 deletions

diff --git a/src/likelihood.c b/src/likelihood.c
index da0614d..786b221 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -387,7 +387,7 @@ static void integrate_path_shower(particle *p, double *x, double *pdf, double T0
 
     *sigma = sqrt(t2-(*time)*(*time));
 
-    if (*mu_direct == 0.0) {
+    if (*mu_direct == 0.0 || *sigma != *sigma) {
         *time = 0.0;
         *sigma = PMT_TTS;
     } else {
@@ -818,7 +818,7 @@ double nll_best(event *ev)
     return kahan_sum(nll,nhit);
 }
 
-double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast, int charge_only)
+double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast, int charge_only, int hit_only)
 {
     /* Returns the negative log likelihood for event `ev` given a particle with
      * id `id`, initial kinetic energy `T0`, position `pos`, direction `dir` and
@@ -953,7 +953,7 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast
 
             /* Skip PMTs which weren't hit when doing the "fast" likelihood
              * calculation. */
-            if (fast && !ev->pmt_hits[i].hit) continue;
+            if (hit_only && !ev->pmt_hits[i].hit) continue;
 
             if (fast) {
                 mu_direct[i][j] = get_total_charge_approx(beta0, v[j].pos, v[j].dir, p, i, smax, theta0, &ts[i][j], &mu_reflected, n_d2o, n_h2o, cos_theta_cerenkov, sin_theta_cerenkov);
@@ -1026,7 +1026,7 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast
 
         /* Skip PMTs which weren't hit when doing the "fast" likelihood
          * calculation. */
-        if (fast && !ev->pmt_hits[i].hit) continue;
+        if (hit_only && !ev->pmt_hits[i].hit) continue;
 
         for (j = 0; j < n; j++) {
             if (fast)
@@ -1048,7 +1048,7 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast
 
         /* Skip PMTs which weren't hit when doing the "fast" likelihood
          * calculation. */
-        if (fast && !ev->pmt_hits[i].hit) continue;
+        if (hit_only && !ev->pmt_hits[i].hit) continue;
 
         log_mu = log(mu[i]);