diff --git a/simulations/airshower_beacon_simulation/ca_period_from_shower.py b/simulations/airshower_beacon_simulation/ca_period_from_shower.py
index 2b361bb..bea9d7e 100755
--- a/simulations/airshower_beacon_simulation/ca_period_from_shower.py
+++ b/simulations/airshower_beacon_simulation/ca_period_from_shower.py
@@ -29,10 +29,17 @@ def find_best_sample_shifts_summing_at_location(test_loc, antennas, allowed_samp
     t_min = 1e9
     t_max = -1e9
     a_maxima = []
+    N_ant = len(antennas)
 
     if dt is None:
         dt = antennas[0].t_AxB[1] - antennas[0].t_AxB[0]
 
+    if not hasattr(plot_iteration_with_shifted_trace, '__len__'):
+        if plot_iteration_with_shifted_trace:
+            plot_iteration_with_shifted_trace = [ plot_iteration_with_shifted_trace ]
+        else:
+            plot_iteration_with_shifted_trace = []
+
     # propagate to test location
     for i, ant in enumerate(antennas):
         aloc = [ant.x, ant.y, ant.z]
@@ -72,9 +79,9 @@ def find_best_sample_shifts_summing_at_location(test_loc, antennas, allowed_samp
             continue
 
         # init figure
-        if i == plot_iteration_with_shifted_trace:
+        if i in plot_iteration_with_shifted_trace:
             fig, ax = plt.subplots()
-            ax.set_title("Traces at ({:.1f},{:.1f},{:.1f})".format(*test_loc))
+            ax.set_title("Traces at ({:.1f},{:.1f},{:.1f}) i={i}/{tot}".format(*test_loc, i=i, tot=N_ant))
             ax.set_xlabel("Time [ns]")
             ax.set_ylabel("Amplitude")
             ax.plot(t_sum, a_sum)
@@ -85,18 +92,38 @@ def find_best_sample_shifts_summing_at_location(test_loc, antennas, allowed_samp
 
             shift_maxima[j] = np.max(augmented_a + a_sum)
 
-            if i == plot_iteration_with_shifted_trace:
-                ax.plot(t_sum, augmented_a, label=f'{shift} shifted')
+            if i in plot_iteration_with_shifted_trace:
+                ax.plot(t_sum, augmented_a, alpha=0.8, label=f'{shift}')
 
         # transform maximum into best_sample_shift
         best_idx = np.argmax(shift_maxima)
+
         best_sample_shifts[i] = allowed_sample_shifts[best_idx]
         a_sum += np.roll(a_int, best_sample_shifts[i])
 
         # cleanup figure
-        if i == plot_iteration_with_shifted_trace:
-            if fig_dir: # note this is a global variable here
-                fig.savefig(path.join(fig_dir, __file__ + '.loc{:.1f}-{:.1f}-{:.1f}'.format(*test_loc) + f'.i{i}{fig_distinguish}.pdf'))
+        if i in plot_iteration_with_shifted_trace:
+            ax.legend( ncol=5)
+            if fig_dir:
+                fname = path.join(fig_dir, __file__ + f'.{fig_distinguish}i{i}' + '.loc{:.1f}-{:.1f}-{:.1f}'.format(*test_loc))
+                if True:
+                    old_xlim = ax.get_xlim()
+
+                    if True: # zoomed on part without peak of this trace
+                        wx = 100
+                        x = max(t_r) - wx
+                        ax.set_xlim(x-wx, x+wx)
+                        fig.savefig(fname + ".zoomed.beacon.pdf")
+
+                    if True: # zoomed on peak of this trace
+                        x = t_r[np.argmax(E_)]
+                        wx = 10 + max(best_sample_shifts) - min(best_sample_shifts)
+                        ax.set_xlim(x-wx, x+wx)
+                        fig.savefig(fname + ".zoomed.peak.pdf")
+
+                    ax.set_xlim(*old_xlim)
+
+                fig.savefig(fname + ".pdf")
             plt.close(fig)
 
     # sort by antenna (undo sorting by maximum)
@@ -259,7 +286,7 @@ if __name__ == "__main__":
             yy.append(y_+yoff)
 
             # Find best k for each antenna
-            shifts, maximum = find_best_sample_shifts_summing_at_location(test_loc, ev.antennas, allowed_sample_shifts, dt=dt, fig_dir=tmp_fig_subdir, plot_iteration_with_shifted_trace=len(ev.antennas)-1, fig_distinguish=f".run{r}")
+            shifts, maximum = find_best_sample_shifts_summing_at_location(test_loc, ev.antennas, allowed_sample_shifts, dt=dt, fig_dir=tmp_fig_subdir, plot_iteration_with_shifted_trace=[ 5, len(ev.antennas)-1], fig_distinguish=f"run{r}.")
 
             # Translate sample shifts back into period multiple k
             ks = np.rint(shifts*f_beacon*dt)