use QLX if QHS is railed

1 files changed, 26 insertions, 9 deletions

diff --git a/src/zdab_utils.c b/src/zdab_utils.c
index 6d05d96..e91d31b 100644
--- a/src/zdab_utils.c
+++ b/src/zdab_utils.c
@@ -140,14 +140,25 @@ int get_event(zebraFile *f, event *ev, zebraBank *bev)
             if (bpmt.qms)
                 ev->pmt_hits[id].flags |= PMT_FLAG_CHARGE;
 
-            /* Currently, the charge model only deals with QHS, so we flag any
-             * hits which have a bad or railed QHS value.
+            /* Determine the best charge to use in the likelihood function (QHS
+             * or QLX).
              *
-             * FIXME: In the future, it would be nice to use the best charge
-             * word (either QHS or QLX) depending on the if QHS is railed or
-             * not, but I need to do more work to see how the QLX values are
-             * normalized and if the existing charge model is good enough. */
-            if (is_mc(&ev_bank))
+             * Technically, the charge model only deals with QHS and so my
+             * original plan was to just flag any PMT hit with a railed or bad
+             * QHS and not include it in the likelihood function, but when
+             * testing I noticed that muons close to the PSUP can accidentally
+             * get reconstructed as electrons inside the AV if I ignore railed
+             * QHS values, so instead we use the QLX values when QHS is railed.
+             *
+             * Although the charge model only deals with QHS, I think it should
+             * be close enough that it's OK to use QLX. Both charges are
+             * normalized in the same way, the only difference would be in
+             * their SPE width.
+             *
+             * If both QHS and QLX are railed or have uncalibrated charges
+             * below 300, we flag the PMT and it's not included in the
+             * likelihood function. */
+            if (is_mc(&ev_bank)) {
                 /* Skip uncalibrated charge check if this is MC. The reason is
                  * that for some reason the uncalibrator in SNOMAN needs the
                  * QSLP tables which aren't available for most runs. According
@@ -163,9 +174,15 @@ int get_event(zebraFile *f, event *ev, zebraBank *bev)
                  * > doesn't work, don't try to use it".
                  *
                  * Therefore, we ignore it for MC. */
-                ;
-            else if(bpmt.pihs >= 4095 || bpmt.pihs < 300)
+                ev->pmt_hits[id].q = bpmt.phs;
+            } else if (bpmt.pihs < 4095 && bpmt.pihs >= 300) {
+                ev->pmt_hits[id].q = bpmt.phs;
+            } else if (bpmt.pilx < 4095 && bpmt.pilx >= 300) {
+                ev->pmt_hits[id].q = bpmt.plx;
+            } else {
+                ev->pmt_hits[id].q = 0.0;
                 ev->pmt_hits[id].flags |= PMT_FLAG_CHARGE;
+            }
 
             if (pmts[id].pmt_type != pmt_types[i]) {
                 get_pmt_type_string(pmts[id].pmt_type,pmt_type_string);