From 275775fba3ff3e7d6bd7d6ce91b00057e7dd5384 Mon Sep 17 00:00:00 2001 From: Eric Teunis de Boone Date: Fri, 23 Dec 2022 11:13:53 +0100 Subject: [PATCH] ZH: prepare code+figure to remove beacon from signal --- .../ca_period_from_shower.py | 66 ++++++++++++++++--- 1 file changed, 58 insertions(+), 8 deletions(-) diff --git a/simulations/airshower_beacon_simulation/ca_period_from_shower.py b/simulations/airshower_beacon_simulation/ca_period_from_shower.py index c772e90..2b361bb 100755 --- a/simulations/airshower_beacon_simulation/ca_period_from_shower.py +++ b/simulations/airshower_beacon_simulation/ca_period_from_shower.py @@ -120,6 +120,7 @@ if __name__ == "__main__": fig_dir = "./figures/periods_from_shower_figures/" fig_subdir = path.join(fig_dir, 'shifts/') + show_plots = False allowed_ks = np.arange(-2, 3, 1) Xref = 400 @@ -137,6 +138,12 @@ if __name__ == "__main__": 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 _, __, antennas = beacon.read_beacon_hdf5(antennas_fname) # Read original REvent @@ -157,17 +164,60 @@ if __name__ == "__main__": allowed_sample_shifts = np.rint(allowed_ks/f_beacon /dt).astype(int) print("Checking:", allowed_ks, ": shifts :", allowed_sample_shifts) - - # Remove time due to true phase - for i, ant in enumerate(ev.antennas): - true_phase = ant.beacon_info[freq_name]['phase'] - true_phase_time = true_phase/(2*np.pi*f_beacon) - - ev.antennas[i].t -= true_phase_time - # Prepare polarisation and passbands rit.set_pol_and_bp(ev, low=0.03, high=0.08) + # Remove time due to true phase + # and optionally remove the beacon + for i, ant in enumerate(ev.antennas): + clock_phase = ant.beacon_info[freq_name]['sigma_phase_mean'] + clock_phase_time = clock_phase/(2*np.pi*f_beacon) + + ev.antennas[i].orig_t = ev.antennas[i].t + ev.antennas[i].t -= clock_phase_time + + if False: # remove beacon from trace + meas_phase = ant.beacon_info[freq_name]['phase'] + f = ant.beacon_info[freq_name]['freq'] + ampl = ant.beacon_info[freq_name]['amplitude'] + calc_beacon = lib.sine_beacon(f, ev.antennas[i].t, amplitude=ampl, phase=meas_phase+clock_phase) + + ev.antennas[i].E_AxB -= calc_beacon + + # Make a figure of the manipulated traces + if i == 2: + orig_beacon_amplifier = ampl/max(ant.beacon) + + fig, ax = plt.subplots() + ax.set_title(f"Signal and Beacon traces Antenna {i}") + ax.set_xlabel("Time [ns]") + ax.set_ylabel("Amplitude [$\\mu V/m$]") + + ax.plot(ant.t_AxB, ant.E_AxB + calc_beacon, alpha=0.6, ls='dashed', label='Signal') # calc_beacon was already removed + 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() + + # save + if fig_dir: + fig.tight_layout() + if True: # zoom + old_xlim = ax.get_xlim() + + wx, x = 100, 0#ant.t_AxB[np.argmax(ant.E_AxB)] + ax.set_xlim(x-wx, x+wx) + + fig.savefig(path.join(fig_dir, __file__+f'.traces.zoomed.A{i}.pdf')) + + ax.set_xlim(*old_xlim) + + fig.savefig(path.join(fig_dir, __file__+f'.traces.A{i}.pdf')) + + if show_plots: + plt.show() + + dXref = atm.distance_to_slant_depth(np.deg2rad(ev.zenith),Xref,0) scale2d = dXref*np.tan(np.deg2rad(2.))