ericteunis/m-thesis-introduction
1
0
Fork 0
mirror of https://gitlab.science.ru.nl/mthesis-edeboone/m-thesis-introduction.git synced 2024-11-13 10:03:32 +01:00
Code Issues Projects Releases Packages Wiki Activity

ZH: Note SNRs in legends of the plots

Browse source
This commit is contained in:
Eric Teunis de Boone 2023-04-12 22:42:59 +02:00
parent 83811bbd1a
commit b20e2f0f03
8 changed files with 67 additions and 22 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

17
airshower_beacon_simulation/ba_measure_beacon_phase.py
View file

@ -59,6 +59,9 @@ if __name__ == "__main__":
print("Antenna file cannot be found, did you try generating a beacon?")
sys.exit(1)
snrs = beacon.read_snr_file(snr_fname)
snr_str = f"$\\langle SNR \\rangle$ = {snrs['mean']: .1e}"
# read in antennas
with h5py.File(antennas_fname, 'a') as fp:
if 'antennas' not in fp.keys():
@ -231,7 +234,7 @@ if __name__ == "__main__":
ax.axvline(0,color='grey',alpha=0.5)
ax.axhline(0,color='grey',alpha=0.5)
ax.legend()
ax.legend(title=snr_str)
if fig_dir:
old_xlims = ax.get_xlim()
@ -271,7 +274,6 @@ if __name__ == "__main__":
# show histogram of found frequencies
if show_plots or fig_dir:
snrs = beacon.read_snr_file(snr_fname)
if True or allow_frequency_fitting:
fig, ax = plt.subplots(figsize=figsize)
@ -279,16 +281,17 @@ if __name__ == "__main__":
ax.set_ylabel("Counts")
ax.axvline(f_beacon, ls='dashed', color='g')
ax.hist(found_data[:,0], bins='sqrt', density=False)
ax.legend(title=snr_str)
if fig_dir:
fig.savefig(path.join(fig_dir, path.basename(__file__) + f".hist_freq.pdf"))
if True:
fig, _ = figlib.fitted_histogram_figure(found_data[:,2], fit_distr=[], mean_snr=snrs['mean'])
fig, _ = figlib.fitted_histogram_figure(found_data[:,2], fit_distr=['rice'])
ax = fig.axes[0]
ax.set_xlabel("Amplitude")
ax.set_ylabel("Counts")
ax.hist(found_data[:,2], bins='sqrt', density=False)
ax.legend()
ax.legend(title=snr_str)
if fig_dir:
fig.savefig(path.join(fig_dir, path.basename(__file__) + f".hist_amp.pdf"))
@ -299,7 +302,7 @@ if __name__ == "__main__":
ax.set_xlabel("Phase")
ax.set_ylabel("#")
ax.hist(noise_data[:,0], bins='sqrt', density=False)
ax.legend()
ax.legend(title=snr_str)
if fig_dir:
fig.savefig(path.join(fig_dir, path.basename(__file__) + f".noise.hist_phase.pdf"))
@ -313,12 +316,12 @@ if __name__ == "__main__":
fig.savefig(path.join(fig_dir, path.basename(__file__) + f".noise.phase_vs_amp.pdf"))
if True:
fig, _ = figlib.fitted_histogram_figure(noise_data[:,1], fit_distr=['rice', 'rayleigh'], mean_snr=snrs['mean'])
fig, _ = figlib.fitted_histogram_figure(noise_data[:,1], fit_distr=['rice', 'rayleigh'])
ax = fig.axes[0]
ax.set_title("Noise Amplitudes")
ax.set_xlabel("Amplitude [a.u.]")
ax.set_ylabel("#")
ax.legend()
ax.legend(title=snr_str)
if fig_dir:
fig.savefig(path.join(fig_dir, path.basename(__file__) + f".noise.hist_amp.pdf"))

16
airshower_beacon_simulation/bb_measure_clock_phase.py
View file

@ -158,6 +158,7 @@ if __name__ == "__main__":
## Histogram
##
snrs = beacon.read_snr_file(snr_fname)
snr_str = f"$\\langle SNR \\rangle$ = {snrs['mean']: .1e}"
fig = figlib.phase_comparison_figure(
loc_c,
@ -166,15 +167,14 @@ if __name__ == "__main__":
f_beacon=f_beacon,
figsize=figsize,
fit_gaussian=plot_residuals,
mean_snr=snrs['mean']
)
if plot_residuals:
fig.suptitle("Difference between Measured and True Clock phases")
else:
fig.suptitle("Comparison Measured and True Clock Phases")
axs = fig.get_axes()
if plot_residuals:
axs[0].set_title("Difference between Measured and True Clock phases")
else:
axs[0].set_title("Comparison Measured and True Clock Phases")
axs[-1].set_xlabel(f'Antenna {title} {color_label}')
#
@ -186,6 +186,10 @@ if __name__ == "__main__":
i=1
axs[i].set_ylabel("Antenna no.")
#
fig.legend(title=snr_str)
# Save figure
if fig_dir:
fig.tight_layout()

6
airshower_beacon_simulation/bc_baseline_phase_deltas.py
View file

@ -34,6 +34,7 @@ if __name__ == "__main__":
fname_dir = args.data_dir
antennas_fname = path.join(fname_dir, beacon.antennas_fname)
time_diffs_fname = 'time_diffs.hdf5' if False else antennas_fname
snr_fname = path.join(fname_dir, beacon.snr_fname)
fig_dir = args.fig_dir # set None to disable saving
@ -92,6 +93,9 @@ if __name__ == "__main__":
##############################
# Compare actual time shifts #
##############################
snrs = beacon.read_snr_file(snr_fname)
snr_str = f"$\\langle SNR \\rangle$ = {snrs['mean']: .1e}"
actual_antenna_clock_phases = { a.name: -2*np.pi*a.attrs['clock_offset']*f_beacon for a in sorted(antennas, key=lambda a: int(a.name)) }
# Compare actual time shifts
@ -130,6 +134,8 @@ if __name__ == "__main__":
axs = fig.get_axes()
axs[0].legend(title=snr_str)
if plot_residuals:
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]")

15
airshower_beacon_simulation/bd_antenna_phase_deltas.py
View file

@ -37,6 +37,7 @@ if __name__ == "__main__":
antennas_fname = path.join(fname_dir, beacon.antennas_fname)
time_diffs_fname = 'time_diffs.hdf5' if False else antennas_fname
fig_dir = args.fig_dir # set None to disable saving
snr_fname = path.join(fname_dir, beacon.snr_fname)
basenames, time_diffs, f_beacons, clock_phase_diffs, k_periods = beacon.read_baseline_time_diffs_hdf5(time_diffs_fname)
@ -61,6 +62,9 @@ if __name__ == "__main__":
for i, b in enumerate(basenames):
idx = (name2idx(b[0]), name2idx(b[1]))
#if idx[1] < idx[0]:
# idx = (idx[1], idx[0])
clock_phase_matrix[(idx[0], idx[1])] = lib.phase_mod(clock_phase_diffs[i])
clock_phase_matrix[(idx[1], idx[0])] = lib.phase_mod(-1*clock_phase_diffs[i])
@ -171,6 +175,9 @@ if __name__ == "__main__":
##############################
# Compare actual time shifts #
##############################
snrs = beacon.read_snr_file(snr_fname)
snr_str = f"$\\langle SNR \\rangle$ = {snrs['mean']: .1e}"
actual_antenna_time_shifts = { a.name: a.attrs['clock_offset'] for a in sorted(antennas, key=lambda a: int(a.name)) }
if True:
@ -199,16 +206,18 @@ if __name__ == "__main__":
figsize=figsize,
hist_kwargs=hist_kwargs,
fit_gaussian=plot_residuals,
fit_ricianphase=plot_residuals,
return_fit_info = True,
)
axs = fig.get_axes()
axs[0].legend(title=snr_str)
if plot_residuals:
fig.suptitle("Difference between Measured and Actual phases (minus global phase)\n for Antenna $i$")
axs[0].set_title("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[0].set_title("Comparison Measured and Actual phases (minus global phase)\n for Antenna $i$")
axs[-1].set_xlabel("Antenna Mean Phase $\\varphi$")
i=1
@ -273,7 +282,7 @@ if __name__ == "__main__":
ax.plot(np.arange(N_base), (actual_baseline_time_shifts+1/(2*f_beacon))%(1/f_beacon) - 1/(2*f_beacon), ls='none', marker='h', label='actual time shifts', alpha=0.8)
ax.plot(np.arange(N_base), measured_baseline_time_diffs, ls='none', alpha=0.8, marker='x', label='calculated')
ax.legend()
ax.legend(title=snr_str)
if fig_dir:
extra_name = ''
if i == 1:

9
airshower_beacon_simulation/ca_period_from_shower.py
View file

@ -190,6 +190,7 @@ if __name__ == "__main__":
antennas_fname = path.join(fname_dir, beacon.antennas_fname)
time_diffs_fname = 'time_diffs.hdf5' if not True else antennas_fname
tx_fname = path.join(fname_dir, beacon.tx_fname)
snr_fname = path.join(fname_dir, beacon.snr_fname)
## This is a file indicating whether the k-finding algorithm was
## stopped early. This happens when the ks do not change between
@ -209,6 +210,9 @@ if __name__ == "__main__":
ev = REvent(args.input_fname)
# .. patch in our antennas
ev.antennas = antennas
# read in snr information
snrs = beacon.read_snr_file(snr_fname)
snr_str = f"$\\langle SNR \\rangle$ = {snrs['mean']: .1e}"
# For now only implement using one freq_name
freq_names = antennas[0].beacon_info.keys()
@ -300,7 +304,7 @@ if __name__ == "__main__":
ax.plot(ant.t_AxB, calc_beacon, alpha=0.6, ls='dashed', label='Calc Beacon')
ax.plot(ant.t_AxB, ant.E_AxB, alpha=0.6, label="Signal - Calc Beacon")
ax.legend()
ax.legend(title=snr_str)
# save
if fig_dir:
@ -431,6 +435,8 @@ if __name__ == "__main__":
sc = axs.scatter(xx/1e3, yy/1e3, c=maxima_per_loc, cmap='Spectral_r', alpha=0.6)
fig.colorbar(sc, ax=axs, label='Max Amplitude [$\\mu V/m$]')
axs.legend(title=snr_str)
# indicate maximum value
idx = np.argmax(maxima_per_loc)
axs.plot(xx[idx]/1e3, yy[idx]/1e3, 'bx', ms=30)
@ -497,6 +503,7 @@ if __name__ == "__main__":
sc = axs.scatter(xx/1e3, yy/1e3, c=maxima, cmap='Spectral_r', alpha=0.6)
fig.colorbar(sc, ax=axs, label='Max Amplitude [$\\mu V/m$]')
axs.legend(title=snr_str)
if fig_dir:
if power_reconstruction:

8
airshower_beacon_simulation/cb_report_measured_antenna_time_offsets.py
View file

@ -32,6 +32,7 @@ if __name__ == "__main__":
fname_dir = args.data_dir
antennas_fname = path.join(fname_dir, beacon.antennas_fname)
time_diffs_fname = 'time_diffs.hdf5' if not True else antennas_fname
snr_fname = path.join(fname_dir, beacon.snr_fname)
# create fig_dir
if fig_dir:
@ -40,6 +41,10 @@ if __name__ == "__main__":
# Read in antennas from file
_, tx, antennas = beacon.read_beacon_hdf5(antennas_fname)
# Read in snr info
snrs = beacon.read_snr_file(snr_fname)
snr_str = f"$\\langle SNR \\rangle$ = {snrs['mean']: .1e}"
# For now only implement using one freq_name
freq_names = antennas[0].beacon_info.keys()
if len(freq_names) > 1:
@ -100,8 +105,9 @@ if __name__ == "__main__":
axs = fig.get_axes()
axs.legend(title=snr_str)
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:

9
airshower_beacon_simulation/da_reconstruction.py
View file

@ -54,6 +54,7 @@ if __name__ == "__main__":
antennas_fname = path.join(fname_dir, beacon.antennas_fname)
pickle_fname = path.join(fname_dir, 'res.pkl')
tx_fname = path.join(fname_dir, beacon.tx_fname)
snr_fname = path.join(fname_dir, beacon.snr_fname)
# create fig_dir
if fig_dir:
@ -68,6 +69,9 @@ if __name__ == "__main__":
ev = REvent(args.input_fname)
# .. patch in our antennas
ev.antennas = antennas
# Read in snr info
snrs = beacon.read_snr_file(snr_fname)
snr_str = f"$\\langle SNR \\rangle$ = {snrs['mean']: .1e}"
# For now only implement using one freq_name
freq_names = antennas[0].beacon_info.keys()
@ -81,6 +85,7 @@ if __name__ == "__main__":
measured_repair_offsets = beacon.read_antenna_clock_repair_offsets(ev.antennas, mode='phases', freq_name=freq_name)
for i, ant in enumerate(ev.antennas):
# t_AxB will be set by the rit.set_pol_and_bp function
ev.antennas[i].orig_t = ev.antennas[i].t
ev.antennas[i].t += measured_repair_offsets[i]
ev.antennas[i].t_AxB += measured_repair_offsets[i]
@ -167,9 +172,9 @@ if __name__ == "__main__":
ax.plot(time, trace, alpha=0.6, label="Signal - Calc Beacon")
if k == 0:
ax.legend()
ax.legend(title=snr_str)
else:
ax.legend(title="Ex")
ax.legend(title="Ex " + snr_str)
# save
if fig_dir:

9
airshower_beacon_simulation/dc_grid_power_time_fixes.py
View file

@ -28,7 +28,7 @@ def save_overlapping_traces_figure(test_location, ev, N_plot = 30, wx=200, title
axs.set_ylabel("Amplitude [$\\mu V/m$]")
if False:
text_loc = (0.02, 0.95)
axs.text(*text_loc, '[' + ', '.join(['{:.2e}'.format(x) for x in test_location]) + ']', ha='left', transform=axs.transAxes)
axs.text(*text_loc, '[' + ', '.join(['{:.1e}'.format(x) for x in test_location]) + ']', ha='left', transform=axs.transAxes)
a_max = [ np.amax(ant.E_AxB) for ant in ev.antennas ]
power_sort_idx = np.argsort(a_max)
@ -123,6 +123,7 @@ if __name__ == "__main__":
antennas_fname = path.join(fname_dir, beacon.antennas_fname)
pickle_fname = path.join(fname_dir, 'res.pkl')
tx_fname = path.join(fname_dir, beacon.tx_fname)
snr_fname = path.join(fname_dir, beacon.snr_fname)
# create fig_dir
if fig_dir:
@ -136,6 +137,9 @@ if __name__ == "__main__":
bak_ants = ev.antennas
# .. patch in our antennas
ev.antennas = antennas
# Read in snr info
snrs = beacon.read_snr_file(snr_fname)
snr_str = f"$\\langle SNR \\rangle$ = {snrs['mean']: .1e}"
##
## Setup grid
@ -288,6 +292,7 @@ if __name__ == "__main__":
fig.suptitle("")
axs.set_title(suptitle +"\n" +plot_titling[case])
#axs.set_aspect('equal', 'datalim')
axs.legend(title=snr_str)
if fig_dir:
fig.tight_layout()
@ -307,7 +312,7 @@ if __name__ == "__main__":
# only add distance horizontally
location = maxp_loc + np.sqrt(dist*1e3)*np.array([1,1,0])
fig = save_overlapping_traces_figure(location, ev, N_plot=N_plot, wx=wx, title_extra = plot_titling[case] + ', ' + scalename + f', x + {dist}km', fname_distinguish=f'{scalename}.{dist}', fig_dir=fig_dir, figsize=figsize)
fig = save_overlapping_traces_figure(location, ev, N_plot=N_plot, wx=wx, title_extra = plot_titling[case] + ', ' + scalename + f', x + {dist}km', fname_distinguish=f'{scalename}.x{dist}', fig_dir=fig_dir, figsize=figsize)
plt.close(fig)