fixes to d2

pysp2/util/normalized_derivative_method.py

@@ -331,6 +331,10 @@ def _to_dataarray(
 
     raise TypeError(f"{name} must be an xarray DataArray or Dataset.")
 
+def _baseline_to_zero(s_event: np.ndarray) -> np.ndarray:
+    """Shift one event's signal so its minimum is 0, matching central_difference's baseline_to_zero."""
+    return s_event - np.nanmin(s_event)
+
 def _moteki_kondo_subset_statistics(
     yk: np.ndarray,
     sk: np.ndarray,
@@ -411,6 +415,7 @@ def mle_tau_moteki_kondo(
     tau_grid: Optional[Union[np.ndarray, xr.DataArray]] = None,
     min_start: int = 15,
     width_metric: str = "fwhm",
+    baseline_to_zero: bool = True,
     config: Optional[MLEConfig] = None,
 ) -> xr.DataArray:
     """
@@ -441,6 +446,8 @@ def mle_tau_moteki_kondo(
         Minimum allowed start index to exclude unusable early samples.
     width_metric : str
         "fwhm" or "fwtm" for defining peak width.
+    baseline_to_zero : bool
+        If True, shift the selected scattering signal so its minimum is zero.
     config : MLEConfig
         Calibration / noise / grid settings.
     """
@@ -565,6 +572,8 @@ def _tau_hat_for_one_event(s_event, y_event):
 
     # Extract event
     s_event = S_std.sel({event_dim: event_index}).values
+    if baseline_to_zero:
+        s_event = _baseline_to_zero(s_event)
     y_event = y_std.sel({event_dim: event_index}).values
 
     tau_hat_1d = _tau_hat_for_one_event(s_event, y_event)
@@ -591,6 +600,7 @@ def compute_d2_moteki_kondo(
     y_sample_dim: Optional[str] = None,
     min_start=15,
     width_metric="fwhm",
+    baseline_to_zero: bool = True,
     config: Optional[MLEConfig] = None,
 ) -> xr.DataArray:
     """
@@ -624,6 +634,8 @@ def compute_d2_moteki_kondo(
         Minimum allowed start index to exclude unusable early samples.
     width_metric : str
         "fwhm" or "fwtm" for defining peak width.
+    baseline_to_zero : bool
+        If True, shift the selected scattering signal so its minimum is zero.
     config : MLEConfig
         Calibration / noise / grid settings.
 
@@ -707,6 +719,8 @@ def compute_d2_moteki_kondo(
     t = np.arange(n_samples) * h
 
     s_event = S_std.sel({event_dim: event_index}).values
+    if baseline_to_zero:
+        s_event = _baseline_to_zero(s_event)
     y_event = y_std.sel({event_dim: event_index}).values
 
     d2_vals = np.full(k_values.size, np.nan)
@@ -751,7 +765,8 @@ def compute_sigma_moteki_kondo(
     y_sample_dim: Optional[str] = None,
     min_start: int = 15,
     width_metric: str = "fwhm",
-    d2_threshold: float = 80000.0,
+    d2_threshold: float = 20.0,
+    baseline_to_zero: bool = True,
     config: Optional[MLEConfig] = None,
 ) -> xr.Dataset:
     """
@@ -768,9 +783,14 @@ def compute_sigma_moteki_kondo(
 
     Notes
     -----
-    The paper applies the sigma estimate after requiring d²(kbest) < 200000.
-    For consistency, this function returns sigma_hat = NaN when the threshold
-    is not met, while still returning diagnostic fields.
+    The paper applies the sigma estimate after requiring d²(kbest) < 20. This is 
+    aligned with Moteki and Kondo (2008). For consistency, this function returns 
+    sigma_hat = NaN when the threshold is not met, while still returning diagnostic fields.
+
+    Parameters
+    ----------
+    baseline_to_zero : bool
+        If True, shift the selected scattering signal so its minimum is zero.
     """
     if config is None:
         raise ValueError("config must be provided.")
@@ -876,6 +896,8 @@ def compute_sigma_moteki_kondo(
 
     # Select the requested event.
     s_event = np.asarray(S_std.sel({event_dim: event_index}).values, dtype=float)
+    if baseline_to_zero:
+        s_event = _baseline_to_zero(s_event)
     y_event = np.asarray(y_std.sel({event_dim: event_index}).values, dtype=float)
 
     if not (np.all(np.isfinite(s_event)) and np.all(np.isfinite(y_event))):

tests/test_ndm.py

@@ -35,9 +35,12 @@ def test_ndm_moteki_kondo():
     h=0.4,           # example: 0.4 microseconds
     sigma_bar= (18.5/ 2.35482 )*0.4,  # example; use your measured average width
     delta_sigma=(1.2/ 2.35482 )*0.4, # example; use your measured width std dev
-    A1=0.37*2.44,
-    A2=(1.6e-2)*2.44**(1/2),
-    A3=6.2e-4,
+    # measured pre-trigger baseline noise sd in raw ADC counts 
+    # A1, A2, and A3 are the coefficients for the polynomial used in the NDM model
+    # A parameters optimized for the PSL dataset
+    A1=37.0,
+    A2=1.56205,
+    A3=0.0008,
 )
 
     ## Test one event ##################################################
@@ -55,7 +58,7 @@ def test_ndm_moteki_kondo():
     tau_val_true = my_binary['Data_ch0'].isel(event_index=event).argmax().item()*0.4
     # Test that the estimated tau for a subset of results is close to the true value for the event
     for i in range(10, 15):
-        np.testing.assert_allclose(tau[i], tau_val_true, atol=0.3)
+        np.testing.assert_allclose(tau[i], tau_val_true, atol=0.4)
 
     d2 = compute_d2_moteki_kondo(
         S=my_binary,
@@ -77,7 +80,7 @@ def test_ndm_moteki_kondo():
     np.testing.assert_allclose(
         tau_best,
         tau_val_true,
-        atol=0.3,  # absolute tolerance = 0.3 microseconds
+        atol=0.4,  # absolute tolerance = 0.4 microseconds
     )
 
     sigma_ds = compute_sigma_moteki_kondo(
@@ -100,7 +103,7 @@ def test_ndm_moteki_kondo():
     np.testing.assert_allclose(
          sigma_ds['sigma_hat'].values,
          sigma_best,
-         atol=0.12,  # absolute tolerance = 1.5 microseconds
+         atol=0.15,  # absolute tolerance = 0.15 microseconds
     )
 
     # Test the normalized irradiance function 
@@ -118,5 +121,5 @@ def test_ndm_moteki_kondo():
         np.testing.assert_allclose(
             (I_norm * np.nanmax(y_scatter_background_shifted))[i],
             y_scatter_background_shifted[i],
-            atol=4500,  # absolute tolerance ~ 10% of the max scattering signal value
+            atol=5000,  # absolute tolerance ~ 11% of the max scattering signal value
         )

tests/test_vis.py

@@ -51,9 +51,12 @@ def test_plot_incident_irradiance():
     h=0.4,           # example: 0.4 microseconds
     sigma_bar= (18.5/ 2.35482 )*0.4,  # example; use your measured average width
     delta_sigma=(1.2/ 2.35482 )*0.4, # example; use your measured width std dev
-    A1=0.37*2.44,
-    A2=(1.6e-2)*2.44**(1/2),
-    A3=6.2e-4,
+    # measured pre-trigger baseline noise sd in raw ADC counts 
+    # A1, A2, and A3 are the coefficients for the polynomial used in the NDM model
+    # A parameters optimized for this dataset in external code (available upon request)
+    A1=37.0,
+    A2=1.56205,
+    A3=0.0008,
 )
 
     tau = mle_tau_moteki_kondo(
@@ -117,9 +120,9 @@ def test_plot_scattering_cross_section():
     h=0.4,           # example: 0.4 microseconds
     sigma_bar= (18.5/ 2.35482 )*0.4,  # example; use your measured average width
     delta_sigma=(1.2/ 2.35482 )*0.4, # example; use your measured width std dev
-    A1=0.37*2.44,
-    A2=(1.6e-2)*2.44**(1/2),
-    A3=6.2e-4,
+    A1=37.0,
+    A2=1.56205,
+    A3=0.0008,
 )
 
     tau = mle_tau_moteki_kondo(
@@ -181,9 +184,9 @@ def test_plot_d2():
     h=0.4,           # example: 0.4 microseconds
     sigma_bar= (18.5/ 2.35482 )*0.4,  # example; use your measured average width
     delta_sigma=(1.2/ 2.35482 )*0.4, # example; use your measured width std dev
-    A1=0.37*2.44,
-    A2=(1.6e-2)*2.44**(1/2),
-    A3=6.2e-4,
+    A1=37.0,
+    A2=1.56205,
+    A3=0.0008,
 )
 
     tau = mle_tau_moteki_kondo(