skip PMTs which weren't hit for the fast likelihood calculation

2 files changed, 23 insertions, 3 deletions

diff --git a/src/fit.c b/src/fit.c
index 4d8d389..3ab0ffe 100644
--- a/src/fit.c
+++ b/src/fit.c
@@ -5253,8 +5253,6 @@ int fit_event(event *ev, double *xopt, double *fmin)
         }
     }
 
-    memcpy(x,xopt,sizeof(x));
-
     /* Reset the lower and upper bounds. */
     lb[0] = Tmin;
     lb[1] = -1000.0;
@@ -5286,6 +5284,12 @@ int fit_event(event *ev, double *xopt, double *fmin)
     nlopt_set_ftol_abs(opt, 1e-5);
     nlopt_set_maxeval(opt, 1000);
 
+    memcpy(x,xopt,sizeof(x));
+
+    nlopt_optimize(opt,x,&fval);
+
+    if (stop) goto end;
+
     do {
         *fmin = fval;
         memcpy(xopt,x,sizeof(x));
diff --git a/src/likelihood.c b/src/likelihood.c
index ccdfc84..5f66a72 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -392,6 +392,10 @@ double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, doubl
     for (i = 0; i < MAX_PMTS; i++) {
         if (pmts[i].pmt_type != PMT_NORMAL) continue;
 
+        /* Skip PMTs which weren't hit when doing the "fast" likelihood
+         * calculation. */
+        if (fast && !ev->pmt_hits[i].hit) continue;
+
         params.i = i;
 
         if (fast) {
@@ -501,13 +505,25 @@ double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, doubl
 
     for (i = 0; i < MAX_PMTS; i++) {
         if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL) continue;
-        mu[i] = mu_direct[i] + mu_indirect + mu_noise;
+
+        /* Skip PMTs which weren't hit when doing the "fast" likelihood
+         * calculation. */
+        if (fast && !ev->pmt_hits[i].hit) continue;
+
+        if (fast)
+            mu[i] = mu_direct[i] + mu_noise;
+        else
+            mu[i] = mu_direct[i] + mu_indirect + mu_noise;
     }
 
     nhit = 0;
     for (i = 0; i < MAX_PMTS; i++) {
         if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL) continue;
 
+        /* Skip PMTs which weren't hit when doing the "fast" likelihood
+         * calculation. */
+        if (fast && !ev->pmt_hits[i].hit) continue;
+
         log_mu = log(mu[i]);
 
         if (ev->pmt_hits[i].hit) {