From daf2209c739d8885ba16a4083ff70f3a8ab28fc0 Mon Sep 17 00:00:00 2001 From: Eric Teunis de Boone Date: Mon, 16 Jan 2023 19:42:21 +0100 Subject: [PATCH 1/5] ZH: new function read_clock_offsets This is moved out from the scripts that already have to modify timing --- .../aa_generate_beacon.py | 29 +++++++++++++++++++ .../ca_period_from_shower.py | 7 ++--- .../da_reconstruction.py | 8 ++--- 3 files changed, 35 insertions(+), 9 deletions(-) diff --git a/simulations/airshower_beacon_simulation/aa_generate_beacon.py b/simulations/airshower_beacon_simulation/aa_generate_beacon.py index d56d232..f3f3321 100755 --- a/simulations/airshower_beacon_simulation/aa_generate_beacon.py +++ b/simulations/airshower_beacon_simulation/aa_generate_beacon.py @@ -20,6 +20,35 @@ tx_fname = 'tx.json' antennas_fname = 'antennas.hdf5' c_light = lib.c_light +def read_antenna_clock_repair_offsets(antennas, mode='all', freq_name=None): + valid_modes = ['orig', 'ks', 'phases', 'all'] + + time_offsets = [] + for i, ant in enumerate(antennas): + _clock_delta = 0 + + # original timing + if mode == 'orig': + _clock_delta = ant.attrs['clock_offset'] + + # phase + if mode in ['all', 'phases']: + clock_phase = ant.beacon_info[freq_name]['sigma_phase_mean'] + f_beacon = ant.beacon_info[freq_name]['freq'] + clock_phase_time = clock_phase/(2*np.pi*f_beacon) + + _clock_delta += clock_phase_time + + # ks + if mode in ['all', 'ks']: + best_k_time = ant.beacon_info[freq_name]['best_k_time'] + + _clock_delta += best_k_time + + time_offsets.append(_clock_delta) + + return time_offsets + def write_tx_file(fname, tx, f_beacon, **kwargs): with open(fname, 'w') as fp: return json.dump( diff --git a/simulations/airshower_beacon_simulation/ca_period_from_shower.py b/simulations/airshower_beacon_simulation/ca_period_from_shower.py index 88654c0..1f0a3c7 100755 --- a/simulations/airshower_beacon_simulation/ca_period_from_shower.py +++ b/simulations/airshower_beacon_simulation/ca_period_from_shower.py @@ -201,14 +201,13 @@ if __name__ == "__main__": # Remove time due to true phase # and optionally remove the beacon + measured_repair_offsets = beacon.read_antenna_clock_repair_offsets(ev.antennas, mode='phases', freq_name=freq_name) for i, ant in enumerate(ev.antennas): - clock_phase = ant.beacon_info[freq_name]['sigma_phase_mean'] - clock_phase_time = clock_phase/(2*np.pi*f_beacon) - ev.antennas[i].orig_t = ev.antennas[i].t_AxB - ev.antennas[i].t_AxB += clock_phase_time + ev.antennas[i].t_AxB += measured_repair_offsets[i] if remove_beacon_from_trace: + clock_phase = measured_repair_offsets[i]*2*np.pi*f_beacon meas_phase = ant.beacon_info[freq_name]['phase'] f = ant.beacon_info[freq_name]['freq'] ampl = ant.beacon_info[freq_name]['amplitude'] diff --git a/simulations/airshower_beacon_simulation/da_reconstruction.py b/simulations/airshower_beacon_simulation/da_reconstruction.py index 6897ef0..3c70a3b 100755 --- a/simulations/airshower_beacon_simulation/da_reconstruction.py +++ b/simulations/airshower_beacon_simulation/da_reconstruction.py @@ -66,12 +66,10 @@ if __name__ == "__main__": f_beacon = ev.antennas[0].beacon_info[freq_name]['freq'] # Repair clock offsets with the measured offsets + measured_repair_offsets = beacon.read_antenna_clock_repair_offsets(ev.antennas, mode='all', freq_name=freq_name) for i, ant in enumerate(ev.antennas): - clock_phase_time = ant.beacon_info[freq_name]['sigma_phase_mean']/(2*np.pi*f_beacon) - best_k_time = ant.beacon_info[freq_name]['best_k_time'] - - total_clock_time = best_k_time + clock_phase_time - ev.antennas[i].t += total_clock_time + # t_AxB will be set by the rit.set_pol_and_bp function + ev.antennas[i].t += measured_repair_offsets[i] N_X, Xlow, Xhigh = 23, 100, 1200 with joblib.parallel_backend("loky"): From 3ba45f4f52407b500c1c0db301ed2ebdd3957072 Mon Sep 17 00:00:00 2001 From: Eric Teunis de Boone Date: Tue, 17 Jan 2023 18:17:12 +0100 Subject: [PATCH 2/5] ZH: renaming phase variables I: phase->beacon_phase --- .../ba_measure_beacon_phase.py | 20 +++++++++---------- .../bb_measure_true_phase.py | 8 ++++---- .../ca_period_from_shower.py | 5 +++-- .../show_beacon_amplitude_antennas.py | 4 ++-- 4 files changed, 19 insertions(+), 18 deletions(-) diff --git a/simulations/airshower_beacon_simulation/ba_measure_beacon_phase.py b/simulations/airshower_beacon_simulation/ba_measure_beacon_phase.py index a1ecd4a..0037897 100755 --- a/simulations/airshower_beacon_simulation/ba_measure_beacon_phase.py +++ b/simulations/airshower_beacon_simulation/ba_measure_beacon_phase.py @@ -136,7 +136,7 @@ if __name__ == "__main__": # Do Fourier Transforms # to find phases and amplitudes if True: - freqs, phases, amps = lib.find_beacon_in_traces( + freqs, beacon_phases, amps = lib.find_beacon_in_traces( test_traces, t_trace, f_beacon_estimate=f_beacon, frequency_fit=allow_frequency_fitting, @@ -147,17 +147,17 @@ if __name__ == "__main__": freqs = [f_beacon] t0 = h5ant.attrs['t0'] - phases = [ 2*np.pi*t0*f_beacon ] + beacon_phases = [ 2*np.pi*t0*f_beacon ] amps = [ 3e-7 ] # choose highest amp idx = 0 - if False and len(phases) > 1: + if False and len(beacon_phases) > 1: #idx = np.argmax(amplitudes, axis=-1) raise NotImplementedError frequency = freqs[idx] - phase = phases[idx] + beacon_phase = beacon_phases[idx] amplitude = amps[idx] orientation = orients[idx] @@ -166,16 +166,16 @@ if __name__ == "__main__": if False: # Subtract phase due to not starting at t=0 # This is already done in beacon_find_traces - phase = lib.phase_mod(phase + corr_phase) + beacon_phase = lib.phase_mod(beacon_phase + corr_phase) # for reporting using plots - found_data[i] = frequency, phase, amplitude + found_data[i] = frequency, beacon_phase, amplitude if (show_plots or fig_dir) and (i == 0 or i == 72 or i == 70): p2t = lambda phase: phase/(2*np.pi*f_beacon) fig, ax = plt.subplots() - ax.set_title(f"Beacon at antenna {h5ant.attrs['name']}\nF:{frequency:.2e}, P:{phase:.4f}, A:{amplitude:.1e}") + ax.set_title(f"Beacon at antenna {h5ant.attrs['name']}\nF:{frequency:.2e}, P:{beacon_phase:.4f}, A:{amplitude:.1e}") ax.set_xlabel("t [ns]") ax.set_ylabel("Amplitude") @@ -191,9 +191,9 @@ if __name__ == "__main__": ax.plot(t_trace - t_0, test_traces[j], marker='.', label='trace '+orients[j]) myt = np.linspace(min(t_trace), max(t_trace), 10*len(t_trace)) - t_0 - ax.plot(myt, lib.sine_beacon(frequency, myt, amplitude=amplitude, t0=0, phase=phase+extra_phase), ls='dotted', label='simulated beacon') + ax.plot(myt, lib.sine_beacon(frequency, myt, amplitude=amplitude, t0=0, phase=beacon_phase+extra_phase), ls='dotted', label='simulated beacon') - ax.axvline( p2t(lib.phase_mod(-1*(phase+extra_phase), low=0)), color='r', ls='dashed', label='$t_\\varphi$') + ax.axvline( p2t(lib.phase_mod(-1*(beacon_phase+extra_phase), low=0)), color='r', ls='dashed', label='$t_\\varphi$') ax.axvline(0,color='grey',alpha=0.5) ax.axhline(0,color='grey',alpha=0.5) @@ -226,7 +226,7 @@ if __name__ == "__main__": h5attrs = h5beacon_freq_info.attrs h5attrs['freq'] = frequency - h5attrs['phase'] = phase + h5attrs['beacon_phase'] = beacon_phase h5attrs['amplitude'] = amplitude h5attrs['orientation'] = orientation diff --git a/simulations/airshower_beacon_simulation/bb_measure_true_phase.py b/simulations/airshower_beacon_simulation/bb_measure_true_phase.py index 80d1817..5a3f2e6 100755 --- a/simulations/airshower_beacon_simulation/bb_measure_true_phase.py +++ b/simulations/airshower_beacon_simulation/bb_measure_true_phase.py @@ -44,20 +44,20 @@ if __name__ == "__main__": # Make sure at least one beacon has been identified if not hasattr(antennas[0], 'beacon_info') or len(antennas[0].beacon_info) == 0: - print(f"No analysed beacon found for {antennas[0].name}, try running the phase analysis script first.") + print(f"No analysed beacon found for {antennas[0].name}, try running the beacon phase analysis script first.") sys.exit(1) # N_beacon_freqs = len(antennas[0].beacon_info) for freq_name in antennas[0].beacon_info.keys(): - measured_phases = np.empty( (len(antennas)) ) + beacon_phases = np.empty( (len(antennas)) ) for i, ant in enumerate(antennas): - measured_phases[i] = ant.beacon_info[freq_name]['phase'] + beacon_phases[i] = ant.beacon_info[freq_name]['beacon_phase'] f_beacon = antennas[0].beacon_info[freq_name]['freq'] - true_phases = lib.remove_antenna_geometry_phase(tx, antennas, f_beacon, measured_phases, c_light=c_light) + true_phases = lib.remove_antenna_geometry_phase(tx, antennas, f_beacon, beacon_phases, c_light=c_light) # Remove the phase from one antenna # this is a free parameter diff --git a/simulations/airshower_beacon_simulation/ca_period_from_shower.py b/simulations/airshower_beacon_simulation/ca_period_from_shower.py index 1f0a3c7..497dc61 100755 --- a/simulations/airshower_beacon_simulation/ca_period_from_shower.py +++ b/simulations/airshower_beacon_simulation/ca_period_from_shower.py @@ -208,10 +208,11 @@ if __name__ == "__main__": if remove_beacon_from_trace: clock_phase = measured_repair_offsets[i]*2*np.pi*f_beacon - meas_phase = ant.beacon_info[freq_name]['phase'] + + beacon_phase = ant.beacon_info[freq_name]['beacon_phase'] f = ant.beacon_info[freq_name]['freq'] ampl = ant.beacon_info[freq_name]['amplitude'] - calc_beacon = lib.sine_beacon(f, ev.antennas[i].t_AxB, amplitude=ampl, phase=meas_phase-clock_phase) + calc_beacon = lib.sine_beacon(f, ev.antennas[i].t_AxB, amplitude=ampl, phase=beacon_phase-clock_phase) ev.antennas[i].E_AxB -= calc_beacon diff --git a/simulations/airshower_beacon_simulation/show_beacon_amplitude_antennas.py b/simulations/airshower_beacon_simulation/show_beacon_amplitude_antennas.py index 4dbd501..2b12b67 100755 --- a/simulations/airshower_beacon_simulation/show_beacon_amplitude_antennas.py +++ b/simulations/airshower_beacon_simulation/show_beacon_amplitude_antennas.py @@ -32,9 +32,9 @@ if __name__ == "__main__": freq_name = list(antennas[0].beacon_info.keys())[0] beacon_frequencies = np.array([ant.beacon_info[freq_name]['freq'] for ant in antennas]) beacon_amplitudes = np.array([ant.beacon_info[freq_name]['amplitude'] for ant in antennas]) - beacon_phases = lib.phase_mod(np.array([ant.beacon_info[freq_name]['phase'] for ant in antennas])) + beacon_phases = lib.phase_mod(np.array([ant.beacon_info[freq_name]['beacon_phase'] for ant in antennas])) - if 'true_phase' in antennas[0].beacon_info[freq_name]: + if False and 'true_phase' in antennas[0].beacon_info[freq_name]: beacon_true_phases = lib.phase_mod(np.array([ant.beacon_info[freq_name]['true_phase'] for ant in antennas])) else: subtitle = " Phases from t0" From 97ebdea54fc356c652e5e8670eff83500286da89 Mon Sep 17 00:00:00 2001 From: Eric Teunis de Boone Date: Thu, 19 Jan 2023 16:24:05 +0100 Subject: [PATCH 3/5] ZH: renaming phase variables II: true_phase->clock_phase --- .../aa_generate_beacon.py | 14 ++++----- .../bb_measure_true_phase.py | 30 +++++++++---------- .../bc_baseline_phase_deltas.py | 20 ++++++------- .../bd_antenna_phase_deltas.py | 6 ++-- .../airshower_beacon_simulation/lib/lib.py | 6 ++-- .../show_beacon_amplitude_antennas.py | 6 ++-- 6 files changed, 41 insertions(+), 41 deletions(-) diff --git a/simulations/airshower_beacon_simulation/aa_generate_beacon.py b/simulations/airshower_beacon_simulation/aa_generate_beacon.py index f3f3321..9eb90eb 100755 --- a/simulations/airshower_beacon_simulation/aa_generate_beacon.py +++ b/simulations/airshower_beacon_simulation/aa_generate_beacon.py @@ -208,13 +208,13 @@ def read_baseline_time_diffs_hdf5(fname): dset = group[dset_name] time_diffs = dset[:,0] f_beacon = dset[:,1] - true_phase_diffs = dset[:,2] + clock_phase_diffs = dset[:,2] k_periods = dset[:,3] - return names, time_diffs, f_beacon, true_phase_diffs, k_periods + return names, time_diffs, f_beacon, clock_phase_diffs, k_periods -def write_baseline_time_diffs_hdf5(fname, baselines, true_phase_diffs, k_periods, f_beacon, time_diffs=None, overwrite=True): +def write_baseline_time_diffs_hdf5(fname, baselines, clock_phase_diffs, k_periods, f_beacon, time_diffs=None, overwrite=True): """ Write a combination of baselines, phase_diff, k_period and f_beacon to file. @@ -225,7 +225,7 @@ def write_baseline_time_diffs_hdf5(fname, baselines, true_phase_diffs, k_periods # this is a single baseline N_baselines = 1 baselines = [baselines] - true_phase_diffs = [true_phase_diffs] + clock_phase_diffs = [clock_phase_diffs] k_periods = [k_periods] f_beacon = np.array([f_beacon]) @@ -237,9 +237,9 @@ def write_baseline_time_diffs_hdf5(fname, baselines, true_phase_diffs, k_periods f_beacon = np.array([f_beacon]*N_baselines) if time_diffs is None: - time_diffs = k_periods/f_beacon + true_phase_diffs/(2*np.pi*f_beacon) + time_diffs = k_periods/f_beacon + clock_phase_diffs/(2*np.pi*f_beacon) - assert len(baselines) == len(true_phase_diffs) == len(k_periods) == len(f_beacon) + assert len(baselines) == len(clock_phase_diffs) == len(k_periods) == len(f_beacon) with h5py.File(fname, 'a') as fp: group_name = 'baseline_time_diffs' @@ -263,7 +263,7 @@ def write_baseline_time_diffs_hdf5(fname, baselines, true_phase_diffs, k_periods base_dset = group.create_dataset(base_dset_name, data=basenames) - data = np.vstack( (time_diffs, f_beacon, true_phase_diffs, k_periods) ).T + data = np.vstack( (time_diffs, f_beacon, clock_phase_diffs, k_periods) ).T dset = group.create_dataset(dset_name, data=data) # }}} vim marker diff --git a/simulations/airshower_beacon_simulation/bb_measure_true_phase.py b/simulations/airshower_beacon_simulation/bb_measure_true_phase.py index 5a3f2e6..59afe8d 100755 --- a/simulations/airshower_beacon_simulation/bb_measure_true_phase.py +++ b/simulations/airshower_beacon_simulation/bb_measure_true_phase.py @@ -57,19 +57,19 @@ if __name__ == "__main__": f_beacon = antennas[0].beacon_info[freq_name]['freq'] - true_phases = lib.remove_antenna_geometry_phase(tx, antennas, f_beacon, beacon_phases, c_light=c_light) + clock_phases = lib.remove_antenna_geometry_phase(tx, antennas, f_beacon, beacon_phases, c_light=c_light) # Remove the phase from one antenna # this is a free parameter # (only required for absolute timing) if remove_absolute_phase_offset_first_antenna or remove_absolute_phase_offset_minimum: if remove_absolute_phase_offset_first_antenna: # just take the first phase - minimum_phase = true_phases[0] + minimum_phase = clock_phases[0] else: # take the minimum - minimum_phase = np.min(true_phases, axis=-1) + minimum_phase = np.min(clock_phases, axis=-1) - true_phases -= minimum_phase - true_phases = lib.phase_mod(true_phases) + clock_phases -= minimum_phase + clock_phases = lib.phase_mod(clock_phases) # Save to antennas in file with h5py.File(antennas_fname, 'a') as fp: @@ -79,7 +79,7 @@ if __name__ == "__main__": h5ant = fp['antennas'][ant.name] h5beacon_freq = h5ant['beacon_info'][freq_name] - h5beacon_freq.attrs['true_phase'] = true_phases[i] + h5beacon_freq.attrs['clock_phase'] = clock_phases[i] # Plot True Phases at their locations if show_plots or fig_dir: @@ -96,7 +96,7 @@ if __name__ == "__main__": #scatter_kwargs['vmax'] = +np.pi color_label='$\\varphi(\\sigma_t)$ [rad]' - sc = ax.scatter(*antenna_locs, c=true_phases, **scatter_kwargs) + sc = ax.scatter(*antenna_locs, c=clock_phases, **scatter_kwargs) fig.colorbar(sc, ax=ax, label=color_label) if False: @@ -114,20 +114,20 @@ if __name__ == "__main__": fig, ax = plt.subplots() fig.suptitle('Clock phase Residuals\nf_beacon={:2.0f}MHz'.format(f_beacon*1e3)) - actual_true_phases = np.array([ -2*np.pi*a.attrs['clock_offset']*f_beacon for a in antennas ]) + actual_clock_phases = np.array([ -2*np.pi*a.attrs['clock_offset']*f_beacon for a in antennas ]) - # Modify actual_true_phases, the same way as true_phases + # Modify actual_clock_phases, the same way as clock_phases # was modified if remove_absolute_phase_offset_first_antenna or remove_absolute_phase_offset_minimum: if remove_absolute_phase_offset_first_antenna: # just take the first phase - minimum_phase = actual_true_phases[0] + minimum_phase = actual_clock_phases[0] else: # take the minimum - minimum_phase = np.min(actual_true_phases, axis=-1) + minimum_phase = np.min(actual_clock_phases, axis=-1) - actual_true_phases -= minimum_phase - actual_true_phases = lib.phase_mod(actual_true_phases) + actual_clock_phases -= minimum_phase + actual_clock_phases = lib.phase_mod(actual_clock_phases) - true_phase_residuals = lib.phase_mod(true_phases - actual_true_phases) + clock_phase_residuals = lib.phase_mod(clock_phases - actual_clock_phases) antenna_locs = list(zip(*[(ant.x, ant.y) for ant in antennas])) ax.set_xlabel('x' if spatial_unit is None else 'x [{}]'.format(spatial_unit)) @@ -136,7 +136,7 @@ if __name__ == "__main__": scatter_kwargs['cmap'] = 'inferno' color_label='$\\Delta\\varphi(\\sigma_t) = \\varphi_{meas} - \\varphi_{true}$ [rad]' - sc = ax.scatter(*antenna_locs, c=true_phase_residuals, **scatter_kwargs) + sc = ax.scatter(*antenna_locs, c=clock_phase_residuals, **scatter_kwargs) fig.colorbar(sc, ax=ax, label=color_label) if fig_dir: diff --git a/simulations/airshower_beacon_simulation/bc_baseline_phase_deltas.py b/simulations/airshower_beacon_simulation/bc_baseline_phase_deltas.py index 725aa4b..ecb1861 100755 --- a/simulations/airshower_beacon_simulation/bc_baseline_phase_deltas.py +++ b/simulations/airshower_beacon_simulation/bc_baseline_phase_deltas.py @@ -70,30 +70,30 @@ if __name__ == "__main__": # Get true phase diffs try: - true_phases = np.array([ant.beacon_info[freq_name]['true_phase'] for ant in base]) - true_phases_diff = lib.phase_mod(lib.phase_mod(true_phases[1]) - lib.phase_mod(true_phases[0])) + clock_phases = np.array([ant.beacon_info[freq_name]['clock_phase'] for ant in base]) + clock_phases_diff = lib.phase_mod(lib.phase_mod(clock_phases[1]) - lib.phase_mod(clock_phases[0])) except IndexError: - # true_phase not determined yet - print(f"Missing true_phases for {freq_name} in baseline {base[0].name},{base[1].name}") - true_phases_diff = np.nan + # clock_phase not determined yet + print(f"Missing clock_phases for {freq_name} in baseline {base[0].name},{base[1].name}") + clock_phases_diff = np.nan # save phase difference with antenna names - phase_diffs[i] = [f_beacon, true_phases_diff] + phase_diffs[i] = [f_beacon, clock_phases_diff] beacon.write_baseline_time_diffs_hdf5(time_diffs_fname, baselines, phase_diffs[:,1], [0]*len(phase_diffs), phase_diffs[:,0]) ############################## # Compare actual time shifts # ############################## - actual_antenna_true_phases = { a.name: -2*np.pi*a.attrs['clock_offset']*f_beacon for a in sorted(antennas, key=lambda a: int(a.name)) } + actual_antenna_clock_phases = { a.name: -2*np.pi*a.attrs['clock_offset']*f_beacon for a in sorted(antennas, key=lambda a: int(a.name)) } # Compare actual time shifts my_phase_diffs = [] for i,b in enumerate(baselines): - actual_true_phase_diff = lib.phase_mod( lib.phase_mod(actual_antenna_true_phases[b[1].name]) - lib.phase_mod(actual_antenna_true_phases[b[0].name])) + actual_clock_phase_diff = lib.phase_mod( lib.phase_mod(actual_antenna_clock_phases[b[1].name]) - lib.phase_mod(actual_antenna_clock_phases[b[0].name])) - this_actual_true_phase_diff = lib.phase_mod( actual_true_phase_diff ) - my_phase_diffs.append(this_actual_true_phase_diff) + this_actual_clock_phase_diff = lib.phase_mod( actual_clock_phase_diff ) + my_phase_diffs.append(this_actual_clock_phase_diff) # Make a plot if True: diff --git a/simulations/airshower_beacon_simulation/bd_antenna_phase_deltas.py b/simulations/airshower_beacon_simulation/bd_antenna_phase_deltas.py index da1da20..809edb7 100755 --- a/simulations/airshower_beacon_simulation/bd_antenna_phase_deltas.py +++ b/simulations/airshower_beacon_simulation/bd_antenna_phase_deltas.py @@ -36,7 +36,7 @@ if __name__ == "__main__": time_diffs_fname = 'time_diffs.hdf5' if False else antennas_fname fig_dir = args.fig_dir # set None to disable saving - basenames, time_diffs, f_beacons, true_phase_diffs, k_periods = beacon.read_baseline_time_diffs_hdf5(time_diffs_fname) + basenames, time_diffs, f_beacons, clock_phase_diffs, k_periods = beacon.read_baseline_time_diffs_hdf5(time_diffs_fname) f_beacon, tx, antennas = beacon.read_beacon_hdf5(antennas_fname) @@ -59,8 +59,8 @@ if __name__ == "__main__": for i, b in enumerate(basenames): idx = (name2idx(b[0]), name2idx(b[1])) - sigma_phase_matrix[(idx[0], idx[1])] = lib.phase_mod(true_phase_diffs[i]) - sigma_phase_matrix[(idx[1], idx[0])] = lib.phase_mod(-1*true_phase_diffs[i]) + sigma_phase_matrix[(idx[0], idx[1])] = lib.phase_mod(clock_phase_diffs[i]) + sigma_phase_matrix[(idx[1], idx[0])] = lib.phase_mod(-1*clock_phase_diffs[i]) mat_kwargs = dict( norm = Normalize(vmin=-np.pi, vmax=+np.pi), diff --git a/simulations/airshower_beacon_simulation/lib/lib.py b/simulations/airshower_beacon_simulation/lib/lib.py index aa54b79..0d5d7db 100644 --- a/simulations/airshower_beacon_simulation/lib/lib.py +++ b/simulations/airshower_beacon_simulation/lib/lib.py @@ -70,16 +70,16 @@ def remove_antenna_geometry_phase(tx, antennas, f_beacon, measured_phases=None, if not hasattr(measured_phases, '__len__'): measured_phases = [measured_phases] - true_phases = np.empty( (len(antennas)) ) + clock_phases = np.empty( (len(antennas)) ) for i, ant in enumerate(antennas): measured_phase = measured_phases[i] geom_time = geometry_time(tx, ant, c_light=c_light) geom_phase = geom_time * 2*np.pi*f_beacon - true_phases[i] = phase_mod(measured_phase - geom_phase) + clock_phases[i] = phase_mod(measured_phase - geom_phase) - return true_phases + return clock_phases """ Fourier """ diff --git a/simulations/airshower_beacon_simulation/show_beacon_amplitude_antennas.py b/simulations/airshower_beacon_simulation/show_beacon_amplitude_antennas.py index 2b12b67..f880761 100755 --- a/simulations/airshower_beacon_simulation/show_beacon_amplitude_antennas.py +++ b/simulations/airshower_beacon_simulation/show_beacon_amplitude_antennas.py @@ -34,8 +34,8 @@ if __name__ == "__main__": beacon_amplitudes = np.array([ant.beacon_info[freq_name]['amplitude'] for ant in antennas]) beacon_phases = lib.phase_mod(np.array([ant.beacon_info[freq_name]['beacon_phase'] for ant in antennas])) - if False and 'true_phase' in antennas[0].beacon_info[freq_name]: - beacon_true_phases = lib.phase_mod(np.array([ant.beacon_info[freq_name]['true_phase'] for ant in antennas])) + if False and 'clock_phase' in antennas[0].beacon_info[freq_name]: + beacon_clock_phases = lib.phase_mod(np.array([ant.beacon_info[freq_name]['clock_phase'] for ant in antennas])) else: subtitle = " Phases from t0" beacon_frequencies = np.array([ f_beacon for ant in antennas ]) @@ -49,7 +49,7 @@ if __name__ == "__main__": colorlabel = '$\\varphi$' sizes = 64*(beacon_amplitudes/np.max(beacon_amplitudes))**2 elif True: # True Phases - vals = beacon_true_phases + vals = beacon_clock_phases colorlabel = '$\\sigma_\\varphi$' plot_phase_field = False plot_tx = False From 6be1bb129fc38cbc9861e4030782300c4ec81f4a Mon Sep 17 00:00:00 2001 From: Eric Teunis de Boone Date: Thu, 19 Jan 2023 17:05:50 +0100 Subject: [PATCH 4/5] ZH: renaming phase variables II: true_phase->clock_phase (file moving) --- simulations/airshower_beacon_simulation/Makefile | 6 +++--- .../{bb_measure_true_phase.py => bb_measure_clock_phase.py} | 0 2 files changed, 3 insertions(+), 3 deletions(-) rename simulations/airshower_beacon_simulation/{bb_measure_true_phase.py => bb_measure_clock_phase.py} (100%) diff --git a/simulations/airshower_beacon_simulation/Makefile b/simulations/airshower_beacon_simulation/Makefile index 7528db2..127768c 100644 --- a/simulations/airshower_beacon_simulation/Makefile +++ b/simulations/airshower_beacon_simulation/Makefile @@ -7,8 +7,8 @@ all: beacon clocks phases findks reconstruct beacon: - ./aa_generate_beacon.py > figures/aa.log - ./ab_modify_clocks.py 0 > figures/ab.log + ./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} clocks: @@ -16,7 +16,7 @@ clocks: phases: ./ba_measure_beacon_phase.py --no-show-plots --fig-dir=${FIG_DIR} - ./bb_measure_true_phase.py --no-show-plots --fig-dir=${FIG_DIR} + ./bb_measure_clock_phase.py --no-show-plots --fig-dir=${FIG_DIR} ./bc_baseline_phase_deltas.py --no-show-plots --fig-dir=${FIG_DIR} ./bd_antenna_phase_deltas.py --no-show-plots --fig-dir=${FIG_DIR} diff --git a/simulations/airshower_beacon_simulation/bb_measure_true_phase.py b/simulations/airshower_beacon_simulation/bb_measure_clock_phase.py similarity index 100% rename from simulations/airshower_beacon_simulation/bb_measure_true_phase.py rename to simulations/airshower_beacon_simulation/bb_measure_clock_phase.py From 75999e6eb361befeb46446882da6d0a48b8e9432 Mon Sep 17 00:00:00 2001 From: Eric Teunis de Boone Date: Thu, 19 Jan 2023 17:13:29 +0100 Subject: [PATCH 5/5] ZH: renaming phase variables III: sigma_phase_*->clock_phase_* --- .../aa_generate_beacon.py | 2 +- .../bd_antenna_phase_deltas.py | 48 +++++++++---------- ...cb_report_measured_antenna_time_offsets.py | 2 +- 3 files changed, 26 insertions(+), 26 deletions(-) diff --git a/simulations/airshower_beacon_simulation/aa_generate_beacon.py b/simulations/airshower_beacon_simulation/aa_generate_beacon.py index 9eb90eb..2b08252 100755 --- a/simulations/airshower_beacon_simulation/aa_generate_beacon.py +++ b/simulations/airshower_beacon_simulation/aa_generate_beacon.py @@ -33,7 +33,7 @@ def read_antenna_clock_repair_offsets(antennas, mode='all', freq_name=None): # phase if mode in ['all', 'phases']: - clock_phase = ant.beacon_info[freq_name]['sigma_phase_mean'] + clock_phase = ant.beacon_info[freq_name]['clock_phase_mean'] f_beacon = ant.beacon_info[freq_name]['freq'] clock_phase_time = clock_phase/(2*np.pi*f_beacon) diff --git a/simulations/airshower_beacon_simulation/bd_antenna_phase_deltas.py b/simulations/airshower_beacon_simulation/bd_antenna_phase_deltas.py index 809edb7..8f45ed7 100755 --- a/simulations/airshower_beacon_simulation/bd_antenna_phase_deltas.py +++ b/simulations/airshower_beacon_simulation/bd_antenna_phase_deltas.py @@ -48,19 +48,19 @@ if __name__ == "__main__": N_ant = len(antennas) # reshape time_diffs into N_ant x N_ant matrix - sigma_phase_matrix = np.full( (N_ant, N_ant), np.nan, dtype=float) + clock_phase_matrix = np.full( (N_ant, N_ant), np.nan, dtype=float) ## set i=i terms to 0 for i in range(N_ant): - sigma_phase_matrix[i,i] = 0 + clock_phase_matrix[i,i] = 0 ## fill matrix name2idx = lambda name: int(name)-1 for i, b in enumerate(basenames): idx = (name2idx(b[0]), name2idx(b[1])) - sigma_phase_matrix[(idx[0], idx[1])] = lib.phase_mod(clock_phase_diffs[i]) - sigma_phase_matrix[(idx[1], idx[0])] = lib.phase_mod(-1*clock_phase_diffs[i]) + clock_phase_matrix[(idx[0], idx[1])] = lib.phase_mod(clock_phase_diffs[i]) + clock_phase_matrix[(idx[1], idx[0])] = lib.phase_mod(-1*clock_phase_diffs[i]) mat_kwargs = dict( norm = Normalize(vmin=-np.pi, vmax=+np.pi), @@ -74,7 +74,7 @@ if __name__ == "__main__": ax.set_ylabel("Antenna no. i") ax.set_xlabel("Antenna no. j") - im = ax.imshow(sigma_phase_matrix, interpolation='none', **mat_kwargs) + im = ax.imshow(clock_phase_matrix, interpolation='none', **mat_kwargs) fig.colorbar(im, ax=ax) if fig_dir: @@ -88,7 +88,7 @@ if __name__ == "__main__": if True: # for each row j subtract the 0,j element from the whole row # and apply phase_mod - first_row = -1*(sigma_phase_matrix[0,:] * np.ones_like(sigma_phase_matrix)).T + first_row = -1*(clock_phase_matrix[0,:] * np.ones_like(clock_phase_matrix)).T # Show subtraction Matrix as figure if True: @@ -105,18 +105,18 @@ if __name__ == "__main__": plt.close(fig) - sigma_phase_matrix = sigma_phase_matrix - first_row - sigma_phase_matrix = lib.phase_mod(sigma_phase_matrix) + clock_phase_matrix = clock_phase_matrix - first_row + clock_phase_matrix = lib.phase_mod(clock_phase_matrix) # Except for the first row, these are all separate measurements # Condense into phase offset per antenna if True: # do not use the first row - my_mask = np.arange(1, len(sigma_phase_matrix), dtype=int) + my_mask = np.arange(1, len(clock_phase_matrix), dtype=int) else: - my_mask = np.arange(0, len(sigma_phase_matrix), dtype=int) + my_mask = np.arange(0, len(clock_phase_matrix), dtype=int) - mean_sigma_phase = np.nanmean(sigma_phase_matrix[my_mask], axis=0) - std_sigma_phase = np.nanstd( sigma_phase_matrix[my_mask], axis=0) + mean_clock_phase = np.nanmean(clock_phase_matrix[my_mask], axis=0) + std_clock_phase = np.nanstd( clock_phase_matrix[my_mask], axis=0) # Show resulting matrix as figure @@ -126,14 +126,14 @@ if __name__ == "__main__": axs[0].set_ylabel("Antenna no. 0") axs[-1].set_xlabel("Antenna no. j") - im = axs[0].imshow(sigma_phase_matrix, interpolation='none', **mat_kwargs) + im = axs[0].imshow(clock_phase_matrix, interpolation='none', **mat_kwargs) fig.colorbar(im, ax=axs) axs[0].set_aspect('auto') - colours = [mat_kwargs['cmap'](mat_kwargs['norm'](x)) for x in mean_sigma_phase] + colours = [mat_kwargs['cmap'](mat_kwargs['norm'](x)) for x in mean_clock_phase] axs[1].set_ylabel("Mean Baseline Phase (0, j)[rad]") - axs[1].errorbar(np.arange(N_ant), mean_sigma_phase, yerr=std_sigma_phase, ls='none') - axs[1].scatter(np.arange(N_ant), mean_sigma_phase, c=colours,s=4) + axs[1].errorbar(np.arange(N_ant), mean_clock_phase, yerr=std_clock_phase, ls='none') + axs[1].scatter(np.arange(N_ant), mean_clock_phase, c=colours,s=4) if fig_dir: fig.savefig(path.join(fig_dir, path.basename(__file__) + f".matrix.modified.pdf")) @@ -157,8 +157,8 @@ if __name__ == "__main__": h5attrs = h5beacon_info[freq_name].attrs idx = name2idx(ant.name) - h5attrs['sigma_phase_mean'] = mean_sigma_phase[idx] - h5attrs['sigma_phase_std'] = std_sigma_phase[idx] + h5attrs['clock_phase_mean'] = mean_clock_phase[idx] + h5attrs['clock_phase_std'] = std_clock_phase[idx] ############################## @@ -171,7 +171,7 @@ if __name__ == "__main__": antenna_names = [int(k)-1 for k,v in actual_antenna_time_shifts.items() ] # Make sure to shift all antennas by a global phase - global_phase_shift = actual_antenna_phase_shifts[0] - mean_sigma_phase[0] + global_phase_shift = actual_antenna_phase_shifts[0] - mean_clock_phase[0] actual_antenna_phase_shifts = lib.phase_mod(actual_antenna_phase_shifts - global_phase_shift ) for i in range(2): @@ -187,7 +187,7 @@ if __name__ == "__main__": secax.set_xlabel('Time $\\Delta\\varphi/(2\\pi f_{beac})$ [ns]') if plot_residuals: - phase_residuals = lib.phase_mod(mean_sigma_phase - actual_antenna_phase_shifts) + phase_residuals = lib.phase_mod(mean_clock_phase - actual_antenna_phase_shifts) fig.suptitle("Difference between Measured and Actual phases (minus global phase)\n for Antenna $i$") axs[-1].set_xlabel("Antenna Phase Residual $\\Delta_\\varphi$") else: @@ -200,7 +200,7 @@ if __name__ == "__main__": if plot_residuals: axs[i].hist(phase_residuals, bins='sqrt', alpha=0.8, color=colors[0]) else: - axs[i].hist(mean_sigma_phase, bins='sqrt', density=False, alpha=0.8, color=colors[0], ls='solid' , histtype='step', label='Measured') + axs[i].hist(mean_clock_phase, bins='sqrt', density=False, alpha=0.8, color=colors[0], ls='solid' , histtype='step', label='Measured') axs[i].hist(actual_antenna_phase_shifts, bins='sqrt', density=False, alpha=0.8, color=colors[1], ls='dashed', histtype='step', label='Actual') @@ -209,7 +209,7 @@ if __name__ == "__main__": if plot_residuals: axs[i].plot(phase_residuals, np.arange(N_ant), alpha=0.6, ls='none', marker='x', color=colors[0]) else: - axs[i].errorbar(mean_sigma_phase, np.arange(N_ant), yerr=std_sigma_phase, marker='4', alpha=0.7, ls='none', color=colors[0], label='Measured') + axs[i].errorbar(mean_clock_phase, np.arange(N_ant), yerr=std_clock_phase, marker='4', alpha=0.7, ls='none', color=colors[0], label='Measured') axs[i].plot(actual_antenna_phase_shifts, antenna_names, ls='none', marker='3', alpha=0.8, color=colors[1], label='Actual') axs[i].legend() @@ -231,10 +231,10 @@ if __name__ == "__main__": actual_baseline_time_shifts.append(actual_baseline_time_shift) - # unpack mean_sigma_phase back into a list of time diffs + # unpack mean_clock_phase back into a list of time diffs measured_baseline_time_diffs = [] for i,b in enumerate(basenames): - phase0, phase1 = mean_sigma_phase[name2idx(b[0])], mean_sigma_phase[name2idx(b[1])] + phase0, phase1 = mean_clock_phase[name2idx(b[0])], mean_clock_phase[name2idx(b[1])] measured_baseline_time_diffs.append(lib.phase_mod(phase1 - phase0)/(2*np.pi*f_beacon)) # Make a plot diff --git a/simulations/airshower_beacon_simulation/cb_report_measured_antenna_time_offsets.py b/simulations/airshower_beacon_simulation/cb_report_measured_antenna_time_offsets.py index 367b521..17edaff 100755 --- a/simulations/airshower_beacon_simulation/cb_report_measured_antenna_time_offsets.py +++ b/simulations/airshower_beacon_simulation/cb_report_measured_antenna_time_offsets.py @@ -51,7 +51,7 @@ if __name__ == "__main__": # TODO: redo matrix sweeping for new timing?? measured_antenna_time_shifts = {} for i, ant in enumerate(antennas): - clock_phase_time = ant.beacon_info[freq_name]['sigma_phase_mean']/(2*np.pi*f_beacon) + clock_phase_time = ant.beacon_info[freq_name]['clock_phase_mean']/(2*np.pi*f_beacon) best_k_time = ant.beacon_info[freq_name]['best_k_time'] total_clock_time = best_k_time + clock_phase_time