ZH: extract code plotting measured and actual phases (or residuals)

simulations/airshower_beacon_simulation/lib/figlib.py

@@ -0,0 +1,65 @@
+import matplotlib.pyplot as plt
+import numpy as np
+
+def phase_comparison_figure(
+        measured_phases,
+        true_phases,
+        plot_residuals=True,
+        f_beacon=None,
+        hist_kwargs={},
+        sc_kwargs={},
+        colors=['blue', 'orange'],
+        legend_on_scatter=True,
+        **fig_kwargs
+        ):
+    """
+    Create a figure comparing measured_phase against true_phase
+    by both plotting the values, and the residuals.
+    """
+    default_fig_kwargs = dict(sharex=True)
+    fig_kwargs = {**default_fig_kwargs, **fig_kwargs}
+
+    do_hist_plot = hist_kwargs is not False
+    do_scatter_plot = sc_kwargs is not False
+
+    fig, axs = plt.subplots(0+do_hist_plot+do_scatter_plot, 1, **fig_kwargs)
+
+    if not hasattr(axs, '__len__'):
+        axs = [axs]
+
+    if f_beacon:
+        phase2time = lambda x: x/(2*np.pi*f_beacon)
+        time2phase = lambda x: 2*np.pi*x*f_beacon
+        secax = axs[0].secondary_xaxis('top', functions=(phase2time, time2phase))
+        secax.set_xlabel('Time $\\varphi/(2\\pi f_{beac})$ [ns]')
+
+    # Histogram
+    if do_hist_plot:
+        i=0
+        default_hist_kwargs = dict(bins='sqrt', density=False, alpha=0.8, histtype='step')
+        hist_kwargs = {**default_hist_kwargs, **hist_kwargs}
+
+        axs[i].set_ylabel("#")
+
+        _counts, _bins, _patches = axs[i].hist(measured_phases, color=colors[0], label='Measured', ls='solid', **hist_kwargs)
+        if not plot_residuals: # also plot the true clock phases
+            axs[i].hist(true_phases, color=colors[1], label='Actual', ls='dashed', **{**hist_kwargs, **dict(bins=_bins)})
+
+    # Scatter plot
+    if do_scatter_plot:
+        i=1
+        default_sc_kwargs = dict(alpha=0.6, ls='none')
+        sc_kwargs = {**default_sc_kwargs, **sc_kwargs}
+
+        axs[i].set_ylabel("Antenna no.")
+        axs[i].plot(measured_phases, np.arange(len(measured_phases)), marker='x' if plot_residuals else '3', color=colors[0], label='Measured', **sc_kwargs)
+
+        if not plot_residuals: # also plot the true clock phases
+            axs[i].plot(true_phases, np.arange(len(true_phases)), marker='4', color=colors[1], label='Actual', **sc_kwargs)
+
+        if not plot_residuals and legend_on_scatter:
+            axs[i].legend()
+
+    fig.tight_layout()
+
+    return fig