Update S2Width Cut

lax/lichens/sciencerun0.py

@@ -433,37 +433,6 @@ def _process(self, df):
         return df
 
 
-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.
-
-    Current version is developed on unblinded Bkg data (paxv6.4.2). It is described in this note:
-    https://xecluster.lngs.infn.it/dokuwiki/doku.php?id=xenon:xenon1t:jacques:s1_single_scatter_cut
-
-    It should be applicable to data regardless of if it is ER or NR.
-
-    Contact: Jacques <jpienaa@purdue.edu>
-    """
-
-    version = 1
-
-    def _process(self, df):
-        s2width = S2Width
-
-        alt_rel_width = df['s2_range_50p_area'] / s2width.s2_width_model(df['alt_s1_interaction_z'])
-        alt_interaction_passes = alt_rel_width < s2width.relative_s2_width_bounds(df.s2.values, kind='high')
-        alt_interaction_passes &= alt_rel_width > s2width.relative_s2_width_bounds(df.s2.values, kind='low')
-
-        df.loc[:, (self.name())] = True ^ alt_interaction_passes
-        return df
-
-
 class S2AreaFractionTop(Lichen):
     """Cuts events with an unusual fraction of S2 on top array.
 
@@ -594,71 +563,76 @@ class S2Threshold(StringLichen):
     string = "200 < s2"
 
 
-class S2Width(ManyLichen):
+class S2Width(Lichen):
+    from scipy.stats import chi2
     """S2 Width cut based on diffusion model
-
     The S2 width cut compares the S2 width to what we could expect based on its depth in the detector. The inputs to
     this are the drift velocity and the diffusion constant. The allowed variation in S2 width is greater at low
     energy (since it is fluctuating statistically) Ref: (arXiv:1102.2865)
-
     It should be applicable to data regardless of if it ER or NR;
     above cS2 = 1e5 pe ERs the acceptance will go down due to track length effects.
+    around S2 = 1e5 pe there are beta-gamma merged peaks from Pb214 that extends the S2 width
+    Tune the diffusion model parameters based on fax data according to note:
+    https://xe1t-wiki.lngs.infn.it/doku.php?id=xenon:xenon1t:sim:notes:tzhu:width_cut_tuning#toy_fax_simulation
+    Contact: Tianyu <tz2263@columbia.edu>, Yuehuan <weiyh@physik.uzh.ch>, Jelle <jaalbers@nikhef.nl>
+    """
+    version = 4
 
-    Tune the diffusion model parameters based on pax v6.4.2 AmBe data according to note:
+    diffusion_constant = 25.26 * ((units.cm)**2) / units.s
+    v_drift = 1.440 * (units.um) / units.ns
+    scg = 23.0  # s2_secondary_sc_gain in pax config
+    scw = 258.41  # s2_secondary_sc_width median
+    SigmaToR50 = 1.349
 
-    xenon:xenon1t:yuehuan:analysis:0sciencerun_s2width_update0#comparison_with_diffusion_model_cut_by_jelle_pax_v642
+    def s2_width_model(self, z):
+        return np.sqrt(- 2 * self.diffusion_constant * z / self.v_drift ** 3)
 
-    Contact: Yuehuan <weiyh@physik.uzh.ch>, Jelle <jaalbers@nikhef.nl>
-    """
-    version = 2
+    def _process(self, df):
+        df.loc[:, 'nElectron'] = np.clip(df['s2'], 0, 5000) / self.scg
+        df.loc[:, 'normWidth'] = (np.square(df['s2_range_50p_area'] / self.SigmaToR50) - np.square(self.scw)) / \
+            np.square(self.s2_width_model(df['z']))
+        df.loc[:, self.name()] = chi2.logpdf(df['normWidth'] * (df['nElectron'] - 1), df['nElectron']) > - 14
+        return df
 
-    def __init__(self):
-        self.lichen_list = [self.S2WidthHigh(),
-                            self.S2WidthLow()]
-
-    @staticmethod
-    def s2_width_model(z):
-        diffusion_constant = 22.8 * ((units.cm)**2) / units.s
-        v_drift = 1.44 * (units.um) / units.ns
-        GausSigmaToR50 = 1.349
-
-        EffectivePar = 1.10795
-        Sigma_0 = 258.41 * units.ns
-        return GausSigmaToR50 * np.sqrt(Sigma_0 ** 2 - EffectivePar * 2 * diffusion_constant * z / v_drift ** 3)
-
-    def subpre(self, df):
-        # relative_s2_width
-        df.loc[:, 'temp'] = df['s2_range_50p_area'] / S2Width.s2_width_model(df['z'])
+    def post(self, df):
+        for temp_column in ['nElectron', 'normWidth']:
+            if temp_column in df.columns:
+                return df.drop(temp_column, 1)
         return df
 
-    @staticmethod
-    def relative_s2_width_bounds(s2, kind='high'):
-        x = 0.5 * np.log10(np.clip(s2, 150, 4500 if kind == 'high' else 2500))
-        if kind == 'high':
-            return 3 - x
-        elif kind == 'low':
-            return -0.9 + x
-        raise ValueError("kind must be high or low")
 
-    class S2WidthHigh(Lichen):
+class S1SingleScatter(Lichen):
+    from scipy.stats import chi2
+    """Requires only one valid interaction between the largest S2, and any S1 recorded before it.
 
-        def pre(self, df):
-            return S2Width.subpre(self, df)
+    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.
 
-        def _process(self, df):
-            df.loc[:, self.name()] = (df.temp <= S2Width.relative_s2_width_bounds(df.s2,
-                                                                                  kind='high'))
-            return df
+    Current version is developed on unblinded Bkg data (paxv6.4.2). It is described in this note:
+    https://xecluster.lngs.infn.it/dokuwiki/doku.php?id=xenon:xenon1t:jacques:s1_single_scatter_cut
 
-    class S2WidthLow(RangeLichen):
+    It should be applicable to data regardless of if it is ER or NR.
 
-        def pre(self, df):
-            return S2Width.subpre(self, df)
+    Contact: Jacques <jpienaa@purdue.edu>
+    """
 
-        def _process(self, df):
-            df.loc[:, self.name()] = (S2Width.relative_s2_width_bounds(df.s2,
-                                                                       kind='low') <= df.temp)
-            return df
+    version = 2
+    s2width = S2Width
+
+    def _process(self, df):
+
+        alt_n_electron = np.clip(df['s2'], 0, 5000) / self.s2width.scg
+        alt_rel_width = (np.square(df['s2_range_50p_area'] / self.s2width.SigmaToR50) - np.square(self.s2width.scw)) / \
+            np.square(self.s2width.s2_width_model(self.s2width, df['alt_s1_interaction_z']))
+
+        alt_interaction_passes = chi2.logpdf(alt_rel_width * (alt_n_electron - 1), alt_n_electron) > - 14
+
+        df.loc[:, (self.name())] = True ^ alt_interaction_passes
+        return df
 
 
 class S1AreaFractionTop(RangeLichen):

lax/lichens/sciencerun1.py

@@ -154,8 +154,6 @@ def pre(self, df):
 
 S1PatternLikelihood = sciencerun0.S1PatternLikelihood
 
-S1SingleScatter = sciencerun0.S1SingleScatter
-
 S2AreaFractionTop = sciencerun0.S2AreaFractionTop
 
 S2SingleScatter = sciencerun0.S2SingleScatter
@@ -177,66 +175,19 @@ class S2PatternLikelihood(StringLichen):
 
 S2Threshold = sciencerun0.S2Threshold
 
-
-class S2Width(ManyLichen):
-    """S2 Width cut based on diffusion model
-    The S2 width cut compares the S2 width to what we could expect based on its depth in the detector. The inputs to
-    this are the drift velocity and the diffusion constant. The allowed variation in S2 width is greater at low
-    energy (since it is fluctuating statistically) Ref: (arXiv:1102.2865)
-    It should be applicable to data regardless of if it ER or NR;
-    above cS2 = 1e5 pe ERs the acceptance will go down due to track length effects.
-    Tune the diffusion model parameters based on pax v6.4.2 AmBe data according to note:
-    xenon:xenon1t:yuehuan:analysis:0sciencerun_s2width_update0#comparison_with_diffusion_model_cut_by_jelle_pax_v642
-    Contact: Yuehuan <weiyh@physik.uzh.ch>, Jelle <jaalbers@nikhef.nl>
-    """
-    version = 3
-
-    def __init__(self):
-        self.lichen_list = [self.S2WidthHigh(),
-                            self.S2WidthLow()]
-
-    @staticmethod
-    def s2_width_model(z):
-        diffusion_constant = 31.73 * ((units.cm)**2) / units.s
-        v_drift = 1.335 * (units.um) / units.ns
-        GausSigmaToR50 = 1.349
-
-        EffectivePar = 0.925
-        Sigma_0 = 229.58 * units.ns
-        return GausSigmaToR50 * np.sqrt(Sigma_0 ** 2 - EffectivePar * 2 * diffusion_constant * z / v_drift ** 3)
-
-    @staticmethod
-    def subpre(self, df):
-        df.loc[:, 'temp'] = df['s2_range_50p_area'] / S2Width.s2_width_model(df['z'])
-        return df
-
-    @staticmethod
-    def relative_s2_width_bounds(s2, kind='high'):
-        x = 0.5 * np.log10(np.clip(s2, 150, 4500 if kind == 'high' else 2500))
-        if kind == 'high':
-            return 3 - x
-        elif kind == 'low':
-            return -0.9 + x
-        raise ValueError("kind must be high or low")
-
-    class S2WidthHigh(Lichen):
-
-        def pre(self, df):
-            return S2Width.subpre(self, df)
-
-        def _process(self, df):
-            df.loc[:, self.name()] = (df.temp <= S2Width.relative_s2_width_bounds(df.s2, kind='high'))
-            return df
-
-    class S2WidthLow(RangeLichen):
-
-        def pre(self, df):
-            return S2Width.subpre(self, df)
-
-        def _process(self, df):
-            df.loc[:, self.name()] = (S2Width.relative_s2_width_bounds(df.s2, kind='low') <= df.temp)
-            return df
-
+class S2Width(sciencerun0.S2Width):
+    from scipy.stats import chi2
+    version = 4
+    diffusion_constant = 29.35 * ((units.cm)**2) / units.s
+    v_drift = 1.335 * (units.um) / units.ns
+    scg = 21.3 # s2_secondary_sc_gain in pax config
+    scw = 229.58  # s2_secondary_sc_width median
+    SigmaToR50 = 1.349
+
+class S1SingleScatter(sciencerun0.S1SingleScatter):
+    from scipy.stats import chi2
+    version = 2
+    s2width = S2Width
 
 S1AreaFractionTop = sciencerun0.S1AreaFractionTop