plotting functions for d2 and scattering cross section

pysp2/util/normalized_derivative_method.py

@@ -1438,4 +1438,113 @@ def plot_scattering_cross_section(
         pad=14,
     )
 
+    return fig
+
+def plot_d2(
+    d2: Union[xr.DataArray, np.ndarray],
+    *,
+    threshold: Optional[float] = None,
+    title: Optional[str] = None,
+    xlabel: str = "k",
+    ylabel: str = r"Statistical distance $d^2$",
+    log_scale: bool = True,
+    highlight_kbest: bool = True,
+) -> plt.Figure:
+    """
+    Plot the Moteki & Kondo statistical distance d^2(k) for one event.
+
+    Parameters
+    ----------
+    d2 : xr.DataArray or np.ndarray
+        Computed d^2 values, typically output from compute_d2_moteki_kondo(...).
+        Expected shape is 1D over k.
+    threshold : float, optional
+        Optional horizontal acceptance threshold line (for example, a chi-square
+        cutoff or your empirical threshold).
+    title : str, optional
+        Custom plot title. If omitted, a default title is used.
+    xlabel : str, default "k"
+        x-axis label.
+    ylabel : str, default r"Statistical distance $d^2$"
+        y-axis label.
+    log_scale : bool, default True
+        If True, use a logarithmic y-axis.
+    highlight_kbest : bool, default True
+        If True, mark the minimum finite d^2 value with a red vertical line.
+
+    Returns
+    -------
+    matplotlib.figure.Figure
+        The created figure.
+    """
+    d2_np = np.asarray(d2.data if isinstance(d2, xr.DataArray) else d2, dtype=float)
+
+    if d2_np.ndim != 1:
+        raise ValueError("d2 must be a 1D array or DataArray.")
+
+    k = (
+         np.asarray(d2["k"].values)
+         if isinstance(d2, xr.DataArray) and "k" in d2.coords
+         else np.arange(d2_np.size)
+     )
+
+    finite = np.isfinite(d2_np)
+    if not np.any(finite):
+        raise ValueError("No finite d2 values available to plot.")
+
+    kbest = int(np.nanargmin(np.where(finite, d2_np, np.nan)))
+
+    plt.rcParams["font.family"] = "Times New Roman"
+    plt.rcParams["mathtext.fontset"] = "stix"
+
+    fig, ax = plt.subplots(figsize=(10, 5))
+
+    ax.plot(
+        k,
+        d2_np,
+        color="blue",
+        marker="o",
+        markersize=4,
+        linewidth=1.2,
+        label=r"$d^2(k)$",
+    )
+
+    if highlight_kbest:
+        ax.axvline(
+            k[kbest],
+            color="red",
+            linestyle="--",
+            linewidth=1.5,
+            label=rf"$k_{{\mathrm{{best}}}}={k[kbest]}$",
+        )
+        ax.scatter(
+            [k[kbest]],
+            [d2_np[kbest]],
+            color="red",
+            s=40,
+            zorder=5,
+        )
+
+    if threshold is not None:
+        ax.axhline(
+            threshold,
+            color="gray",
+            linestyle=":",
+            linewidth=1.5,
+            label=rf"Threshold = {threshold:g}",
+        )
+
+    if log_scale:
+        ax.set_yscale("log")
+
+    ax.set_xlabel(xlabel)
+    ax.set_ylabel(ylabel)
+    ax.grid(True, alpha=0.3, which="both" if log_scale else "major")
+
+    if title is None:
+        title = r"Moteki & Kondo Statistical Distance $d^2(k)$"
+    ax.set_title(title, pad=14)
+
+    ax.legend(loc="best", fontsize=10)
+
     return fig

tests/test_ndm.py

@@ -33,8 +33,8 @@ def test_ndm_moteki_kondo():
 
     cfg = MLEConfig(
     h=0.4,           # example: 0.4 microseconds
-    sigma_bar= 18.5*0.4,  # example; use your measured average width
-    delta_sigma=1.2*0.4, # example; use your measured width std dev
+    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,

tests/test_vis.py

@@ -10,6 +10,7 @@
 from pysp2.util.normalized_derivative_method import compute_sigma_moteki_kondo
 from pysp2.util.normalized_derivative_method import plot_incident_irradiance
 from pysp2.util.normalized_derivative_method import plot_scattering_cross_section
+from pysp2.util.normalized_derivative_method import plot_d2
 from pysp2.vis.plot_wave import plot_wave
 
 matplotlib.use("Agg")
@@ -48,8 +49,8 @@ def test_plot_incident_irradiance():
 
     cfg = MLEConfig(
     h=0.4,           # example: 0.4 microseconds
-    sigma_bar= 18.5*0.4,  # example; use your measured average width
-    delta_sigma=1.2*0.4, # example; use your measured width std dev
+    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,
@@ -114,8 +115,8 @@ def test_plot_scattering_cross_section():
     # print("dSdt dimensions:", dSdt.dims)
     cfg = MLEConfig(
     h=0.4,           # example: 0.4 microseconds
-    sigma_bar= 18.5*0.4,  # example; use your measured average width
-    delta_sigma=1.2*0.4, # example; use your measured width std dev
+    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,
@@ -168,4 +169,55 @@ def test_plot_scattering_cross_section():
         time_units="us"
         )
 
+    return fig
+
+@pytest.mark.mpl_image_compare(tolerance=10)
+def test_plot_d2():
+
+    my_binary = pysp2.util.gaussian_fit(my_sp2b, my_ini, parallel=False, baseline_to_zero=True)
+    dSdt = pysp2.util.central_difference(my_binary, normalize=True, baseline_to_zero=True)
+
+    cfg = MLEConfig(
+    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,
+)
+
+    tau = mle_tau_moteki_kondo(
+        S=my_binary,
+        norm_deriv=dSdt,
+        p=13,
+        ch="Data_ch0",
+        event_index=event,
+        min_start=15,
+        width_metric="fwtm",
+        config=cfg,
+    )
+
+    d2 = compute_d2_moteki_kondo(
+        S=my_binary,
+        norm_deriv=dSdt,
+        tau_hat=tau,
+        p=13,
+        ch="Data_ch0",
+        event_index=event,
+        min_start=15,
+        width_metric="fwtm",
+        config=cfg,
+    )
+
+    # Test the plotting function for channel 0 
+    fig = plot_d2(
+        d2=d2,
+        threshold=None,
+        title=None,
+        xlabel="k",
+        ylabel=r"Statistical distance $d^2$",
+        log_scale=True,
+        highlight_kbest=True,
+    )
+
     return fig