update nll_best() so that we don't get negative psi values

This commit updates the time estimate in nll_best() to take into account the
fact that for PMTs which get hit by a lot of photons the time distribution we
expect is the first order statistic of the overall photon hit time
distribution.

1 files changed, 12 insertions, 2 deletions

diff --git a/src/likelihood.c b/src/likelihood.c
index dfe580e..b02aeef 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -1084,7 +1084,7 @@ double nll_best(event *ev)
 
         if (!ev->pmt_hits[i].hit) {
             /* If the PMT wasn't hit we assume the expected number of PE just
-             * come from noise hits. */
+             * comes from noise hits. */
             mu[i] = mu_noise;
             continue;
         }
@@ -1095,7 +1095,17 @@ double nll_best(event *ev)
          *
          */
         mu[i] = get_most_likely_mean_pe(ev->pmt_hits[i].qhs);
-        ts[i] = ev->pmt_hits[i].t;
+        /* We want to estimate the mean time which is most likely to produce a
+         * first order statistic of the actual hit time given there were
+         * approximately mu[i] PE. As far as I know there are no closed form
+         * solutions for the first order statistic of a Gaussian, but there is
+         * an approximate form in the paper "Expected Normal Order Statistics
+         * (Exact and Approximate)" by J. P. Royston which is what I use here.
+         *
+         * I found this via the stack exchange question here:
+         * https://stats.stackexchange.com/questions/9001/approximate-order-statistics-for-normal-random-variables. */
+        double alpha = 0.375;
+        ts[i] = ev->pmt_hits[i].t - gsl_cdf_ugaussian_Pinv((1-alpha)/(mu[i]-2*alpha+1))*PMT_TTS;
     }
 
     ts_sigma = PMT_TTS;