Pulse: snr plot: indicate masking

simulations/11_pulsed_timing.py

@@ -5,7 +5,7 @@ from lib import util
 from scipy import signal, interpolate, stats
 import matplotlib.pyplot as plt
 import numpy as np
-from itertools import zip_longest
+from itertools import zip_longest, pairwise
 import h5py
 from copy import deepcopy
 
@@ -232,6 +232,7 @@ if __name__ == "__main__":
     h5_cache_fname = f'11_pulsed_timing.hdf5'
 
     time_accuracies = np.zeros(len(snr_factors))
+    mask_counts = np.zeros(len(snr_factors))
     for k, snr_sigma_factor in tqdm(enumerate(snr_factors)):
         # Read in cached time residuals
         if True:
@@ -448,65 +449,90 @@ if __name__ == "__main__":
 
         # Make a plot of the time residuals
         if N_residuals > 1:
-            time_accuracies[k] = np.std(time_residuals[:N_residuals])
+            time_residuals = time_residuals[:N_residuals]
 
-            hist_kwargs = dict(bins='sqrt', density=False, alpha=0.8, histtype='step')
+            for i in range(1 + cut_wrong_peak_matches):
-            fig, ax = plt.subplots()
+                mask_count = 0
-            ax.set_title(
-                    "Template Correlation Lag finding"
-                    + f"\n template dt: {template_dt*1e3: .1e}ps"
-                    + f"; antenna dt: {antenna_dt: .1e}ns"
-                    + f"; noise_factor: {noise_sigma_factor: .1e}"
-                    )
-            ax.set_xlabel("Time Residual [ns]")
-            ax.set_ylabel("#")
 
-            counts, bins, _patches = ax.hist(time_residuals, **hist_kwargs)
+                if i==1: # if cut_wrong_peak_matches:
-            if True: # fit gaussian to histogram
+                    wrong_peak_condition = lambda t_res: abs(t_res) > antenna_dt*4
-                min_x = min(time_residuals)
-                max_x = max(time_residuals)
 
-                dx = bins[1] - bins[0]
+                    mask = wrong_peak_condition(time_residuals)
-                scale = len(time_residuals) * dx
 
-                xs = np.linspace(min_x, max_x)
+                    mask_count = np.count_nonzero(mask)
 
-                # do the fit
+                    print("Masking {} residuals".format(mask_count))
-                name = "Norm"
+                    time_residuals = time_residuals[~mask]
-                param_names = [ "$\\mu$", "$\\sigma$" ]
-                distr_func = stats.norm
 
-                label = name +"(" + ','.join(param_names) + ')'
+                    if not mask_count:
+                        print("Continuing")
+                        continue
 
-                # plot
+                time_accuracies[k] = np.std(time_residuals)
-                fit_params = distr_func.fit(time_residuals)
+                mask_counts[k] = mask_count
-                fit_ys = scale * distr_func.pdf(xs, *fit_params)
+
-                ax.plot(xs, fit_ys, label=label)
+                hist_kwargs = dict(bins='sqrt', density=False, alpha=0.8, histtype='step')
+                fig, ax = plt.subplots()
+                ax.set_title(
+                        "Template Correlation Lag finding"
+                        + f"\n template dt: {template_dt: .1e}ns"
+                        + f"; antenna dt: {antenna_dt: .1e}ns"
+                        + ";" if not mask_count else "\n"
+                        + f"snr_factor: {snr_sigma_factor: .1e}"
+                        + "" if not mask_count else f"; N_masked: {mask_count}"
+                        )
+                ax.set_xlabel("Time Residual [ns]")
+                ax.set_ylabel("#")
+
+                counts, bins, _patches = ax.hist(time_residuals, **hist_kwargs)
+                if True: # fit gaussian to histogram
+                    min_x = min(time_residuals)
+                    max_x = max(time_residuals)
+
+                    dx = bins[1] - bins[0]
+                    scale = len(time_residuals) * dx
+
+                    xs = np.linspace(min_x, max_x)
+
+                    # do the fit
+                    name = "Norm"
+                    param_names = [ "$\\mu$", "$\\sigma$" ]
+                    distr_func = stats.norm
+
+                    label = name +"(" + ','.join(param_names) + ')'
+
+                    # plot
+                    fit_params = distr_func.fit(time_residuals)
+                    fit_ys = scale * distr_func.pdf(xs, *fit_params)
+                    ax.plot(xs, fit_ys, label=label)
+
+                    # chisq
+                    ct = np.diff(distr_func.cdf(bins, *fit_params))*np.sum(counts)
+                    if True:
+                        ct *= np.sum(counts)/np.sum(ct)
+                    c2t = stats.chisquare(counts, ct, ddof=len(fit_params))
+                    chisq_strs = [
+                                f"$\\chi^2$/dof = {c2t[0]: .2g}/{len(fit_params)}"
+                                ]
+
+                    # text on plot
+                    text_str = "\n".join(
+                        [label]
+                        +
+                        [ f"{param} = {value: .2e}" for param, value in zip_longest(param_names, fit_params, fillvalue='?') ]
+                        +
+                        chisq_strs
+                        )
+
+                    ax.text( *(0.02, 0.95), text_str, fontsize=12, ha='left', va='top', transform=ax.transAxes)
+
+                if mask_count:
+                    fig.savefig(f"figures/11_time_residual_hist_tdt{template_dt:0.1e}_n{snr_sigma_factor:.1e}_masked.pdf")
+                else:
+                    fig.savefig(f"figures/11_time_residual_hist_tdt{template_dt:0.1e}_n{snr_sigma_factor:.1e}.pdf")
 
-                # chisq
-                ct = np.diff(distr_func.cdf(bins, *fit_params))*np.sum(counts)
                 if True:
-                    ct *= np.sum(counts)/np.sum(ct)
+                    plt.close(fig)
-                c2t = stats.chisquare(counts, ct, ddof=len(fit_params))
-                chisq_strs = [
-                            f"$\\chi^2$/dof = {c2t[0]: .2g}/{len(fit_params)}"
-                            ]
-
-                # text on plot
-                text_str = "\n".join(
-                    [label]
-                    +
-                    [ f"{param} = {value: .2e}" for param, value in zip_longest(param_names, fit_params, fillvalue='?') ]
-                    +
-                    chisq_strs
-                    )
-
-                ax.text( *(0.02, 0.95), text_str, fontsize=12, ha='left', va='top', transform=ax.transAxes)
-
-            fig.savefig(f"figures/11_time_residual_hist_tdt{template_dt:0.1e}_n{noise_sigma_factor: .1e}.pdf")
-
-            if True:
-                plt.close(fig)
 
     # SNR time accuracy plot
     if True:
@@ -526,7 +552,17 @@ if __name__ == "__main__":
             ax.set_yscale('log')
 
         # plot the values
-        ax.plot(np.asarray(snr_factors), time_accuracies, ls='none', marker='o')
+        l = None
+        for j, mask_threshold in enumerate(pairwise([np.inf, 250, 50, 1, 0])):
+            kwargs = dict(
+                    ls='none',
+                    marker=['^', 'v','8', 'o',][j],
+                    color=None if l is None else l[0].get_color(),
+                    )
+            mask = mask_counts >= mask_threshold[1]
+            mask &= mask_counts < mask_threshold[0]
+
+            l = ax.plot(snr_factors[mask], time_accuracies[mask], **kwargs)
 
         if True: # limit y-axis to 1e0
             ax.set_ylim([None, 1e1])