From 496ed103da12d50349fc229d9633ff1ead1df066 Mon Sep 17 00:00:00 2001 From: Eric Teunis de Boone Date: Tue, 24 Jan 2023 10:58:52 +0100 Subject: [PATCH 1/2] ZH: power on grid for different time fixes The script shows how the power changes after incoporating the various clock offset contributions --- .../airshower_beacon_simulation/Makefile | 2 +- .../dc_grid_power_time_fixes.py | 154 ++++++++++++++++++ 2 files changed, 155 insertions(+), 1 deletion(-) create mode 100755 simulations/airshower_beacon_simulation/dc_grid_power_time_fixes.py diff --git a/simulations/airshower_beacon_simulation/Makefile b/simulations/airshower_beacon_simulation/Makefile index 127768c..7734097 100644 --- a/simulations/airshower_beacon_simulation/Makefile +++ b/simulations/airshower_beacon_simulation/Makefile @@ -27,7 +27,7 @@ findks: reconstruct: ./da_reconstruction.py --no-show-plots --fig-dir=${FIG_DIR} ./db_longitudinal_figure.py --no-show-plots --fig-dir=${FIG_DIR} - + ./dc_grid_power_time_fixes.py --no-show-plots --fig-dir=${FIG_DIR} dist-clean: rm -f ZH_airshower/antennas.hdf5 rm -f ZH_airshower/clocks.csv diff --git a/simulations/airshower_beacon_simulation/dc_grid_power_time_fixes.py b/simulations/airshower_beacon_simulation/dc_grid_power_time_fixes.py new file mode 100755 index 0000000..87b3913 --- /dev/null +++ b/simulations/airshower_beacon_simulation/dc_grid_power_time_fixes.py @@ -0,0 +1,154 @@ +#!/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 + + +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() + + 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() + + wanted_cases = args.figures + + if not wanted_cases or 'all' in wanted_cases: + wanted_cases = valid_cases + + fname = "ZH_airshower/mysim.sry" + + fig_dir = args.fig_dir + show_plots = args.show_plots + + remove_beacon_from_traces = True + + #### + fname_dir = path.dirname(fname) + antennas_fname = path.join(fname_dir, beacon.antennas_fname) + pickle_fname = path.join(fname_dir, 'res.pkl') + + # 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) + # Read original REvent + ev = REvent(fname) + # .. patch in our antennas + ev.antennas = antennas + + rit.set_pol_and_bp(ev) + + ## + ## 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, + } + + # 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 ] + + 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. + if remove_beacon_from_traces: + 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) + + # Only need to manipulate E_AxB + ev.antennas[i].E_AxB -= calc_beacon + + 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 + + # Measure power on grid + 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) + + fig, axs = rit.slice_figure(ev, X, xx, yy, p, mode='sp') + + 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')) From a7c66bb6da5d74724fcb72c4ee92fa0c152d8088 Mon Sep 17 00:00:00 2001 From: Eric Teunis de Boone Date: Mon, 30 Jan 2023 10:24:38 +0100 Subject: [PATCH 2/2] ZH: include trace windowing for power on grid script --- .../aa_generate_beacon.py | 2 +- .../dc_grid_power_time_fixes.py | 86 +++++++++++++++++-- 2 files changed, 79 insertions(+), 9 deletions(-) diff --git a/simulations/airshower_beacon_simulation/aa_generate_beacon.py b/simulations/airshower_beacon_simulation/aa_generate_beacon.py index 2b08252..46e9c41 100755 --- a/simulations/airshower_beacon_simulation/aa_generate_beacon.py +++ b/simulations/airshower_beacon_simulation/aa_generate_beacon.py @@ -29,7 +29,7 @@ def read_antenna_clock_repair_offsets(antennas, mode='all', freq_name=None): # original timing if mode == 'orig': - _clock_delta = ant.attrs['clock_offset'] + _clock_delta = -1*ant.attrs['clock_offset'] # phase if mode in ['all', 'phases']: 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 87b3913..54035b4 100755 --- a/simulations/airshower_beacon_simulation/dc_grid_power_time_fixes.py +++ b/simulations/airshower_beacon_simulation/dc_grid_power_time_fixes.py @@ -50,11 +50,13 @@ if __name__ == "__main__": show_plots = args.show_plots remove_beacon_from_traces = True + apply_signal_window_from_max = True #### fname_dir = path.dirname(fname) 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: @@ -62,13 +64,13 @@ if __name__ == "__main__": # 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(fname) + bak_ants = ev.antennas # .. patch in our antennas ev.antennas = antennas - rit.set_pol_and_bp(ev) - ## ## Setup grid ## @@ -86,10 +88,6 @@ if __name__ == "__main__": 'scale02d': scale02d, } - # 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 ] - plot_titling = { 'no_offset': "no clock offset", 'repair_none': "unrepaired clock offset", @@ -108,16 +106,56 @@ if __name__ == "__main__": # # 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) - # Only need to manipulate E_AxB + # 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 + + # 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 + + for i, ant in enumerate(ev.antennas): + max_idx = np.argmax(ant.E_AxB) + + 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] + + # 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 ] + + ## Apply polarisation and bandpass filter + rit.set_pol_and_bp(ev) + with joblib.parallel_backend("loky"): for case in wanted_cases: print(f"Starting {case} figure") @@ -133,9 +171,41 @@ if __name__ == "__main__": 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 = 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(bak_ants[i].t[0], ev.antennas[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() + 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(power_sort_idx): + if i > N_plot: + break + + alpha = max(0.7, 1/len(a_)) + + axs.plot(t_[idx], a_[idx], color='r', alpha=alpha) + + 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) + # Measure power on grid for scalename, scale in scales.items(): wx, wy = scale, scale