add a goodness of fit parameter psi to the fit

1 files changed, 106 insertions, 5 deletions

diff --git a/src/likelihood.c b/src/likelihood.c
index 0e78800..64f28d4 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -720,7 +720,104 @@ static double getKineticEnergy(double x, void *p)
     return particle_get_energy(x, (particle *) p);
 }
 
-double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast)
+double nll_best(event *ev)
+{
+    /* Returns the negative log likelihood of the "best hypothesis" for the event `ev`.
+     *
+     * By "best hypothesis" I mean we restrict the model to assign a single
+     * mean number of PE for each PMT in the event assuming that the number of
+     * PE hitting each PMT is poisson distributed. In addition, the model picks
+     * a single time for the light to arrive with the PDF being given by a sum
+     * of dark noise and a Gaussian around this time with a width equal to the
+     * single PE transit time spread.
+     *
+     * This calculation is intended to be used as a goodness of fit test for
+     * the fitter by computing a likelihood ratio between the best fit
+     * hypothesis and this ideal case. See Chapter 9 in Jayne's "Probability
+     * Theory: The Logic of Science" for more details. */
+    size_t i, j, nhit, maxj;
+    static double logp[MAX_PE], nll[MAX_PMTS], mu[MAX_PMTS], ts[MAX_PMTS], ts_shower, ts_sigma, mu_shower;
+    double log_mu, max_logp, mu_noise, mu_indirect_total, min_ratio;
+
+    mu_noise = DARK_RATE*GTVALID*1e-9;
+
+    /* Compute the "best" number of expected PE for each PMT. */
+    for (i = 0; i < MAX_PMTS; i++) {
+        if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL) continue;
+
+        if (!ev->pmt_hits[i].hit) {
+            /* If the PMT wasn't hit we assume the expected number of PE just
+             * come from noise hits. */
+            mu[i] = mu_noise;
+            continue;
+        }
+
+        /* Technically we should be computing:
+         *
+         *     mu[i] = max(p(q|mu))
+         *
+         * where by max I mean we find the expected number of PE which
+         * maximizes p(q|mu). However, I don't know of any easy analytical
+         * solution to this. So instead we just compute the number of PE which
+         * maximizes p(q|n). */
+        for (j = 1; j < MAX_PE; j++) {
+            logp[j] = get_log_pq(ev->pmt_hits[i].qhs,j);
+
+            if (j == 1 || logp[j] > max_logp) {
+                maxj = j;
+                max_logp = logp[j];
+            }
+        }
+
+        mu[i] = maxj;
+        ts[i] = ev->pmt_hits[i].t;
+    }
+
+    mu_shower = 0;
+    ts_shower = 0;
+    ts_sigma = PMT_TTS;
+
+    mu_indirect_total = 0.0;
+
+    min_ratio = MIN_RATIO;
+
+    /* Now, we actually compute the negative log likelihood for the best
+     * hypothesis.
+     *
+     * Currently this calculation is identical to the one in nll() so it should
+     * probably be made into a function. */
+    nhit = 0;
+    for (i = 0; i < MAX_PMTS; i++) {
+        if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL) continue;
+
+        log_mu = log(mu[i]);
+
+        if (ev->pmt_hits[i].hit) {
+            for (j = 1; j < MAX_PE; j++) {
+                logp[j] = get_log_pq(ev->pmt_hits[i].qhs,j) - mu[i] + j*log_mu - lnfact(j) + log_pt(ev->pmt_hits[i].t, j, mu_noise, mu_indirect_total, &mu[i], &mu_shower, 1, &ts[i], &ts_shower, ts[i], PMT_TTS, &ts_sigma);
+
+                if (j == 1 || logp[j] > max_logp) max_logp = logp[j];
+
+                if (logp[j] - max_logp < min_ratio*ln(10)) {
+                    j++;
+                    break;
+                }
+            }
+
+            nll[nhit++] = -logsumexp(logp+1, j-1);
+        } else {
+            logp[0] = -mu[i];
+            for (j = 1; j < MAX_PE_NO_HIT; j++) {
+                logp[j] = get_log_pmiss(j) - mu[i] + j*log_mu - lnfact(j);
+            }
+            nll[nhit++] = -logsumexp(logp, MAX_PE_NO_HIT);
+        }
+    }
+
+    return kahan_sum(nll,nhit);
+}
+
+double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast, int charge_only)
 {
     /* Returns the negative log likelihood for event `ev` given a particle with
      * id `id`, initial kinetic energy `T0`, position `pos`, direction `dir` and
@@ -959,10 +1056,14 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast
 
         if (ev->pmt_hits[i].hit) {
             for (j = 1; j < MAX_PE; j++) {
-                if (fast)
-                    logp[j] = get_log_pq(ev->pmt_hits[i].qhs,j) - mu[i] + j*log_mu - lnfact(j) + log_pt_fast;
-                else
-                    logp[j] = get_log_pq(ev->pmt_hits[i].qhs,j) - mu[i] + j*log_mu - lnfact(j) + log_pt(ev->pmt_hits[i].t, j, mu_noise, mu_indirect_total, &mu_direct[i][0], &mu_shower[i][0], n, &ts[i][0], &ts_shower[i][0], ts[i][0], PMT_TTS, &ts_sigma[i][0]);
+                logp[j] = get_log_pq(ev->pmt_hits[i].qhs,j) - mu[i] + j*log_mu - lnfact(j);
+
+                if (!charge_only) {
+                    if (fast)
+                        logp[j] += log_pt_fast;
+                    else
+                        logp[j] += log_pt(ev->pmt_hits[i].t, j, mu_noise, mu_indirect_total, &mu_direct[i][0], &mu_shower[i][0], n, &ts[i][0], &ts_shower[i][0], ts[i][0], PMT_TTS, &ts_sigma[i][0]);
+                }
 
                 if (j == 1 || logp[j] > max_logp) max_logp = logp[j];