diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/likelihood.c | 23 | ||||
| -rw-r--r-- | src/likelihood.h | 4 |
2 files changed, 24 insertions, 3 deletions
diff --git a/src/likelihood.c b/src/likelihood.c index f77cab6..6f0e5d9 100644 --- a/src/likelihood.c +++ b/src/likelihood.c @@ -1285,7 +1285,7 @@ double nll_best(event *ev) double nll(event *ev, vertex *v, size_t n, double dx, int ns, const int fast, int charge_only, int hit_only) { size_t i, j, k, nhit; - static double logp[MAX_PE], nll[MAX_PMTS]; + static double logp[MAX_PE], nll[MAX_PMTS], tmp[2]; double range, theta0, E0, p0, beta0, smax, log_mu, max_logp; particle *p; double cos_theta_cerenkov, sin_theta_cerenkov; @@ -1552,7 +1552,15 @@ double nll(event *ev, vertex *v, size_t n, double dx, int ns, const int fast, in } } - nll[nhit++] = -logsumexp(logp+1, j-1); + /* Now, we include the probability that the channel is miscalibrated as: + * + * P(q,t) = P(q,t|calibrated)P(calibrated) + P(q,t|miscalibrated)P(miscalibrated) + * + * For P(q,t|miscalibrated) we just assume a uniform distribution + * over all 4096 possible values. */ + tmp[0] = logsumexp(logp+1, j-1) + log(1-P_MISCALIBRATION); + tmp[1] = log(1/4096.0) + log(1/4096.0) + log(P_MISCALIBRATION); + nll[nhit++] = -logsumexp(tmp, 2); } else { logp[0] = -mu_sum[i]; if (fast) { @@ -1562,7 +1570,16 @@ double nll(event *ev, vertex *v, size_t n, double dx, int ns, const int fast, in for (j = 1; j < MAX_PE_NO_HIT; j++) { logp[j] = get_log_pmiss(j) - mu_sum[i] + j*log_mu - lnfact(j); } - nll[nhit++] = -logsumexp(logp, MAX_PE_NO_HIT); + + /* Now, we include the probability that the channel is miscalibrated as: + * + * P(not hit) = P(not hit|calibrated)P(calibrated) + P(not hit|miscalibrated)P(miscalibrated) + * + * where by miscalibrated here we really mean a PMT that wasn't + * read out or included in the event for whatever reason. */ + tmp[0] = logsumexp(logp, MAX_PE_NO_HIT) + log(1-P_MISCALIBRATION); + tmp[1] = log(P_MISCALIBRATION); + nll[nhit++] = -logsumexp(tmp, 2); } } diff --git a/src/likelihood.h b/src/likelihood.h index 304fd70..ddd8397 100644 --- a/src/likelihood.h +++ b/src/likelihood.h @@ -22,6 +22,10 @@ #include <gsl/gsl_interp.h> #include <gsl/gsl_spline.h> +/* Probability that a channel is miscalibrated and/or not working even though + * it's marked as OK in DQXX, etc. */ +#define P_MISCALIBRATION 1e-4 + #define EPSILON 1e-10 #define PSUP_REFLECTION_TIME 80.0 |