ericteunis/m-thesis-introduction
1
0
Fork 0
mirror of https://gitlab.science.ru.nl/mthesis-edeboone/m-thesis-introduction.git synced 2024-12-22 11:33:32 +01:00
Code Issues Projects Releases Packages Wiki Activity

ZH: prepare code+figure to remove beacon from signal

Browse source
This commit is contained in:
Eric Teunis de Boone 2022-12-23 11:13:53 +01:00
parent 7851a13065
commit 275775fba3
1 changed files with 58 additions and 8 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

66
simulations/airshower_beacon_simulation/ca_period_from_shower.py
View file

@ -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.))