ZH: let view_beacon script plot both raw and filtered traces

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

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]

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)
+    for i_fig in range(2):
-    if not plot_ft_amplitude:
+        name_dist=''
-        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')
 
-        axs[2].set_ylabel("Phase")
+        if i_fig == 1: #read in the raw_traces
-        axs[2].set_xlabel('f [GHz]')
+            _, __, antennas = beacon.read_beacon_hdf5(antennas_fname, traces_key='prefiltered_traces')
-        axs[2].set_ylim(-np.pi,+np.pi)
+            name_dist='.raw'
 
-    colorlist = []
-    for i in idx:
-        ant = antennas[i]
 
-        n_samples = len(ant.t)
+        fig1, axs = plt.subplots(1+plot_ft_amplitude*2, figsize=figsize)
-        samplerate = (ant.t[-1] - ant.t[0])/n_samples
+        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 = {}
+        if True:
-        for p in args.polarisations:
+            axs[0].set_xlim(-250, 250)
-            pattr = 'E'+str(p)
-            if p == 'b':
-                pattr = 'beacon'
-            elif p == 'AxB':
-                pattr = 'E_AxB'
 
-            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()):
+            axs[2].set_ylabel("Phase")
-            l = axs[0].plot(ant.t, trace, label=f"$E_{{{direction}}}$ {ant.name}", alpha=0.7)
+            axs[2].set_xlabel('f [GHz]')
+            axs[2].set_ylim(-np.pi,+np.pi)
 
-            #if False and j == 0 and 't0' in ant.attrs:
+        colorlist = []
-            #    axs[0].axvline(ant.attrs['t0'], color=l[0].get_color(), alpha=0.5)
+        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:
+            axs[0].axvline(ant.t[0], color='k', alpha=0.5)
-                continue
 
-            freqs = ft.fftfreq(n_samples, 1/samplerate)[:n_samples//2]
+            mydict = {}
-            fft = 2*ft.fft(trace)[:n_samples//2]/n_samples
+            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:
+            for j, (direction, trace) in enumerate(mydict.items()):
-                cft = lib.direct_fourier_transform(f_beacon, ant.t, trace)
+                l = axs[0].plot(ant.t, trace, label=f"$E_{{{direction}}}$ {ant.name}", alpha=0.7)
-                amp = 2*len(ant.t) * (cft[0]**2 + cft[1]**2)
 
-                #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)
+                colorlist.append(l[0].get_color())
-                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:
+                if not plot_ft_amplitude:
-        fig1.legend(loc='center right', ncol=min(2, len(idx)))
+                    continue
-    else:
+
-        axs[0].legend(loc='upper right', ncol=min(3, len(idx)))
+                freqs = ft.fftfreq(n_samples, 1/samplerate)[:n_samples//2]
-    if fig_dir:
+                fft = 2*ft.fft(trace)[:n_samples//2]/n_samples
-        fig1.savefig(path.join(fig_dir, path.basename(__file__) + f".trace.pdf"))
+
+                #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)