diff options
Diffstat (limited to 'src/likelihood.c')
| -rw-r--r-- | src/likelihood.c | 14 |
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; |