ZH: add gaussian noise to traces when beaconing

Note that this is a single noise realisation that is added to the three traces.
It will be ~3 times stronger for E_AxB
This commit is contained in:
Eric Teunis de Boone 2023-01-10 12:08:37 +01:00
parent 0447df4f43
commit 265cb16ce2

View file

@ -82,6 +82,8 @@ def Antenna_from_h5ant(h5ant, traces_key='traces', raise_exception=True, read_Ax
ant.Ez = deepcopy(h5ant[traces_key][3]) ant.Ez = deepcopy(h5ant[traces_key][3])
if len(h5ant[traces_key]) > 4: if len(h5ant[traces_key]) > 4:
ant.beacon = deepcopy(h5ant[traces_key][4]) ant.beacon = deepcopy(h5ant[traces_key][4])
if len(h5ant[traces_key]) > 5:
ant.noise = deepcopy(h5ant[traces_key][5])
# E_AxB # E_AxB
if read_AxB and 'E_AxB' in h5ant: if read_AxB and 'E_AxB' in h5ant:
@ -240,6 +242,8 @@ def write_baseline_time_diffs_hdf5(fname, baselines, true_phase_diffs, k_periods
if __name__ == "__main__": if __name__ == "__main__":
from os import path from os import path
rng = np.random.default_rng()
fname = "ZH_airshower/mysim.sry" fname = "ZH_airshower/mysim.sry"
# Transmitter # Transmitter
@ -276,6 +280,9 @@ if __name__ == "__main__":
beacon_amplitudes *= ampl beacon_amplitudes *= ampl
# Noise properties
noise_sigma = 1e-4 # set to None to ignore
# Disable block_filter # Disable block_filter
if False: if False:
block_filter = lambda x, dt, low, high: x block_filter = lambda x, dt, low, high: x
@ -294,6 +301,7 @@ if __name__ == "__main__":
print("Beacon amplitude at tx [muV/m]:", beacon_amplitudes) print("Beacon amplitude at tx [muV/m]:", beacon_amplitudes)
print("Tx location:", [tx.x, tx.y, tx.z]) print("Tx location:", [tx.x, tx.y, tx.z])
print("Noise sigma:", noise_sigma)
# read in antennas # read in antennas
ev = REvent(fname) ev = REvent(fname)
@ -328,13 +336,20 @@ if __name__ == "__main__":
print(f"Modified trace lengths by {pre_N},{after_N-N_samples}") print(f"Modified trace lengths by {pre_N},{after_N-N_samples}")
beacon = lib.beacon_from(tx, antenna, f_beacon, antenna.t, c_light=c_light, radiate_rsq=beacon_radiate_rsq) beacon = lib.beacon_from(tx, antenna, f_beacon, antenna.t, c_light=c_light, radiate_rsq=beacon_radiate_rsq)
traces = np.array([antenna.Ex, antenna.Ey, antenna.Ez, beacon])
# noise realisation
noise_realisation = 0
if noise_sigma is not None:
noise_realisation = rng.normal(0, noise_sigma, size=len(beacon))
# Collect all data to be saved (with the first 3 values the E fields)
traces = np.array([antenna.Ex, antenna.Ey, antenna.Ez, beacon, noise_realisation])
# add to relevant polarisation # add to relevant polarisation
# and apply block filter # and apply block filter
dt = antenna.t[1] - antenna.t[0] dt = antenna.t[1] - antenna.t[0]
for j, amp in enumerate(beacon_amplitudes): for j, amp in enumerate(beacon_amplitudes):
traces[j] = block_filter(traces[j] + amp*beacon, dt, low_bp, high_bp) traces[j] = block_filter(traces[j] + amp*beacon + noise_realisation, dt, low_bp, high_bp)
append_antenna_hdf5( antennas_fname, antenna, traces, name='traces', prepend_time=True) append_antenna_hdf5( antennas_fname, antenna, traces, name='traces', prepend_time=True)