add a data cleaning cut to tag incoming muons

This commit adds a data cleaning cut to tag incoming muons by looking for early
OWL hits. It also significantly updates the flasher cut to catch more flashers.
In particular, the flasher cut now does the following:

- loops over *all* paddle cards with at least 4 hits instead of just the paddle
  cards with the most hits
- uses QLX to look for charge outliers in the paddle card
- fixes a few bugs (for example, uninitialized values in the charge array)
- adds a check to to see if the given slot is early with respect to all PMTs
  within 4 meters to catch the case where the flashing channel is missing from
  the event

4 files changed, 314 insertions, 84 deletions

diff --git a/src/dc.c b/src/dc.c
index 5056fbc..9f427a7 100644
--- a/src/dc.c
+++ b/src/dc.c
@@ -24,6 +24,98 @@
 #include <gsl/gsl_sort.h>
 #include "zdab_utils.h"
 #include "zebra.h"
+#include <gsl/gsl_statistics_double.h>
+
+/* Returns 1 if the event is tagged as an incoming muon. The goal of this cut
+ * is to tag external muons *without* tagging muons which start within the PSUP
+ * and then exit. In order to achieve that we use the ECA time of the OWL
+ * tubes.
+ *
+ * Ideally we would have the full calibrated time of the OWL tubes, but there
+ * was never any PCA for the OWL tubes. I also discussed this with Josh and he
+ * said that the cables were never changed for the OWL tubes so we should at
+ * least be able to rely on them being consistent.
+ *
+ * Since the average PCA correction is of order +/- 5 ns and the round trip
+ * time for light from the edge of the PSUP to the cavity and back should be
+ * longer than 15 ns (4 meters), the ECA time should be enough to distinguish
+ * incoming vs outgoing muons.
+ *
+ * This function tags an event as an incoming muon if the following criteria
+ * are met:
+ *
+ * - calibrated nhit >= 150
+ * - at least 5 calibrated OWL hits
+ * - at least 1 OWL hit with a QHS >= 20.0
+ * - at least 2 OWL hits with ECA times earlier than the 10th percentile of the
+ *   PSUP PMT times
+ *
+ * Idea: Should we also require at least one of the early tubes to have a high
+ * charge? */
+int is_muon(event *ev)
+{
+    size_t i;
+    double ept[MAX_PMTS];
+    size_t nhit, nhit_owl, nhit_owl_early;
+    double early_time, max_owl_qhs;
+
+    nhit = 0;
+    nhit_owl = 0;
+    for (i = 0; i < MAX_PMTS; i++) {
+        /* Require good calibrations.
+         *
+         * Note: Typically I would do this with the following check:
+         *
+         *     ev->pmt_hits[i].pf & (KPF_NO_CAL | KPF_BAD_CAL)
+         *
+         * However, Javi mentioned in an email that for very early hits which
+         * are in the TAC curl region, the calibration processor will flag
+         * them, set the time to -9999.0, and set the KPF_BAD_CAL bit. However,
+         * in this case for some reason the ECA time still looks reasonable.
+         * Therefore, since we only use the ECA time here we instead just make
+         * sure that the ECA time is something reasonable. */
+        if (!ev->pmt_hits[i].hit || ev->pmt_hits[i].ept <= -100) continue;
+
+        switch (pmts[i].pmt_type) {
+        case PMT_OWL:
+            if (nhit_owl == 0) {
+                max_owl_qhs = ev->pmt_hits[i].qhs;
+            } else if (ev->pmt_hits[i].qhs > max_owl_qhs) {
+                max_owl_qhs = ev->pmt_hits[i].qhs;
+            }
+            nhit_owl += 1;
+            break;
+        case PMT_NORMAL:
+            if (!ev->pmt_hits[i].flags) ept[nhit++] = ev->pmt_hits[i].ept;
+            break;
+        }
+    }
+
+    /* Require at least 150 normal hits. */
+    if (nhit < MUON_MIN_NHIT) return 0;
+
+    /* Require at least 5 OWL hits. */
+    if (nhit_owl < MUON_MIN_OWL_NHIT) return 0;
+
+    /* Require at least 1 OWL tube to have a QHS >= 20.0. */
+    if (max_owl_qhs < MUON_OWL_QHS) return 0;
+
+    /* Sort the times. */
+    gsl_sort(ept,1,nhit);
+
+    /* Calculate the time to define early OWL hits. */
+    early_time = gsl_stats_quantile_from_sorted_data(ept,1,nhit,MUON_OWL_TIME_PERCENTILE);
+
+    nhit_owl_early = 0;
+    for (i = 0; i < MAX_PMTS; i++) {
+        /* Require good calibrations. */
+        if (!ev->pmt_hits[i].hit || ev->pmt_hits[i].ept <= -100) continue;
+
+        if (pmts[i].pmt_type == PMT_OWL && ev->pmt_hits[i].ept < early_time) nhit_owl_early += 1;
+    }
+
+    return nhit_owl_early >= 2;
+}
 
 /* Returns 1 if the event fails the junk cut. The definition for the junk cut
  * comes from the SNOMAN companion:
@@ -31,11 +123,14 @@
  *     This will remove orphans, ECA data and events containing the same PMT
  *     more than once.
  *
+ * I also added a check that will flag the event if it contains more than
+ * MAX_PMTS PMT hits (which shouldn't be possible).
+ *
  * Note: This function may return -1 if there is an error traversing the PMT
  * banks. */
 int junk_cut(zebraFile *f, zebraBank *bmast, zebraBank *bev)
 {
-    size_t i;
+    int i;
     int ids[MAX_PMTS];
     int nhit;
     static int pmt_links[] = {KEV_PMT,KEV_OWL,KEV_LG,KEV_FECD,KEV_BUTT,KEV_NECK};
@@ -68,7 +163,14 @@ int junk_cut(zebraFile *f, zebraBank *bmast, zebraBank *bev)
 
         while (1) {
             unpack_pmt(b.data, &bpmt);
+
             ids[nhit++] = bpmt.pin;
+
+            if (nhit >= MAX_PMTS) {
+                fprintf(stderr, "event with more than %i PMT hits\n", MAX_PMTS);
+                return 1;
+            }
+
             if (!b.next) break;
 
             rv = zebra_get_bank(f,&b,b.next);
@@ -240,14 +342,17 @@ int is_neck_event(event *ev)
     double avg_ept_psup;
     int nhit;
 
-    /* First, we compute the average ECA time for all normal PMTs. */
+    /* First, we compute the average ECA time for all the PSUP PMTs with z <= 0. */
 
     nhit = 0;
     avg_ept_psup = 0.0;
     for (i = 0; i < MAX_PMTS; i++) {
-        if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL) continue;
+        if (!ev->pmt_hits[i].hit || ev->pmt_hits[i].ept <= -100) continue;
 
-        if (ev->pmt_hits[i].hit) {
+        /* Require good calibrations. */
+        if (ev->pmt_hits[i].pf & (KPF_NO_CAL | KPF_BAD_CAL)) continue;
+
+        if (ev->pmt_hits[i].hit && pmts[i].pos[2] < 0) {
             avg_ept_psup += ev->pmt_hits[i].ept;
             nhit += 1;
         }
@@ -263,7 +368,10 @@ int is_neck_event(event *ev)
     n = 0;
     high_charge = 0;
     for (i = 0; i < MAX_PMTS; i++) {
-        if (pmts[i].pmt_type == PMT_NECK && ev->pmt_hits[i].hit) {
+        /* Require good calibrations. */
+        if (!ev->pmt_hits[i].hit || ev->pmt_hits[i].ept <= -100) continue;
+
+        if (pmts[i].pmt_type == PMT_NECK) {
             n += 1;
 
             /* FIXME: The SNOMAN companion says "Early if defined by the neck
@@ -280,6 +388,81 @@ int is_neck_event(event *ev)
     return 0;
 }
 
+/* Returns 1 if the average time for the hits in the same card as
+ * `flasher_pmt_id` is at least 10 ns earlier than all the hits within 4
+ * meters of the potential flasher.
+ *
+ * This cut is designed to help flag flashers in the is_flasher() function for
+ * events in which the flashing channel may be missing from the event, so there
+ * is no high charge channel to flag it. */
+int is_slot_early(event *ev, int flasher_pmt_id)
+{
+    int i;
+    double t_slot, t_nearby;
+    int n_slot, n_nearby;
+    double flasher_pos[3];
+    double pmt_dir[3];
+    double distance;
+
+    COPY(flasher_pos,pmts[flasher_pmt_id].pos);
+
+    n_slot = 0;
+    n_nearby = 0;
+    t_slot = 0.0;
+    t_nearby = 0.0;
+    for (i = 0; i < MAX_PMTS; i++) {
+        if (!ev->pmt_hits[i].hit) continue;
+
+        /* Require good calibrations.
+         *
+         * Note: For some reason PMT hits in the flasher slot sometimes are
+         * marked as having bad calibrations and get a time of -9999. I still
+         * don't understand why this is, so for now we don't check the
+         * calibration bits and instead just check to see that the ECA time is
+         * reasonable.
+         *
+         * Here is what Javi said via email:
+         *
+         * If SNO calibrations are the same than SNO+ (and I'm pretty sure they
+         * are, at last for this matter) there is a hit-level time validation,
+         * i.e., if the TAC value is determined to be in the curl region of the
+         * time slope (high TAC = very early hits), they get flagged and a
+         * calibration is not provided (which is compatible with a flashing
+         * channel). If this is the case, you can fall back to raw TACs and QHS
+         * when calibration is not available. */
+        if (ev->pmt_hits[i].ept < -100) continue;
+
+        if (i/32 == flasher_pmt_id/32) {
+            /* This hit is in the same slot as the potential flasher. */
+            t_slot += ev->pmt_hits[i].ept;
+            n_slot += 1;
+            continue;
+        }
+
+        /* Calculate the distance from the current channel to the high charge
+         * channel. */
+        SUB(pmt_dir,pmts[i].pos,flasher_pos);
+        distance = NORM(pmt_dir);
+
+        if (distance < 400.0) {
+            t_nearby += ev->pmt_hits[i].ept;
+            n_nearby += 1;
+            continue;
+        }
+    }
+
+    /* If there are no PMTs within 4 meters of the slot, it's almost certainly
+     * a flasher. */
+    if (n_slot > 4 && n_nearby == 0) return 1;
+
+    if (n_slot == 0 || n_nearby == 0) return 0;
+
+    t_slot /= n_slot;
+    t_nearby /= n_nearby;
+
+    return t_slot < t_nearby + 10.0;
+}
+
 /* Returns 1 if the event is a flasher and 0 if it isn't. The definition of
  * a flasher event comes from
  * http://detector.sno.laurentian.ca/~detector/private/snoop/doc/eventIDs.html,
@@ -317,18 +500,19 @@ int is_neck_event(event *ev)
 int is_flasher(event *ev)
 {
     int hits_pc[1280] = {0};
-    int qhs, qhl;
+    int qhs, qhl, qlx;
     int qhl_pc[8] = {0};
     int qhs_pc[8] = {0};
+    int qlx_pc[8] = {0};
     double t_pc[8] = {0};
+    int channel_pc[8] = {0};
     size_t crate, card, channel, id;
     int i, j;
-    size_t index[1280], index_qhs[8], index_qhl[8];
+    size_t index[1280], index_qhs[8], index_qhl[8], index_qlx[8];
     double t;
     int flasher_pc;
     double flasher_pos[3];
     double pmt_dir[3];
-    int flasher = 0;
     double distance;
     int nhit, nhit_late, nhit_far;
 
@@ -358,12 +542,12 @@ int is_flasher(event *ev)
     for (j = 1279; j >= 0; j--) {
 	id = index[j];
 
-        /* If this paddle card doesn't have as many hits as the most hit paddle
-         * card and we haven't found a flasher yet, then it's not a flasher. */
-        if (hits_pc[id] < hits_pc[index[1279]]) return 0;
+        /* Make sure this paddle card has at least 4 hits. */
+        if (hits_pc[id] < 4) return 0;
 
+        nhit = 0;
 	for (i = 0; i < MAX_PMTS; i++) {
-            if (!ev->pmt_hits[i].hit) continue;
+            if (!ev->pmt_hits[i].hit || pmts[i].pmt_type != PMT_NORMAL) continue;
 
             crate = i/512;
             card = (i % 512)/32;
@@ -376,101 +560,122 @@ int is_flasher(event *ev)
             /* Get the uncalibrated QHS and QHL charges. */
 	    qhs = ev->pmt_hits[i].qihs;
 	    qhl = ev->pmt_hits[i].qihl;
+	    qlx = ev->pmt_hits[i].qilx;
 
             if (qhs < 300) {
                 /* QHS values below 300 are set to 4095. */
-                qhs_pc[channel % 8] = 4095;
+                qhs_pc[nhit] = 4095;
             } else {
-                qhs_pc[channel % 8] = qhs;
+                qhs_pc[nhit] = qhs;
             }
 
             if (qhl < 300) {
                 /* QHL values below 300 are set to 4095. */
-                qhl_pc[channel % 8] = 4095;
+                qhl_pc[nhit] = 4095;
             } else {
-                qhl_pc[channel % 8] = qhl;
+                qhl_pc[nhit] = qhl;
             }
 
-            t_pc[channel % 8] = ev->pmt_hits[i].t;
-	}
+            if (qlx < 300) {
+                /* QLX values below 300 are set to 4095. */
+                qlx_pc[nhit] = 4095;
+            } else {
+                qlx_pc[nhit] = qlx;
+            }
 
-        /* Sort the QHS and QHL values by size. */
-	argsort(qhs_pc, LEN(qhs_pc), index_qhs);
-	argsort(qhl_pc, LEN(qhl_pc), index_qhl);
-
-        /* Check if this paddle card has QHS or QHL values > 4000, and the QHS
-         * and QHL values for the other channels in the card are less than 700. */
-        if (qhs_pc[index_qhs[7]] > 4000 || qhl_pc[index_qhl[7]] > 4000) {
-	    if (qhs_pc[index_qhs[6]] < 800 && qhl_pc[index_qhl[6]] < 800) {
-                /* Flasher! */
-                flasher = 1;
-                break;
-	    }
+            t_pc[nhit] = ev->pmt_hits[i].t;
+
+            channel_pc[nhit] = channel;
+
+            nhit += 1;
 	}
-    }
 
-    /* If we didn't find a flasher, return 0. */
-    if (!flasher) return 0;
-
-    flasher_pc = index[j];
-    crate = index[j]/64;
-    card = (index[j] % 64)/4;
-    if (qhs_pc[index_qhs[7]] > 4000) {
-        channel = (index[j] % 4)*8 + index_qhs[7];
-        t = t_pc[index_qhs[7]];
-    } else {
-        channel = (index[j] % 4)*8 + index_qhl[7];
-        t = t_pc[index_qhl[7]];
-    }
-    COPY(flasher_pos,pmts[crate*512 + card*32 + channel].pos);
+        if (nhit < 2) {
+            fprintf(stderr, "only 2 hits on %zu/%zu/PC%zu but expected at least 4!\n",
+                    id/64, (id % 64)/4, id % 4);
+            continue;
+        }
 
-    /* Check that 70% of the regular PMTs in the event which are not in the
-     * flasher slot fired more than 50 ns after the high charge channel and are
-     * more than 12 meters away from the high charge channel. */
-    nhit = 0;
-    /* Number of PMTs which fired 50 ns after the high charge channel. */
-    nhit_late = 0;
-    /* Number of PMTs which are 12 meters or more away from the high charge
-     * channel. */
-    nhit_far = 0;
-    for (i = 0; i < MAX_PMTS; i++) {
-        /* Skip PMTs which weren't hit. */
-        if (!ev->pmt_hits[i].hit) continue;
+        /* Sort the QHS and QHL values by size. */
+	argsort(qhs_pc, nhit, index_qhs);
+	argsort(qhl_pc, nhit, index_qhl);
+	argsort(qlx_pc, nhit, index_qlx);
 
-        /* Only count regular PMTs without any flags. */
-        if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL) continue;
+        flasher_pc = index[j];
+        crate = index[j]/64;
+        card = (index[j] % 64)/4;
+        channel = channel_pc[index_qhs[nhit-1]];
 
-	crate = i/512;
-	card = (i % 512)/32;
-	channel = i % 32;
+        /* Check if this paddle card has QHS, QHL, or QLX values > 4000, and the QHS,
+         * QHL, and QLX values for the other channels in the card are less than 700.
+         *
+         * Alternatively, if there isn't an obvious charge outlier, it's still
+         * possible that the flasher PMT hit is missing from the event or that
+         * the QHS/QHL/QLX values were just outside of the threshold. In this
+         * case, we look to see if the TAC for the proposed slot is early with
+         * respect to neighboring PMTs. */
+        if ((qhs_pc[index_qhs[nhit-1]] < 4000 || qhs_pc[index_qhs[nhit-2]] > 800) &&
+	    (qhl_pc[index_qhl[nhit-1]] < 4000 || qhl_pc[index_qhl[nhit-2]] > 800) &&
+	    (qlx_pc[index_qlx[nhit-1]] < 4000 || qlx_pc[index_qlx[nhit-2]] > 800) &&
+            !is_slot_early(ev, crate*512 + card*32 + channel)) continue;
+
+        if (qhs_pc[index_qhs[nhit-1]] > 4000) {
+            channel = channel_pc[index_qhs[nhit-1]];
+            t = t_pc[index_qhs[nhit-1]];
+        } else if (qhl_pc[index_qhl[nhit-1]] > 4000) {
+            channel = channel_pc[index_qhl[nhit-1]];
+            t = t_pc[index_qhl[nhit-1]];
+        } else {
+            channel = channel_pc[index_qlx[nhit-1]];
+            t = t_pc[index_qlx[nhit-1]];
+        }
+        COPY(flasher_pos,pmts[crate*512 + card*32 + channel].pos);
+
+        /* Check that 70% of the regular PMTs in the event which are not in the
+         * flasher slot fired more than 50 ns after the high charge channel and are
+         * more than 12 meters away from the high charge channel. */
+        nhit = 0;
+        /* Number of PMTs which fired 50 ns after the high charge channel. */
+        nhit_late = 0;
+        /* Number of PMTs which are 12 meters or more away from the high charge
+         * channel. */
+        nhit_far = 0;
+        for (i = 0; i < MAX_PMTS; i++) {
+            /* Skip PMTs which weren't hit. */
+            if (!ev->pmt_hits[i].hit) continue;
 
-	id = crate*64 + card*4 + channel/8;
+            /* Only count regular PMTs without any flags. */
+            if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL) continue;
 
-        /* Skip PMTs in the same paddle card as the high charge channel. */
-        if (id == flasher_pc) continue;
+            crate = i/512;
+            card = (i % 512)/32;
+            channel = i % 32;
 
-        nhit += 1;
+            id = crate*64 + card*4 + channel/8;
 
-        /* Calculate the distance from the current channel to the high charge
-         * channel. */
-        SUB(pmt_dir,pmts[i].pos,flasher_pos);
-        distance = NORM(pmt_dir);
+            /* Skip PMTs in the same card as the high charge channel. */
+            if (id/4 == flasher_pc/4) continue;
 
-        /* If this channel is 12 meters or more away, increment nhit_far. */
-        if (distance >= 1200.0) nhit_far += 1;
+            nhit += 1;
 
-        /* If this channel fired more than 50 ns after the high charge channel,
-         * increment nhit_late. */
-        if (ev->pmt_hits[i].t > t + 50.0) nhit_late += 1;
-    }
+            /* Calculate the distance from the current channel to the high charge
+             * channel. */
+            SUB(pmt_dir,pmts[i].pos,flasher_pos);
+            distance = NORM(pmt_dir);
+
+            /* If this channel is 12 meters or more away, increment nhit_far. */
+            if (distance >= 1200.0) nhit_far += 1;
 
-    /* If less than 70% of the regular PMTs fired within 50 ns of the high
-     * charge channel, then it's not a flasher. */
-    if (nhit_late < nhit*0.7) return 0;
+            /* If this channel fired more than 50 ns after the high charge channel,
+             * increment nhit_late. */
+            if (ev->pmt_hits[i].t > t + 50.0) nhit_late += 1;
+        }
 
-    /* If less than 70% of the regular PMTs were closer than 12 meters away
-     * from the high charge channel, then it's not a flasher. */
-    if (nhit_far < nhit*0.7) return 0;
+        /* If at least 70% of the regular PMTs fired within 50 ns and were at
+         * least 12 meters away from the high charge channel, then it's a
+         * flasher! */
+        if (nhit_late >= 7*nhit/10 && nhit_far >= 7*nhit/10) return 1;
+    }
 
-    return 1;
+    return 0;
 }
diff --git a/src/dc.h b/src/dc.h
index 9df82c8..89a7670 100644
--- a/src/dc.h
+++ b/src/dc.h
@@ -20,9 +20,21 @@
 #include "event.h"
 #include "zebra.h"
 
+/* Minimum number of normal PMTs which must be hit to be tagged as an incoming
+ * muon. */
+#define MUON_MIN_NHIT 150
+/* Minimum number of OWL PMTs which must be hit to be tagged as an incoming
+ * muon. */
+#define MUON_MIN_OWL_NHIT 5
+/* Minimum value for the maximum OWL QHS value. */
+#define MUON_OWL_QHS 20.0
+/* Percentile of PSUP PMT times to define early OWL hits. */
+#define MUON_OWL_TIME_PERCENTILE 0.1
+
 /* QvNHIT ratio threshold. */
 #define QRATIO_THRESHOLD 0.25
 
+int is_muon(event *ev);
 int junk_cut(zebraFile *f, zebraBank *bmast, zebraBank *bev);
 int crate_isotropy(event *ev);
 int qvnhit(event *ev);
diff --git a/src/fit.c b/src/fit.c
index 40f67f5..e4e054f 100644
--- a/src/fit.c
+++ b/src/fit.c
@@ -6172,7 +6172,7 @@ skip_mc:
                 last_run = ev.run;
             }
 
-            rv = get_event(f,&ev,&b);
+            if (get_event(f,&ev,&b)) goto skip_event;
 
             nhit = get_nhit(&ev);
 
@@ -6181,6 +6181,11 @@ skip_mc:
                 fprintf(fout, "        gtid: %i\n", ev.gtid);
                 fprintf(fout, "        nhit: %zu\n", nhit);
                 fprintf(fout, "        is_flasher: %i\n", is_flasher(&ev));
+                fprintf(fout, "        is_muon: %i\n", is_muon(&ev));
+                fprintf(fout, "        is_junk: %i\n", junk_cut(f, &bmast, &b));
+                fprintf(fout, "        is_crate_isotropy: %i\n", crate_isotropy(&ev));
+                fprintf(fout, "        is_qvnhit: %i\n", qvnhit(&ev));
+                fprintf(fout, "        is_neck_event: %i\n", is_neck_event(&ev));
             }
 
             if (nhit < min_nhit) goto skip_event;
diff --git a/src/zdab_utils.c b/src/zdab_utils.c
index ba76b63..6dffe0f 100644
--- a/src/zdab_utils.c
+++ b/src/zdab_utils.c
@@ -78,6 +78,11 @@ int get_event(zebraFile *f, event *ev, zebraBank *bev)
             channel = bpmt.pin % 32;
             id = crate*512 + card*32 + channel;
 
+            if (id >= MAX_PMTS) {
+                fprintf(stderr, "PMT hit from %i/%i/%i\n", crate, card, channel);
+                return -1;
+            }
+
             if (ev->pmt_hits[id].hit) {
                 fprintf(stderr, "%i/%i/%i is in the PMT bank twice!\n", crate, card, channel);
             }
@@ -105,6 +110,9 @@ int get_event(zebraFile *f, event *ev, zebraBank *bev)
                 fprintf(stderr, "%i/%i/%i has PMT type %s but expected %s based on bank\n", crate, card, channel, pmt_type_string, pmt_names[i]);
             }
 
+            /* For now we assume SNOMAN has it correct. */
+            pmts[id].pmt_type = pmt_types[i];
+
             if (!b.next) break;
 
             rv = zebra_get_bank(f,&b,b.next);