ZH: total time_diffs saved to file

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]
+        time_diffs = dset[:,0]
-        true_phase_diffs = dset[:,1]
+        f_beacon = dset[:,1]
-        k_periods = dset[:,2]
+        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)
 
 

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
+        # read f_beacon from the first antenna
-        # how many samples do we need to shift
+        f_beacon = base[0].beacon_info[freq_name]['freq']
-        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
 
-        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, t_rest = np.divmod(delta_t_coherence, 1/f_beacon)
-        k_period, 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
+            # true_phase not determined yet
-            # or 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 ]:
+        if show_plots and (i in [ 0, 1 ] or k_period > 3):
-            print('i',i,'k[T]',k_period, 'rest[ns]',rest, 'T[ns]',1/f_beacon)
+            # 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.set_title(
-            ax.plot(base[0].t, traces[0], label=f'Reference {base[0].name}', alpha=0.5)
+                ", ".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, 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 + 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='Beacon only + Periods', alpha=0.6)
+            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])

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) )
+    N_samples = int( len(ref_x) )
-    ks = np.arange(0, max_k, step=k_step)
+    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)