initial commit of likelihood fit for muons

This commit contains code to fit for the energy, position, and direction of
muons in the SNO detector. Currently, we read events from SNOMAN zebra files
and fill an event struct containing the PMT hits and fit it with the Nelder
Mead simplex algorithm from GSL.

I've also added code to read in ZEBRA title bank files to read in the DQXX
files for a specific run. Any problems with channels in the DQCH and DQCR banks
are flagged in the event struct by masking in a bit in the flags variable and
these PMT hits are not included in the likelihood calculation.

The likelihood for an event is calculated by integrating along the particle
track for each PMT and computing the expected number of PE. The charge
likelihood is then calculated by looping over all possible number of PE and
computing:

    P(q|n)*P(n|mu)

where q is the calibrated QHS charge, n is the number of PE, and mu is the
expected number of photoelectrons. The latter is calculated assuming the
distribution of PE at a given PMT follows a Poisson distribution (which I think
should be correct given the track, but is probably not perfect for tracks which
scatter a lot).

The time part of the likelihood is calculated by integrating over the track for
each PMT and calculating the average time at which the PMT is hit. We then
assume the PDF for the photons to arrive is approximately a delta function and
compute the first order statistic for a given time to compute the probability
that the first photon arrived at a given time. So far I've only tested this
with single tracks but the method was designed to be easy to use when you are
fitting for multiple particles.

1 files changed, 246 insertions, 0 deletions

diff --git a/sno_charge.c b/sno_charge.c
new file mode 100644
index 0000000..faf7826
--- /dev/null
+++ b/sno_charge.c
@@ -0,0 +1,246 @@
+#include "sno_charge.h"
+#include <gsl/gsl_sf_gamma.h>
+#include <math.h>
+#include <gsl/gsl_spline.h>
+#include <gsl/gsl_integration.h>
+#include <unistd.h>
+#include <stdio.h>
+#include "misc.h"
+
+#define MAX_PE 10
+
+double qlo = -1.0;
+double qhi = 100.0;
+size_t nq = 10000;
+static gsl_spline *splines[MAX_PE];
+static gsl_interp_accel *acc;
+
+static double qmean, qstd;
+
+static double pmiss[MAX_PE];
+
+static int initialized = 0;
+
+/* Parameters for the single PE charge distribution. These numbers come from
+ * mc_generator.dat in SNOMAN 5.0294. */
+
+/* Polya parameter. */
+static double M = 8.1497;
+/* Amplitude of 1st Polya. */
+static double NG1 = 0.27235e-1;
+/* Slope of Polya Gain 1 vs. high-half point. */
+static double SLPG1 = 0.71726;
+/* Offset of Polya Gain 1 vs. high-half point. */
+static double OFFG1 = 1.2525;
+/* Amplitude of 2nd Polya. */
+static double NG2 = 0.7944e-3;
+/* Slope of Polya Gain 2 vs. high-half point. */
+static double SLPG2 = 0.71428;
+/* Offset of Polya Gain 2 vs. high-half point. */
+static double OFFG2 = 118.21;
+/* Amplitude of exponential noise peak. */
+static double NEXP = 0.81228e-1;
+/* Falloff parameter of exponential noise peak. */
+static double Q0 = 20.116;
+/* Mean of high-half point distribution. */
+static double MEAN_HIPT = 46.0;
+/* Width of noise before discriminator in ADC counts. */
+static double QNOISE = 0.61;
+/* Width of smearing after discriminator in ADC counts. */
+static double QSMEAR_ADC = 3.61;
+/* Mean of simulated threshold distribution. */
+static double MEAN_THRESH = 8.5;
+
+double spe_pol2exp(double q)
+{
+    /* Returns the probability distribution to get a charge `q` on a PMT with a
+     * high-half point `hhp` from a single photoelectron.
+     *
+     * We assume that the charge has been normalized by the high-half point.
+     *
+     * This function models the single PE charge distribution as a double Polya
+     * with an exponential noise term. */
+    double qpe1, qpe2, funpol1, funpol2, gmratio, g1, g2, norm, hhp;
+
+    if (q < 0) return 0.0;
+
+    hhp = MEAN_HIPT;
+
+    q *= MEAN_HIPT;
+
+    g1 = OFFG1 + SLPG1*hhp;
+    g2 = OFFG2 + SLPG2*hhp;
+
+    qpe1 = q/g1;
+    qpe2 = q/g2;
+
+    funpol1 = pow(M*qpe1,M-1)*exp(-M*qpe1);
+    funpol2 = pow(M*qpe2,M-1)*exp(-M*qpe2);
+
+    gmratio = M/gsl_sf_gamma(M);
+
+    norm = (NG1*g1 + NG2*g2 + NEXP*Q0)/MEAN_HIPT;
+
+    return (gmratio*(NG1*funpol1 + NG2*funpol2) + NEXP*exp(-q/Q0))/norm;
+}
+
+double pq(double q, int n)
+{
+    if (!initialized) {
+        fprintf(stderr, "charge interpolation hasn't been initialized!\n");
+        exit(1);
+    }
+
+    if (n > MAX_PE) {
+        /* Assume the distribution is gaussian by the central limit theorem. */
+        return norm(q,n*qmean,qstd);
+    }
+
+    if (q < qlo || q > qhi) return 0.0;
+
+    return gsl_spline_eval(splines[n-1], q, acc);
+}
+
+double get_pmiss(int n)
+{
+    if (!initialized) {
+        fprintf(stderr, "charge interpolation hasn't been initialized!\n");
+        exit(1);
+    }
+
+    if (n == 0) {
+        return 1.0;
+    } else if (n > MAX_PE) {
+        /* Assume the distribution is gaussian by the central limit theorem. */
+        return 0.0;
+    }
+
+    return pmiss[n-1];
+}
+
+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);
+}
+
+static double gsl_charge2(double x, void *params)
+{
+    int n = (int) ((double *) params)[0];
+
+    return pq(x,n);
+}
+
+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);
+}
+
+static double sno_charge1(double q, void *params)
+{
+    return q*spe_pol2exp(q);
+}
+
+static double sno_charge2(double q, void *params)
+{
+    return q*q*spe_pol2exp(q);
+}
+
+void init_charge(void)
+{
+    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;
+
+    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);
+    }
+
+    w = gsl_integration_cquad_workspace_alloc(100);
+
+    F.function = &sno_charge1;
+    F.params = NULL;
+
+    gsl_integration_cquad(&F, 0, qhi, 0, 1e-9, w, &result, &error, &nevals);
+    q = result;
+
+    F.function = &sno_charge2;
+    F.params = NULL;
+
+    gsl_integration_cquad(&F, 0, qhi, 0, 1e-9, w, &result, &error, &nevals);
+    q2 = result;
+
+    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);
+    for (j = 0; j < nq; j++) {
+        y[j] = spe_pol2exp(x[j]);
+    }
+    gsl_spline_init(splines[0],x,y,nq);
+
+    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;
+                continue;
+            }
+            gsl_integration_cquad(&F, 0, x[j], 0, 1e-4, w, &result, &error, &nevals);
+            y[j] = result;
+        }
+        gsl_spline_init(splines[i-1],x,y,nq);
+    }
+
+    /* Integrate the charge distribution before smearing to figure out the
+     * probability that the discriminator didn't fire.
+     *
+     * Note: Technically there is some smearing before the discriminator, but
+     * it is very small, so we ignore it. */
+    F.function = &gsl_charge2;
+    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] = result;
+    }
+
+    F.function = &gsl_smear;
+
+    for (i = 1; i <= MAX_PE; i++) {
+        for (j = 0; j < nq; j++) {
+            params[0] = x[j];
+            params[1] = i;
+            gsl_integration_cquad(&F, x[j]-QSMEAR_ADC*10/MEAN_HIPT, x[j]+QSMEAR_ADC*10/MEAN_HIPT, 0, 1e-4, w, &result, &error, &nevals);
+            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);
+    }
+
+    gsl_integration_cquad_workspace_free(w);
+    free(x);
+    free(y);
+}