From 72a98d3f0906d840d7e3e3bd5fb3ebe22733c50b Mon Sep 17 00:00:00 2001 From: Eric Teunis de Boone Date: Mon, 9 Jan 2023 13:39:36 +0100 Subject: [PATCH] ZH: remove global time shift when comparing measured and actual time shift --- ...cb_report_measured_antenna_time_offsets.py | 88 ++++++++++++++----- 1 file changed, 64 insertions(+), 24 deletions(-) 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 d4ea78b..0f89f0d 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 @@ -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()