ZH: fit gaussian to antenna clock phase figures

simulations/airshower_beacon_simulation/lib/figlib.py

@@ -1,6 +1,8 @@
 import matplotlib.pyplot as plt
 import numpy as np
 
+from scipy.stats import norm
+
 def phase_comparison_figure(
         measured_phases,
         true_phases,
@@ -11,6 +13,7 @@ def phase_comparison_figure(
         colors=['blue', 'orange'],
         legend_on_scatter=True,
         secondary_axis='time',
+        fit_gaussian=False,
         **fig_kwargs
         ):
     """
@@ -53,6 +56,24 @@ def phase_comparison_figure(
         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)})
 
+        if fit_gaussian:# Fit a gaussian to the histogram
+            gauss_kwargs = dict(color=colors[0], ls='dotted')
+            xs = np.linspace(min(_bins), max(_bins))
+
+            mean, std = norm.fit(measured_phases)
+            dx = _bins[1] - _bins[0]
+            scale = len(measured_phases)*dx
+            fit_ys = norm.pdf(xs, mean, std) * scale
+
+            axs[i].axvline(mean, ls='dotted', color='k', lw=2)
+
+            axs[i].plot( xs, fit_ys, label='fit', **gauss_kwargs)
+
+            # put mean, sigma and chisq on plot
+            text_str = f"$\\mu$ = {mean: .2e}\n$\\sigma$ = {std: .2e}"
+            text_kwargs = dict( transform=axs[i].transAxes, fontsize=14, verticalalignment='top')
+            axs[i].text(0.02, 0.95, text_str, **text_kwargs)
+
     # Scatter plot
     if do_scatter_plot:
         i=1