diff --git a/simulations/airshower_beacon_simulation/dc_grid_power_time_fixes.py b/simulations/airshower_beacon_simulation/dc_grid_power_time_fixes.py index d1e460d..1d834ed 100755 --- a/simulations/airshower_beacon_simulation/dc_grid_power_time_fixes.py +++ b/simulations/airshower_beacon_simulation/dc_grid_power_time_fixes.py @@ -19,6 +19,65 @@ import lib from lib import rit +def save_overlapping_traces_figure(test_location, ev, N_plot = 30, wx=200, title_extra=None, fname_distinguish='', fig_dir=None, **fig_kwargs): + P, t_, a_, a_sum, t_sum = rit.pow_and_time(test_location, ev, dt=1) + + fig, axs = plt.subplots(**fig_kwargs) + axs.set_title("Antenna traces" + (("\n" + title_extra) if title_extra is not None else '') ) + axs.set_xlabel("Time [ns]") + axs.set_ylabel("Amplitude [$\\mu V/m$]") + if False: + text_loc = (0.02, 0.95) + axs.text(*text_loc, '[' + ', '.join(['{:.2e}'.format(x) for x in test_location]) + ']', ha='left', transform=axs.transAxes) + + a_max = [ np.amax(ant.E_AxB) for ant in ev.antennas ] + power_sort_idx = np.argsort(a_max) + + for i, idx in enumerate(reversed(power_sort_idx)): + if i >= N_plot: + break + + alpha = max(0.4, 1/N_plot) + + axs.plot(t_[idx], a_[idx], color='r', alpha=alpha, lw=2) + + if fig_dir: + if fname_distinguish: + fname_distinguish = "." + fname_distinguish + + fig.tight_layout() + + fig.savefig(path.join(fig_dir, path.basename(__file__) + f'{fname_distinguish}.trace_overlap.{case}.pdf')) + fig.savefig(path.join(fig_dir, path.basename(__file__) + f'{fname_distinguish}.trace_overlap.{case}.png'), transparent=True) + + # Take center between t_low and t_high + if True: + orig_xlims = axs.get_xlim() + + if not True: # t_high and t_low from strongest signal + t_low = np.min(t_[power_sort_idx[-1]]) + t_high = np.max(t_[power_sort_idx[-1]]) + else: # take t_high and t_low from plotted signals + a = [np.min(t_[idx]) for idx in power_sort_idx[-N_plot:]] + t_low = np.nanmin(a) + b = [np.max(t_[idx]) for idx in power_sort_idx[-N_plot:]] + t_high = np.nanmax(b) + + if False: + axs.plot(a, [0]*N_plot, 'gx', ms=10) + axs.plot(b, [0]*N_plot, 'b+', ms=10) + + center_x = (t_high - t_low)/2 + t_low + + low_xlim = max(orig_xlims[0], center_x - wx) + high_xlim = min(orig_xlims[1], center_x + wx) + + axs.set_xlim(low_xlim, high_xlim) + fig.savefig(path.join(fig_dir, path.basename(__file__) + f'{fname_distinguish}.trace_overlap.zoomed.{case}.pdf')) + fig.savefig(path.join(fig_dir, path.basename(__file__) + f'{fname_distinguish}.trace_overlap.zoomed.{case}.png'), transparent=True) + + return fig + if __name__ == "__main__": valid_cases = ['no_offset', 'repair_none', 'repair_phases', 'repair_all'] @@ -95,6 +154,9 @@ if __name__ == "__main__": 'scale02d': scale02d, } + N_plot = 30 + trace_zoom_wx = 100 + plot_titling = { 'no_offset': "no clock offset", 'repair_none': "unrepaired clock offset", @@ -175,6 +237,9 @@ if __name__ == "__main__": backup_antenna_t = [ ant.t for ant in ev.antennas ] backup_antenna_t_AxB = [ ant.t_AxB for ant in ev.antennas ] + fig = save_overlapping_traces_figure([0,0,0], ev, N_plot=1, wx=trace_zoom_wx, title_extra = plot_titling[case], fname_distinguish=f'single', fig_dir=fig_dir, figsize=figsize) + plt.close(fig) + with joblib.parallel_backend("loky"): for case in wanted_cases: print(f"Starting {case} figure") @@ -195,73 +260,29 @@ if __name__ == "__main__": if i == 0: # Specifically compare the times - print(bak_ants[i].t[0], ev.antennas[i].t[0], ev.antennas[i].t[0], ev.antennas[i].attrs['clock_offset'], measured_offsets[i]) + print("backup time, time with measured_offset, true clock offset, measured clock offset") + print(bak_ants[i].t[0], ev.antennas[i].t[0], ev.antennas[i].attrs['clock_offset'], measured_offsets[i]) # # Plot overlapping traces at 0,0,0 # - if True: - P, t_, a_, a_sum, t_sum = rit.pow_and_time([0,0,0], ev, dt=1) - - fig, axs = plt.subplots(figsize=figsize) - axs.set_title("Antenna traces" + "\n" + plot_titling[case]) - axs.set_xlabel("Time [ns]") - axs.set_ylabel("Amplitude [$\\mu V/m$]") - - a_max = [ np.amax(ant.E_AxB) for ant in ev.antennas ] - power_sort_idx = np.argsort(a_max) - - N_plot = 30 - for i, idx in enumerate(reversed(power_sort_idx)): - if i > N_plot: - break - - alpha = max(0.4, 1/len(a_)) - - axs.plot(t_[idx], a_[idx], color='r', alpha=alpha, lw=2) - - if fig_dir: - fig.tight_layout() - - fig.savefig(path.join(fig_dir, path.basename(__file__) + f'.trace_overlap.{case}.pdf')) - fig.savefig(path.join(fig_dir, path.basename(__file__) + f'.trace_overlap.{case}.png'), transparent=True) - - # Take center between t_low and t_high - if True: - orig_xlims = axs.get_xlim() - - if not True: # t_high and t_low from strongest signal - t_low = np.min(t_[power_sort_idx[-1]]) - t_high = np.max(t_[power_sort_idx[-1]]) - else: # take t_high and t_low from plotted signals - a = [np.min(t_[idx]) for idx in power_sort_idx[-N_plot:]] - axs.plot(a, [0]*N_plot, 'gx', ms=10) - t_low = np.nanmin(a) - b = [np.max(t_[idx]) for idx in power_sort_idx[-N_plot:]] - axs.plot(b, [0]*N_plot, 'b+', ms=10) - t_high = np.nanmax(b) - - center_x = (t_high - t_low)/2 + t_low - wx = 200 - - low_xlim = max(orig_xlims[0], center_x - wx) - high_xlim = min(orig_xlims[1], center_x + wx) - - axs.set_xlim(low_xlim, high_xlim) - fig.savefig(path.join(fig_dir, path.basename(__file__) + f'.trace_overlap.zoomed.{case}.pdf')) - fig.savefig(path.join(fig_dir, path.basename(__file__) + f'.trace_overlap.zoomed.{case}.png'), transparent=True) - - - if True: - continue + fig = save_overlapping_traces_figure([0,0,0], ev, N_plot=N_plot, wx=trace_zoom_wx, title_extra = plot_titling[case], fname_distinguish=f'{case}.0', fig_dir=fig_dir, figsize=figsize) + plt.close(fig) # Measure power on grid + # and plot overlapping traces at position with highest power for scalename, scale in scales.items(): wx, wy = scale, scale print(f"Starting grid measurement for figure {case} with {scalename}") - xx, yy, p, maxp = rit.shower_plane_slice(ev, X=X, Nx=Nx, Ny=Nx, wx=wx, wy=wy) + xx, yy, p, maxp_loc = rit.shower_plane_slice(ev, X=X, Nx=Nx, Ny=Nx, wx=wx, wy=wy, zgr=zgr) - fig, axs = rit.slice_figure(ev, X, xx, yy, p, mode='sp') + fig, axs = rit.slice_figure(ev, X, xx, yy, p, mode='sp', scatter_kwargs=dict( + vmax=1e5, + vmin=0, + s=150, + cmap='inferno', +# edgecolor='black', + )) suptitle = fig._suptitle.get_text() fig.suptitle("") @@ -271,7 +292,24 @@ if __name__ == "__main__": if fig_dir: fig.tight_layout() fig.savefig(path.join(fig_dir, path.basename(__file__) + f'.X{X}.{case}.{scalename}.pdf')) + plt.close(fig) + + # + # Plot overlapping traces at highest power of each scale + # + fig = save_overlapping_traces_figure(maxp_loc, ev, N_plot=N_plot, wx=trace_zoom_wx, title_extra = plot_titling[case] + ', ' + scalename + ' best', fname_distinguish=scalename+'.best', fig_dir=fig_dir, figsize=figsize) + + # + # and plot overlapping traces at two other locations + # + if True: + for dist in [ 0.5, 5, 10, 50, 100]: + # only add distance horizontally + location = maxp_loc + np.sqrt(dist*1e3)*np.array([1,1,0]) + + fig = save_overlapping_traces_figure(location, ev, N_plot=N_plot, wx=wx, title_extra = plot_titling[case] + ', ' + scalename + f', x + {dist}km', fname_distinguish=f'{scalename}.{dist}', fig_dir=fig_dir, figsize=figsize) + plt.close(fig) + if args.show_plots: plt.show() -