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

diff --git a/zdab_utils.h b/zdab_utils.h
new file mode 100644
index 0000000..c5fc2c0
--- /dev/null
+++ b/zdab_utils.h
@@ -0,0 +1,202 @@
+#ifndef ZDAB_UTILS
+#define ZDAB_UTILS
+
+#include "Record_Info.h"
+#include <stdint.h>
+#include <stdlib.h> /* for size_t */
+
+// the builder won't put out events with NHIT > 10000
+// (note that these are possible due to hardware problems)
+// but XSNOED can write an event with up to 10240 channels
+#define MAX_NHIT 10240
+
+#ifdef SWAP_BYTES
+#define SWAP_INT32(a,b) swap_int32((int32_t *)(a),(b))
+#define SWAP_INT16(a,b) swap_int16((int16_t *)(a),(b))
+#else
+#define SWAP_INT32(a,b)
+#define SWAP_INT16(a,b)
+#endif
+
+/* Bit masks for the PMT flags field in the PMT bank. */
+#define KPF_DIS      0x00000001
+#define KPF_NOI      0x00000004
+#define KPF_PRP      0x00000008
+#define KPF_AFP      0x00000010
+#define KPF_FECD     0x00000020
+#define KPF_FT_MASK1 0x00000040
+#define KPF_XT       0x04000000
+#define KPF_NO_CAL   0x08000000
+#define KPF_BAD_CAL  0x10000000
+
+typedef struct EVBank {
+    /* Run number. */
+    uint32_t run;
+    /* Event number. */
+    uint32_t evn;
+    /* Data type. */
+    uint32_t dtp;
+    /* Julian date. */
+    uint32_t jdy;
+    /* Universal time (seconds). */
+    uint32_t ut1;
+    /* Universal time (nanoseconds). */
+    uint32_t ut2;
+    /* Date (format: yyyymmdd) */
+    uint32_t dte;
+    /* Time (format: hhmmsscc - cc is centisec). */
+    uint32_t hmsc;
+    /* Global trigger time in nsec (first word). */
+    uint32_t gtr1;
+    /* Second word. */
+    uint32_t gtr2;
+    /* Number of PMTs that fired. */
+    uint32_t npm;
+    /* Integrated charge in the event. */
+    uint32_t nph;
+    /* Sub-run number, or -1 if unknown. Only defined for run number >= 10614. */
+    uint32_t sub_run;
+    /* Packed word. */
+    uint32_t mc_pck;
+    /* ZDAB input/output record type, e.g. PMT, NCD, CMA, RUN, etc. */
+    uint32_t rec;
+    /* ZDAB_PMT format number (x10). */
+    uint32_t vpck;
+    /* Global trigger ID. */
+    uint32_t gtr_id;
+    /* Trigger word. */
+    uint32_t trg_type;
+    /* Digitized peak of analog sum. */
+    uint32_t peak;
+    /* Analog sum derivative. */
+    uint32_t diff;
+    /* Digitized integral of analog sum. */
+    uint32_t integral;
+    /* Trigger error bits plus spares. */
+    uint32_t err;
+    /* Data splitter blindness word. */
+    uint32_t data_set;
+    uint32_t spare1[3];
+    uint32_t ncd_status;
+    /* Number of multiplexer global records. */
+    uint32_t num_muxg;
+    /* Number of multiplexer records. */
+    uint32_t num_mux;
+    /* Number of scope records. */
+    uint32_t num_scope;
+    uint32_t spare2[5];
+    /* Upper 24 bits of flock */
+    uint32_t ncd_clk_up;
+    /* Lower clock register (32 bits). */
+    uint32_t ncd_clk_lw;
+    /* Latch register id. */
+    uint32_t ncd_reg;
+    /* GTID. */
+    uint32_t ncd_gtid;
+    /* Sync clear error. */
+    uint32_t ncd_sync;
+    uint32_t spare3[10];
+} EVBank;
+
+/* PMT Bank struct. This is a struct which mimics the PMT data bank in SNOMAN.
+ * The fields and comments were taken from the SNOMAN companion. */
+typedef struct PMTBank {
+    /* Tube number. For FECD, which has no associated tube number, this word is
+     * a copy of PIN. */
+    uint32_t pn;
+    /* Set of 1-bit flags. All undefined bits are set to zero. The following
+     * bits are defined:
+     *
+     * 0  KPF_DIS       Discard hit for all processing including fitting.
+     * 2  KPF_NOI       Noise.
+     * 3  KPF_PRP       Pre-pulse.
+     * 4  KPF_AFP       After-pulse.
+     * 5  KPF_FECD      Front-End Card Data i.e. non-PMT channel.
+     *                  Should only be set on the FECD bank.
+     * 6  KPF_FT_MASK+i Discarded by fitter method i (i = 1,2,..).
+     *                  For description of fitter methods see the FTx bank.
+     * 26 KPF_XT        Crosstalk.
+     * 27 KPF_NO_CAL    No calibration.
+     * 28 KPF_BAD_CAL   Bad calibration. */
+    uint32_t pf;
+    /* Time in nano-secs relative to event T0 - see EV. */
+    float pt;
+    /* Integrated charge. */
+    float phl;
+    /* Short-time integrated charge. */
+    float phs;
+    /* Low-gain integrated charge. */
+    float plx;
+    /* Time with the PMT jitter removed. */
+    float pt0;
+    /* Set of n-bit flags for the ith hit in list. All undefined bits are set
+     * to zero. The following bits are defined:
+     *
+     * Bit Field  DAQ Meaning
+     * 0-3        NCELL     (CMOS Cell address 0 to 15)
+     * 4-5        FLAG1      Bit 4 = CGT ES16.
+     *                       Bit 5 = CGT ES24.
+     * 6-9        FLAG2      Bit 6 = Missed Count
+     *                       Bit 7 = NC/CC Flag
+     *                       Bit 8 = LGISELECT (1=long sample)
+     *                       Bit 9 = CMOS ES16. */
+    uint32_t pif;
+    /* Uncalibrated times. */
+    float pit;
+    /* Uncalibrated high-gain, long int. charge. */
+    float pihl;
+    /* Uncalibrated high-gain, short int. charge. */
+    float pihs;
+    /* Uncalibrated log-gain, long int. charge. */
+    float pilx;
+    /* Uncalibrated unjittered time. */
+    float pit0;
+    /* PMT channel cell number. */
+    uint32_t cell;
+    /* CCC DAQ Circuit Number
+     *     = 1024*card + 32*crate + channel.
+     * Before uncalibration it holds the tube no. */
+    uint32_t pin;
+    /* TSLH of PMT. */
+    uint32_t tslh;
+    /* HCA information. */
+    uint32_t hca;
+    /* ECA validation status word. */
+    uint32_t eca_val;
+    /* PCA validation status word. */
+    uint32_t pca_val;
+    /* ANXX validation status word. */
+    uint32_t anxx;
+    /* ECA calibrated time (nsec). */
+    uint32_t ept;
+    /* ECA calibrated QHL (pedestal subtracted). */
+    uint32_t ehl;
+    /* ECA calibrated QHS (pedestal subtracted). */
+    uint32_t ehs;
+    /* ECA calibrated QLX (pedestal subtracted). */
+    uint32_t elx;
+    /* Non-walk corrected PMT time. */
+    uint32_t pt1;
+    /* Multiphoton PCA time. */
+    uint32_t ptm;
+    /* Multiphoton PCA PMT transit time RMS. */
+    uint32_t ptms;
+    /* Best charge (either QHS or QLX). */
+    uint32_t qm;
+    /* Best charge status word. */
+    uint32_t qms;
+    /* Charge correction for rate-dependent shifting. */
+    uint32_t qrc;
+} PMTBank;
+
+void unpack_ev(uint32_t *data, EVBank *b);
+void unpack_pmt(uint32_t *data, PMTBank *b);
+
+int isOrphan(aPmtEventRecord *pmtRecord);
+void swap_int32(int32_t *val_pt, int count);
+void swap_int16(int16_t *val_pt, int count);
+int swap_PmtRecord(aPmtEventRecord *aPmtRecord, size_t size);
+void swap_TrigRecord(struct TriggerInfo *aTrigRecord);
+void swap_RunRecord(struct RunRecord *aRunRecord);
+
+#endif