#!/usr/bin/env python3 # vim: fdm=indent ts=4 """ Do a reconstruction of airshower after correcting for the clock offsets. """ import matplotlib.pyplot as plt import numpy as np from os import path from earsim import REvent from atmocal import AtmoCal import aa_generate_beacon as beacon import lib from lib import rit if __name__ == "__main__": import sys import os import matplotlib if os.name == 'posix' and "DISPLAY" not in os.environ: matplotlib.use('Agg') atm = AtmoCal() fname = "ZH_airshower/mysim.sry" fig_dir = "./figures/" fig_subdir = path.join(fig_dir, 'reconstruction') show_plots = False #### 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 # Read in antennas from file _, tx, antennas = beacon.read_beacon_hdf5(antennas_fname) # Read original REvent ev = REvent(fname) # .. patch in our antennas ev.antennas = antennas # 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) best_k_time = ant.beacon_info[freq_name]['best_k_time'] 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) pickle.dump(res,fp) if show_plots(): plt.show()