From 81b3502de53715af47b290c91de930966deae17b Mon Sep 17 00:00:00 2001
From: Eric Teunis de Boone <ericteunis@deboone.nl>
Date: Mon, 24 Apr 2023 17:52:19 +0200
Subject: [PATCH] Pulse: fit gaussian

---
 simulations/11_pulsed_timing.py | 60 +++++++++++++++++++++++++++++----
 1 file changed, 54 insertions(+), 6 deletions(-)

diff --git a/simulations/11_pulsed_timing.py b/simulations/11_pulsed_timing.py
index e67c8ba..ae4eb38 100755
--- a/simulations/11_pulsed_timing.py
+++ b/simulations/11_pulsed_timing.py
@@ -2,9 +2,10 @@
 
 
 from lib import util
-from scipy import signal, interpolate
+from scipy import signal, interpolate, stats
 import matplotlib.pyplot as plt
 import numpy as np
+from itertools import zip_longest
 
 try:
     from tqdm import tqdm
@@ -203,8 +204,6 @@ if __name__ == "__main__":
 
             antenna.peak_sample = antenna.peak_time/antenna.dt
             antenna_true_signal = antenna.signal
-
-
         true_time_offset = antenna.peak_time - template.peak_time
 
         if do_plots:
@@ -296,6 +295,9 @@ if __name__ == "__main__":
         best_time_lag = best_sample_lag * lag_dt
         time_residuals[j] = best_time_lag - true_time_offset
 
+        if not do_plots:
+            continue
+
         if do_plots and axs2:
             axs2[-1].axvline(best_time_lag, color='r', alpha=0.5, linewidth=2)
             axs2[-1].axvline(true_time_offset, color='g', alpha=0.5, linewidth=2)
@@ -368,13 +370,59 @@ if __name__ == "__main__":
 
     # Make a plot of the time residuals
     if len(time_residuals) > 1:
+
+        hist_kwargs = dict(bins='sqrt', density=False, alpha=0.8, histtype='step')
         fig, ax = plt.subplots()
-        ax.set_title("Template Correlation Lag finding")
+        ax.set_title(
+                "Template Correlation Lag finding"
+                + f"\n template dt: {template.dt*1e3: .1e}ps"
+                + f"; antenna dt: {antenna.dt: .1e}ns"
+                + f"; noise_factor: {noise_sigma_factor: .1e}"
+                )
         ax.set_xlabel("Time Residual [ns]")
         ax.set_ylabel("#")
-        ax.hist(time_residuals, bins='sqrt', density=False)
 
-        ax.legend(title=f"template dt: {template.dt: .1e}ns\nantenna dt: {antenna.dt: .1e}ns")
+        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)
+
+            dx = bins[1] - bins[0]
+            scale = len(time_residuals) * dx
+
+            xs = np.linspace(min_x, max_x)
+
+            # do the fit
+            name = "Norm"
+            param_names = [ "$\\mu$", "$\\sigma$" ]
+            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)
 
         fig.savefig("figures/11_time_residual_hist.pdf")