Pulse: snr plot multiple template_dt curves

simulations/11_pulsed_timing.py

@@ -401,11 +401,14 @@ if __name__ == "__main__":
         matplotlib.use('Agg')
 
     bp_freq = (30e-3, 80e-3) # GHz
-    template_dt = 5e-2 # ns
     interp_template_dt = 5e-5 # ns
     template_length = 200 # ns
 
+    antenna_dt = 2 # ns
+    antenna_timelength = 1024 # ns
+
     N_residuals = 50*3 if len(sys.argv) < 2 else int(sys.argv[1])
+    template_dts = np.array([antenna_dt, 5e-1, 5e-2]) # ns
     snr_factors = np.concatenate( # 1/noise_amplitude factor
             (
                 #[0.25, 0.5, 0.75],
@@ -415,8 +418,6 @@ if __name__ == "__main__":
             ),
             axis=None, dtype=float)
 
-    antenna_dt = 2 # ns
-    antenna_timelength = 1024 # ns
 
     cut_wrong_peak_matches = True
     normalise_noise = False
@@ -454,112 +455,113 @@ if __name__ == "__main__":
         if True:
             plt.close(fig)
 
-    #
-    # Create the template
-    # This is sampled at a lower samplerate than the interpolation template
-    #
-    template, _ = create_template(dt=template_dt, timelength=template_length, bp_freq=bp_freq, name='Template')
 
     #
     # Find time accuracies as a function of signal strength
     #
-    time_accuracies = np.zeros(len(snr_factors))
+    time_accuracies = np.zeros((len(template_dts), len(snr_factors)))
-    mask_counts = np.zeros(len(snr_factors))
+    mask_counts = np.zeros_like(time_accuracies)
-    for k, snr_sigma_factor in tqdm(enumerate(snr_factors)):
+    for l, template_dt in tqdm(enumerate(template_dts)):
 
-        time_residuals = get_time_residuals_for_template(
+        # Create the template
-                N_residuals, template, interpolation_template=interp_template,
+        # This is sampled at a lower samplerate than the interpolation template
-                antenna_dt=antenna_dt, antenna_timelength=antenna_timelength,
+        template, _ = create_template(dt=template_dt, timelength=template_length, bp_freq=bp_freq, name='Template')
-                snr_sigma_factor=snr_sigma_factor, bp_freq=bp_freq, normalise_noise=normalise_noise,
-                h5_cache_fname=h5_cache_fname, rng=rng, tqdm=tqdm)
 
-        print()# separating tqdm
+        for k, snr_sigma_factor in tqdm(enumerate(snr_factors)):
-        print()# separating tqdm
 
-        # Make a plot of the time residuals
+            # get the time residuals
-        if N_residuals > 1:
+            time_residuals = get_time_residuals_for_template(
-            for i in range(1 + cut_wrong_peak_matches):
+                    N_residuals, template, interpolation_template=interp_template,
-                mask_count = 0
+                    antenna_dt=antenna_dt, antenna_timelength=antenna_timelength,
+                    snr_sigma_factor=snr_sigma_factor, bp_freq=bp_freq, normalise_noise=normalise_noise,
+                    h5_cache_fname=h5_cache_fname, rng=rng, tqdm=tqdm)
 
-                if i==1: # if cut_wrong_peak_matches:
+            print()# separating tqdm
-                    wrong_peak_condition = lambda t_res: abs(t_res) > antenna_dt*4
+            print()# separating tqdm
 
-                    mask = wrong_peak_condition(time_residuals)
+            # Make a plot of the time residuals
+            if N_residuals > 1:
+                for i in range(1 + cut_wrong_peak_matches):
+                    mask_count = 0
 
-                    mask_count = np.count_nonzero(mask)
+                    if i==1: # if cut_wrong_peak_matches:
+                        wrong_peak_condition = lambda t_res: abs(t_res) > antenna_dt*4
 
-                    print("Masking {} residuals".format(mask_count))
+                        mask = wrong_peak_condition(time_residuals)
-                    time_residuals = time_residuals[~mask]
 
-                    if not mask_count:
+                        mask_count = np.count_nonzero(mask)
-                        print("Continuing")
-                        continue
 
-                time_accuracies[k] = np.std(time_residuals)
+                        print("Masking {} residuals".format(mask_count))
-                mask_counts[k] = mask_count
+                        time_residuals = time_residuals[~mask]
 
-                hist_kwargs = dict(bins='sqrt', density=False, alpha=0.8, histtype='step')
+                        if not mask_count:
-                fig, ax = plt.subplots()
+                            continue
-                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)
+                    time_accuracies[l, k] = np.std(time_residuals)
-                if True: # fit gaussian to histogram
+                    mask_counts[l, k] = mask_count
-                    min_x = min(time_residuals)
-                    max_x = max(time_residuals)
 
-                    dx = bins[1] - bins[0]
+                    hist_kwargs = dict(bins='sqrt', density=False, alpha=0.8, histtype='step')
-                    scale = len(time_residuals) * dx
+                    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("#")
 
-                    xs = np.linspace(min_x, max_x)
+                    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)
 
-                    # do the fit
+                        dx = bins[1] - bins[0]
-                    name = "Norm"
+                        scale = len(time_residuals) * dx
-                    param_names = [ "$\\mu$", "$\\sigma$" ]
-                    distr_func = stats.norm
 
-                    label = name +"(" + ','.join(param_names) + ')'
+                        xs = np.linspace(min_x, max_x)
 
-                    # plot
+                        # do the fit
-                    fit_params = distr_func.fit(time_residuals)
+                        name = "Norm"
-                    fit_ys = scale * distr_func.pdf(xs, *fit_params)
+                        param_names = [ "$\\mu$", "$\\sigma$" ]
-                    ax.plot(xs, fit_ys, label=label)
+                        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)
-
-                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")
-
-                if True:
-                    plt.close(fig)
 
     # SNR time accuracy plot
     if True:
@@ -567,10 +569,11 @@ if __name__ == "__main__":
         ax.set_title(f"Template matching SNR vs time accuracy")
         ax.set_xlabel("Signal to Noise Factor")
         ax.set_ylabel("Time Accuracy [ns]")
+        ax.grid()
 
         ax.legend(title="\n".join([
             f"N={N_residuals}",
-            f"template_dt={template_dt:0.1e}ns",
+            #f"template_dt={template_dt:0.1e}ns",
             f"antenna_dt={antenna_dt:0.1e}ns",
             ]))
 
@@ -578,28 +581,36 @@ if __name__ == "__main__":
             ax.set_xscale('log')
             ax.set_yscale('log')
 
-        # plot the values
+        # plot the values per template_dt slice
-        l = None
+        template_dt_colors = [None]*len(template_dts)
-        for j, mask_threshold in enumerate(pairwise([np.inf, 250, 50, 1, 0])):
+        for k, template_dt in enumerate(template_dts):
-            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)
+            # indicate masking values
+            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 template_dt_colors[k] is None else template_dt_colors[k]
+                        )
+                mask = mask_counts[k] >= mask_threshold[1]
+                mask &= mask_counts[k] < mask_threshold[0]
 
-        if True: # limit y-axis to 1e1
+                l = ax.plot(snr_factors[mask], time_accuracies[k][mask], **kwargs)
-            ax.set_ylim([None, 1e1])
+                template_dt_colors[k] = l[0].get_color()
+
+            # indicate threshold
+            if True:
+                ax.axhline(template_dt/np.sqrt(12), ls='--', alpha=0.7, color=template_dt_colors[k], label=f'Template dt:{template_dt:0.1e}ns')
 
 
         # Set horizontal line at 1 ns
         if True:
             ax.axhline(1, ls='--', alpha=0.8, color='g')
-            ax.grid()
+
-            ax.axhline(template_dt/np.sqrt(12), ls='--', alpha=0.7, color='b')
+        ax.legend()
+
+        if True: # limit y-axis to 1e1
+            ax.set_ylim([None, 1e1])
 
         fig.tight_layout()
         fig.savefig(f"figures/11_time_res_vs_snr_tdt{template_dt:0.1e}.pdf")