From 39049252039b2b429cb0f837ac77a3877ebc0e63 Mon Sep 17 00:00:00 2001 From: Eric Teunis de Boone Date: Fri, 20 Jan 2023 12:27:50 +0100 Subject: [PATCH] ZH: remove beacon from traces before reconstruction --- .../da_reconstruction.py | 27 ++++++++++++++++++- 1 file changed, 26 insertions(+), 1 deletion(-) diff --git a/simulations/airshower_beacon_simulation/da_reconstruction.py b/simulations/airshower_beacon_simulation/da_reconstruction.py index 3c70a3b..e41ff14 100755 --- a/simulations/airshower_beacon_simulation/da_reconstruction.py +++ b/simulations/airshower_beacon_simulation/da_reconstruction.py @@ -39,10 +39,13 @@ if __name__ == "__main__": fig_subdir = path.join(fig_dir, 'reconstruction') show_plots = args.show_plots + remove_beacon_from_traces = True + #### fname_dir = path.dirname(fname) 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) # create fig_dir if fig_dir: @@ -52,6 +55,7 @@ if __name__ == "__main__": # Read in antennas from file _, tx, antennas = beacon.read_beacon_hdf5(antennas_fname) + _, __, txdata = beacon.read_tx_file(tx_fname) # Read original REvent ev = REvent(fname) # .. patch in our antennas @@ -66,11 +70,32 @@ if __name__ == "__main__": f_beacon = ev.antennas[0].beacon_info[freq_name]['freq'] # Repair clock offsets with the measured offsets - measured_repair_offsets = beacon.read_antenna_clock_repair_offsets(ev.antennas, mode='all', freq_name=freq_name) + 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].t += measured_repair_offsets[i] + # .. and remove the beacon from the traces + if remove_beacon_from_traces: + clock_phase = measured_repair_offsets[i]*2*np.pi*f_beacon + + beacon_phase = ant.beacon_info[freq_name]['beacon_phase'] + f = ant.beacon_info[freq_name]['freq'] + ampl_AxB = ant.beacon_info[freq_name]['amplitude'] + calc_beacon = lib.sine_beacon(f, ev.antennas[i].t, amplitude=ampl_AxB, phase=beacon_phase-clock_phase) + + tx_amps = txdata['amplitudes'] + tx_amps_sum = np.sum(tx_amps) + + # Split up contribution to the various polarisations + for j, amp in enumerate(tx_amps): + if j == 0: + ev.antennas[i].Ex -= amp*(1/tx_amps_sum)*calc_beacon + elif j == 1: + ev.antennas[i].Ey -= amp*(1/tx_amps_sum)*calc_beacon + elif j == 2: + ev.antennas[i].Ez -= amp*(1/tx_amps_sum)*calc_beacon + N_X, Xlow, Xhigh = 23, 100, 1200 with joblib.parallel_backend("loky"): res = rit.reconstruction(ev, outfile=fig_subdir+'/fig.pdf', slice_outdir=fig_subdir+'/', Xlow=Xlow, N_X=N_X, Xhigh=Xhigh)