add a bunch of SNO data quality cuts

This commit adds the following data quality cuts used in SNOMAN:

- neck
- qvnhit
- crate isotropy
- junk

Still need to test these.

5 files changed, 299 insertions, 0 deletions

diff --git a/src/dc.c b/src/dc.c
index 9d7ee8e..5056fbc 100644
--- a/src/dc.c
+++ b/src/dc.c
@@ -14,12 +14,271 @@
  * this program. If not, see <https://www.gnu.org/licenses/>.
  */
 
+#include "dc.h"
 #include "event.h"
 #include "sort.h"
 #include <stddef.h> /* for size_t */
 #include "misc.h"
 #include "pmt.h"
 #include "vector.h"
+#include <gsl/gsl_sort.h>
+#include "zdab_utils.h"
+#include "zebra.h"
+
+/* Returns 1 if the event fails the junk cut. The definition for the junk cut
+ * comes from the SNOMAN companion:
+ *
+ *     This will remove orphans, ECA data and events containing the same PMT
+ *     more than once.
+ *
+ * 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 ids[MAX_PMTS];
+    int nhit;
+    static int pmt_links[] = {KEV_PMT,KEV_OWL,KEV_LG,KEV_FECD,KEV_BUTT,KEV_NECK};
+    static char *pmt_names[] = {"PMT","OWL","LG","FECD","BUTT","NECK"};
+    size_t index[MAX_PMTS];
+    PMTBank bpmt;
+    zebraBank b;
+    int rv;
+
+    /* Fail orphans. */
+    if (bev->status & KEV_ORPHAN) return 1;
+
+    /* Fail any events with the ECA bit set. */
+    if (bmast->status & KMAST_ECA) return 1;
+
+    /* Next we check to see if any PMTs got hit more than once. To do this we
+     * construct an array of the PMT PIN numbers, sort it, and then check for
+     * duplicates. */
+
+    nhit = 0;
+    for (i = 0; i < LEN(pmt_links); i++) {
+        if (bev->links[pmt_links[i]-1] == 0) continue;
+
+        rv = zebra_get_bank(f,&b,bev->links[pmt_links[i]-1]);
+
+        if (rv) {
+            fprintf(stderr, "error getting %s bank: %s\n", pmt_names[i], zebra_err);
+            return -1;
+        }
+
+        while (1) {
+            unpack_pmt(b.data, &bpmt);
+            ids[nhit++] = bpmt.pin;
+            if (!b.next) break;
+
+            rv = zebra_get_bank(f,&b,b.next);
+
+            if (rv) {
+                fprintf(stderr, "error getting %s bank: %s\n", pmt_names[i], zebra_err);
+                return -1;
+            }
+        }
+    }
+
+    argsort(ids,nhit,index);
+
+    for (i = 0; i < nhit - 1; i++) {
+        /* Check if we have duplicate PMT hits. */
+        if (ids[i] == ids[i+1]) return 1;
+    }
+
+    return 0;
+}
+
+/* Returns 1 if the event fails the crate isotropy test. The definition of
+ * this cut comes from the SNOMAN companion:
+ *
+ *     This tags events where more than 70% of the hits occur on one crate and
+ *     more than 80% of those hits occur on two adjacent cards, including a
+ *     wrap around to count cards 0 and 15 together. 
+ *
+ */
+int crate_isotropy(event *ev)
+{
+    size_t i;
+    int crate_nhit[20];
+    int card_nhit[16];
+    int nhit;
+    int crate, card;
+    int max_crate;
+
+    for (i = 0; i < LEN(crate_nhit); i++) crate_nhit[i] = 0;
+
+    /* Loop over all the hits and add one to the crate_nhit array for each hit
+     * in a crate. */
+
+    nhit = 0;
+    for (i = 0; i < MAX_PMTS; i++) {
+        if (!ev->pmt_hits[i].hit) continue;
+
+        crate = i/512;
+
+        crate_nhit[crate] += 1;
+
+        nhit += 1;
+    }
+
+    /* Check if any crate has more than 70% of the hits. */
+    max_crate = -1;
+    for (i = 0; i < LEN(crate_nhit); i++) {
+        if (crate_nhit[i] > 7*nhit/10) max_crate = i;
+    }
+
+    /* If no crates had more than 70% of the hits, return 0. */
+    if (max_crate == -1) return 0;
+
+    for (i = 0; i < LEN(card_nhit); i++) card_nhit[i] = 0;
+
+    /* Now we count up the number of hits in adjacent cards. We store the
+     * number of hits in the card_nhit array as follows: Hits from card 15 or
+     * card 0 will be added to card_nhit[0], hits from cards 0 or 1 will go
+     * into card_nhit[1], etc. */
+
+    nhit = 0;
+    for (i = 0; i < MAX_PMTS; i++) {
+        if (!ev->pmt_hits[i].hit) continue;
+
+        crate = i/512;
+
+        if (crate != max_crate) continue;
+
+        card = (i % 512)/32;
+
+        card_nhit[card] += 1;
+        card_nhit[(card + 1) % LEN(card_nhit)] += 1;
+
+        nhit += 1;
+    }
+
+    /* Check to see if 80% of the hits in this crate came from two adjacent
+     * cards. */
+    for (i = 0; i < LEN(card_nhit); i++) {
+        if (card_nhit[i] > 8*nhit/10) return 1;
+    }
+
+    return 0;
+}
+
+/* Returns 1 if the event is tagged by the QvNHIT cut. The definition of
+ * this cut comes from the SNOMAN companion:
+ *
+ *     Using pedestal subtracted charge and a hardwired gain of 32.3 first
+ *     throw away the 10% of the tubes that have the largest charge and then
+ *     divide this by 0.9* NHIT. Tag the event if this ratio is below 0.25.
+ *
+ * I copied the logic from the SNOMAN file flt_q_nhit_cut.for. */
+int qvnhit(event *ev)
+{
+    size_t i;
+    double ehl[MAX_PMTS];
+    int nhit, max;
+    double sum, qhl_ratio;
+
+    nhit = 0;
+    for (i = 0; i < MAX_PMTS; i++) {
+        if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL) continue;
+
+        /* Require good calibrations. */
+        if (ev->pmt_hits[i].pf & (KPF_NO_CAL | KPF_BAD_CAL)) continue;
+
+        /* Remove unphysical charge hits. */
+        if (ev->pmt_hits[i].qhl <= -100) continue;
+
+        /* FIXME: SNOMAN code checks for pmt terminator. I don't know where
+         * that info is stored. */
+
+        /* Note: We use QHL instead of QHS here because that's how SNOMAN did
+         * it. In the SNOMAN code there is a comment which says:
+         *
+         *     we use QHL since it is a better measure of pickup (long integrate)
+         *
+         */
+        if (ev->pmt_hits[i].hit) {
+            ehl[nhit++] = ev->pmt_hits[i].ehl/32.3;
+        }
+    }
+
+    gsl_sort(ehl,1,nhit);
+
+    max = nhit*9/10;
+
+    /* Check that we have enough calibrated tubes. */
+    if (max < 5) return 0;
+
+    sum = 0.0;
+    for (i = 0; i < max; i++) {
+        sum += ehl[i];
+    }
+
+    qhl_ratio = sum/max;
+
+    return qhl_ratio < QRATIO_THRESHOLD;
+}
+
+/* Returns 1 if the event is classified as a neck event. The definition of
+ * a neck event comes from the SNOMAN companion:
+ *
+ *     This cuts events containing neck tubes. It requires that either both
+ *     tubes in the neck fire, or that one of those tubes fires and it has a
+ *     high charge and is early. High charge is defined by a neck tube having a
+ *     pedestal subtracted charge greater than 70 or less than -110. Early if
+ *     defined by the neck tube having an ECA time 70ns or more before the
+ *     average ECA time of the PSUP PMTS with z les than 0. After the cable
+ *     changes to the neck tubes this time difference changes to 15ns. 
+ *
+ */
+int is_neck_event(event *ev)
+{
+    size_t i;
+    int n;
+    int high_charge;
+    double avg_ept_psup;
+    int nhit;
+
+    /* First, we compute the average ECA time for all normal PMTs. */
+
+    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) {
+            avg_ept_psup += ev->pmt_hits[i].ept;
+            nhit += 1;
+        }
+    }
+
+    if (nhit < 1) return 0;
+
+    avg_ept_psup /= nhit;
+
+    /* Now, we check if two or more neck tubes fired *or* a single neck tube
+     * fired with a high charge and 70 ns before the average time of the normal
+     * PMTs. */
+    n = 0;
+    high_charge = 0;
+    for (i = 0; i < MAX_PMTS; i++) {
+        if (pmts[i].pmt_type == PMT_NECK && ev->pmt_hits[i].hit) {
+            n += 1;
+
+            /* FIXME: The SNOMAN companion says "Early if defined by the neck
+             * tube having an ECA time 70ns or more before the average ECA time
+             * of the PSUP PMTS with z les than 0." When does this change take
+             * place? */
+            if ((ev->pmt_hits[i].ehs > 70 || ev->pmt_hits[i].ehs < -110) &&
+                (ev->pmt_hits[i].ept < avg_ept_psup - 70.0)) high_charge = 1;
+        }
+    }
+
+    if (n >= 2 || high_charge) return 1;
+
+    return 0;
+}
 
 /* Returns 1 if the event is a flasher and 0 if it isn't. The definition of
  * a flasher event comes from
diff --git a/src/dc.h b/src/dc.h
index 503fb25..9df82c8 100644
--- a/src/dc.h
+++ b/src/dc.h
@@ -18,7 +18,15 @@
 #define DC_H
 
 #include "event.h"
+#include "zebra.h"
 
+/* QvNHIT ratio threshold. */
+#define QRATIO_THRESHOLD 0.25
+
+int junk_cut(zebraFile *f, zebraBank *bmast, zebraBank *bev);
+int crate_isotropy(event *ev);
+int qvnhit(event *ev);
+int is_neck_event(event *ev);
 int is_flasher(event *ev);
 
 #endif
diff --git a/src/event.h b/src/event.h
index a561a1f..0f1efce 100644
--- a/src/event.h
+++ b/src/event.h
@@ -28,7 +28,11 @@
 
 typedef struct pmt_hit {
     int hit;
+    float ept;
     float t;
+    float ehl;
+    float ehs;
+    float elx;
     float qhl;
     float qhs;
     float qlx;
@@ -36,6 +40,7 @@ typedef struct pmt_hit {
     uint16_t qihs;
     uint16_t qilx;
     int flags;
+    int pf;
 } pmt_hit;
 
 typedef struct event {
diff --git a/src/zdab_utils.c b/src/zdab_utils.c
index f0651d5..ba76b63 100644
--- a/src/zdab_utils.c
+++ b/src/zdab_utils.c
@@ -77,14 +77,24 @@ int get_event(zebraFile *f, event *ev, zebraBank *bev)
             crate = (bpmt.pin % 1024)/32;
             channel = bpmt.pin % 32;
             id = crate*512 + card*32 + channel;
+
+            if (ev->pmt_hits[id].hit) {
+                fprintf(stderr, "%i/%i/%i is in the PMT bank twice!\n", crate, card, channel);
+            }
+
             ev->pmt_hits[id].hit = 1;
+            ev->pmt_hits[id].ept = bpmt.ept;
             ev->pmt_hits[id].t = bpmt.pt;
             ev->pmt_hits[id].qihl = bpmt.pihl;
             ev->pmt_hits[id].qihs = bpmt.pihs;
             ev->pmt_hits[id].qilx = bpmt.pilx;
+            ev->pmt_hits[id].ehl = bpmt.ehl;
+            ev->pmt_hits[id].ehs = bpmt.ehs;
+            ev->pmt_hits[id].elx = bpmt.elx;
             ev->pmt_hits[id].qhl = bpmt.phl;
             ev->pmt_hits[id].qhs = bpmt.phs;
             ev->pmt_hits[id].qlx = bpmt.plx;
+            ev->pmt_hits[id].pf = bpmt.pf;
             /* Clear the PMT_FLAG_DIS bit. */
             ev->pmt_hits[id].flags &= ~PMT_FLAG_DIS;
             if (bpmt.pf & KPF_DIS)
diff --git a/src/zdab_utils.h b/src/zdab_utils.h
index a2804e1..1c138b9 100644
--- a/src/zdab_utils.h
+++ b/src/zdab_utils.h
@@ -50,6 +50,12 @@
 #define KMC_MCVX 2
 #define KMC_MCGN 1
 
+/* Status bits for the EV bank. */
+#define KEV_ORPHAN    0    /* Event contains orphan data. */
+#define KEV_TSLH      1    /* Event contains TSLH information. */
+#define KEV_HCA       2    /* Event contains HCA information. */
+#define KEV_MC_GE_ERR 3    /* Event produced from MC that had geometry errors suppressed. */
+
 /* Links for the EV bank. */
 #define KEV_NPFA     32
 #define KEV_NECL     31
@@ -84,6 +90,17 @@
 #define KEV_PT       2
 #define KEV_FT       1
 
+/* Status bits for the MAST bank. */
+
+#define KMAST_PMT    1    /* PMT data. */
+#define KMAST_PERM   2    /* Permanent data of any kind e.g. Run Header. */
+#define KMAST_NCD    3    /* NCD data. */
+#define KMAST_CMA    4    /* CMA data. */
+#define KMAST_ECA    5    /* ECA data. */
+#define KMAST_SUBRUN 6    /* Sub-run data, i.e. permanent bank data that is duplicated
+                             at the start of the second and subsequent sub runs to
+                             ensure the context is maintained across file boundaries. */
+
 /* Links for the MAST bank. */
 #define KMAST_RLOG 25 /* Runlog data. */
 #define KMAST_RLAI 24 /* Runlog analysis data. */