From 929c6c7748ce402291648b314a1f4cd16e6a17e7 Mon Sep 17 00:00:00 2001 From: Eric Teunis de Boone Date: Mon, 19 Dec 2022 21:34:49 +0100 Subject: [PATCH] ZH: update baseline phase script --- .../bc_baseline_phase_deltas.py | 83 ++++++++++--------- 1 file changed, 46 insertions(+), 37 deletions(-) diff --git a/simulations/airshower_beacon_simulation/bc_baseline_phase_deltas.py b/simulations/airshower_beacon_simulation/bc_baseline_phase_deltas.py index afe2f61..40ea6d3 100755 --- a/simulations/airshower_beacon_simulation/bc_baseline_phase_deltas.py +++ b/simulations/airshower_beacon_simulation/bc_baseline_phase_deltas.py @@ -21,7 +21,7 @@ if __name__ == "__main__": fname = "ZH_airshower/mysim.sry" c_light = 3e8*1e-9 - show_plots = not True + show_plots = True ref_ant_id = None # leave None for all baselines #### @@ -66,7 +66,6 @@ if __name__ == "__main__": # Get true phase diffs try: - true_phases = np.array([ant.beacon_info[freq_name]['true_phase'] for ant in base]) true_phases_diff = lib.phase_mod(lib.phase_mod(true_phases[1]) - lib.phase_mod(true_phases[0])) except IndexError: @@ -79,59 +78,69 @@ if __name__ == "__main__": beacon.write_baseline_time_diffs_hdf5(time_diffs_fname, baselines, phase_diffs[:,1], [0]*len(phase_diffs), phase_diffs[:,0]) - # Read actual phases from antenna hdf5 - actual_antenna_measured_phases = { a.name: 2*np.pi*a.attrs['clock_offset']*f_beacon for a in antennas } + ############################## + # Compare actual time shifts # + ############################## + actual_antenna_true_phases = { a.name: -2*np.pi*a.attrs['clock_offset']*f_beacon for a in antennas } # Compare actual time shifts my_phase_diffs = [] for i,b in enumerate(baselines): - actual_phase_measured_diff = lib.phase_mod( lib.phase_mod(actual_antenna_measured_phases[b[1].name]) - lib.phase_mod(actual_antenna_measured_phases[b[0].name])) + actual_true_phase_diff = lib.phase_mod( lib.phase_mod(actual_antenna_true_phases[b[1].name]) - lib.phase_mod(actual_antenna_true_phases[b[0].name])) - # remove phase due to time delay from transmitter difference - tds = np.array([ lib.geometry_time(tx, ant, c_light=c_light) for ant in base]) - delta_td = tds[1] - tds[0] - delta_td_phase = lib.phase_mod(delta_td*2*np.pi*f_beacon) - - actual_phase_diff = lib.phase_mod( actual_phase_measured_diff - delta_td_phase) - my_phase_diffs.append(actual_phase_diff) + this_actual_true_phase_diff = lib.phase_mod( actual_true_phase_diff ) + my_phase_diffs.append(this_actual_true_phase_diff) # Make a plot if True: N_base = len(baselines) N_ant = len(antennas) - phase_residuals = lib.phase_mod(phase_diffs[:,1] - my_phase_diffs) + for i in range(2): + plot_residuals = i == 1 + colors = ['blue', 'orange'] - fig, axs = plt.subplots(2, 1, sharex=True) - axs[0].set_title("Measured phase difference - Actual phase difference") - axs[0].set_xlabel("Phase $\\Delta\\varphi = \\varphi_{meas} - \\varphi_{true}$") - axs[0].tick_params(bottom=True, labelbottom=True) - #axs[1].tick_params(top=True, labeltop=True, bottom=False, labelbottom=False) + fig, axs = plt.subplots(2, 1, sharex=True) - if True: - forward = lambda x: x/(2*np.pi*f_beacon) - inverse = lambda x: 2*np.pi*x*f_beacon - secax = axs[0].secondary_xaxis('top', functions=(forward, inverse)) - secax.set_xlabel('Time $\\Delta\\varphi/(2\\pi f_{beac})$ [ns]') + if True: + forward = lambda x: x/(2*np.pi*f_beacon) + inverse = lambda x: 2*np.pi*x*f_beacon + secax = axs[0].secondary_xaxis('top', functions=(forward, inverse)) + secax.set_xlabel('Time $\\Delta\\varphi/(2\\pi f_{beac})$ [ns]') - i=0 - axs[i].set_ylabel("#") - axs[i].hist(phase_residuals, bins='sqrt', density=False) + if plot_residuals: + phase_residuals = lib.phase_mod(phase_diffs[:,1] - my_phase_diffs) - i=1 - axs[i].set_ylabel("Baseline\n combination #") - if not True: - axs[i].plot(my_phase_diffs, np.arange(N_base), ls='none', marker='+', label='actual time shifts') - l = axs[i].plot(phase_diffs[:,1], np.arange(N_base), ls='none', marker='x', label='calculated') + axs[0].set_title("Difference between Measured and Actual phase difference for Baseline (i,j)") + axs[1].set_xlabel("Baseline Phase Residual $\\varphi_{ij_{meas}} - \\varphi_{ij_{true}}$ [rad]") + else: + axs[0].set_title("Comparison Measured and Actual phase difference for Baseline (i,j)") + axs[1].set_xlabel("Baseline Phase $\\varphi_{ij}$ [rad]") - axs[i].legend() + i=0 + axs[i].set_ylabel("#") + if plot_residuals: + axs[i].hist(phase_residuals, bins='sqrt', density=False, alpha=0.8, color=colors[0]) + else: + axs[i].hist(phase_diffs[:,1], bins='sqrt', density=False, alpha=0.8, color=colors[0], ls='solid' , histtype='step', label='Measured') + axs[i].hist(my_phase_diffs, bins='sqrt', density=False, alpha=0.8, color=colors[1], ls='dashed', histtype='step', label='Actual') - else: - axs[i].plot(phase_residuals, np.arange(N_base), ls='none', marker='x') - fig.tight_layout() + i=1 + axs[i].set_ylabel("Baseline no.") + if not plot_residuals: + axs[i].plot(phase_diffs[:,1], np.arange(N_base), alpha=0.8, color=colors[0], ls='none', marker='x', label='calculated') + axs[i].plot(my_phase_diffs, np.arange(N_base), alpha=0.8, color=colors[1], ls='none', marker='+', label='actual time shifts') - if fig_dir: - fig.savefig(path.join(fig_dir, __file__ + f".F{freq_name}.pdf")) + axs[i].legend() + else: + axs[i].plot(phase_residuals, np.arange(N_base), alpha=0.6, ls='none', marker='x', color=colors[0]) + fig.tight_layout() + + if fig_dir: + extra_name = "measured" + if plot_residuals: + extra_name = "differences" + fig.savefig(path.join(fig_dir, __file__ + f".{extra_name}.F{freq_name}.pdf")) if show_plots: plt.show()