From 83f1e08ded354daf0f71cc7ec67f8de23bf69ad2 Mon Sep 17 00:00:00 2001
From: tlatorre <tlatorre@uchicago.edu>
Date: Tue, 27 Aug 2019 14:16:11 -0500
Subject: add the in time channel data cleaning cut from SNO

---
 src/dc.c | 43 +++++++++++++++++++++++++++++++++++++++++++
 src/dc.h |  8 ++++++++
 2 files changed, 51 insertions(+)

(limited to 'src')

diff --git a/src/dc.c b/src/dc.c
index ac062c6..7b5b1d7 100644
--- a/src/dc.c
+++ b/src/dc.c
@@ -29,6 +29,47 @@
 
 #define MAX_PAIRS 10000
 
+/* Returns 1 if the event is tagged as failing the ITC cut. The definition for
+ * the ITC cut comes from the SNOMAN companion:
+ *
+ *     This tags events having fewer than 60% of the tubes within 93 ns
+ *     coincidence ("in-time"). The number of in-time hits is the maximum
+ *     number in any position of a 93 ns sliding window applied to the PMT time
+ *     spectrum of the event.
+ *
+ */
+int is_itc(event *ev)
+{
+    size_t i, j;
+    double t_array[MAX_PMTS];
+    int nhit, n, maxn;
+    double start;
+
+    nhit = 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;
+
+        t_array[nhit++] = ev->pmt_hits[i].t;
+    }
+
+    /* Sort the times. */
+    gsl_sort(t_array,1,nhit);
+
+    maxn = 0;
+    for (i = 0; i < nhit; i++) {
+        start = t_array[i];
+        n = 1;
+        for (j = i+1; j < nhit; j++) {
+            if (t_array[j] < start + ITC_TIME_WINDOW) n += 1;
+            else break;
+        }
+        if (n > maxn) maxn = n;
+    }
+
+    return maxn/(double) nhit < ITC_TIME_FRACTION;
+}
+
 /* Returns 1 if the event is tagged as failing the FTS cut. The definition for
  * the FTS cut comes from the SNOMAN companion:
  *
@@ -153,6 +194,8 @@ uint32_t get_dc_word(event *ev, zebraFile *f, zebraBank *bmast, zebraBank *bev)
         status |= DC_OWL_TRIGGER;
     if (is_fts(ev))
         status |= DC_FTS;
+    if (is_itc(ev))
+        status |= DC_ITC;
 
     return status;
 }
diff --git a/src/dc.h b/src/dc.h
index 6569e0f..fdd7b0b 100644
--- a/src/dc.h
+++ b/src/dc.h
@@ -32,6 +32,13 @@
 #define DC_OWL            0x80
 #define DC_OWL_TRIGGER    0x100
 #define DC_FTS            0x200
+#define DC_ITC            0x400
+
+/* Length of the sliding window used in the ITC cut (ns). */
+#define ITC_TIME_WINDOW 93.0
+/* Fraction of hits which must be in a sliding window of 93 ns to pass the ITC
+ * cut. */
+#define ITC_TIME_FRACTION 0.6
 
 /* Height of PMTs which are considered to be at the "bottom" of the detector
  * for the neck cut. Comes from the SNOOP definition of a neck event in SNO. */
@@ -64,6 +71,7 @@
 /* QvNHIT ratio threshold. */
 #define QRATIO_THRESHOLD 0.25
 
+int is_itc(event *ev);
 int is_fts(event *ev);
 int is_owl_trigger(event *ev);
 int is_owl(event *ev);
-- 
cgit