add a goodness of fit parameter psi to the fit

3 files changed, 140 insertions, 10 deletions

diff --git a/src/fit.c b/src/fit.c
index e4fb97f..8a68105 100644
--- a/src/fit.c
+++ b/src/fit.c
@@ -43,6 +43,7 @@ typedef struct fitParams {
     int fast;
     int print;
     int id;
+    int charge_only;
 } fitParams;
 
 int particles[] = {
@@ -4975,7 +4976,7 @@ static struct startingParameters {
     { 800.0, 800.0, 800.0},
 };
 
-static double nopt_nll(unsigned int n, const double *x, double *grad, void *params)
+static double nlopt_nll(unsigned int n, const double *x, double *grad, void *params)
 {
     fitParams *fpars = (fitParams *) params;
     double theta, phi;
@@ -5007,7 +5008,7 @@ static double nopt_nll(unsigned int n, const double *x, double *grad, void *para
     v.n = 1;
 
     gettimeofday(&tv_start, NULL);
-    fval = nll(fpars->ev, &v, 1, fpars->dx, fpars->dx_shower, fpars->fast);
+    fval = nll(fpars->ev, &v, 1, fpars->dx, fpars->dx_shower, fpars->fast, fpars->charge_only);
     gettimeofday(&tv_stop, NULL);
 
     long long elapsed = (tv_stop.tv_sec - tv_start.tv_sec)*1000 + (tv_stop.tv_usec - tv_start.tv_usec)/1000;
@@ -5118,7 +5119,7 @@ int fit_event(event *ev, double *xopt, double *fmin, int id)
     struct timeval tv_start, tv_stop;
 
     opt = nlopt_create(NLOPT_LN_BOBYQA, 9);
-    nlopt_set_min_objective(opt,nopt_nll,&fpars);
+    nlopt_set_min_objective(opt,nlopt_nll,&fpars);
 
     /* Make a guess as to the energy. Right now we just use a simple
      * approximation that the muon produces approximately 6 photons/MeV.
@@ -5219,6 +5220,7 @@ int fit_event(event *ev, double *xopt, double *fmin, int id)
     fpars.dx_shower = 10.0;
     fpars.fast = 1;
     fpars.print = 0;
+    fpars.charge_only = 0;
     nlopt_set_ftol_abs(opt, 1.0);
     nlopt_set_maxeval(opt, 20);
 
@@ -5313,6 +5315,7 @@ int fit_event(event *ev, double *xopt, double *fmin, int id)
     fpars.dx_shower = 10.0;
     fpars.fast = 0;
     fpars.print = 1;
+    fpars.charge_only = 0;
     nlopt_set_ftol_abs(opt, 1e-5);
     nlopt_set_maxeval(opt, 1000);
 
@@ -5408,6 +5411,25 @@ end:
     return rv;
 }
 
+size_t get_nhit(event *ev)
+{
+    /* Returns the number of PMT hits in event `ev`.
+     *
+     * Note: Only hits on normal PMTs which aren't flagged are counted. */
+    size_t i, nhit;
+
+    nhit = 0;
+    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) continue;
+
+        nhit++;
+    }
+
+    return nhit;
+}
+
 int main(int argc, char **argv)
 {
     int i;
@@ -5425,6 +5447,8 @@ int main(int argc, char **argv)
     FILE *fout = NULL;
     int skip_second_event = 0;
     struct timeval tv_start, tv_stop;
+    long long elapsed;
+    double fmin_best;
 
     for (i = 1; i < argc; i++) {
         if (strlen(argv[i]) >= 2 && !strncmp(argv[i], "--", 2)) {
@@ -5581,6 +5605,7 @@ int main(int argc, char **argv)
             if (fout) {
                 if (first_ev) fprintf(fout, "    ev:\n");
                 fprintf(fout, "      - gtid: %i\n", ev.gtid);
+                fprintf(fout, "        nhit: %zu\n", get_nhit(&ev));
                 fprintf(fout, "        fit:\n");
             }
 
@@ -5592,7 +5617,9 @@ int main(int argc, char **argv)
                 }
                 gettimeofday(&tv_stop, NULL);
 
-                long long elapsed = (tv_stop.tv_sec - tv_start.tv_sec)*1000 + (tv_stop.tv_usec - tv_start.tv_usec)/1000;
+                elapsed = (tv_stop.tv_sec - tv_start.tv_sec)*1000 + (tv_stop.tv_usec - tv_start.tv_usec)/1000;
+
+                fmin_best = nll_best(&ev);
 
                 if (fout) {
                     fprintf(fout, "          %i:\n", particles[i]);
@@ -5607,6 +5634,7 @@ int main(int argc, char **argv)
                     fprintf(fout, "            z2: %.2f\n", xopt[8]);
                     fprintf(fout, "            fmin: %.2f\n", fmin);
                     fprintf(fout, "            time: %lld\n", elapsed);
+                    fprintf(fout, "            psi: %.2f\n", fmin-fmin_best);
                     fflush(fout);
                 }
             }
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];
 
diff --git a/src/likelihood.h b/src/likelihood.h
index 69d6896..1ef2736 100644
--- a/src/likelihood.h
+++ b/src/likelihood.h
@@ -83,6 +83,7 @@ double particle_get_energy(double x, particle *p);
 void particle_free(particle *p);
 double time_pdf(double t, double mu_noise, double mu_indirect, double *mu_direct, double *mu_shower, size_t n, double *ts, double *ts_shower, double tmean, double sigma, double *ts_sigma);
 double time_cdf(double t, double mu_noise, double mu_indirect, double *mu_direct, double *mu_shower, size_t n, double *ts, double *ts_shower, double tmean, double sigma, double *ts_sigma);
-double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast);
+double nll_best(event *ev);
+double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast, int charge_only);
 
 #endif