update sno_charge.c

This commit adds lots of comments to sno_charge.c and makes a couple of other
changes:

- use interp1d() instead of the GSL interpolation routines
- increase MAX_PE to 100

I increased MAX_PE because I determined that it had a rather large impact on
the likelihood function for 500 MeV electrons. This unfortunately slows down
the initialization by a lot. I think I could speed this up by convolving the
single PE charge distribution with a gaussian *before* convolving the charge
distributions to compute the charge distributions for multiple PE.

6 files changed, 81 insertions, 68 deletions

diff --git a/src/fit.c b/src/fit.c
index 23eab74..0634116 100644
--- a/src/fit.c
+++ b/src/fit.c
@@ -5619,7 +5619,6 @@ int main(int argc, char **argv)
 end:
 
     free_interpolation();
-    free_charge();
     pmt_response_free();
 
     db_free(db);
@@ -5632,7 +5631,6 @@ end:
 
 err:
     free_interpolation();
-    free_charge();
     pmt_response_free();
 
     db_free(db);
diff --git a/src/likelihood.c b/src/likelihood.c
index 40261f9..93740a0 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -944,9 +944,9 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t
         if (ev->pmt_hits[i].hit) {
             for (j = 1; j < MAX_PE; j++) {
                 if (fast)
-                    logp[j] = log_pq(ev->pmt_hits[i].qhs,j) - mu[i] + j*log_mu - lnfact(j) + log_pt_fast;
+                    logp[j] = get_log_pq(ev->pmt_hits[i].qhs,j) - mu[i] + j*log_mu - lnfact(j) + log_pt_fast;
                 else
-                    logp[j] = log_pq(ev->pmt_hits[i].qhs,j) - mu[i] + j*log_mu - lnfact(j) + log_pt(ev->pmt_hits[i].t, j, mu_noise, mu_indirect, &mu_direct[i], &mu_shower[i], 1, &ts[i], &ts_shower[i], ts[i], PMT_TTS, &ts_sigma[i]);
+                    logp[j] = get_log_pq(ev->pmt_hits[i].qhs,j) - mu[i] + j*log_mu - lnfact(j) + log_pt(ev->pmt_hits[i].t, j, mu_noise, mu_indirect, &mu_direct[i], &mu_shower[i], 1, &ts[i], &ts_shower[i], ts[i], PMT_TTS, &ts_sigma[i]);
 
                 if (j == 1 || logp[j] > max_logp) max_logp = logp[j];
 
diff --git a/src/sno_charge.c b/src/sno_charge.c
index a019134..de0f2b3 100644
--- a/src/sno_charge.c
+++ b/src/sno_charge.c
@@ -7,20 +7,40 @@
 #include <stdio.h>
 #include "misc.h"
 
-#define MAX_PE 10
+/* Maximum number of PE to compute charge PDFs for.
+ *
+ * For PE greater than this we compute the PDF using a Gaussian distribution. */
+#define MAX_PE 100
 
+/* Minimum charge to compute the PDF at (in units of PE). */
 double qlo = -1.0;
+/* Maximum charge to compute the PDF at (in units of PE). */
 double qhi = 100.0;
-size_t nq = 10000;
-static gsl_spline *splines[MAX_PE];
-static gsl_interp_accel *acc;
-static gsl_spline *log_splines[MAX_PE];
-static gsl_interp_accel *log_acc;
-
+/* Number of points in the charge PDF. */
+#define nq 10000
+
+/* Charge PDFs. */
+static double x[nq];
+static double pq[MAX_PE][nq];
+/* Log of the charge PDFs used for calculations in the likelihood. */
+static double log_pq[MAX_PE][nq];
+
+/* Mean and standard deviation of the single PE charge distribution.
+ *
+ * This is used to evaluate the probability of a high charge where we don't
+ * have a precomputed PDF. We assume the central limit theorem and just compute
+ * the probability density as a Gaussian with mean n*qmean and standard
+ * deviation sqrt(n)*qstd where `n` is the number of PE. */
 static double qmean, qstd;
 
-static double pmiss[MAX_PE];
+/* Lookup table for the log of the probability of the discriminator not firing
+ * as a function of the number of PE which hit the tube. We compute this in
+ * init_charge() by integrating the charge PDF up to the mean threshold for all
+ * channels. */
+static double log_pmiss[MAX_PE];
 
+/* Keep track of whether init_charge() has been called so we can exit if a
+ * function is called before being initialized. */
 static int initialized = 0;
 
 /* Parameters for the single PE charge distribution. These numbers come from
@@ -86,8 +106,12 @@ double spe_pol2exp(double q)
     return (gmratio*(NG1*funpol1 + NG2*funpol2) + NEXP*exp(-q/Q0))/norm;
 }
 
-double log_pq(double q, int n)
+double get_log_pq(double q, int n)
 {
+    /* Returns the log of the probability of observing a charge `q` given that
+     * `n` PE were generated.
+     *
+     * `q` should be in units of PE. */
     if (!initialized) {
         fprintf(stderr, "charge interpolation hasn't been initialized!\n");
         exit(1);
@@ -100,11 +124,15 @@ double log_pq(double q, int n)
 
     if (q < qlo || q > qhi) return -INFINITY;
 
-    return interp1d(q,log_splines[n-1]->x,log_splines[n-1]->y,log_splines[n-1]->size);
+    return interp1d(q,x,log_pq[n-1],nq);
 }
 
-double pq(double q, int n)
+double get_pq(double q, int n)
 {
+    /* Returns the probability of observing a charge `q` given that `n` PE were
+     * generated.
+     *
+     * `q` should be in units of PE. */
     if (!initialized) {
         fprintf(stderr, "charge interpolation hasn't been initialized!\n");
         exit(1);
@@ -117,23 +145,29 @@ double pq(double q, int n)
 
     if (q < qlo || q > qhi) return 0.0;
 
-    return gsl_spline_eval(splines[n-1], q, acc);
+    return interp1d(q,x,pq[n-1],nq);
 }
 
 double get_log_pmiss(int n)
 {
+    /* Returns the log of the probability that the channel discriminator didn't
+     * fire 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, so the log is zero. */
         return 0.0;
     } else if (n > MAX_PE) {
+        /* For n > MAX_PE we assume there is a 100% chance that the
+         * discriminator fires. */
         return -INFINITY;
     }
 
-    return pmiss[n-1];
+    return log_pmiss[n-1];
 }
 
 static double gsl_charge(double x, void *params)
@@ -141,14 +175,14 @@ static double gsl_charge(double x, void *params)
     double q = ((double *) params)[0];
     int n = (int) ((double *) params)[1];
 
-    return pq(x,1)*pq(q-x,n-1);
+    return get_pq(x,1)*get_pq(q-x,n-1);
 }
 
 static double gsl_charge2(double x, void *params)
 {
     int n = (int) ((double *) params)[0];
 
-    return pq(x,n);
+    return get_pq(x,n);
 }
 
 static double gsl_smear(double x, void *params)
@@ -156,7 +190,7 @@ static double gsl_smear(double x, void *params)
     double q = ((double *) params)[0];
     int n = (int) ((double *) params)[1];
 
-    return pq(x,n)*norm(q-x,0.0,QSMEAR_ADC/MEAN_HIPT);
+    return get_pq(x,n)*norm(q-x,0.0,QSMEAR_ADC/MEAN_HIPT);
 }
 
 static double sno_charge1(double q, void *params)
@@ -171,20 +205,18 @@ static double sno_charge2(double q, void *params)
 
 void init_charge(void)
 {
+    /* Compute the charge PDFs for up to MAX_PE PE. */
     int i, j;
     gsl_integration_cquad_workspace *w;
     double result, error;
     size_t nevals;
-    double *x, *y;
     double params[2];
     gsl_function F;
     double q, q2;
+    double y[nq];
 
     initialized = 1;
 
-    x = malloc(nq*sizeof(double));
-    y = malloc(nq*sizeof(double));
-
     for (i = 0; i < nq; i++) {
         x[i] = qlo + (qhi-qlo)*i/(nq-1);
     }
@@ -194,42 +226,42 @@ void init_charge(void)
     F.function = &sno_charge1;
     F.params = NULL;
 
+    /* Compute the average single PE charge. */
     gsl_integration_cquad(&F, 0, qhi, 0, 1e-9, w, &result, &error, &nevals);
     q = result;
 
     F.function = &sno_charge2;
     F.params = NULL;
 
+    /* Compute the expected value of the single PE charge squared. */
     gsl_integration_cquad(&F, 0, qhi, 0, 1e-9, w, &result, &error, &nevals);
     q2 = result;
 
+    /* Compute the mean and standard deviation of the single PE charge PDF. */
     qmean = q;
     qstd = sqrt(q2 - q*q);
 
-    F.function = &gsl_charge;
-    F.params = &params;
-
-    acc = gsl_interp_accel_alloc();
-
-    splines[0] = gsl_spline_alloc(gsl_interp_linear,nq);
+    /* Compute the single PE charge PDF. */
     for (j = 0; j < nq; j++) {
-        y[j] = spe_pol2exp(x[j]);
+        pq[0][j] = spe_pol2exp(x[j]);
     }
-    gsl_spline_init(splines[0],x,y,nq);
 
+    F.function = &gsl_charge;
+    F.params = &params;
+
+    /* Compute the charge PDFs for `n` PE by convolving the `n-1` charge PDF
+     * with the single PE charge PDF. */
     for (i = 2; i <= MAX_PE; i++) {
-        splines[i-1] = gsl_spline_alloc(gsl_interp_linear,nq);
         for (j = 0; j < nq; j++) {
             params[0] = x[j];
             params[1] = i;
             if (x[j] < 0) {
-                y[j] = 0.0;
+                pq[i-1][j] = 0.0;
                 continue;
             }
             gsl_integration_cquad(&F, 0, x[j], 0, 1e-4, w, &result, &error, &nevals);
-            y[j] = result;
+            pq[i-1][j] = result;
         }
-        gsl_spline_init(splines[i-1],x,y,nq);
     }
 
     /* Integrate the charge distribution before smearing to figure out the
@@ -241,46 +273,34 @@ void init_charge(void)
     for (i = 1; i <= MAX_PE; i++) {
         params[0] = i;
         gsl_integration_cquad(&F, 0, MEAN_THRESH/MEAN_HIPT, 0, 1e-9, w, &result, &error, &nevals);
-        pmiss[i-1] = log(result);
+        log_pmiss[i-1] = log(result);
     }
 
     F.function = &gsl_smear;
 
+    /* Convolve the charge PDFs with a gaussian distribution to represent
+     * smearing from the electronics. */
     for (i = 1; i <= MAX_PE; i++) {
         for (j = 0; j < nq; j++) {
             params[0] = x[j];
             params[1] = i;
             gsl_integration_cquad(&F, fmax(0.0,x[j]-QSMEAR_ADC*10/MEAN_HIPT), x[j]+QSMEAR_ADC*10/MEAN_HIPT, 0, 1e-4, w, &result, &error, &nevals);
+            /* Store the smeared charge PDF in a temporary array since we need
+             * to keep the original PDF for the next iteration of the loop. */
             y[j] = result;
         }
-        gsl_spline_free(splines[i-1]);
-        splines[i-1] = gsl_spline_alloc(gsl_interp_linear,nq);
-        gsl_spline_init(splines[i-1],x,y,nq);
+        for (j = 0; j < nq; j++) {
+            pq[i-1][j] = y[j];
+        }
     }
 
     gsl_integration_cquad_workspace_free(w);
 
+    /* Compute the log of the charge distributions to speed up the likelihood
+     * calculation. */
     for (i = 1; i <= MAX_PE; i++) {
         for (j = 0; j < nq; j++) {
-            y[j] = log(splines[i-1]->y[j]);
+            log_pq[i-1][j] = log(pq[i-1][j]);
         }
-        log_splines[i-1] = gsl_spline_alloc(gsl_interp_linear,nq);
-        gsl_spline_init(log_splines[i-1],x,y,nq);
     }
-
-    free(x);
-    free(y);
-}
-
-void free_charge(void)
-{
-    size_t i;
-
-    for (i = 0; i < MAX_PE; i++) {
-        gsl_spline_free(splines[i]);
-    }
-
-    gsl_interp_accel_free(acc);
-
-    initialized = 0;
 }
diff --git a/src/sno_charge.h b/src/sno_charge.h
index 50b018e..6f513f9 100644
--- a/src/sno_charge.h
+++ b/src/sno_charge.h
@@ -2,10 +2,9 @@
 #define SNO_CHARGE_H
 
 void init_charge(void);
-double log_pq(double q, int n);
-double pq(double q, int n);
+double get_log_pq(double q, int n);
+double get_pq(double q, int n);
 double get_log_pmiss(int n);
 double spe_pol2exp(double q);
-void free_charge(void);
 
 #endif
diff --git a/src/test-charge.c b/src/test-charge.c
index 2fe199b..2373beb 100644
--- a/src/test-charge.c
+++ b/src/test-charge.c
@@ -70,7 +70,7 @@ int main(int argc, char **argv)
 
     for (i = 1; i <= n; i++) {
         for (j = 0; j < nq; j++) {
-            fprintf(pipe, "%.10f %.10f\n", x[j], pq(x[j],i));
+            fprintf(pipe, "%.10f %.10f\n", x[j], get_pq(x[j],i));
         }
         fprintf(pipe, "\n\n");
     }
@@ -82,7 +82,5 @@ int main(int argc, char **argv)
 
     free(x);
 
-    free_charge();
-
     return 0;
 }
diff --git a/src/test.c b/src/test.c
index a3f98a9..8bd0688 100644
--- a/src/test.c
+++ b/src/test.c
@@ -404,7 +404,7 @@ int test_solid_angle(char *err)
 static double sno_charge(double q, void *params)
 {
     int n = ((int *) params)[0];
-    return pq(q,n);
+    return get_pq(q,n);
 }
 
 int test_sno_charge_integral(char *err)
@@ -431,8 +431,6 @@ int test_sno_charge_integral(char *err)
         goto err;
     }
 
-    free_charge();
-
     gsl_integration_cquad_workspace_free(w);
 
     return 0;