add a function to tag flasher events

7 files changed, 321 insertions, 22 deletions

diff --git a/src/Makefile b/src/Makefile
index 4a366d7..214bf3f 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -1,7 +1,7 @@
 release_hdr := $(shell sh -c './mkreleasehdr.sh')
 
-CFLAGS=-fdiagnostics-color -O2 -Wall -g -DSWAP_BYTES
-LDLIBS=-fdiagnostics-color -lm -lgsl -lgslcblas -lnlopt_cxx -lstdc++
+CFLAGS=-fsanitize=address -fdiagnostics-color -O2 -Wall -g -DSWAP_BYTES
+LDLIBS=-fsanitize=address -fdiagnostics-color -lm -lgsl -lgslcblas -lnlopt_cxx -lstdc++
 
 all: test test-vector test-likelihood fit test-charge test-path calculate-csda-range test-time-pdf test-zebra test-find-peaks
 
@@ -28,7 +28,7 @@ test-time-pdf: test-time-pdf.o zebra.o likelihood.o pmt.o vector.o misc.o muon.o
 
 test-zebra: test-zebra.o zebra.o pack2b.o
 
-fit: fit.o zebra.o likelihood.o pmt.o vector.o misc.o muon.o optics.o quantum_efficiency.o solid_angle.o pdg.o scattering.o zdab_utils.o pack2b.o sno_charge.o db.o dqxx.o dict.o siphash.o path.o pmt_response.o release.o electron.o proton.o find_peaks.o quad.o
+fit: fit.o zebra.o likelihood.o pmt.o vector.o misc.o muon.o optics.o quantum_efficiency.o solid_angle.o pdg.o scattering.o zdab_utils.o pack2b.o sno_charge.o db.o dqxx.o dict.o siphash.o path.o pmt_response.o release.o electron.o proton.o find_peaks.o quad.o dc.o sort.o
 
 test-find-peaks: test-find-peaks.o zebra.o likelihood.o pmt.o vector.o misc.o muon.o optics.o quantum_efficiency.o solid_angle.o pdg.o scattering.o zdab_utils.o pack2b.o sno_charge.o db.o dqxx.o dict.o siphash.o path.o pmt_response.o release.o electron.o proton.o find_peaks.o
 
diff --git a/src/dc.c b/src/dc.c
new file mode 100644
index 0000000..897a645
--- /dev/null
+++ b/src/dc.c
@@ -0,0 +1,201 @@
+#include "event.h"
+#include "sort.h"
+#include <stddef.h> /* for size_t */
+#include "misc.h"
+#include "pmt.h"
+#include "vector.h"
+
+/* Returns 1 if the event is a flasher and 0 if it isn't. The definition of
+ * a flasher event comes from
+ * http://detector.sno.laurentian.ca/~detector/private/snoop/doc/eventIDs.html,
+ * with a few small changes.
+ *
+ * Here we define a flasher as an event with the following criteria:
+ *
+ * - 31 <= nhit <= 1000
+ *
+ * - Look for the paddle card with the most PMTs hit. This paddle card should
+ *   have 4 or more hits.
+ *
+ * - If such a paddle card exists, look for one channel in the card that hosts
+ *   this paddle card which has either QHS or QHL values (uncalib'd) > 4000,
+ *   whereas the QHS and QHL values for the other channels in the card are <
+ *   800. Note that QHS and QHL values below 300 are sent above 4096.
+ *
+ *   Note: This is changed from the original definition. We look for all the
+ *   other channels to have a QHL and QHS below 800 instead of 700 since the
+ *   latter was failing to tag many obvious flashers.
+ *
+ * - At least 70% of the (regular) pmts in the event which are not in the
+ *   flasher slot should have fired more than 50 ns after the high charge
+ *   channel.
+ *
+ * - At least 70% of the (regular) pmts in the event which are not in the
+ *   flasher slot should be at a distance >= 12.0 m from the high charge pmt.
+ *
+ * The one criteria we *don't* use that is on the webpage is the AMB cut:
+ *
+ * - AMB Int >= 200
+ *
+ * since I can't seem to figure out where that is stored in the data structure.
+ */
+int is_flasher(event *ev)
+{
+    int hits_pc[1280] = {0};
+    int qhs, qhl;
+    int qhl_pc[8] = {0};
+    int qhs_pc[8] = {0};
+    double t_pc[8] = {0};
+    size_t crate, card, channel, id;
+    int i, j;
+    size_t index[1280], index_qhs[8], index_qhl[8];
+    double t;
+    int flasher_pc;
+    double flasher_pos[3];
+    double pmt_dir[3];
+    int flasher = 0;
+    double distance;
+    int nhit, nhit_late, nhit_far;
+
+    /* Flasher event must have an nhit between 31 and 1000. */
+    if (ev->nhit < 31 || 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;
+
+	crate = i/512;
+	card = (i % 512)/32;
+	channel = i % 32;
+
+	id = crate*64 + card*4 + channel/8;
+
+	hits_pc[id] += 1;
+    }
+
+    argsort(hits_pc, LEN(hits_pc), index);
+
+    /* The paddle card with the most hits must have >= 4 hits. */
+    if (hits_pc[index[1279]] < 4) return 0;
+
+    /* Loop over the most hit paddle cards in order from most hits to least
+     * hits. We do this because it's possible that more than one paddle card
+     * has the same number of channels hit, and so we need to test all of them. */
+    for (j = 1279; j >= 0; j--) {
+	id = index[j];
+
+        /* If this paddle card doesn't have as many hits as the most hit paddle
+         * card and we haven't found a flasher yet, then it's not a flasher. */
+        if (hits_pc[id] < hits_pc[index[1279]]) return 0;
+
+	for (i = 0; i < MAX_PMTS; i++) {
+            if (!ev->pmt_hits[i].hit) continue;
+
+            crate = i/512;
+            card = (i % 512)/32;
+            channel = i % 32;
+
+	    id = crate*64 + card*4 + channel/8;
+
+	    if (id != index[j]) continue;
+
+            /* Get the uncalibrated QHS and QHL charges. */
+	    qhs = ev->pmt_hits[i].qihs;
+	    qhl = ev->pmt_hits[i].qihl;
+
+            if (qhs < 300) {
+                /* QHS values below 300 are set to 4095. */
+                qhs_pc[channel % 8] = 4095;
+            } else {
+                qhs_pc[channel % 8] = qhs;
+            }
+
+            if (qhl < 300) {
+                /* QHL values below 300 are set to 4095. */
+                qhl_pc[channel % 8] = 4095;
+            } else {
+                qhl_pc[channel % 8] = qhl;
+            }
+
+            t_pc[channel % 8] = ev->pmt_hits[i].t;
+	}
+
+        /* Sort the QHS and QHL values by size. */
+	argsort(qhs_pc, LEN(qhs_pc), index_qhs);
+	argsort(qhl_pc, LEN(qhl_pc), index_qhl);
+
+        /* Check if this paddle card has QHS or QHL values > 4000, and the QHS
+         * and QHL values for the other channels in the card are less than 700. */
+        if (qhs_pc[index_qhs[7]] > 4000 || qhl_pc[index_qhl[7]] > 4000) {
+	    if (qhs_pc[index_qhs[6]] < 800 && qhl_pc[index_qhl[6]] < 800) {
+                /* Flasher! */
+                flasher = 1;
+                break;
+	    }
+	}
+    }
+
+    /* If we didn't find a flasher, return 0. */
+    if (!flasher) return 0;
+
+    flasher_pc = index[j];
+    crate = index[j]/64;
+    card = (index[j] % 64)/4;
+    if (qhs_pc[index_qhs[7]] > 4000) {
+        channel = (index[j] % 4)*8 + index_qhs[7];
+        t = t_pc[index_qhs[7]];
+    } else {
+        channel = (index[j] % 4)*8 + index_qhl[7];
+        t = t_pc[index_qhl[7]];
+    }
+    COPY(flasher_pos,pmts[crate*512 + card*32 + channel].pos);
+
+    /* Check that 70% of the regular PMTs in the event which are not in the
+     * flasher slot fired more than 50 ns after the high charge channel and are
+     * more than 12 meters away from the high charge channel. */
+    nhit = 0;
+    /* Number of PMTs which fired 50 ns after the high charge channel. */
+    nhit_late = 0;
+    /* Number of PMTs which are 12 meters or more away from the high charge
+     * channel. */
+    nhit_far = 0;
+    for (i = 0; i < MAX_PMTS; i++) {
+        /* Skip PMTs which weren't hit. */
+        if (!ev->pmt_hits[i].hit) continue;
+
+        /* Only count regular PMTs without any flags. */
+        if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL) continue;
+
+	crate = i/512;
+	card = (i % 512)/32;
+	channel = i % 32;
+
+	id = crate*64 + card*4 + channel/8;
+
+        /* Skip PMTs in the same paddle card as the high charge channel. */
+        if (id == flasher_pc) continue;
+
+        nhit += 1;
+
+        /* Calculate the distance from the current channel to the high charge
+         * channel. */
+        SUB(pmt_dir,pmts[i].pos,flasher_pos);
+        distance = NORM(pmt_dir);
+
+        /* If this channel is 12 meters or more away, increment nhit_far. */
+        if (distance >= 1200.0) nhit_far += 1;
+
+        /* If this channel fired more than 50 ns after the high charge channel,
+         * increment nhit_late. */
+        if (ev->pmt_hits[i].t > t + 50.0) nhit_late += 1;
+    }
+
+    /* If less than 70% of the regular PMTs fired within 50 ns of the high
+     * charge channel, then it's not a flasher. */
+    if (nhit_late < nhit*0.7) return 0;
+
+    /* If less than 70% of the regular PMTs were closer than 12 meters away
+     * from the high charge channel, then it's not a flasher. */
+    if (nhit_far < nhit*0.7) return 0;
+
+    return 1;
+}
diff --git a/src/dc.h b/src/dc.h
new file mode 100644
index 0000000..34ddf09
--- /dev/null
+++ b/src/dc.h
@@ -0,0 +1,8 @@
+#ifndef DC_H
+#define DC_H
+
+#include "event.h"
+
+int is_flasher(event *ev);
+
+#endif
diff --git a/src/event.h b/src/event.h
index b530347..64e8df1 100644
--- a/src/event.h
+++ b/src/event.h
@@ -12,12 +12,16 @@ typedef struct pmt_hit {
     float qhl;
     float qhs;
     float qlx;
+    uint16_t qihl;
+    uint16_t qihs;
+    uint16_t qilx;
     int flags;
 } pmt_hit;
 
 typedef struct event {
     int run;
     uint32_t gtid;
+    int nhit;
     pmt_hit pmt_hits[MAX_PMTS];
 } event;
 
diff --git a/src/fit.c b/src/fit.c
index f3ae311..ef22d44 100644
--- a/src/fit.c
+++ b/src/fit.c
@@ -5776,6 +5776,25 @@ void sigint_handler(int dummy)
     stop = 1;
 }
 
+size_t get_nhit(event *ev)
+{
+    /* Returns the number of PMT hits in event `ev`.
+     *
+     * Note: Only hits on normal PMTs which aren't flagged are counted. */
+    size_t i, nhit;
+
+    nhit = 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) continue;
+
+        nhit++;
+    }
+
+    return nhit;
+}
+
 int get_event(zebraFile *f, event *ev, zebraBank *bev)
 {
     /* Read all the PMT banks from the zebra file and update `ev`.
@@ -5810,6 +5829,9 @@ int get_event(zebraFile *f, event *ev, zebraBank *bev)
         id = crate*512 + card*32 + channel;
         ev->pmt_hits[id].hit = 1;
         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].qhl = bpmt.phl;
         ev->pmt_hits[id].qhs = bpmt.phs;
         ev->pmt_hits[id].qlx = bpmt.plx;
@@ -5825,26 +5847,9 @@ int get_event(zebraFile *f, event *ev, zebraBank *bev)
         }
     }
 
-    return 0;
-}
-
-size_t get_nhit(event *ev)
-{
-    /* Returns the number of PMT hits in event `ev`.
-     *
-     * Note: Only hits on normal PMTs which aren't flagged are counted. */
-    size_t i, nhit;
+    ev->nhit = get_nhit(ev);
 
-    nhit = 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) continue;
-
-        nhit++;
-    }
-
-    return nhit;
+    return 0;
 }
 
 void sprintf_particle_string(int *id, size_t n, char *str)
diff --git a/src/sort.c b/src/sort.c
new file mode 100644
index 0000000..788c757
--- /dev/null
+++ b/src/sort.c
@@ -0,0 +1,72 @@
+#include <stdlib.h>
+#include "sort.h"
+
+typedef struct fargsort_element {
+    double value;
+    size_t index;
+} fargsort_element;
+
+typedef struct argsort_element {
+    int value;
+    size_t index;
+} argsort_element;
+
+int fargsort_compare(const void *a, const void *b)
+{
+    const double da = ((fargsort_element *) a)->value;
+    const double db = ((fargsort_element *) b)->value;
+
+    return (da > db) - (da < db);
+}
+
+int argsort_compare(const void *a, const void *b)
+{
+    const int da = ((argsort_element *) a)->value;
+    const int db = ((argsort_element *) b)->value;
+
+    return (da > db) - (da < db);
+}
+
+void fargsort(double *base, size_t num, size_t *index_array)
+{
+    /* Returns the indices that would sort the array `base`. */
+    size_t i;
+    fargsort_element *array;
+
+    array = (fargsort_element *) malloc(num*sizeof(fargsort_element));
+
+    for (i = 0; i < num; i++) {
+        array[i].value = base[i];
+        array[i].index = i;
+    }
+
+    qsort(array, num, sizeof(fargsort_element), fargsort_compare);
+
+    for (i = 0; i < num; i++) {
+        index_array[i] = array[i].index;
+    }
+
+    free(array);
+}
+
+void argsort(int *base, size_t num, size_t *index_array)
+{
+    /* Returns the indices that would sort the array `base`. */
+    size_t i;
+    argsort_element *array;
+
+    array = (argsort_element *) malloc(num*sizeof(argsort_element));
+
+    for (i = 0; i < num; i++) {
+        array[i].value = base[i];
+        array[i].index = i;
+    }
+
+    qsort(array, num, sizeof(argsort_element), argsort_compare);
+
+    for (i = 0; i < num; i++) {
+        index_array[i] = array[i].index;
+    }
+
+    free(array);
+}
diff --git a/src/sort.h b/src/sort.h
new file mode 100644
index 0000000..4c4b734
--- /dev/null
+++ b/src/sort.h
@@ -0,0 +1,9 @@
+#ifndef SORT_H
+#define SORT_H
+
+#include <stddef.h> /* for size_t */
+
+void fargsort(double *base, size_t num, size_t *index_array);
+void argsort(int *base, size_t num, size_t *index_array);
+
+#endif