From 5fca4a8bdaeada7752926fd0294a8867daeeeb3e Mon Sep 17 00:00:00 2001 From: Eric Teunis de Boone Date: Wed, 1 Feb 2023 13:51:43 +0100 Subject: [PATCH] ZH: let view_beacon script plot both raw and filtered traces --- .../airshower_beacon_simulation/Makefile | 2 +- .../aa_generate_beacon.py | 2 +- .../view_beaconed_antenna.py | 120 ++++++++++-------- 3 files changed, 71 insertions(+), 53 deletions(-) diff --git a/simulations/airshower_beacon_simulation/Makefile b/simulations/airshower_beacon_simulation/Makefile index 29f3e3e..d7e3fb1 100644 --- a/simulations/airshower_beacon_simulation/Makefile +++ b/simulations/airshower_beacon_simulation/Makefile @@ -10,7 +10,7 @@ beacon: ./aa_generate_beacon.py | tee figures/aa.log ./ab_modify_clocks.py 0 | tee figures/ab.log ./ac_show_signal_to_noise.py --no-show-plots --fig-dir=${FIG_DIR} - ./view_beaconed_antenna.py 72 -p x -p y -p z -p AxB --no-show-plots --fig-dir=${FIG_DIR} + ./view_beaconed_antenna.py 72 -p x -p y -p z --no-show-plots --fig-dir=${FIG_DIR} clocks: ./ab_modify_clocks.py 15 --gaussian diff --git a/simulations/airshower_beacon_simulation/aa_generate_beacon.py b/simulations/airshower_beacon_simulation/aa_generate_beacon.py index 68b3222..0f77fd2 100755 --- a/simulations/airshower_beacon_simulation/aa_generate_beacon.py +++ b/simulations/airshower_beacon_simulation/aa_generate_beacon.py @@ -400,7 +400,7 @@ if __name__ == "__main__": for j, amp in enumerate(beacon_amplitudes): traces[j] = traces[j] + amp*beacon + noise_realisation - append_antenna_hdf5( antennas_fname, antenna, traces, name='raw_traces', prepend_time=True) + append_antenna_hdf5( antennas_fname, antenna, traces, name='prefiltered_traces', prepend_time=True) # .. and apply block_filter to every trace dt = antenna.t[1] - antenna.t[0] diff --git a/simulations/airshower_beacon_simulation/view_beaconed_antenna.py b/simulations/airshower_beacon_simulation/view_beaconed_antenna.py index c9ca252..06990f8 100755 --- a/simulations/airshower_beacon_simulation/view_beaconed_antenna.py +++ b/simulations/airshower_beacon_simulation/view_beaconed_antenna.py @@ -20,7 +20,7 @@ if __name__ == "__main__": from scriptlib import MyArgumentParser parser = MyArgumentParser() - parser.add_argument('ant_idx', default=[72], nargs='*', help='Antenna Indices') + parser.add_argument('ant_idx', default=[72], nargs='*', type=int, help='Antenna Indices') parser.add_argument('-p', '--polarisations', choices=['x', 'y', 'z', 'b', 'AxB'], action='append', help='Default: x,y,z') parser.add_argument('--geom', action='store_true', help='Make a figure containg the geometry from tx to antenna(s)') parser.add_argument('--ft', action='store_true', help='Add FT strenghts of antenna traces') @@ -28,6 +28,7 @@ if __name__ == "__main__": args = parser.parse_args() fname = "ZH_airshower/mysim.sry" + figsize = (9,6) plot_ft_amplitude = args.ft plot_geometry = args.geom @@ -57,71 +58,88 @@ if __name__ == "__main__": idx = [ i for i, ant in enumerate(antennas) if int(ant.name) in names ] - fig1, axs = plt.subplots(1+plot_ft_amplitude*2) - if not plot_ft_amplitude: - axs = [axs] - axs[0].set_xlabel('t [ns]') - axs[0].set_ylabel('[$\mu$V/m]') - if plot_ft_amplitude: - axs[1].set_xlabel('f [GHz]') - axs[1].set_ylabel('Power') + for i_fig in range(2): + name_dist='' - axs[2].set_ylabel("Phase") - axs[2].set_xlabel('f [GHz]') - axs[2].set_ylim(-np.pi,+np.pi) + if i_fig == 1: #read in the raw_traces + _, __, antennas = beacon.read_beacon_hdf5(antennas_fname, traces_key='prefiltered_traces') + name_dist='.raw' - colorlist = [] - for i in idx: - ant = antennas[i] - n_samples = len(ant.t) - samplerate = (ant.t[-1] - ant.t[0])/n_samples + fig1, axs = plt.subplots(1+plot_ft_amplitude*2, figsize=figsize) + if not plot_ft_amplitude: + axs = [axs] + axs[0].set_xlabel('t [ns]') + axs[0].set_ylabel('[$\mu$V/m]') - axs[0].axvline(ant.t[0], color='k', alpha=0.5) + if i_fig == 1: + axs[0].set_title("UnFiltered traces") + else: + axs[0].set_title("Filtered traces") - mydict = {} - for p in args.polarisations: - pattr = 'E'+str(p) - if p == 'b': - pattr = 'beacon' - elif p == 'AxB': - pattr = 'E_AxB' + if True: + axs[0].set_xlim(-250, 250) - mydict[p] = getattr(ant, pattr) + if plot_ft_amplitude: + axs[1].set_xlabel('f [GHz]') + axs[1].set_ylabel('Power') - for j, (direction, trace) in enumerate(mydict.items()): - l = axs[0].plot(ant.t, trace, label=f"$E_{{{direction}}}$ {ant.name}", alpha=0.7) + axs[2].set_ylabel("Phase") + axs[2].set_xlabel('f [GHz]') + axs[2].set_ylim(-np.pi,+np.pi) - #if False and j == 0 and 't0' in ant.attrs: - # axs[0].axvline(ant.attrs['t0'], color=l[0].get_color(), alpha=0.5) + colorlist = [] + for i in idx: + ant = antennas[i] - colorlist.append(l[0].get_color()) + n_samples = len(ant.t) + samplerate = (ant.t[-1] - ant.t[0])/n_samples - if not plot_ft_amplitude: - continue + axs[0].axvline(ant.t[0], color='k', alpha=0.5) - freqs = ft.fftfreq(n_samples, 1/samplerate)[:n_samples//2] - fft = 2*ft.fft(trace)[:n_samples//2]/n_samples + mydict = {} + for p in args.polarisations: + pattr = 'E'+str(p) + if p == 'b': + pattr = 'beacon' + elif p == 'AxB': + pattr = 'E_AxB' - #axs[1].plot(freqs, np.abs(fft)**2, color=l[0].get_color()) + mydict[p] = getattr(ant, pattr) - if True: - cft = lib.direct_fourier_transform(f_beacon, ant.t, trace) - amp = 2*len(ant.t) * (cft[0]**2 + cft[1]**2) + for j, (direction, trace) in enumerate(mydict.items()): + l = axs[0].plot(ant.t, trace, label=f"$E_{{{direction}}}$ {ant.name}", alpha=0.7) - #axs[0].axhline(amp, color=l[0].get_color()) + #if False and j == 0 and 't0' in ant.attrs: + # axs[0].axvline(ant.attrs['t0'], color=l[0].get_color(), alpha=0.5) - print(amp) - phase = np.arctan2(cft[0],cft[1]) - axs[1].plot(f_beacon, amp, color=l[0].get_color(), marker='3', alpha=0.8, ms=30) - axs[2].plot(f_beacon, phase, color=l[0].get_color(), marker='3', alpha=0.8, ms=30) + colorlist.append(l[0].get_color()) - if plot_ft_amplitude: - fig1.legend(loc='center right', ncol=min(2, len(idx))) - else: - axs[0].legend(loc='upper right', ncol=min(3, len(idx))) - if fig_dir: - fig1.savefig(path.join(fig_dir, path.basename(__file__) + f".trace.pdf")) + if not plot_ft_amplitude: + continue + + freqs = ft.fftfreq(n_samples, 1/samplerate)[:n_samples//2] + fft = 2*ft.fft(trace)[:n_samples//2]/n_samples + + #axs[1].plot(freqs, np.abs(fft)**2, color=l[0].get_color()) + + if True: + cft = lib.direct_fourier_transform(f_beacon, ant.t, trace) + amp = 2*len(ant.t) * (cft[0]**2 + cft[1]**2) + + #axs[0].axhline(amp, color=l[0].get_color()) + + print(amp) + phase = np.arctan2(cft[0],cft[1]) + axs[1].plot(f_beacon, amp, color=l[0].get_color(), marker='3', alpha=0.8, ms=30) + axs[2].plot(f_beacon, phase, color=l[0].get_color(), marker='3', alpha=0.8, ms=30) + + if plot_ft_amplitude: + fig1.legend(loc='center right', ncol=min(2, len(idx))) + else: + axs[0].legend(loc='upper right', ncol=min(3, len(idx))) + if fig_dir: + fig1.savefig(path.join(fig_dir, path.basename(__file__) + f".trace{name_dist}.pdf")) if plot_geometry: if len(mydict) == 1: @@ -130,7 +148,7 @@ if __name__ == "__main__": # only take the colour belonging to mydict[0] geom_colorlist = [ colorlist[len(mydict)*(i)] for i in range(len(colorlist)//len(mydict)) ] - fig2, axs2 = plt.subplots(1) + fig2, axs2 = plt.subplots(1, figsize=figsize) plot_antenna_geometry(antennas, ax=axs2, plot_max_values=False, color='grey', plot_names=False) plot_antenna_geometry([ antennas[i] for i in idx], ax=axs2, colors=geom_colorlist, plot_max_values=False)