small update to the charge likelihood calculation

This commit updates the charge likelihood calculation to calculate:

    P(hit,q|n) = P(q|hit,n)*P(hit|n)

This has almost no effect on the fit results, but is technically correct.

1 files changed, 35 insertions, 0 deletions

diff --git a/src/sno_charge.c b/src/sno_charge.c
index f78da38..4b51954 100644
--- a/src/sno_charge.c
+++ b/src/sno_charge.c
@@ -54,6 +54,7 @@ static double qmean, qstd;
  * init_charge() by integrating the charge PDF up to the mean threshold for all
  * channels. */
 static double log_pmiss[MAX_PE];
+static double log_phit[MAX_PE];
 
 /* Keep track of whether init_charge() has been called so we can exit if a
  * function is called before being initialized. */
@@ -164,6 +165,28 @@ double get_pq(double q, int n)
     return interp1d(q,x,pq[n-1],nq);
 }
 
+double get_log_phit(int n)
+{
+    /* Returns the log of the probability that the channel discriminator fired
+     * given `n` PE were generated. */
+    if (!initialized) {
+        fprintf(stderr, "charge interpolation hasn't been initialized!\n");
+        exit(1);
+    }
+
+    if (n == 0) {
+        /* If no PE were generated, there is a 100% chance that the
+         * discriminator didn't fire. */
+        return -INFINITY;
+    } else if (n > MAX_PE) {
+        /* For n > MAX_PE we assume there is a 100% chance that the
+         * discriminator fires. */
+        return 0.0;
+    }
+
+    return log_phit[n-1];
+}
+
 double get_log_pmiss(int n)
 {
     /* Returns the log of the probability that the channel discriminator didn't
@@ -290,6 +313,18 @@ void init_charge(void)
         params[0] = i;
         gsl_integration_cquad(&F, 0, MEAN_THRESH/MEAN_HIPT, 0, 1e-9, w, &result, &error, &nevals);
         log_pmiss[i-1] = log(result);
+        log_phit[i-1] = log(1.0-result);
+    }
+
+    /* Set the PDF for the charge to 0 before the discriminator level. */
+    for (i = 1; i <= MAX_PE; i++) {
+        for (j = 0; j < nq; j++) {
+            if (x[j] < MEAN_THRESH/MEAN_HIPT) pq[i-1][j] = 0.0;
+        }
+        double norm = trapz(&pq[i-1][0],x[1]-x[0],nq);
+        for (j = 0; j < nq; j++) {
+            pq[i-1][j] /= norm;
+        }
     }
 
     F.function = &gsl_smear;