diff --git a/simulations/airshower_beacon_simulation/cc_power_along_axis_reconstruction.py b/simulations/airshower_beacon_simulation/cc_power_along_axis_reconstruction.py index 0e4399b..a8eb731 100755 --- a/simulations/airshower_beacon_simulation/cc_power_along_axis_reconstruction.py +++ b/simulations/airshower_beacon_simulation/cc_power_along_axis_reconstruction.py @@ -7,8 +7,10 @@ clock offsets. """ import matplotlib.pyplot as plt +from mpl_toolkits.mplot3d import Axes3D # required for projection='3d' on old matplotliblib versions import numpy as np from os import path +import pickle from earsim import REvent from atmocal import AtmoCal @@ -31,12 +33,19 @@ if __name__ == "__main__": fig_dir = "./figures/" fig_subdir = path.join(fig_dir, 'reconstruction') show_plots = False + pickle_fname = path.join(fig_subdir, 'res.pkl') #### 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 + # create fig_dir + if fig_dir: + os.makedirs(fig_dir, exist_ok=True) + if fig_subdir: + os.makedirs(fig_subdir, exist_ok=True) + # Read in antennas from file _, tx, antennas = beacon.read_beacon_hdf5(antennas_fname) # Read original REvent @@ -44,6 +53,14 @@ if __name__ == "__main__": # .. patch in our antennas ev.antennas = antennas + # For now only implement using one freq_name + freq_names = antennas[0].beacon_info.keys() + if len(freq_names) > 1: + raise NotImplementedError + + freq_name = next(iter(freq_names)) + f_beacon = ev.antennas[0].beacon_info[freq_name]['freq'] + # Repair clock offsets with the measured offsets for i, ant in enumerate(ev.antennas): clock_phase_time = ant.beacon_info[freq_name]['sigma_phase_mean']/(2*np.pi*f_beacon) @@ -52,12 +69,12 @@ if __name__ == "__main__": total_clock_time = best_k_time + clock_phase_time ev.antennas[i].t += total_clock_time - ## Start a pickle - with open(pickle_fname, 'wb') as fp: - N_X, Xlow, Xhigh = 23, 100, 1200 - res = rit.reconstruction(ev, outfile=fig_subdir+'/fig.pdf', slice_outdir=fig_subdir, Xlow=Xlow, N_X=N_X, Xhigh=Xhigh) + N_X, Xlow, Xhigh = 23, 100, 1200 + res = rit.reconstruction(ev, outfile=fig_subdir+'/fig.pdf', slice_outdir=fig_subdir+'/', Xlow=Xlow, N_X=N_X, Xhigh=Xhigh) + ## Save a pickle + with open(pickle_fname, 'wb') as fp: pickle.dump(res,fp) - if show_plots(): + if show_plots: plt.show()