From 8fc7d4bc0c2c3aab7ed298faeb4540423c2263fd Mon Sep 17 00:00:00 2001 From: Eric Teunis de Boone Date: Tue, 7 Feb 2023 17:58:45 +0100 Subject: [PATCH] ZH: employ phase_comparison_figure function --- .../bc_baseline_phase_deltas.py | 47 +++++++++---------- .../bd_antenna_phase_deltas.py | 42 +++++++---------- ...cb_report_measured_antenna_time_offsets.py | 47 ++++++++----------- .../airshower_beacon_simulation/lib/figlib.py | 14 ++++-- 4 files changed, 71 insertions(+), 79 deletions(-) diff --git a/simulations/airshower_beacon_simulation/bc_baseline_phase_deltas.py b/simulations/airshower_beacon_simulation/bc_baseline_phase_deltas.py index a3c067f..f9d7307 100755 --- a/simulations/airshower_beacon_simulation/bc_baseline_phase_deltas.py +++ b/simulations/airshower_beacon_simulation/bc_baseline_phase_deltas.py @@ -8,6 +8,7 @@ import numpy as np import aa_generate_beacon as beacon import lib +from lib import figlib if __name__ == "__main__": @@ -108,45 +109,41 @@ if __name__ == "__main__": for i in range(2): plot_residuals = i == 1 - colors = ['blue', 'orange'] - fig, axs = plt.subplots(2, 1, sharex=True, figsize=figsize) + true_phases = my_phase_diffs + measured_phases = phase_diffs[:,1] - 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=(phase2time, time2phase)) - secax.set_xlabel('Time $\\Delta\\varphi/(2\\pi f_{beac})$ [ns]') + hist_kwargs = {} + if plot_residuals: + measured_phases = lib.phase_mod(measured_phases - true_phases) + hist_kwargs['histtype'] = 'stepfilled' + + fig = figlib.phase_comparison_figure( + measured_phases, + true_phases, + plot_residuals=plot_residuals, + f_beacon=f_beacon, + figsize=figsize, + hist_kwargs=hist_kwargs, + ) + + axs = fig.get_axes() if plot_residuals: - phase_residuals = lib.phase_mod(phase_diffs[:,1] - my_phase_diffs) - fig.suptitle("Difference between Measured and Actual phase difference\n for Baselines (i,j" + (')' if not ref_ant_id else '='+str([ int(a.name) for a in ref_ants])+')')) axs[-1].set_xlabel("Baseline Phase Residual $\\Delta\\varphi_{ij_{meas}} - \\Delta\\varphi_{ij_{true}}$ [rad]") else: fig.suptitle("Comparison Measured and Actual phase difference\n for Baselines (i,j" + (')' if not ref_ant_id else '='+str([ int(a.name) for a in ref_ants])+')')) axs[-1].set_xlabel("Baseline Phase $\\Delta\\varphi_{ij}$ [rad]") - + # 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') - + secax = axs[i].child_axes[0] + secax.set_xlabel('Time $\\Delta\\varphi/(2\\pi f_{beac})$ [ns]') + # i=1 axs[i].set_ylabel("Baseline no.") - if plot_residuals: - axs[i].plot(phase_residuals, np.arange(N_base), alpha=0.6, ls='none', marker='x', color=colors[0]) - else: - 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') - - axs[i].legend() - fig.tight_layout() if fig_dir: extra_name = "measured" diff --git a/simulations/airshower_beacon_simulation/bd_antenna_phase_deltas.py b/simulations/airshower_beacon_simulation/bd_antenna_phase_deltas.py index 6e71f6e..afc0f8f 100755 --- a/simulations/airshower_beacon_simulation/bd_antenna_phase_deltas.py +++ b/simulations/airshower_beacon_simulation/bd_antenna_phase_deltas.py @@ -10,6 +10,7 @@ import numpy as np import aa_generate_beacon as beacon import lib +from lib import figlib if __name__ == "__main__": @@ -176,43 +177,36 @@ if __name__ == "__main__": for i in range(2): plot_residuals = i == 1 - colors = ['blue', 'orange'] + true_phases = actual_antenna_phase_shifts + measured_phases = mean_clock_phase - fig, axs = plt.subplots(2, 1, sharex=True, figsize=figsize) + hist_kwargs = {} + if plot_residuals: + measured_phases = lib.phase_mod(measured_phases - actual_antenna_phase_shifts) + hist_kwargs['histtype'] = 'stepfilled' - 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=(phase2time, time2phase)) - secax.set_xlabel('Time $\\Delta\\varphi/(2\\pi f_{beac})$ [ns]') + fig = figlib.phase_comparison_figure( + measured_phases, + true_phases, + plot_residuals=plot_residuals, + f_beacon=f_beacon, + figsize=figsize, + hist_kwargs=hist_kwargs, + ) + + axs = fig.get_axes() if plot_residuals: - phase_residuals = lib.phase_mod(mean_clock_phase - actual_antenna_phase_shifts) fig.suptitle("Difference between Measured and Actual phases (minus global phase)\n for Antenna $i$") axs[-1].set_xlabel("Antenna Mean Phase Residual $\\Delta_\\varphi$") else: fig.suptitle("Comparison Measured and Actual phases (minus global phase)\n for Antenna $i$") axs[-1].set_xlabel("Antenna Mean Phase $\\varphi$") - - i=0 - axs[i].set_ylabel("#") - if plot_residuals: - axs[i].hist(phase_residuals, bins='sqrt', alpha=0.8, color=colors[0]) - else: - axs[i].hist(mean_clock_phase, bins='sqrt', density=False, alpha=0.8, color=colors[0], ls='solid' , histtype='step', label='Measured') - axs[i].hist(actual_antenna_phase_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(phase_residuals, np.arange(N_ant), alpha=0.6, ls='none', marker='x', color=colors[0]) - else: - axs[i].errorbar(mean_clock_phase, np.arange(N_ant), yerr=std_clock_phase, marker='4', alpha=0.7, ls='none', color=colors[0], label='Measured') - axs[i].plot(actual_antenna_phase_shifts, antenna_names, ls='none', marker='3', alpha=0.8, color=colors[1], label='Actual') + #axs[i].errorbar(mean_clock_phase, np.arange(N_ant), yerr=std_clock_phase, marker='4', alpha=0.7, ls='none', color=colors[0], label='Measured') - axs[i].legend() fig.tight_layout() if fig_dir: 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 547216f..7fbf1f0 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 @@ -10,6 +10,7 @@ import numpy as np from os import path import aa_generate_beacon as beacon +from lib import figlib if __name__ == "__main__": import sys @@ -76,15 +77,27 @@ if __name__ == "__main__": for i in range(2): plot_residuals = i == 1 - colors = ['blue', 'orange'] - fig, axs = plt.subplots(2, 1, sharex=True, figsize=figsize) - 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]') + true_phases = actual_time_shifts + measured_phases = measured_time_shifts + + hist_kwargs = {} + if plot_residuals: + measured_phases = measured_phases - true_phases + hist_kwargs['histtype'] = 'stepfilled' + + fig = figlib.phase_comparison_figure( + measured_phases, + true_phases, + plot_residuals=plot_residuals, + f_beacon=f_beacon, + figsize=figsize, + hist_kwargs=hist_kwargs, + secondary_axis='phase', + ) + + axs = fig.get_axes() if plot_residuals: time_shift_residuals = measured_time_shifts - actual_time_shifts @@ -94,26 +107,6 @@ if __name__ == "__main__": 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: diff --git a/simulations/airshower_beacon_simulation/lib/figlib.py b/simulations/airshower_beacon_simulation/lib/figlib.py index bdd8bde..6904ff6 100644 --- a/simulations/airshower_beacon_simulation/lib/figlib.py +++ b/simulations/airshower_beacon_simulation/lib/figlib.py @@ -10,6 +10,7 @@ def phase_comparison_figure( sc_kwargs={}, colors=['blue', 'orange'], legend_on_scatter=True, + secondary_axis='time', **fig_kwargs ): """ @@ -27,11 +28,18 @@ def phase_comparison_figure( if not hasattr(axs, '__len__'): axs = [axs] - if f_beacon: + if f_beacon and secondary_axis in ['phase', 'time']: 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=(phase2time, time2phase)) - secax.set_xlabel('Time $\\varphi/(2\\pi f_{beac})$ [ns]') + + if secondary_axis == 'time': + functions = (phase2time, time2phase) + label = 'Time $\\varphi/(2\\pi f_{beac})$ [ns]' + else: + functions = (time2phase, phase2time) + label = 'Phase $2\\pi t f_{beac}$ [rad]' + + secax = axs[0].secondary_xaxis('top', functions=functions) # Histogram if do_hist_plot: