From 8da9d55c568cee9b23df94385d3f5289acf2a30c Mon Sep 17 00:00:00 2001 From: Eric Teunis de Boone Date: Thu, 24 Nov 2022 14:47:29 +0100 Subject: [PATCH] ZH: total time_diffs saved to file --- .../aa_generate_beacon.py | 20 ++++--- .../bc_beacon_periods.py | 58 +++++++++++-------- .../airshower_beacon_simulation/lib/lib.py | 16 ++++- 3 files changed, 60 insertions(+), 34 deletions(-) diff --git a/simulations/airshower_beacon_simulation/aa_generate_beacon.py b/simulations/airshower_beacon_simulation/aa_generate_beacon.py index 05b1669..7d1b7cb 100755 --- a/simulations/airshower_beacon_simulation/aa_generate_beacon.py +++ b/simulations/airshower_beacon_simulation/aa_generate_beacon.py @@ -167,14 +167,15 @@ def read_baseline_time_diffs_hdf5(fname): names = group[base_dset_name][:] dset = group[dset_name] - f_beacon = dset[:,0] - true_phase_diffs = dset[:,1] - k_periods = dset[:,2] + time_diffs = dset[:,0] + f_beacon = dset[:,1] + true_phase_diffs = dset[:,2] + k_periods = dset[:,3] - return names, f_beacon, true_phase_diffs, k_periods + return names, time_diffs, f_beacon, true_phase_diffs, k_periods -def write_baseline_time_diffs_hdf5(fname, baselines, true_phase_diffs, k_periods, f_beacon, overwrite=True): +def write_baseline_time_diffs_hdf5(fname, baselines, true_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. @@ -187,14 +188,17 @@ def write_baseline_time_diffs_hdf5(fname, baselines, true_phase_diffs, k_periods baselines = [baselines] true_phase_diffs = [true_phase_diffs] k_periods = [k_periods] - f_beacon = [f_beacon] + f_beacon = np.array([f_beacon]) else: N_baselines = len(baselines) # Expand the f_beacon list if not hasattr(f_beacon, '__len__'): - f_beacon = [f_beacon]*N_baselines + 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) assert len(baselines) == len(true_phase_diffs) == len(k_periods) == len(f_beacon) @@ -220,7 +224,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( (f_beacon, true_phase_diffs, k_periods) ).T + data = np.vstack( (time_diffs, f_beacon, true_phase_diffs, k_periods) ).T dset = group.create_dataset(dset_name, data=data) diff --git a/simulations/airshower_beacon_simulation/bc_beacon_periods.py b/simulations/airshower_beacon_simulation/bc_beacon_periods.py index c1db54a..f2ff709 100755 --- a/simulations/airshower_beacon_simulation/bc_beacon_periods.py +++ b/simulations/airshower_beacon_simulation/bc_beacon_periods.py @@ -31,7 +31,8 @@ if __name__ == "__main__": baselines = list(combinations(antennas,2)) # use ref_ant else: - ref_ant = antennas[ref_ant_idx] + ref_ant = antennas[ref_ant_id] + print(f"Doing all baselines with {ref_ant.name}") baselines = list(zip_longest([], antennas, fillvalue=ref_ant)) freq_names = antennas[0].beacon_info.keys() @@ -47,38 +48,41 @@ if __name__ == "__main__": # which traces to keep track of traces = [ base[0].E_AxB, base[1].E_AxB ] - sampling_dt = (base[1].t[1] - base[1].t[0]) # ns - # how many samples do we need to shift - ks, maxima = lib.coherence_sum_maxima(traces[0], traces[1]) - max_idx = np.argmax(maxima) - best_k = ks[max_idx] - delta_t_coherence = sampling_dt*best_k # ns + # read f_beacon from the first antenna + f_beacon = base[0].beacon_info[freq_name]['freq'] - print('A1:', base[0].name, 'A2:', base[1].name, "K:", best_k, '= [ns]', delta_t_coherence) + # how many samples do we need to shift + sample_shifts, maxima = lib.coherence_sum_maxima(traces[0], traces[1], periodic=False) + best_sample_shift = sample_shifts[np.argmax(maxima)] + + # turn sample_shift into time + sampling_dt = (base[1].t[1] - base[1].t[0]) # ns + delta_t_coherence = sampling_dt*best_sample_shift # ns # get the amount of periods to move - f_beacon = base[0].beacon_info[freq_name]['freq'] - k_period, rest = np.divmod(delta_t_coherence, 1/f_beacon) + k_period, t_rest = np.divmod(delta_t_coherence, 1/f_beacon) # always keep the reference before traces[1] - if rest < 0: + if t_rest < 0: # np.divmod already does this k_period -= 1 + t_rest = 1/f_beacon + t_rest # 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(true_phases[0] - true_phases[1]) except IndexError: - # freq_name not in beacon_info - # or true_phase not determined yet + # 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 # save k_period with antenna names time_diffs[i] = [true_phases_diff, k_period, f_beacon] # Plotting for one or two iterations - if show_plots and i in [ 0, 1 ]: - print('i',i,'k[T]',k_period, 'rest[ns]',rest, 'T[ns]',1/f_beacon) + if show_plots and (i in [ 0, 1 ] or k_period > 3): + # More than three periods is quite much so report it + print('i',i,'k[T]',k_period, 'rest[ns]',t_rest, 'T[ns]',1/f_beacon, 'dT_coher[ns]', delta_t_coherence) # Show correlation maxima plot if not True: @@ -97,21 +101,29 @@ if __name__ == "__main__": true_phases_diff = 0 delta_t_beacon = true_phases_diff/(2*np.pi*f_beacon) - print("t0[ns]", delta_t_antennas, "t_beacon[ns]", delta_t_beacon, "phase", true_phases_diff) fig, ax = plt.subplots() - ax.set_xlabel('t') - ax.plot(base[0].t, traces[0], label=f'Reference {base[0].name}', alpha=0.5) + ax.set_title( + ", ".join([ + f"$\\Delta$t0 [ns] : {delta_t_antennas:.2f}", + f"$\\Delta$t_beacon [ns]: {delta_t_beacon:.2f}", + f"$\\Delta\\sigma_\\varphi$: {true_phases_diff:.4f}", + f"", + ]) + ) + ax.set_xlabel('Sampling t [ns]') + ax.set_ylabel('Amplitude [a.u.]') + ax.plot(base[0].t, traces[0], label=f'Reference: {base[0].name}', alpha=0.5) # plot vertical lines indicating f_beacon min_t, max_t = base[0].t[0], base[0].t[-1] N_lines = int( (max_t - min_t)*f_beacon) +1 for i, t in enumerate(np.arange(N_lines)/f_beacon): - ax.axvline( min_t + t, color='k', alpha=0.5, label=None if i!=0 else 'P_beacon') + ax.axvline( min_t + t, color='k', alpha=0.3) - ax.plot(base[1].t + delta_t_antennas, traces[1], label=f'Original {base[1].name} (t0 removed)', alpha=0.4, marker='+', ms=5) - ax.plot(base[1].t + delta_t_antennas + k_period/f_beacon + rest, traces[1], label='Coherence', alpha=0.3, marker='x', ms=5) - ax.plot(base[1].t + delta_t_antennas + k_period/f_beacon + delta_t_beacon, traces[1], label='Beacon only + Periods', alpha=0.6) + ax.plot(base[1].t + delta_t_antennas, traces[1], label=f'Original: {base[1].name} (t0 removed)', alpha=0.4, marker='+', ms=5) + ax.plot(base[1].t + delta_t_antennas + k_period/f_beacon + t_rest, traces[1], label='Coherence', alpha=0.3, marker='x', ms=5) + ax.plot(base[1].t + delta_t_antennas + k_period/f_beacon + delta_t_beacon, traces[1], label=f'$\\Delta t_\\varphi$ + $k={k_period:.0f}$ Periods', alpha=0.6) - ax.legend() + ax.legend(fancybox=True, framealpha=0.5) # Save integer periods to antennas beacon.write_baseline_time_diffs_hdf5(antennas_fname, baselines, time_diffs[:,0], time_diffs[:,1], time_diffs[:,2]) diff --git a/simulations/airshower_beacon_simulation/lib/lib.py b/simulations/airshower_beacon_simulation/lib/lib.py index 8d07687..7441686 100644 --- a/simulations/airshower_beacon_simulation/lib/lib.py +++ b/simulations/airshower_beacon_simulation/lib/lib.py @@ -200,16 +200,26 @@ def find_beacon_in_traces( return frequencies, phases, amplitudes -def coherence_sum_maxima(ref_x, y, k_step=1): +def coherence_sum_maxima(ref_x, y, k_step=1, k_start=0, k_end=None, periodic=True): """ Use the maximum of a coherent sum to determine the best number of samples to move """ - max_k = int( len(ref_x) ) - ks = np.arange(0, max_k, step=k_step) + N_samples = int( len(ref_x) ) + k_end = N_samples if k_end is None or k_end > max_k else k_end + + ks = np.arange(k_start, k_end, step=k_step) maxima = np.empty(len(ks)) + if periodic is False: + # prepend zeros + N_zeros = N_samples + preshift = 0 # only required for testing purposes + + ref_x = np.pad(ref_x, (N_zeros-0,0), 'constant') + y = np.pad(y, (N_zeros-preshift,preshift), 'constant') + for i,k in enumerate(ks, 0): augmented_y = np.roll(y, k) maxima[i] = max(ref_x + augmented_y)