add best_uncal_q to the pmt_hit struct and use it in the muon cut

This commit adds a field to the pmt_hit struct called best_uncal_q which
represents the best ECA calibrated charge (in units of QHS counts above
pedestals). This is then used in the muon data cleaning cut.

1 files changed, 30 insertions, 0 deletions

diff --git a/src/zdab_utils.c b/src/zdab_utils.c
index 9238d9d..c4b04ea 100644
--- a/src/zdab_utils.c
+++ b/src/zdab_utils.c
@@ -42,6 +42,13 @@ static double TAC_PER_NS = 8.9;
 /* Average y offset of the TAC slope. */
 static double TAC_OFFSET = 300.0;
 
+/* Average gain between QLX and QHS.
+ *
+ * Note: This is not an accurate number it's just something I got from the SNO
+ * document "SNO Electronic Calibration Constants" and confirmed by looking at
+ * a plot of QHS vs QLX. */
+static double QLX_TO_QHS = 12.0;
+
 size_t get_nhit(event *ev)
 {
     /* Returns the number of PMT hits in event `ev`.
@@ -72,6 +79,26 @@ int is_mc(EVBank *ev)
     return ev->dtp/10 == 2;
 }
 
+/* Returns the "best" ECA calibrated charge (in units of EHS).
+ *
+ * This is sort of modelled after the SNOMAN routine cal_best_q() but here we
+ * do it for uncalibrated charges (in counts above pedestal). One issue here is
+ * that QHS and QLX have different gains, so to deal with that we just assume a
+ * constant gain and multiply any QLX values by QLX_TO_QHS.
+ *
+ * This best uncalibrated charge is mainly useful for data cleaning cuts. */
+double cal_best_q(float pihs, float pilx, float ehs, float elx)
+{
+    if (pilx < 300 || pilx > 4000 || pihs < 300 || pihs > 4000) {
+        /* QHS or QLX is railed, so use QLX. */
+        if (pilx < 300)
+            return 4095.0*QLX_TO_QHS;
+        return elx*QLX_TO_QHS;
+    }
+
+    return ehs;
+}
+
 int get_event(zebraFile *f, event *ev, zebraBank *bev)
 {
     /* Read all the PMT banks from the zebra file and update `ev`.
@@ -136,6 +163,7 @@ int get_event(zebraFile *f, event *ev, zebraBank *bev)
             ev->pmt_hits[id].ehl = bpmt.ehl;
             ev->pmt_hits[id].ehs = bpmt.ehs;
             ev->pmt_hits[id].elx = bpmt.elx;
+            ev->pmt_hits[id].best_uncal_q = cal_best_q(bpmt.pihs, bpmt.pilx, bpmt.ehs, bpmt.elx);
             ev->pmt_hits[id].qhl = bpmt.phl;
             ev->pmt_hits[id].qhs = bpmt.phs;
             ev->pmt_hits[id].qlx = bpmt.plx;
@@ -237,6 +265,8 @@ int get_event(zebraFile *f, event *ev, zebraBank *bev)
                 ev->pmt_hits[id].ehl = ev->pmt_hits[id].qihl - MEAN_PEDESTAL;
                 ev->pmt_hits[id].ehs = ev->pmt_hits[id].qihs - MEAN_PEDESTAL;
                 ev->pmt_hits[id].elx = ev->pmt_hits[id].qilx - MEAN_PEDESTAL;
+
+                ev->pmt_hits[id].best_uncal_q = cal_best_q(ev->pmt_hits[id].qihs, ev->pmt_hits[id].qilx, ev->pmt_hits[id].ehs, ev->pmt_hits[id].elx);
             } else if (bpmt.pihs < 4095 && bpmt.pihs >= 300) {
                 ev->pmt_hits[id].q = bpmt.phs;
             } else if (bpmt.pilx < 4095 && bpmt.pilx >= 300) {