add a bunch more data cleaning cuts from SNO

This commit adds the following data cleaning cuts:

- FTS
- OWL trigger
- OWL
- ESUM

3 files changed, 130 insertions, 0 deletions

diff --git a/src/Record_Info.h b/src/Record_Info.h
index 91e060f..f3241c6 100644
--- a/src/Record_Info.h
+++ b/src/Record_Info.h
@@ -214,6 +214,7 @@ typedef struct TriggerInfo {
 #define TRIG_SPECIAL_RAW        0x00800000
 #define TRIG_NCD_MUX            0x01000000
 #define TRIG_SOFT_GT            0x02000000
+#define TRIG_MISS               0x04000000
 
 // -------------------------------------------------------------------------------------
 //   EPED record
diff --git a/src/dc.c b/src/dc.c
index 230e373..807d5cc 100644
--- a/src/dc.c
+++ b/src/dc.c
@@ -25,6 +25,108 @@
 #include "zdab_utils.h"
 #include "zebra.h"
 #include <gsl/gsl_statistics_double.h>
+#include "Record_Info.h"
+
+#define MAX_PAIRS 10000
+
+/* Returns 1 if the event is tagged as failing the FTS cut. The definition for
+ * the FTS cut comes from the SNOMAN companion:
+ *
+ *     This tags events where the median time difference of hit PMT pairs
+ *     within 3 m of each other is greater than 6.8 ns. Pairs with large time
+ *     differences (> 25 ns) are not included in the median and at least 15
+ *     pairs satisfying the distance and time difference criteria are required
+ *     before the event can be tagged.
+ *
+ * I used the SNOMAN implementation in flt_fts_cut.for for reference.
+ *
+ */
+int is_fts(event *ev)
+{
+    size_t i, j;
+    double pmt_dir[3], distance, dt;
+    double deltat_array[MAX_PAIRS];
+    size_t count;
+    double dtmedian;
+
+    count = 0;
+    for (i = 0; i < MAX_PMTS; i++) {
+        /* Only loop over hit PMTs with good calibrations. */
+        if (!ev->pmt_hits[i].hit || ev->pmt_hits[i].flags) continue;
+
+        for (j = 0; j < MAX_PMTS; j++) {
+            /* Only loop over hit PMTs with good calibrations. */
+            if (i == j || !ev->pmt_hits[j].hit || ev->pmt_hits[j].flags) continue;
+
+            /* Calculate distance between PMTs. */
+            SUB(pmt_dir,pmts[i].pos,pmts[j].pos);
+            distance = NORM(pmt_dir);
+
+            /* Get the absolute value of the time difference. */
+            dt = fabs(ev->pmt_hits[i].t - ev->pmt_hits[j].t);
+
+            if (dt <= 25.0 && distance < 300.0) {
+                deltat_array[count++] = dt;
+                if (count >= MAX_PAIRS) goto end;
+            }
+        }
+    }
+
+end:
+
+    if (count > FTS_COUNT_THRESH) {
+        gsl_sort(deltat_array,1,count);
+        dtmedian = gsl_stats_median_from_sorted_data(deltat_array,1,count);
+        return dtmedian < FTS_MEDIAN_CUT;
+    }
+
+    return 0;
+}
+
+/* Returns 1 if the event is tagged as an OWL trigger event. The definition for
+ * the OWL trigger cut comes from the SNOMAN companion:
+ *
+ *     This tags events where the OWLESUMHI trigger fires.
+ *
+ */
+int is_owl_trigger(event *ev)
+{
+    return (ev->trigger_type & TRIG_OWLE_HI) != 0;
+}
+
+/* Returns 1 if the event is tagged as an OWL event. The definition for the OWL
+ * cut comes from the SNOMAN companion:
+ *
+ *     This tags events where the number of OWLs + BUTTs is 3 or more.
+ *
+ */
+int is_owl(event *ev)
+{
+    size_t i;
+    size_t nhit;
+
+    nhit = 0;
+    for (i = 0; i < MAX_PMTS; i++) {
+        if (ev->pmt_hits[i].hit && (pmts[i].pmt_type == PMT_OWL || pmts[i].pmt_type == PMT_BUTT))
+            nhit += 1;
+    }
+
+    return nhit >= 3;
+}
+
+/* Returns 1 if the event fails the ESUM cut. The definition for the ESUM cut
+ * comes from the SNOMAN companion:
+ *
+ *     Tags an event if the ESUMLO AND ESUMHI trigger bits are set, AND NONE of
+ *     the following trigger bits are set: (NHIT100LO, NHIT100MED,
+ *     NHIT100HI,NHIT20, NHIT20LB, OWLN, OWLELO, OWLEHI, PULSE_GT, PRESCALE,
+ *     MISS).
+ *
+ */
+int is_esum(event *ev)
+{
+    return (ev->trigger_type & (TRIG_ESUM_LO | TRIG_ESUM_HI | TRIG_NHIT_100_LO | TRIG_NHIT_100_MED | TRIG_NHIT_100_HI | TRIG_NHIT_20 | TRIG_NHIT_20_LB | TRIG_OWLN | TRIG_OWLE_LO | TRIG_OWLE_HI | TRIG_PULSE_GT | TRIG_PRESCALE | TRIG_MISS)) == (TRIG_ESUM_LO | TRIG_ESUM_HI);
+}
 
 /* Returns the data cleaning bitmask for the event `ev`. */
 uint32_t get_dc_word(event *ev, zebraFile *f, zebraBank *bmast, zebraBank *bev)
@@ -43,6 +145,14 @@ uint32_t get_dc_word(event *ev, zebraFile *f, zebraBank *bmast, zebraBank *bev)
         status |= DC_NECK;
     if (is_flasher(ev))
         status |= DC_FLASHER;
+    if (is_esum(ev))
+        status |= DC_ESUM;
+    if (is_owl(ev))
+        status |= DC_OWL;
+    if (is_owl_trigger(ev))
+        status |= DC_OWL_TRIGGER;
+    if (is_fts(ev))
+        status |= DC_FTS;
 
     return status;
 }
diff --git a/src/dc.h b/src/dc.h
index 4635806..f613b93 100644
--- a/src/dc.h
+++ b/src/dc.h
@@ -28,6 +28,21 @@
 #define DC_QVNHIT         0x8
 #define DC_NECK           0x10
 #define DC_FLASHER        0x20
+#define DC_ESUM           0x40
+#define DC_OWL            0x80
+#define DC_OWL_TRIGGER    0x100
+#define DC_FTS            0x200
+
+/* Delta-t threshold for PMT pairs in the FTS cut (ns). */
+#define FTS_DT_THRESH 25.0
+/* Distance threshold for PMT pairs in the FTS cut (cm). */
+#define FTS_DIST_THRESH 300.0
+/* Minimum number of PMT pairs required for the FTS cut. */
+#define FTS_COUNT_THRESH 15
+/* For the FTS cut the median time difference between PMT pairs within 3 meters
+ * is computed and if this difference is greater than FTS_MEDIAN_CUT the event
+ * doesn't pass. */
+#define FTS_MEDIAN_CUT 6.8
 
 /* Minimum number of normal PMTs which must be hit to be tagged as an incoming
  * muon. */
@@ -43,6 +58,10 @@
 /* QvNHIT ratio threshold. */
 #define QRATIO_THRESHOLD 0.25
 
+int is_fts(event *ev);
+int is_owl_trigger(event *ev);
+int is_owl(event *ev);
+int is_esum(event *ev);
 uint32_t get_dc_word(event *ev, zebraFile *f, zebraBank *bmast, zebraBank *bev);
 int is_muon(event *ev);
 int junk_cut(zebraFile *f, zebraBank *bmast, zebraBank *bev);