add a first draft of a data cleaning cut to detect breakdown events

1 files changed, 43 insertions, 0 deletions

diff --git a/src/dc.c b/src/dc.c
index 5f2d4c6..6195eac 100644
--- a/src/dc.c
+++ b/src/dc.c
@@ -196,6 +196,8 @@ uint32_t get_dc_word(event *ev, zebraFile *f, zebraBank *bmast, zebraBank *bev)
         status |= DC_FTS;
     if (is_itc(ev))
         status |= DC_ITC;
+    if (is_breakdown(ev))
+        status |= DC_BREAKDOWN;
 
     return status;
 }
@@ -604,6 +606,47 @@ int is_neck_event(event *ev)
     return 0;
 }
 
+/* Returns 1 if the event is classified as a breakdown event. The definition
+ * for a breakdown event is:
+ *
+ *   - nhit > 1000
+ *
+ *   - a single crate has at least 256 hits and the median TAC for that crate
+ *     is more than 500 from the median TAC for all the other crates with at
+ *     least 20 hits
+ *
+ */
+int is_breakdown(event *ev)
+{
+    int i;
+    double tac[20][512];
+    int nhit[20] = {0};
+    double median_tac[20];
+    size_t index[20];
+
+    /* Breakdown event must have an nhit greater than 1000. */
+    if (ev->nhit < 1000) return 0;
+
+    for (i = 0; i < MAX_PMTS; i++) {
+        if (!ev->pmt_hits[i].hit || pmts[i].pmt_type != PMT_NORMAL) continue;
+
+        tac[i/512][nhit[i/512]++] = ev->pmt_hits[i].tac;
+    }
+
+    for (i = 0; i < 19; i++) {
+        if (nhit[i] <= MIN_NHIT_CRATE) {
+            median_tac[i] = 0.0;
+            continue;
+        }
+        gsl_sort(&tac[i][0],1,nhit[i]);
+        median_tac[i] = gsl_stats_median_from_sorted_data(&tac[i][0],1,nhit[i]);
+    }
+
+    fargsort(median_tac,19,index);
+
+    return (nhit[index[18]] >= MIN_NHIT_BREAKDOWN && median_tac[index[18]] > median_tac[index[17]] + 500);
+}
+
 /* 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.