diff --git a/simulations/airshower_beacon_simulation/bc_baseline_phase_deltas.py b/simulations/airshower_beacon_simulation/bc_baseline_phase_deltas.py
index 6bf9032..ce1783b 100755
--- a/simulations/airshower_beacon_simulation/bc_baseline_phase_deltas.py
+++ b/simulations/airshower_beacon_simulation/bc_baseline_phase_deltas.py
@@ -2,7 +2,7 @@
 # vim: fdm=indent ts=4
 
 import h5py
-from itertools import combinations, zip_longest
+from itertools import combinations, product
 import matplotlib.pyplot as plt
 import numpy as np
 
@@ -22,7 +22,7 @@ if __name__ == "__main__":
     fname = "ZH_airshower/mysim.sry"
     c_light = 3e8*1e-9
     show_plots = True
-    ref_ant_id = 50 # leave None for all baselines
+    ref_ant_id = None if not True else [50] # leave None for all baselines
 
     ####
     fname_dir = path.dirname(fname)
@@ -44,9 +44,9 @@ if __name__ == "__main__":
         baselines = list(combinations(antennas,2))
     # use ref_ant
     else:
-        ref_ant = antennas[ref_ant_id]
-        print(f"Doing all baselines with {ref_ant.name}")
-        baselines = list(zip_longest([], antennas, fillvalue=ref_ant))
+        ref_ants = [antennas[i] for i in ref_ant_id]
+        print("Doing all baselines with {}".format([int(a.name) for a in ref_ants]))
+        baselines = list(product(ref_ants, antennas))
 
     # For now, only one beacon_frequency is supported
     freq_names = antennas[0].beacon_info.keys()
@@ -111,11 +111,12 @@ if __name__ == "__main__":
             if plot_residuals:
                 phase_residuals = lib.phase_mod(phase_diffs[:,1] - my_phase_diffs)
 
-                axs[0].set_title("Difference between Measured and Actual phase difference\n for Baseline (i,j" + ( '='+str(ref_ant_id) if ref_ant_id is not None else None) + ')')
-                axs[1].set_xlabel("Baseline Phase Residual $\\varphi_{ij_{meas}} - \\varphi_{ij_{true}}$ [rad]")
+                fig.suptitle("Difference between Measured and Actual phase difference\n for Baselines (i,j" + (')' if ref_ant_id is None else '='+str([ int(a.name) for a in ref_ants])+')'))
+                axs[-1].set_xlabel("Baseline Phase Residual $\\varphi_{ij_{meas}} - \\varphi_{ij_{true}}$ [rad]")
             else:
-                axs[0].set_title("Comparison Measured and Actual phase difference\n for Baseline (i,j" + ( '='+str(ref_ant_id) if ref_ant_id is not None else None) + ')')
-                axs[1].set_xlabel("Baseline Phase $\\varphi_{ij}$ [rad]")
+                fig.suptitle("Comparison Measured and Actual phase difference\n for Baselines (i,j"  + (')' if ref_ant_id is None else '='+str([ int(a.name) for a in ref_ants])+')'))
+                axs[-1].set_xlabel("Baseline Phase $\\varphi_{ij}$ [rad]")
+
 
             i=0
             axs[i].set_ylabel("#")
@@ -125,15 +126,16 @@ if __name__ == "__main__":
                 axs[i].hist(phase_diffs[:,1], bins='sqrt', density=False, alpha=0.8, color=colors[0], ls='solid' , histtype='step', label='Measured')
                 axs[i].hist(my_phase_diffs,   bins='sqrt', density=False, alpha=0.8, color=colors[1], ls='dashed', histtype='step', label='Actual')
 
+
             i=1
             axs[i].set_ylabel("Baseline no.")
-            if not plot_residuals:
+            if plot_residuals:
+                axs[i].plot(phase_residuals, np.arange(N_base), alpha=0.6, ls='none', marker='x', color=colors[0])
+            else:
                 axs[i].plot(phase_diffs[:,1], np.arange(N_base), alpha=0.8, color=colors[0], ls='none', marker='x', label='calculated')
                 axs[i].plot(my_phase_diffs,   np.arange(N_base), alpha=0.8, color=colors[1], ls='none', marker='+', label='actual time shifts')
 
                 axs[i].legend()
-            else:
-                axs[i].plot(phase_residuals, np.arange(N_base), alpha=0.6, ls='none', marker='x', color=colors[0])
             fig.tight_layout()
 
             if fig_dir: