diff --git a/simulations/airshower_beacon_simulation/bb_measure_clock_phase.py b/simulations/airshower_beacon_simulation/bb_measure_clock_phase.py
index b2efc53..1b66d77 100755
--- a/simulations/airshower_beacon_simulation/bb_measure_clock_phase.py
+++ b/simulations/airshower_beacon_simulation/bb_measure_clock_phase.py
@@ -161,7 +161,8 @@ if __name__ == "__main__":
                         actual_clock_phases,
                         plot_residuals=plot_residuals,
                         f_beacon=f_beacon,
-                        figsize=figsize
+                        figsize=figsize,
+                        fit_gaussian=plot_residuals,
                         )
 
                 if plot_residuals:
diff --git a/simulations/airshower_beacon_simulation/bc_baseline_phase_deltas.py b/simulations/airshower_beacon_simulation/bc_baseline_phase_deltas.py
index f9d7307..0175df1 100755
--- a/simulations/airshower_beacon_simulation/bc_baseline_phase_deltas.py
+++ b/simulations/airshower_beacon_simulation/bc_baseline_phase_deltas.py
@@ -125,6 +125,7 @@ if __name__ == "__main__":
                     f_beacon=f_beacon,
                     figsize=figsize,
                     hist_kwargs=hist_kwargs,
+                    fit_gaussian=plot_residuals,
                     )
 
             axs = fig.get_axes()
diff --git a/simulations/airshower_beacon_simulation/bd_antenna_phase_deltas.py b/simulations/airshower_beacon_simulation/bd_antenna_phase_deltas.py
index 22bbd96..897a719 100755
--- a/simulations/airshower_beacon_simulation/bd_antenna_phase_deltas.py
+++ b/simulations/airshower_beacon_simulation/bd_antenna_phase_deltas.py
@@ -191,6 +191,7 @@ if __name__ == "__main__":
                     f_beacon=f_beacon,
                     figsize=figsize,
                     hist_kwargs=hist_kwargs,
+                    fit_gaussian=plot_residuals,
                     )
 
             axs = fig.get_axes()
diff --git a/simulations/airshower_beacon_simulation/cb_report_measured_antenna_time_offsets.py b/simulations/airshower_beacon_simulation/cb_report_measured_antenna_time_offsets.py
index 7fbf1f0..2871a03 100755
--- a/simulations/airshower_beacon_simulation/cb_report_measured_antenna_time_offsets.py
+++ b/simulations/airshower_beacon_simulation/cb_report_measured_antenna_time_offsets.py
@@ -95,6 +95,7 @@ if __name__ == "__main__":
                     figsize=figsize,
                     hist_kwargs=hist_kwargs,
                     secondary_axis='phase',
+                    fit_gaussian=True,
                     )
 
             axs = fig.get_axes()
diff --git a/simulations/airshower_beacon_simulation/lib/figlib.py b/simulations/airshower_beacon_simulation/lib/figlib.py
index 6904ff6..f18e10d 100644
--- a/simulations/airshower_beacon_simulation/lib/figlib.py
+++ b/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