add probability of a miscalibrated channel to the likelihood

This commit adds the probability that a channel is miscalibrated and/or doesn't
make it into the event to the likelihood. This was added because I noticed when
looking at the likelihood for one very high energy event that there was a
single PMT that should have been hit that wasn't in the event and which was not
marked as bad in DQXX.

I did some testing and the addition of this term does not seem to significantly
affect that atmospheric MC or the psi values for flashers. One unexpected
improvement is that it seems that external muons are more likely to correctly
reconstruct at the PSUP with this change. I haven't determined the exact cause
but I suspect it's because there is some mismodelling of the likelihood for
muons near the edge of the detector when they exit and that adding this term
allows the likelihood to ignore these PMT hits.

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