mirror of
https://gitlab.science.ru.nl/mthesis-edeboone/m-thesis-introduction.git
synced 2024-11-14 02:23:32 +01:00
Eric Teunis de Boone
7a7647f231
It now produces trace view aswell as power on grid for different positions while looking for the best location
315 lines
12 KiB
Python
Executable file
315 lines
12 KiB
Python
Executable file
#!/usr/bin/env python3
|
|
# vim: fdm=indent ts=4
|
|
|
|
"""
|
|
Show how the Power changes when incorporating the
|
|
various clock offsets by plotting on a grid.
|
|
"""
|
|
|
|
import matplotlib.pyplot as plt
|
|
from mpl_toolkits.mplot3d import Axes3D # required for projection='3d' on old matplotliblib versions
|
|
import numpy as np
|
|
from os import path
|
|
import joblib
|
|
|
|
from earsim import REvent
|
|
from atmocal import AtmoCal
|
|
import aa_generate_beacon as beacon
|
|
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']
|
|
|
|
import sys
|
|
import os
|
|
import matplotlib
|
|
if os.name == 'posix' and "DISPLAY" not in os.environ:
|
|
matplotlib.use('Agg')
|
|
|
|
atm = AtmoCal()
|
|
|
|
from scriptlib import MyArgumentParser
|
|
parser = MyArgumentParser()
|
|
parser.add_argument('--input-fname', type=str, default=None, help='Path to mysim.sry, either directory or path. If empty it takes DATA_DIR and appends mysim.sry. (Default: %(default)s)')
|
|
|
|
group = parser.add_argument_group('figures')
|
|
for case in valid_cases:
|
|
group.add_argument('--'+case.replace('_','-'), dest='figures', action='append_const', const=case)
|
|
|
|
args = parser.parse_args()
|
|
|
|
if not args.input_fname:
|
|
args.input_fname = args.data_dir
|
|
|
|
if path.isdir(args.input_fname):
|
|
args.input_fname = path.join(args.input_fname, "mysim.sry")
|
|
|
|
wanted_cases = args.figures
|
|
|
|
if not wanted_cases or 'all' in wanted_cases:
|
|
wanted_cases = valid_cases
|
|
|
|
figsize = (12,8)
|
|
|
|
fig_dir = args.fig_dir
|
|
show_plots = args.show_plots
|
|
|
|
remove_beacon_from_traces = True
|
|
apply_signal_window_from_max = True
|
|
|
|
####
|
|
fname_dir = args.data_dir
|
|
antennas_fname = path.join(fname_dir, beacon.antennas_fname)
|
|
pickle_fname = path.join(fname_dir, 'res.pkl')
|
|
tx_fname = path.join(fname_dir, beacon.tx_fname)
|
|
|
|
# create fig_dir
|
|
if fig_dir:
|
|
os.makedirs(fig_dir, exist_ok=True)
|
|
|
|
# Read in antennas from file
|
|
_, tx, antennas = beacon.read_beacon_hdf5(antennas_fname)
|
|
_, __, txdata = beacon.read_tx_file(tx_fname)
|
|
# Read original REvent
|
|
ev = REvent(args.input_fname)
|
|
bak_ants = ev.antennas
|
|
# .. patch in our antennas
|
|
ev.antennas = antennas
|
|
|
|
##
|
|
## Setup grid
|
|
##
|
|
X = 400
|
|
zgr = 0 #not exact
|
|
dXref = atm.distance_to_slant_depth(np.deg2rad(ev.zenith),750,zgr+ev.core[2])
|
|
scale2d = dXref*np.tan(np.deg2rad(2.))
|
|
scale4d = dXref*np.tan(np.deg2rad(4.))
|
|
scale02d = dXref*np.tan(np.deg2rad(0.2))
|
|
|
|
Nx, Ny = 21, 21
|
|
scales = {
|
|
'scale2d': scale2d,
|
|
'scale4d': scale4d,
|
|
'scale02d': scale02d,
|
|
}
|
|
|
|
N_plot = 30
|
|
trace_zoom_wx = 100
|
|
|
|
plot_titling = {
|
|
'no_offset': "no clock offset",
|
|
'repair_none': "unrepaired clock offset",
|
|
'repair_phases': "phase resolved clock offsets repaired",
|
|
'repair_all': "final measured clock offsets repaired"
|
|
}
|
|
|
|
# For now only implement using one freq_name
|
|
freq_names = ev.antennas[0].beacon_info.keys()
|
|
if len(freq_names) > 1:
|
|
raise NotImplementedError
|
|
|
|
freq_name = next(iter(freq_names))
|
|
|
|
# Pre remove the beacon from the traces
|
|
#
|
|
# We need to remove it here, so we do not shoot ourselves in
|
|
# the foot when changing to the various clock offsets.
|
|
#
|
|
# Note that the bandpass filter is applied only after E_AxB is
|
|
# reconstructed so we have to manipulate the original traces.
|
|
if remove_beacon_from_traces:
|
|
tx_amps = txdata['amplitudes']
|
|
tx_amps_sum = np.sum(tx_amps)
|
|
|
|
for i, ant in enumerate(ev.antennas):
|
|
beacon_phase = ant.beacon_info[freq_name]['beacon_phase']
|
|
f = ant.beacon_info[freq_name]['freq']
|
|
ampl_AxB = ant.beacon_info[freq_name]['amplitude']
|
|
calc_beacon = lib.sine_beacon(f, ev.antennas[i].t, amplitude=ampl_AxB, phase=beacon_phase)
|
|
|
|
# Split up contribution to the various polarisations
|
|
for j, amp in enumerate(tx_amps):
|
|
if j == 0:
|
|
ev.antennas[i].Ex -= amp*(1/tx_amps_sum)*calc_beacon
|
|
elif j == 1:
|
|
ev.antennas[i].Ey -= amp*(1/tx_amps_sum)*calc_beacon
|
|
elif j == 2:
|
|
ev.antennas[i].Ez -= amp*(1/tx_amps_sum)*calc_beacon
|
|
|
|
# Subtract the beacon from E_AxB
|
|
ev.antennas[i].E_AxB -= calc_beacon
|
|
|
|
# Slice the traces to a small part around the peak
|
|
if apply_signal_window_from_max:
|
|
N_pre, N_post = 250, 250 # TODO: make this configurable
|
|
|
|
for i, ant in enumerate(ev.antennas):
|
|
# Get max idx from all the traces
|
|
# and select the strongest
|
|
max_idx = []
|
|
maxs = []
|
|
|
|
for trace in [ant.Ex, ant.Ey, ant.Ez]:
|
|
idx = np.argmax(np.abs(trace))
|
|
max_idx.append(idx)
|
|
maxs.append( np.abs(trace[idx]) )
|
|
|
|
idx = np.argmax(maxs)
|
|
max_idx = max_idx[idx]
|
|
|
|
|
|
low_idx = max(0, max_idx-N_pre)
|
|
high_idx = min(len(ant.t), max_idx+N_post)
|
|
|
|
ev.antennas[i].t = ant.t[low_idx:high_idx]
|
|
ev.antennas[i].t_AxB = ant.t_AxB[low_idx:high_idx]
|
|
|
|
ev.antennas[i].Ex = ant.Ex[low_idx:high_idx]
|
|
ev.antennas[i].Ey = ant.Ey[low_idx:high_idx]
|
|
ev.antennas[i].Ez = ant.Ez[low_idx:high_idx]
|
|
ev.antennas[i].E_AxB = ant.E_AxB[low_idx:high_idx]
|
|
|
|
## Apply polarisation and bandpass filter
|
|
rit.set_pol_and_bp(ev)
|
|
|
|
# backup antenna times
|
|
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")
|
|
# Repair clock offsets with the measured offsets
|
|
transl_modes = {'no_offset':'orig', 'repair_phases':'phases', 'repair_all':'all'}
|
|
|
|
if case in transl_modes:
|
|
transl_mode = transl_modes[case]
|
|
measured_offsets = beacon.read_antenna_clock_repair_offsets(antennas, mode=transl_mode, freq_name=freq_name)
|
|
else:
|
|
measured_offsets = [0]*len(ev.antennas)
|
|
|
|
for i, ant in enumerate(ev.antennas):
|
|
total_clock_offset = measured_offsets[i]
|
|
|
|
ev.antennas[i].t = backup_antenna_t[i] + total_clock_offset
|
|
ev.antennas[i].t_AxB = backup_antenna_t_AxB[i] + total_clock_offset
|
|
|
|
if i == 0:
|
|
# Specifically compare the times
|
|
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
|
|
#
|
|
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_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', scatter_kwargs=dict(
|
|
vmax=1e5,
|
|
vmin=0,
|
|
s=150,
|
|
cmap='inferno',
|
|
# edgecolor='black',
|
|
))
|
|
|
|
suptitle = fig._suptitle.get_text()
|
|
fig.suptitle("")
|
|
axs.set_title(suptitle +"\n" +plot_titling[case])
|
|
#axs.set_aspect('equal', 'datalim')
|
|
|
|
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()
|