ZH: remove global time shift when comparing measured and actual time shift

simulations/airshower_beacon_simulation/cb_report_measured_antenna_time_offsets.py

@@ -5,14 +5,11 @@
 Report best time offset per frequency for each antenna
 """
 
-import h5py
 import matplotlib.pyplot as plt
 import numpy as np
 from os import path
-from scipy.interpolate import interp1d
 
 import aa_generate_beacon as beacon
-import lib
 
 if __name__ == "__main__":
     import sys
@@ -23,9 +20,8 @@ if __name__ == "__main__":
 
     fname = "ZH_airshower/mysim.sry"
 
-    fig_dir = "./figures/periods_from_shower_figures/"
-    fig_subdir = path.join(fig_dir, 'shifts/')
-    show_plots = False
+    fig_dir = "./figures" # set None to disable saving
+    show_plots = not False
 
     ####
     fname_dir = path.dirname(fname)
@@ -35,8 +31,6 @@ if __name__ == "__main__":
     # create fig_dir
     if fig_dir:
         os.makedirs(fig_dir, exist_ok=True)
-        if fig_subdir:
-            os.makedirs(fig_subdir, exist_ok=True)
 
     # Read in antennas from file
     _, tx, antennas = beacon.read_beacon_hdf5(antennas_fname)
@@ -50,31 +44,77 @@ if __name__ == "__main__":
     f_beacon = antennas[0].beacon_info[freq_name]['freq']
 
 
+    # TODO: redo matrix sweeping for new timing??
+    measured_antenna_time_shifts = {}
     for i, ant in enumerate(antennas):
         clock_phase_time = ant.beacon_info[freq_name]['sigma_phase_mean']/(2*np.pi*f_beacon)
-
         best_k_time = ant.beacon_info[freq_name]['best_k_time']
-        best_k = ant.beacon_info[freq_name]['best_k_time']
 
         total_clock_time = best_k_time + clock_phase_time
+        measured_antenna_time_shifts[ant.name] = -1*total_clock_time
 
-        actual_clock_time = ant.attrs['clock_offset']
+    ###
+    # Compare actual vs measured time shifts
+    ###
+    actual_antenna_time_shifts = { a.name: a.attrs['clock_offset'] for a in sorted(antennas, key=lambda a: int(a.name)) }
 
-        # filter k == 0
-        if best_k == 0:
-            continue
+    N_ant = len(antennas)
 
-        # Report to stdout
-        print("Timing of antenna", ant.name)
-        print(" + k:", best_k, ";(-)", best_k_time, 'ns')
-        print(" + phase_time:", clock_phase_time, 'ns')
-        print("==============================")
-        print("Sum: (-) ", total_clock_time, 'ns')
-        print("Actual:", actual_clock_time, 'ns')
-        print("Diff (A - - S):", actual_clock_time - -total_clock_time, 'ns')
-        print()
+    if True:
+        # keep dataset in the same ordering
+        antenna_names = [int(k)-1 for k,v in actual_antenna_time_shifts.items()]
+        actual_time_shifts =   np.array([ v for k,v in actual_antenna_time_shifts.items()])
+        measured_time_shifts = np.array([ measured_antenna_time_shifts[k] for k,v in actual_antenna_time_shifts.items() ])
+
+        # remove global shift
+        global_shift = actual_time_shifts[0] - measured_time_shifts[0]
+        actual_time_shifts -= global_shift
+
+        for i in range(2):
+            plot_residuals = i == 1
+            colors = ['blue', 'orange']
+
+            fig, axs = plt.subplots(2, 1, sharex=True)
+
+            if True:
+                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=(time2phase, phase2time))
+                secax.set_xlabel('Phase $2\\pi t f_{beac}$ [rad]')
+
+            if plot_residuals:
+                time_shift_residuals = measured_time_shifts - actual_time_shifts
+                fig.suptitle("Difference between Measured and Actual clock offsets")
+                axs[-1].set_xlabel("Antenna Time Offset Residual $\\Delta_t$ [ns]")
+            else:
+                fig.suptitle("Comparison Measured and Actual clock offset")
+                axs[-1].set_xlabel("Antenna Time Offset $t_c = \\left(\\frac{\\Delta\\varphi}{2\\pi} + k\\right) / f_{beac}$ [ns]")
+
+            i=0
+            axs[i].set_ylabel("#")
+            if plot_residuals:
+                axs[i].hist(time_shift_residuals, bins='sqrt', alpha=0.8, color=colors[0])
+            else:
+                axs[i].hist(measured_time_shifts, bins='sqrt', density=False, alpha=0.8, color=colors[0], ls='solid' , histtype='step', label='Measured')
+                axs[i].hist(actual_time_shifts,   bins='sqrt', density=False, alpha=0.8, color=colors[1], ls='dashed', histtype='step', label='Actual')
 
 
+            i=1
+            axs[i].set_ylabel("Antenna no.")
+            if plot_residuals:
+                axs[i].plot(time_shift_residuals, np.arange(N_ant), alpha=0.6, ls='none', marker='x', color=colors[0])
+            else:
+                axs[i].errorbar(measured_time_shifts, np.arange(N_ant), yerr=None, marker='4', alpha=0.7, ls='none', color=colors[0], label='Measured')
+                axs[i].plot(actual_time_shifts,  antenna_names, ls='none', marker='3', alpha=0.8, color=colors[1], label='Actual')
 
+                axs[i].legend()
+            fig.tight_layout()
 
+            if fig_dir:
+                extra_name = "comparison"
+                if plot_residuals:
+                    extra_name = "residuals"
+                fig.savefig(path.join(fig_dir, path.basename(__file__) + f".time.{extra_name}.pdf"))
 
+    if show_plots:
+        plt.show()