From e8ff82812086863ef377ac6778e7aa7dbd93198c Mon Sep 17 00:00:00 2001 From: Eric Teunis de Boone Date: Wed, 21 Dec 2022 17:24:44 +0100 Subject: [PATCH] ZH: improve antenna matrix solving --- .../bc_baseline_phase_deltas.py | 10 +- .../bd_antenna_phase_deltas.py | 152 +++++++++++++++--- 2 files changed, 138 insertions(+), 24 deletions(-) diff --git a/simulations/airshower_beacon_simulation/bc_baseline_phase_deltas.py b/simulations/airshower_beacon_simulation/bc_baseline_phase_deltas.py index 40ea6d3..6bf9032 100755 --- a/simulations/airshower_beacon_simulation/bc_baseline_phase_deltas.py +++ b/simulations/airshower_beacon_simulation/bc_baseline_phase_deltas.py @@ -22,12 +22,12 @@ if __name__ == "__main__": fname = "ZH_airshower/mysim.sry" c_light = 3e8*1e-9 show_plots = True - ref_ant_id = None # leave None for all baselines + ref_ant_id = 50 # leave None for all baselines #### fname_dir = path.dirname(fname) antennas_fname = path.join(fname_dir, beacon.antennas_fname) - time_diffs_fname = 'time_diffs.hdf5' if not True else antennas_fname + time_diffs_fname = 'time_diffs.hdf5' if False else antennas_fname fig_dir = "./figures" # set None to disable saving @@ -111,10 +111,10 @@ if __name__ == "__main__": if plot_residuals: phase_residuals = lib.phase_mod(phase_diffs[:,1] - my_phase_diffs) - axs[0].set_title("Difference between Measured and Actual phase difference for Baseline (i,j)") + axs[0].set_title("Difference between Measured and Actual phase difference\n for Baseline (i,j" + ( '='+str(ref_ant_id) if ref_ant_id is not None else None) + ')') 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[0].set_title("Comparison Measured and Actual phase difference\n for Baseline (i,j" + ( '='+str(ref_ant_id) if ref_ant_id is not None else None) + ')') axs[1].set_xlabel("Baseline Phase $\\varphi_{ij}$ [rad]") i=0 @@ -139,7 +139,7 @@ if __name__ == "__main__": if fig_dir: extra_name = "measured" if plot_residuals: - extra_name = "differences" + extra_name = "residuals" fig.savefig(path.join(fig_dir, __file__ + f".{extra_name}.F{freq_name}.pdf")) if show_plots: diff --git a/simulations/airshower_beacon_simulation/bd_antenna_phase_deltas.py b/simulations/airshower_beacon_simulation/bd_antenna_phase_deltas.py index d4824f5..d484c8c 100755 --- a/simulations/airshower_beacon_simulation/bd_antenna_phase_deltas.py +++ b/simulations/airshower_beacon_simulation/bd_antenna_phase_deltas.py @@ -14,6 +14,11 @@ if __name__ == "__main__": from os import path import sys + import os + import matplotlib + if os.name == 'posix' and "DISPLAY" not in os.environ: + matplotlib.use('Agg') + fname = "ZH_airshower/mysim.sry" show_plots = True @@ -22,35 +27,144 @@ if __name__ == "__main__": #### fname_dir = path.dirname(fname) antennas_fname = path.join(fname_dir, beacon.antennas_fname) - time_diffs_fname = antennas_fname + time_diffs_fname = 'time_diffs.hdf5' if False else antennas_fname + fig_dir = "./figures" # set None to disable saving - basenames, time_diffs, f_beacon, true_phase_diffs, k_periods = beacon.read_baseline_time_diffs_hdf5(time_diffs_fname) + basenames, time_diffs, f_beacons, true_phase_diffs, k_periods = beacon.read_baseline_time_diffs_hdf5(time_diffs_fname) f_beacon, tx, antennas = beacon.read_beacon_hdf5(antennas_fname) + # TODO: allow multiple frequencies + if (f_beacon != f_beacons).any(): + raise NotImplementedError + + N_base = len(basenames) + N_ant = len(antennas) + + # reshape time_diffs into N_ant x N_ant array + sigma_phase_matrix = np.full( (N_ant, N_ant), np.nan, dtype=float) + #sigma_phase_matrix = np.arange(0, N_ant**2, dtype=float).reshape( (N_ant,N_ant)) + + name2idx = lambda name: int(name)-1 + for i, b in enumerate(basenames): + idx = (name2idx(b[0]), name2idx(b[1])) + + if idx[0] == idx[1]: + pass + + sigma_phase_matrix[(idx[0], idx[1])] = true_phase_diffs[i] + sigma_phase_matrix[(idx[1], idx[0])] = true_phase_diffs[i] + + # for each row j subtract the 0,j element from the whole row + # and apply phase_mod + first_row = sigma_phase_matrix[0] + + sigma_phase_matrix = sigma_phase_matrix - first_row[:,np.newaxis] + sigma_phase_matrix = lib.phase_mod(sigma_phase_matrix) + + + # Except for the first row, these are all separate measurements + # Condense into phase offset per antenna + if True: # do not use the first row + my_mask = np.arange(1, len(sigma_phase_matrix), dtype=int) + else: + my_mask = np.arange(0, len(sigma_phase_matrix), dtype=int) + + mean_sigma_phase = np.nanmean(sigma_phase_matrix[my_mask], axis=0) + std_sigma_phase = np.nanstd( sigma_phase_matrix[my_mask], axis=0) + + # write into antenna hdf5 + with h5py.File(antennas_fname, 'a') as fp: + group = fp['antennas'] + + freq_name = None + for i, ant in enumerate(antennas): + h5ant = group[ant.name] + + h5beacon_info = h5ant['beacon_info'] + + # find out freq_name + if freq_name is None: + freq_name = [ k for k in h5beacon_info.keys() if np.isclose(h5beacon_info[k].attrs['freq'], f_beacon)][0] + + h5attrs = h5beacon_info[freq_name].attrs + + idx = name2idx(ant.name) + h5attrs['sigma_phase_mean'] = mean_sigma_phase[idx] + h5attrs['sigma_phase_std'] = std_sigma_phase[idx] + + + ############################## + # Compare actual time shifts # + ############################## antenna_time_shifts = { a.name: a.attrs['clock_offset'] for a in antennas } + if True: + # show means and std + fig, axs = plt.subplots(1,2, sharey=True) + axs[0].set_title("") + axs[0].set_xlabel("Antenna no.") + axs[0].set_ylabel("Antenna Phase $\\Delta_\\varphi$") + axs[1].set_xlabel("#") + if True: + forward = lambda x: x/(2*np.pi*f_beacon) + inverse = lambda x: 2*np.pi*x*f_beacon + secax = axs[1].secondary_yaxis('right', functions=(forward, inverse)) + secax.set_ylabel('Time $\\Delta\\varphi/(2\\pi f_{beac})$ [ns]') - my_time_shifts = [] + l = axs[0].errorbar(np.arange(N_ant), mean_sigma_phase, yerr=std_sigma_phase, marker='.', alpha=0.7) + + axs[1].hist(mean_sigma_phase, bins='sqrt', density=False, orientation='horizontal', color=l[0].get_color(), histtype='step') + + # Actual time shifts + if True: + actual_phase_shifts = [ -1*lib.phase_mod(2*np.pi*f_beacon*v) for k,v in antenna_time_shifts.items() ] + antenna_names = [int(k) for k,v in antenna_time_shifts.items() ] + + # make sure to keep the same offset + if True: + phase_offset = mean_sigma_phase[0] - actual_phase_shifts[0] + + actual_phase_shifts += phase_offset + + l = axs[0].plot(antenna_names, actual_phase_shifts, ls='none', marker='3', alpha=0.8) + + if True: + axs[1].hist(actual_phase_shifts, bins='sqrt', density=False, orientation='horizontal', ls='dashed', color=l[0].get_color(), histtype='step') + + fig.tight_layout() + if fig_dir: + fig.savefig(path.join(fig_dir, __file__ + f".residuals.pdf")) + + ########################## + ########################## + ########################## + + actual_time_shifts = [] for i,b in enumerate(basenames): - actual_time_shift = antenna_time_shifts[b[0]] - antenna_time_shifts[b[1]] - my_time_shifts.append(actual_time_shift) + actual_time_shift = lib.phase_mod(lib.phase_mod(antenna_time_shifts[b[0]]*2*np.pi*f_beacon) - lib.phase_mod(antenna_time_shifts[b[1]]*2*np.pi*f_beacon)) - print( - f'B({b[0]}, {b[1]}):', - time_diffs[i], - k_periods[i], - 'A', actual_time_shift - ) + actual_time_shifts.append(actual_time_shift) + + # unpack mean_sigma_phase back into a list of time diffs + measured_time_diffs = [] + for i,b in enumerate(basenames): + time0, time1 = mean_sigma_phase[name2idx(b[0])], mean_sigma_phase[name2idx(b[1])] + measured_time_diffs.append(time1 - time0) # Make a plot - N = len(basenames) - fig, ax = plt.subplots() - ax.set_xlabel("Baseline combo") - ax.set_ylabel("[ns]") - ax.plot(np.arange(N), my_time_shifts, marker='+', label='actual time shifts') - ax.plot(np.arange(N), time_diffs, marker='x', label='calculated') + if True: + fig, ax = plt.subplots() + ax.set_xlabel("Baseline no.") + ax.set_ylabel("$\\Delta t$[ns]") + if True: # indicate single beacon period span + ax.plot((-1, -1), (0, 1/f_beacon), marker='3', ms=10, label='1/f_beacon') + ax.plot(np.arange(N_base), actual_time_shifts, marker='+', label='actual time shifts') + ax.plot(np.arange(N_base), measured_time_diffs, marker='x', label='calculated') - ax.legend() + ax.legend() + if fig_dir: + fig.savefig(path.join(fig_dir, __file__ + f".calculated_shifts.pdf")) - plt.show() + if show_plots: + plt.show()