diff --git a/hax/treemakers/posrec.py b/hax/treemakers/posrec.py
index 4263096..4dd2ec4 100644
--- a/hax/treemakers/posrec.py
+++ b/hax/treemakers/posrec.py
@@ -34,7 +34,7 @@ class PositionReconstruction(TreeMaker):
        - s1_area_upper_injection_fraction: s1 area fraction near Rn220 injection points (near PMT 131)
        - s1_area_lower_injection_fraction: s1 area fraction near Rn220 injection points (near PMT 243)
     """
-    __version__ = '0.16'
+    __version__ = '0.17'
     extra_branches = ['peaks.area_per_channel[260]',
                       'peaks.hits_per_channel[260]',
                       'peaks.n_saturated_per_channel[260]',
@@ -175,10 +175,10 @@ def extract_data(self, event):
                                                                             self.nn_tensorflow,
                                                                             self.loaded_tfnn_model_json_name)
 
-        e = event.event_number
+        event_num = event.event_number
 
         try:
-            i = self.indices.index(e)
+            event_index = self.indices.index(event_num)
         except Exception:
             return event_data
 
@@ -189,8 +189,8 @@ def extract_data(self, event):
         # Position reconstruction based on NN from TensorFlow
         s2apc = np.array(list(s2.area_per_channel))
         s2apc_clean = []
-        for i, s2_t in enumerate(s2apc):
-            if i not in self.list_bad_pmts and i < self.ntop_pmts:
+        for ipmt, s2_t in enumerate(s2apc):
+            if ipmt not in self.list_bad_pmts and ipmt < self.ntop_pmts:
                 s2apc_clean.append(s2_t)
         s2apc_clean = np.asarray(s2apc_clean)
         s2apc_clean_norm = s2apc_clean / s2apc_clean.sum()
@@ -218,7 +218,7 @@ def extract_data(self, event):
             size_top = s1.area * s1.area_fraction_top
             size_tot = s1.area
 
-        aft = self.s1_aft_map.get_value(self.x[i], self.y[i], self.z[i])
+        aft = self.s1_aft_map.get_value(self.x[event_index], self.y[event_index], self.z[event_index])
         event_data['s1_area_fraction_top_probability_hax'] = binom_test_pax(size_top, size_tot, aft)
 
         # Now do s1_pattern_fit
@@ -232,38 +232,40 @@ def extract_data(self, event):
                 confused_s1_channels.append(a)
 
         try:
+
+            # Create PMT array of booleans for use in likelihood calculation
+            is_pmt_in = np.ones(len(self.tpc_channels), dtype=bool)  # Default True
+            is_pmt_in[confused_s1_channels] = False  # Ignore saturated channels
+
             event_data['s1_pattern_fit_hax'] = self.pattern_fitter.compute_gof(
-                (self.x[i], self.y[i], self.z[i]),
+                (self.x[event_index], self.y[event_index], self.z[event_index]),
                 apc[self.tpc_channels],
-                pmt_selection=np.setdiff1d(self.tpc_channels,
-                                           confused_s1_channels),
+                pmt_selection=is_pmt_in,
                 statistic=self.statistic)
 
             event_data['s1_pattern_fit_hits_hax'] = self.pattern_fitter.compute_gof(
-                (self.x[i], self.y[i], self.z[i]),
+                (self.x[event_index], self.y[event_index], self.z[event_index]),
                 hpc[self.tpc_channels],
-                pmt_selection=np.setdiff1d(self.tpc_channels,
-                                           confused_s1_channels),
+                pmt_selection=is_pmt_in,
                 statistic=self.statistic)
 
-            pmts_bottom = np.setdiff1d(self.tpc_channels, confused_s1_channels)
-            pmts_bottom[self.pax_config['DEFAULT']['channels_top']] = 0
+            # Switch to bottom PMTs only
+            is_pmt_in[self.pax_config['DEFAULT']['channels_top']] = False
 
             event_data['s1_pattern_fit_bottom_hax'] = self.pattern_fitter.compute_gof(
-                (self.x[i], self.y[i], self.z[i]),
+                (self.x[event_index], self.y[event_index], self.z[event_index]),
                 apc[self.tpc_channels],
-                pmt_selection=pmts_bottom,
+                pmt_selection=is_pmt_in,
                 statistic=self.statistic)
 
             event_data['s1_pattern_fit_bottom_hits_hax'] = self.pattern_fitter.compute_gof(
-                (self.x[i], self.y[i], self.z[i]),
+                (self.x[event_index], self.y[event_index], self.z[event_index]),
                 hpc[self.tpc_channels],
-                pmt_selection=pmts_bottom,
+                pmt_selection=is_pmt_in,
                 statistic=self.statistic)
 
-        except exceptions.CoordinateOutOfRangeException as E:
+        except exceptions.CoordinateOutOfRangeException as _:
             # pax does this too. happens when event out of TPC (usually z)
-            log.error(E)
             return event_data
 
         return event_data