Pulse: show residuals in SNR vs accuracy plot: Finished

simulations/11_pulsed_timing.py

@@ -625,8 +625,8 @@ if __name__ == "__main__":
     # SNR time accuracy plot
     #
     if True:
-        threshold_markers = ['^', 'v', '8', 'o']
+        threshold_markers = ['^', 'v', '8', 'X'] # make sure to have filled markers here
-        mask_thresholds = [np.inf, N_residuals*0.5, N_residuals*0.1, 1, 0]
+        mask_thresholds = np.array([np.inf, N_residuals*0.5, N_residuals*0.1, 1, 0])
 
         fig, ax = plt.subplots()
         ax.set_title(f"Template matching SNR vs time accuracy")
@@ -640,30 +640,66 @@ if __name__ == "__main__":
             f"antenna_dt={antenna_dt:0.1e}ns",
             ]))
 
-        if True:
-            ax.set_xscale('log')
-            ax.set_yscale('log')
-
         # plot the values per template_dt slice
         template_dt_colors = [None]*len(template_dts)
-        for k, template_dt in enumerate(template_dts):
+        for a, template_dt in enumerate(template_dts):
+            for k, snr_sigma_factor in enumerate(snr_factors):
+                time_residuals, snrs, valid_mask = time_residuals_data[a][k]
 
-            # indicate masking values
+                valid_mask = np.array(valid_mask, dtype=bool)
-            for j, mask_threshold in enumerate(pairwise(mask_thresholds)):
+
+                mean_residual = np.mean(time_residuals[valid_mask])
+                time_accuracy = np.std(time_residuals[valid_mask])
+
+                residual_mean_deviation = np.sqrt( (time_residuals - mean_residual)**2 )
+
+                scatter_kwargs = dict(
+                        ls='none',
+                        marker='.',
+                        alpha=0.3,
+                        zorder=1.8,
+                )
+
+                y_values = residual_mean_deviation
+
+                # snr_sigma_factor is a factor 2 too low
+                snr_sigma_factor *= 2
+
+                # plot all invalid datapoints
+                if True:
+                    ax.plot(snrs[~valid_mask], y_values[~valid_mask], color='grey', **scatter_kwargs)
+
+                # plot valid datapoints
+                if True:
+                    if template_dt_colors[a] is not None:
+                        scatter_kwargs['color'] = template_dt_colors[a]
+
+                    l = ax.plot(snrs[valid_mask], y_values[valid_mask], **scatter_kwargs)
+
+                    template_dt_colors[a] = l[0].get_color()
+
+                masked_count = np.count_nonzero(~valid_mask)
+
+                # plot accuracy indicating masking counts
                 kwargs = dict(
                         ls='none',
-                        marker=threshold_markers[j],
+                        color= None if template_dt_colors[a] is None else template_dt_colors[a],
-                        color= None if template_dt_colors[k] is None else template_dt_colors[k]
+                        marker=threshold_markers[np.argmin( masked_count <= mask_thresholds)-1],
-                        )
+                        ms=10,
-                mask = mask_counts[k] >= mask_threshold[1]
+                        markeredgecolor='white',
-                mask &= mask_counts[k] < mask_threshold[0]
+                        markeredgewidth=1,
+                    )
 
-                l = ax.plot(snr_factors[mask], time_accuracies[k][mask], **kwargs)
+                #l = ax.plot(snr_sigma_factor, np.sqrt(np.mean(y_values[valid_mask])**2), **{**kwargs, **dict(ms=50)})
-                template_dt_colors[k] = l[0].get_color()
+                l = ax.plot(snr_sigma_factor, np.std(time_residuals[valid_mask]), **kwargs)
 
-            # indicate threshold
+                # set color if not yet set
+                template_dt_colors[a] = l[0].get_color()
+
+
+            # indicate boxcar 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')
+                ax.axhline(template_dt/np.sqrt(12), ls='--', alpha=0.7, color=template_dt_colors[a], label=f'Template dt:{template_dt:0.1e}ns')
 
 
         # Set horizontal line at 1 ns
@@ -672,8 +708,32 @@ if __name__ == "__main__":
 
         ax.legend()
 
-        if True: # limit y-axis to 1e1
+        fig.tight_layout()
-            ax.set_ylim([None, 1e1])
+        fig.savefig(f"figures/11_time_res_vs_snr_full_linear.pdf")
+
+        # logscaling
+        if True:
+            ax.set_xscale('log')
+            ax.set_yscale('log')
+
+            # limit y-axis upper limit to 1e1
+            if True:
+                this_lim = 1e1
+                ax.set_ylim([None, this_lim])
+
+            # require y-axis lower limit to be at least 1e-1
+            if True:
+                this_lim = 1e-1
+                low_ylims = ax.get_ylim()[0]
+                if low_ylims >= this_lim:
+                    ax.set_ylim([this_lim, None])
+
+            # .. but keep it above 1e-3
+            if True:
+                this_lim = 1e-3
+                low_ylims = ax.get_ylim()[0]
+                if low_ylims <= this_lim:
+                    ax.set_ylim([this_lim, None])
 
         if True: # require y-axis lower limit to be at least 1e-1
             low_ylims = ax.get_ylim()[0]