s1singlescatter update

lax/lichens/postsr1.py

@@ -173,3 +173,57 @@ def _process(self, df):
         df.loc[:, self.name()] = (np.nan_to_num(df.largest_s2_before_main_s2_area) < self.cutline(df.cs1)) | \
                                  (df.cs1 < self.min_s1)
         return df
+
+
+class S1SingleScatter(Lichen):
+    """
+    Added more requirement for largest_other_s1 to be identified as problematic. We found many largest_other_s1 are
+    from AP after S1, thus we required the largest hit area in one PMT shouldn't not exceed certain amount of the
+    largest_other_s1. The function is a bit different from S1 MAX PMT, as the current one seems not strict at
+    very low energy (~5 PE S1).
+
+    Requires Extended minitrees.
+
+    The cut is investigated for both SR1 and SR2 here:
+    https://xe1t-wiki.lngs.infn.it/doku.php?id=xenon:xenon1t:low_energy_er:s1singlescattercut
+
+    It should be applicable to data regardless whether it is ER or NR.
+    Contacts:
+    Jacques Pienaar, <[email protected]>
+    Joran Angevaare, <[email protected]>
+    Jingqiang Ye, <[email protected]>
+    """
+
+    version = 6
+    s2width = S2Width
+    alt_s1_coincidence_threshold = 3
+
+    @staticmethod
+    def largest_area_threshold(s1):
+        """
+        threshold of largest hit area in a single PMT for largest_other_s1, similar to S1 MAX PMT cut.
+        """
+        return np.minimum(0.052 * s1 + 4.15, 0.6 * s1 - 0.5)
+
+    def _process(self, df):
+        df.loc[:, self.name()] = True  # Default is True
+        mask = (df.alt_s1_interaction_drift_time > self.s2width.DriftTimeFromGate) & (
+                df.alt_s1_tight_coincidence >= self.alt_s1_coincidence_threshold)
+
+        # S2 width cut for alternate S1 - main S2 interaction
+        alt_n_electron = np.clip(df.loc[mask, 's2'], 0, 5000) / self.s2width.scg
+
+        alt_rel_width = np.square(df.loc[mask,
+                                         's2_range_50p_area'] / self.s2width.SigmaToR50) - np.square(self.s2width.scw)
+        alt_rel_width /= np.square(self.s2width.s2_width_model(self.s2width,
+                                                               df.loc[mask, 'alt_s1_interaction_drift_time']))
+
+        alt_interaction_passes = stats.chi2.logpdf(
+            alt_rel_width * (alt_n_electron - 1), alt_n_electron) > - 20
+
+        alt_pmt_passes = df.loc[mask, 'alt_s1_largest_hit_area'] < \
+                         self.largest_area_threshold(df.loc[mask, 'largest_other_s1'])
+
+        df.loc[mask, (self.name())] = True ^ (alt_interaction_passes & alt_pmt_passes)
+
+        return df

lax/lichens/sciencerun2.py

@@ -8,8 +8,6 @@
 from lax.lichens import postsr1
 DATA_DIR = sr1.DATA_DIR
 
-from scipy.stats import chi2
-
 ##
 # Combination cut packages
 ##
@@ -291,50 +289,8 @@ def _process(self, df):
         return df
 
 # S1 single scatter
-# Contact: Joran
-class S1SingleScatter(Lichen):
-    """Requires only one valid interaction between the largest S2, and any S1 recorded before it.
-
-    The S1 cut checks that any possible secondary S1s recorded in a waveform, could not have also
-    produced a valid interaction with the primary S2. To check whether an interaction between the
-    second largest S1 and the largest S2 is valid, we use the S2Width cut. If the event would pass
-    the S2Width cut, a valid second interaction exists, and we may have mis-identified which S1 to
-    pair with the primary S2. Therefore we cut this event. If it fails the S2Width cut the event is
-    not removed.
-
-    Requires Extended minitrees. 
-
-    The cut is investigated for SR2 here:
-    https://xe1t-wiki.lngs.infn.it/doku.php?id=xenon:xenon1t:analysis:sciencerun2:s1singlescattercut
-
-    It should be applicable to data regardless whether it is ER or NR.
-    Contacts: Jacques Pienaar, <[email protected]>
-    (for SR2) Joran Angevaare, <[email protected]>  
-    """
-
-    version = 5
-    s2width = S2Width
-    alt_s1_coincidence_threshold = 3
-
-    def _process(self, df):
-        df.loc[:, self.name()] = True  # Default is True
-        mask = (df.alt_s1_interaction_drift_time > self.s2width.DriftTimeFromGate) & (
-                df.alt_s1_tight_coincidence >= self.alt_s1_coincidence_threshold)      
-
-        # S2 width cut for alternate S1 - main S2 interaction
-        alt_n_electron = np.clip(df.loc[mask, 's2'], 0, 5000) / self.s2width.scg
-
-        alt_rel_width = np.square(df.loc[mask,
-                's2_range_50p_area'] / self.s2width.SigmaToR50) - np.square(self.s2width.scw)
-        alt_rel_width /= np.square(self.s2width.s2_width_model(self.s2width,
-                df.loc[mask, 'alt_s1_interaction_drift_time']))
-
-        alt_interaction_passes = chi2.logpdf(
-                alt_rel_width * (alt_n_electron - 1), alt_n_electron) > - 20
-
-        df.loc[mask, (self.name())] = True ^ alt_interaction_passes
-
-        return df
+# Contact: Joran, Jingqiang
+S1SingleScatter = postsr1.S1SingleScatter
 
 
 ##