update muon cut

This commit updates the muon cut to require at least 1 OWL hit which has a high
charge and is early relative to nearby normal PMTs. This replaces the previous
cut criteria which required 2 OWL hits to be early relative to the 10%
percentile of all the normal PMT hits.

This new cut should target external muons better although I still need to do
some testing. In the future I'd like to optimize the distance at which PMTs are
considered nearby. Currently I just use 3 meters which seemed like a good value
based on some initial testing.

2 files changed, 42 insertions, 47 deletions

diff --git a/src/dc.c b/src/dc.c
index fdc89d2..87e6480 100644
--- a/src/dc.c
+++ b/src/dc.c
@@ -222,22 +222,25 @@ uint32_t get_dc_word(event *ev, zebraFile *f, zebraBank *bmast, zebraBank *bev)
  *
  * - calibrated nhit >= 150
  * - at least 5 calibrated OWL hits
- * - at least 1 OWL hit with a QHS >= 20.0
- * - at least 2 OWL hits with ECA times earlier than the 10th percentile of the
- *   PSUP PMT times
+ * - at least 1 OWL hit which has a high charge and is early relative to nearby
+ *   normal PMTs
  *
- * Idea: Should we also require at least one of the early tubes to have a high
- * charge? */
+ */
 int is_muon(event *ev)
 {
-    size_t i;
-    double ept[MAX_PMTS];
-    size_t nhit, nhit_owl, nhit_owl_early;
-    double early_time, max_owl_ehs;
+    size_t i, j;
+    double t_nearby[MAX_PMTS], q_nearby[MAX_PMTS];
+    int n_nearby;
+    size_t nhit_owl, nhit_owl_early;
+    double pmt_dir[3];
+    double distance;
+    double t_median, q_median;
+
+    /* Require at least 150 normal hits. */
+    if (ev->nhit < MUON_MIN_NHIT) return 0;
 
-    nhit = 0;
     nhit_owl = 0;
-    max_owl_ehs = 0.0;
+    nhit_owl_early = 0;
     for (i = 0; i < MAX_PMTS; i++) {
         /* Require good calibrations.
          *
@@ -251,47 +254,40 @@ int is_muon(event *ev)
          * in this case for some reason the ECA time still looks reasonable.
          * Therefore, since we only use the ECA time here we instead just make
          * sure that the ECA time is something reasonable. */
-        if (!ev->pmt_hits[i].hit || ev->pmt_hits[i].ept <= -100) continue;
+        if (!ev->pmt_hits[i].hit || pmts[i].pmt_type != PMT_OWL || ev->pmt_hits[i].ept <= -100) continue;
 
-        switch (pmts[i].pmt_type) {
-        case PMT_OWL:
-            if (nhit_owl == 0) {
-                max_owl_ehs = ev->pmt_hits[i].ehs;
-            } else if (ev->pmt_hits[i].ehs > max_owl_ehs) {
-                max_owl_ehs = ev->pmt_hits[i].ehs;
-            }
-            nhit_owl += 1;
-            break;
-        case PMT_NORMAL:
-            if (!ev->pmt_hits[i].flags) ept[nhit++] = ev->pmt_hits[i].ept;
-            break;
-        }
-    }
+        n_nearby = 0;
+        for (j = 0; j < MAX_PMTS; j++) {
+            if (!ev->pmt_hits[j].hit || pmts[j].pmt_type != PMT_NORMAL || ev->pmt_hits[j].ept <= -100) continue;
 
-    /* Require at least 150 normal hits. */
-    if (nhit < MUON_MIN_NHIT) return 0;
-
-    /* Require at least 5 OWL hits. */
-    if (nhit_owl < MUON_MIN_OWL_NHIT) return 0;
+            SUB(pmt_dir,pmts[i].pos,pmts[j].pos);
+            distance = NORM(pmt_dir);
 
-    /* Require at least 1 OWL tube to have a QHS >= 20.0. */
-    if (max_owl_ehs < MUON_OWL_QHS) return 0;
+            if (distance < MUON_OWL_NEARBY_DISTANCE) {
+                t_nearby[n_nearby] = ev->pmt_hits[j].ept;
+                q_nearby[n_nearby] = ev->pmt_hits[j].ehs;
+                n_nearby += 1;
+            }
+        }
 
-    /* Sort the times. */
-    gsl_sort(ept,1,nhit);
+        /* Not sure what to do if there are no nearby PMTs. */
+        if (n_nearby == 0) continue;
 
-    /* Calculate the time to define early OWL hits. */
-    early_time = gsl_stats_quantile_from_sorted_data(ept,1,nhit,MUON_OWL_TIME_PERCENTILE);
+        gsl_sort(t_nearby,1,n_nearby);
+        gsl_sort(q_nearby,1,n_nearby);
+        t_median = gsl_stats_median_from_sorted_data(t_nearby,1,n_nearby);
+        q_median = gsl_stats_median_from_sorted_data(q_nearby,1,n_nearby);
 
-    nhit_owl_early = 0;
-    for (i = 0; i < MAX_PMTS; i++) {
-        /* Require good calibrations. */
-        if (!ev->pmt_hits[i].hit || ev->pmt_hits[i].ept <= -100) continue;
+        if (ev->pmt_hits[i].ehs/q_median > 1.0 && ev->pmt_hits[i].ept < t_median)
+            nhit_owl_early += 1;
 
-        if (pmts[i].pmt_type == PMT_OWL && ev->pmt_hits[i].ept < early_time) nhit_owl_early += 1;
+        nhit_owl += 1;
     }
 
-    return nhit_owl_early >= 2;
+    /* Require at least 5 OWL hits. */
+    if (nhit_owl < MUON_MIN_OWL_NHIT) return 0;
+
+    return nhit_owl_early >= 1;
 }
 
 /* Returns 1 if the event fails the junk cut. The definition for the junk cut
diff --git a/src/dc.h b/src/dc.h
index c0c1791..e999235 100644
--- a/src/dc.h
+++ b/src/dc.h
@@ -71,10 +71,9 @@
 /* Minimum number of OWL PMTs which must be hit to be tagged as an incoming
  * muon. */
 #define MUON_MIN_OWL_NHIT 5
-/* Minimum value for the maximum OWL QHS value. */
-#define MUON_OWL_QHS 20.0
-/* Percentile of PSUP PMT times to define early OWL hits. */
-#define MUON_OWL_TIME_PERCENTILE 0.1
+/* Maximum distance between OWL and normal PMTs when computing the median time
+ * and charge. */
+#define MUON_OWL_NEARBY_DISTANCE 300.0
 
 /* QvNHIT ratio threshold. */
 #define QRATIO_THRESHOLD 0.25