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

4 files changed, 69 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.
diff --git a/src/dc.h b/src/dc.h
index fdd7b0b..c0c1791 100644
--- a/src/dc.h
+++ b/src/dc.h
@@ -33,6 +33,14 @@
 #define DC_OWL_TRIGGER    0x100
 #define DC_FTS            0x200
 #define DC_ITC            0x400
+#define DC_BREAKDOWN      0x800
+
+/* Minimum number of hits required for a crate to be considered the source of a
+ * breakdown. */
+#define MIN_NHIT_BREAKDOWN 256
+/* Minimum number of hits for a crate to calculate the median TAC in the
+ * breakdown cut. */
+#define MIN_NHIT_CRATE 20
 
 /* Length of the sliding window used in the ITC cut (ns). */
 #define ITC_TIME_WINDOW 93.0
@@ -83,6 +91,7 @@ int crate_isotropy(event *ev);
 int qvnhit(event *ev);
 double get_neck_tube_cut_time_diff(event *ev);
 int is_neck_event(event *ev);
+int is_breakdown(event *ev);
 int is_flasher(event *ev);
 
 #endif
diff --git a/src/event.h b/src/event.h
index ddb8be3..c07d648 100644
--- a/src/event.h
+++ b/src/event.h
@@ -30,6 +30,8 @@
 typedef struct pmt_hit {
     /* Set to 1 if the PMT was hit. */
     int hit;
+    /* Uncalibrated time. */
+    float tac;
     /* ECA calibrated time (ns). */
     float ept;
     /* Time in nano-secs relative to event T0 (ns). */
diff --git a/src/zdab_utils.c b/src/zdab_utils.c
index acb2fd9..9238d9d 100644
--- a/src/zdab_utils.c
+++ b/src/zdab_utils.c
@@ -37,6 +37,11 @@ static double MEAN_HIPT = 46.0;
  * to "uncalibrate" the Monte Carlo. */
 static double MEAN_PEDESTAL = 600.0;
 
+/* Average value for the slope of the TAC slope. */
+static double TAC_PER_NS = 8.9;
+/* Average y offset of the TAC slope. */
+static double TAC_OFFSET = 300.0;
+
 size_t get_nhit(event *ev)
 {
     /* Returns the number of PMT hits in event `ev`.
@@ -123,6 +128,7 @@ int get_event(zebraFile *f, event *ev, zebraBank *bev)
             }
 
             ev->pmt_hits[id].hit = 1;
+            ev->pmt_hits[id].tac = bpmt.pit;
             ev->pmt_hits[id].ept = bpmt.ept;
             ev->pmt_hits[id].qihl = bpmt.pihl;
             ev->pmt_hits[id].qihs = bpmt.pihs;
@@ -142,6 +148,15 @@ int get_event(zebraFile *f, event *ev, zebraBank *bev)
                 /* Use standard time if this is MC since the multiphoton time
                  * isn't calculated for MC. */
                 ev->pmt_hits[id].t = bpmt.pt;
+
+                /* Do a *very* simple conversion between time and TAC by
+                 * assuming a linear relationship. Note: This ignores *lots* of
+                 * effects like the TAC curl, any non-linearities, PCA offsets,
+                 * charge walk, etc. */
+                ev->pmt_hits[id].tac = fmax(fmin(bpmt.pt*TAC_PER_NS + TAC_OFFSET,4095.0),0.0);
+
+                /* Assume the ECA calibrated time and the ECA + PCA no walk
+                 * calibrated time are just equal to the calibrated time. */
                 ev->pmt_hits[id].ept = bpmt.pt;
                 ev->pmt_hits[id].pt1 = bpmt.pt;
             } else {