ZH: save figures to dir

This commit is contained in:
Eric Teunis de Boone 2022-12-05 17:48:58 +01:00
parent 2ef87343f5
commit 31392f4fa0
6 changed files with 63 additions and 13 deletions

View file

@ -16,7 +16,8 @@ import lib
if __name__ == "__main__": if __name__ == "__main__":
from os import path from os import path
import sys import sys
import matplotlib
matplotlib.use('Agg')
f_beacon_band = (49e-3,55e-3) #GHz f_beacon_band = (49e-3,55e-3) #GHz
allow_frequency_fitting = False allow_frequency_fitting = False
@ -35,6 +36,8 @@ if __name__ == "__main__":
fname_dir = path.dirname(fname) fname_dir = path.dirname(fname)
antennas_fname = path.join(fname_dir, beacon.antennas_fname) antennas_fname = path.join(fname_dir, beacon.antennas_fname)
fig_dir = "./figures" # set None to disable saving
if not path.isfile(antennas_fname): if not path.isfile(antennas_fname):
print("Antenna file cannot be found, did you try generating a beacon?") print("Antenna file cannot be found, did you try generating a beacon?")
sys.exit(1) sys.exit(1)
@ -140,7 +143,7 @@ if __name__ == "__main__":
# for reporting using plots # for reporting using plots
found_data[i] = frequency, phase, amplitude found_data[i] = frequency, phase, amplitude
if show_plots and (i == 60 or i == 72): if (show_plots or fig_dir) and (i == 60 or i == 72):
fig, ax = plt.subplots() fig, ax = plt.subplots()
ax.set_title(f"Beacon at antenna {h5ant.attrs['name']}\nF:{frequency:.2e}, P:{phase:.4f}, A:{amplitude:.1e}") ax.set_title(f"Beacon at antenna {h5ant.attrs['name']}\nF:{frequency:.2e}, P:{phase:.4f}, A:{amplitude:.1e}")
ax.set_xlabel("t [ns]") ax.set_xlabel("t [ns]")
@ -154,6 +157,9 @@ if __name__ == "__main__":
ax.legend() ax.legend()
if fig_dir:
fig.savefig(path.join(fig_dir, __file__ + f".A{h5ant.attrs['name']}.pdf"))
# save to file # save to file
h5beacon_info = h5ant.require_group('beacon_info') h5beacon_info = h5ant.require_group('beacon_info')
@ -178,18 +184,23 @@ if __name__ == "__main__":
print("Beacon Phases, Amplitudes and Frequencies written to", antennas_fname) print("Beacon Phases, Amplitudes and Frequencies written to", antennas_fname)
# show histogram of found frequencies # show histogram of found frequencies
if show_plots: if show_plots or fig_dir:
if True or allow_frequency_fitting: if True or allow_frequency_fitting:
fig, ax = plt.subplots() fig, ax = plt.subplots()
ax.set_xlabel("Frequency") ax.set_xlabel("Frequency")
ax.set_ylabel("Counts") ax.set_ylabel("Counts")
ax.axvline(f_beacon, ls='dashed', color='g') ax.axvline(f_beacon, ls='dashed', color='g')
ax.hist(found_data[:,0], bins='sqrt', density=False) ax.hist(found_data[:,0], bins='sqrt', density=False)
if fig_dir:
fig.savefig(path.join(fig_dir, __file__ + f".hist_freq.pdf"))
if True: if True:
fig, ax = plt.subplots() fig, ax = plt.subplots()
ax.set_xlabel("Amplitudes") ax.set_xlabel("Amplitudes")
ax.set_ylabel("Counts") ax.set_ylabel("Counts")
ax.hist(found_data[:,2], bins='sqrt', density=False) ax.hist(found_data[:,2], bins='sqrt', density=False)
if fig_dir:
fig.savefig(path.join(fig_dir, __file__ + f".hist_amp.pdf"))
if show_plots:
plt.show() plt.show()

View file

@ -37,6 +37,9 @@ if __name__ == "__main__":
from os import path from os import path
import sys import sys
import matplotlib
matplotlib.use('Agg')
fname = "ZH_airshower/mysim.sry" fname = "ZH_airshower/mysim.sry"
c_light = 3e8*1e-9 # m/ns c_light = 3e8*1e-9 # m/ns
show_plots = True show_plots = True
@ -48,6 +51,8 @@ if __name__ == "__main__":
fname_dir = path.dirname(fname) fname_dir = path.dirname(fname)
antennas_fname = path.join(fname_dir, beacon.antennas_fname) antennas_fname = path.join(fname_dir, beacon.antennas_fname)
fig_dir = "./figures" # set None to disable saving
if not path.isfile(antennas_fname): if not path.isfile(antennas_fname):
print("Antenna file cannot be found, did you try generating a beacon?") print("Antenna file cannot be found, did you try generating a beacon?")
sys.exit(1) sys.exit(1)
@ -89,7 +94,7 @@ if __name__ == "__main__":
h5beacon_freq.attrs['true_phase'] = true_phases[i] h5beacon_freq.attrs['true_phase'] = true_phases[i]
# Plot True Phases at their locations # Plot True Phases at their locations
if show_plots: if show_plots or fig_dir:
fig, ax = plt.subplots() fig, ax = plt.subplots()
spatial_unit=None spatial_unit=None
fig.suptitle('True phases\nf_beacon= {:2.0f}MHz'.format(f_beacon*1e3)) fig.suptitle('True phases\nf_beacon= {:2.0f}MHz'.format(f_beacon*1e3))
@ -106,4 +111,10 @@ if __name__ == "__main__":
sc = ax.scatter(*antenna_locs, c=true_phases, **scatter_kwargs) sc = ax.scatter(*antenna_locs, c=true_phases, **scatter_kwargs)
fig.colorbar(sc, ax=ax, label=color_label) fig.colorbar(sc, ax=ax, label=color_label)
if fig_dir:
fig.savefig(path.join(fig_dir, __file__ + f".F{freq_name}.pdf"))
print(f"True phases written to", antennas_fname)
if show_plots:
plt.show() plt.show()

View file

@ -14,6 +14,9 @@ if __name__ == "__main__":
from os import path from os import path
import sys import sys
import matplotlib
matplotlib.use('Agg')
fname = "ZH_airshower/mysim.sry" fname = "ZH_airshower/mysim.sry"
c_light = 3e8*1e-9 c_light = 3e8*1e-9
show_plots = not True show_plots = not True
@ -24,6 +27,8 @@ if __name__ == "__main__":
antennas_fname = path.join(fname_dir, beacon.antennas_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 not True else antennas_fname
fig_dir = "./figures" # set None to disable saving
if not path.isfile(antennas_fname): if not path.isfile(antennas_fname):
print("Antenna file cannot be found, did you try generating a beacon?") print("Antenna file cannot be found, did you try generating a beacon?")
sys.exit(1) sys.exit(1)
@ -51,7 +56,7 @@ if __name__ == "__main__":
# Get phase difference per baseline # Get phase difference per baseline
phase_diffs = np.empty( (len(baselines), 2) ) phase_diffs = np.empty( (len(baselines), 2) )
for i, base in enumerate(baselines): for i, base in enumerate(baselines):
if i%50==0: if i%100==0:
print(i, "out of", len(baselines)) print(i, "out of", len(baselines))
# read f_beacon from the first antenna # read f_beacon from the first antenna
@ -117,4 +122,8 @@ if __name__ == "__main__":
else: else:
axs[i].plot(phase_residuals, np.arange(N_base), ls='none', marker='x') axs[i].plot(phase_residuals, np.arange(N_base), ls='none', marker='x')
if fig_dir:
fig.savefig(path.join(fig_dir, __file__ + f".F{freq_name}.pdf"))
if show_plots:
plt.show() plt.show()

View file

@ -10,6 +10,7 @@ import h5py
from itertools import combinations, zip_longest, product from itertools import combinations, zip_longest, product
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import numpy as np import numpy as np
from os import path
from scipy.interpolate import interp1d from scipy.interpolate import interp1d
from earsim import REvent from earsim import REvent
@ -18,7 +19,7 @@ import aa_generate_beacon as beacon
import lib import lib
from lib import rit from lib import rit
def find_best_sample_shifts_summing_at_location(test_loc, antennas, allowed_sample_shifts, dt=None, sample_shift_first_trace=0, plot_iteration_with_shifted_trace=None): def find_best_sample_shifts_summing_at_location(test_loc, antennas, allowed_sample_shifts, dt=None, sample_shift_first_trace=0, plot_iteration_with_shifted_trace=None, fig_dir=None, fig_distinguish=None):
""" """
Find the best sample_shift for each antenna by summing the antenna traces Find the best sample_shift for each antenna by summing the antenna traces
and seeing how to get the best alignment. and seeing how to get the best alignment.
@ -63,6 +64,7 @@ def find_best_sample_shifts_summing_at_location(test_loc, antennas, allowed_samp
best_sample_shifts[i] = sample_shift_first_trace best_sample_shifts[i] = sample_shift_first_trace
continue continue
# init figure
if i == plot_iteration_with_shifted_trace: if i == plot_iteration_with_shifted_trace:
fig, ax = plt.subplots() fig, ax = plt.subplots()
ax.set_title("Traces at ({:.1f},{:.1f},{:.1f})".format(*test_loc)) ax.set_title("Traces at ({:.1f},{:.1f},{:.1f})".format(*test_loc))
@ -84,18 +86,26 @@ def find_best_sample_shifts_summing_at_location(test_loc, antennas, allowed_samp
best_sample_shifts[i] = allowed_sample_shifts[best_idx] best_sample_shifts[i] = allowed_sample_shifts[best_idx]
a_sum += np.roll(a_int, best_sample_shifts[i]) a_sum += np.roll(a_int, best_sample_shifts[i])
# cleanup figure
if i == plot_iteration_with_shifted_trace:
if fig_dir: # note this is a global variable here
fig.savefig(path.join(fig_dir, __file__ + '.loc{:.1f}-{:.1f}-{:.1f}'.format(*test_loc) + f'.i{i}{fig_distinguish}.pdf'))
plt.close(fig)
# Return ks # Return ks
return best_sample_shifts, np.max(a_sum) return best_sample_shifts, np.max(a_sum)
if __name__ == "__main__": if __name__ == "__main__":
from os import path
import sys import sys
import matplotlib
matplotlib.use('Agg')
atm = AtmoCal() atm = AtmoCal()
fname = "ZH_airshower/mysim.sry" fname = "ZH_airshower/mysim.sry"
savepath = "./periods_from_shower_figures/" fig_dir = "./figures/periods_from_shower_figures/"
fig_subdir = path.join(fig_dir, 'shifts/')
allowed_ks = np.arange(-2, 3, 1) allowed_ks = np.arange(-2, 3, 1)
Xref = 400 Xref = 400
@ -173,15 +183,19 @@ if __name__ == "__main__":
xx = [] xx = []
yy = [] yy = []
N_loc = len(maxima_per_loc) N_loc = len(maxima_per_loc)
for i, (x_, y_) in enumerate(product(x,y)): for i, (x_, y_) in enumerate(product(x,y)):
tmp_fig_subdir = None
if i % 10 ==0: if i % 10 ==0:
print(f"Testing location {i} out of {N_loc}") print(f"Testing location {i} out of {N_loc}")
tmp_fig_subdir = fig_subdir
test_loc = (x_+xoff)* ev.uAxB + (y_+yoff)*ev.uAxAxB + dXref *ev.uA test_loc = (x_+xoff)* ev.uAxB + (y_+yoff)*ev.uAxAxB + dXref *ev.uA
xx.append(x_+xoff) xx.append(x_+xoff)
yy.append(y_+yoff) yy.append(y_+yoff)
# Find best k for each antenna # Find best k for each antenna
shifts, maximum = find_best_sample_shifts_summing_at_location(test_loc, ev.antennas, allowed_sample_shifts, dt=dt) shifts, maximum = find_best_sample_shifts_summing_at_location(test_loc, ev.antennas, allowed_sample_shifts, dt=dt, fig_dir=tmp_fig_subdir, plot_iteration_with_shifted_trace=len(ev.antennas)-1, fig_distinguish=f".run{r}")
# Translate sample shifts back into period multiple k # Translate sample shifts back into period multiple k
ks = np.rint(shifts*f_beacon*dt) ks = np.rint(shifts*f_beacon*dt)
@ -194,7 +208,7 @@ if __name__ == "__main__":
locs = list(zip(xx, yy)) locs = list(zip(xx, yy))
## Save maxima to file ## Save maxima to file
np.savetxt(savepath + __file__+f'.maxima.run{r}.txt', np.column_stack((locs, maxima_per_loc, ks_per_loc)) ) np.savetxt(fig_dir + __file__+f'.maxima.run{r}.txt', np.column_stack((locs, maxima_per_loc, ks_per_loc)) )
if True: #plot maximum at test locations if True: #plot maximum at test locations
fig, axs = plt.subplots() fig, axs = plt.subplots()
@ -204,11 +218,12 @@ if __name__ == "__main__":
sc = axs.scatter(xx/1e3, yy/1e3, c=maxima_per_loc, cmap='Spectral_r', alpha=0.6) sc = axs.scatter(xx/1e3, yy/1e3, c=maxima_per_loc, cmap='Spectral_r', alpha=0.6)
fig.colorbar(sc, ax=axs) fig.colorbar(sc, ax=axs)
fig.savefig(savepath + __file__+f'.maxima.run{r}.pdf') if fig_dir:
fig.savefig(path.join(fig_dir, __file__+f'.maxima.run{r}.pdf'))
## Save ks to file ## Save ks to file
best_idx = np.argmax(maxima_per_loc) best_idx = np.argmax(maxima_per_loc)
np.savetxt(savepath + __file__+f'.bestk.run{r}.txt', ks_per_loc[best_idx] ) np.savetxt(fig_dir + __file__+f'.bestk.run{r}.txt', ks_per_loc[best_idx] )
print('Best k:', ks_per_loc[best_idx]) print('Best k:', ks_per_loc[best_idx])
# Abort if no improvement # Abort if no improvement

View file

@ -0,0 +1,2 @@
*
!.gitignore

View file

@ -0,0 +1,2 @@
*
!.gitignore