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, 256 insertions, 0 deletions

diff --git a/db.c b/db.c
new file mode 100644
index 0000000..9160cbf
--- /dev/null
+++ b/db.c
@@ -0,0 +1,256 @@
+#include "db.h"
+#include <stdio.h> /* for fopen(), etc. */
+#include <errno.h> /* for strerror(), etc. */
+#include <string.h> /* for strncmp(), etc. */
+#include <stdint.h> /* for uint32_t */
+#include <stdlib.h> /* for atoi() */
+#include "dict.h"
+
+char db_err[256];
+
+static uint64_t dbHash(const void *key)
+{
+    return dictGenHashFunction(key,8);
+}
+
+static int dbKeyCompare(void *privdata, const void *key1, const void *key2)
+{
+    return memcmp(key1,key2,8) == 0;
+}
+
+static void dbFree(void *privdata, void *val)
+{
+    free(val);
+}
+
+static dictType titleDB = {
+    dbHash,
+    NULL,
+    NULL,
+    dbKeyCompare,
+    dbFree,
+    dbFree
+};
+
+/* Add a bank to the database.
+ *
+ * Returns 0 on success, -1 on error. */
+int add_bank(dict *db, const char name[4], uint32_t id, dbval *data)
+{
+    uint32_t *buf = malloc(sizeof(uint32_t)*2);
+
+    memcpy(buf,name,4);
+    buf[1] = id;
+
+    if (dictAdd(db, buf, data) != DICT_OK) {
+        sprintf(db_err, "failed to add bank to database!\n");
+        goto err;
+    }
+
+    return 0;
+
+err:
+    free(buf);
+    return -1;
+}
+
+/* Get a bank from the database.
+ *
+ * Returns a pointer to the first value in the bank, or NULL if the bank
+ * doesn't exist. */
+dbval *get_bank(dict *db, const char name[4], uint32_t id)
+{
+    uint32_t buf[2];
+
+    memcpy(buf,name,4);
+    buf[1] = id;
+
+    return dictFetchValue(db, buf);
+}
+
+/* Function to iterate over the fields in a title bank text file. Works sort of
+ * like strtok(). Example:
+ *
+ *     char *item = iter_field(line);
+ *     while (item) {
+ *          // do something with item
+ *          item = iter_field(NULL);
+ *     }
+ */
+static char *iter_field(char *str)
+{
+    static char *ptr;
+    static char buf[81];
+
+    if (!str) str = ptr;
+
+    ptr = buf;
+
+    while (*str != '\x0' && ptr < (buf + 81)) {
+        if (ptr == (buf + 2) && buf[0] == '#' && buf[1] == '.') {
+            ptr = buf;
+            while (*str != '\x0' && *str++ != '#');
+            continue;
+        }
+
+        if (*str == ' ' || *str == '\n') {
+            str++;
+            if (ptr == buf) continue;
+            break;
+        }
+
+        *ptr++ = *str++;
+    }
+
+    if (ptr == buf) return NULL;
+
+    *ptr = '\x0';
+    ptr = str;
+    return buf;
+}
+
+/* Create a new database. */
+dict *db_init(void)
+{
+    dict *db = dictCreate(&titleDB, NULL);
+    return db;
+}
+
+/* Free a database. */
+void db_free(dict *db)
+{
+    dictRelease(db);
+}
+
+/* Load a title bank file into the database `db`.
+ *
+ * Returns 0 on success, -1 on error. */
+int load_file(dict *db, const char *filename)
+{
+    int i, index;
+    char *item, *s;
+    char line[256];
+    char idh[4];
+    uint32_t idn;
+    int mul;
+    int64_t value;
+    double float_value;
+    dbval *buf;
+    int buf_size;
+
+    buf_size = 10;
+
+    buf = malloc(buf_size*sizeof(dbval));
+
+    if (!buf) {
+        strcpy(db_err,strerror(errno));
+        return -1;
+    }
+
+    FILE *f = fopen(filename, "r");
+
+    if (!f) {
+        sprintf(db_err, "unable to open file '%s': %s", filename, strerror(errno));
+        return -1;
+    }
+
+    index = 0;
+    while (fgets(line, sizeof(line), f)) {
+        if (!strncmp(line,"*---",4)) {
+            /* Comment. */
+            continue;
+        } else if (!strncmp(line,"*.--",4)) {
+            /* Comment. */
+            continue;
+        } else if (!strncmp(line, "*LOG",4)) {
+            /* Control line which we don't care about. */
+            continue;
+        } else if (!strncmp(line, "*PRI",4)) {
+            /* Control line which we don't care about. */
+            continue;
+        } else if (!strncmp(line, "*US",3)) {
+            /* Control line which we don't care about. */
+            continue;
+        } else if (!strncmp(line, "*KI",3)) {
+            /* Control line which we don't care about. */
+            continue;
+        } else if (!strncmp(line, "*ANYWAY",7)) {
+            /* Control line which we don't care about. */
+            continue;
+        } else if (!strncmp(line, "*FIN",4)) {
+            /* Control line which we don't care about. */
+            continue;
+        } else if (!strncmp(line, "*DO",3)) {
+            if (index > 0) {
+                /* Save previous bank. */
+                buf = realloc(buf,sizeof(dbval)*index);
+
+                if (add_bank(db, idh, idn, buf) == DICT_ERR) {
+                    sprintf(db_err, "bank '%.4s' already exists in the database", idh);
+                    goto err;
+                }
+
+                buf_size = 10;
+
+                buf = malloc(buf_size*sizeof(dbval));
+
+                index = 0;
+            }
+
+            item = iter_field(line);
+            item = iter_field(NULL);
+            strncpy(idh,item,4);
+            item = iter_field(NULL);
+            idn = atoi(item);
+
+            continue;
+        }
+
+        item = iter_field(line);
+        mul = 1;
+        while (item) {
+            if ((s = strchr(item,'*'))) {
+                *s = '\0';
+                mul = atoi(item);
+                item = s+1;
+            }
+
+            if (buf_size < (index+1)*sizeof(dbval)) {
+                buf_size *= 2;
+                buf = realloc(buf, buf_size*sizeof(dbval));
+            }
+
+            if (!strncmp(item,"#x",2)) {
+                /* Hexadecimal input. */
+                value = strtol(item+2,NULL,16);
+                for (i = 0; i < mul; i++) buf[index++].u32 = value;
+            } else if (strchr(item,'E') || strchr(item,'.')) {
+                /* Floating point input. */
+                float_value = strtod(item,NULL);
+                for (i = 0; i < mul; i++) buf[index++].f = float_value;
+            } else {
+                /* Assume it's an integer. */
+                value = strtol(item,NULL,10);
+                for (i = 0; i < mul; i++) buf[index++].u32 = value;
+            }
+            item = iter_field(NULL);
+        }
+    }
+
+    buf = realloc(buf,sizeof(dbval)*index);
+
+    if (add_bank(db, idh, idn, buf) == DICT_ERR) {
+        sprintf(db_err, "bank '%.4s' already exists in the database", idh);
+        goto err;
+    }
+
+    fclose(f);
+
+    return 0;
+
+err:
+
+    if (buf) free(buf);
+    fclose(f);
+    return -1;
+}