diff --git a/magicctapipe/io/io.py b/magicctapipe/io/io.py
index 133ef46c..4912bea7 100644
--- a/magicctapipe/io/io.py
+++ b/magicctapipe/io/io.py
@@ -1541,17 +1541,34 @@ def make_time_select(
         comb_array = np.array(np.meshgrid(t_lst_all, t_magic_all)).reshape(2, -1)
         time_offsets = comb_array[0] - comb_array[1]
         time_window = 500e-9
-        b = np.array(
-            [np.sum(np.abs(time_offsets - T) < time_window) for T in time_offsets]
-        )
-        n_coincident = np.max(b)
-        time_offset_best = np.mean(time_offsets[b == np.max(b)])
+
+        # original method should have a complexity of len (time_offsets)^2
+        # b = np.array(
+        #    [np.sum(np.abs(time_offsets - T) < time_window) for T in time_offsets]
+        # )
+        # n_coincident = np.max(b)
+        # time_offset_best = np.mean(time_offsets[b == np.max(b)])
+
+        # method 2
+        # should have complexity of
+        # len(time_offsets) * ln (len(time_offsets)) # sorting
+        # len(time_offsets) * n_coincident # searching for highest n_coincident
+        offsets_sorted = np.sort(time_offsets)
+        diffs = np.diff(offsets_sorted)
+        n_coincident = 1
+        while np.min(diffs) <= 2 * time_window:
+            n_coincident += 1
+            diffs = offsets_sorted[n_coincident:] - offsets_sorted[:-n_coincident]
+            # e.g. n_coin=3, time diffs 1-3, 2-4, 3-5, ...
+        pos = np.argmin(diffs)
+        # median seem to work a bit better than the mean
+        # (and more consistent with the original method
+        time_offset_best = np.median(offsets_sorted[pos : pos + n_coincident])
     else:
         t_magic_ave = np.mean(data_magic_["trigger_time"].values)
         event_id_magic_ave = np.mean(data_magic_["event_id"].values[N_start:N_end])
         n_coincident = 0
         time_offset_best = 0
-
     return t_magic_ave, event_id_magic_ave, time_offset_best, n_coincident
 
 
diff --git a/magicctapipe/scripts/lst1_magic/coincident_brokenGPS.py b/magicctapipe/scripts/lst1_magic/coincident_brokenGPS.py
index 00e9c8b3..d3725dc9 100644
--- a/magicctapipe/scripts/lst1_magic/coincident_brokenGPS.py
+++ b/magicctapipe/scripts/lst1_magic/coincident_brokenGPS.py
@@ -369,8 +369,7 @@ def print_gti(df1, df2):
         # PLEASE REWRITE THIS according to your analysis environment (e.g., SLURM)
         output = subprocess.run(
             [
-                "python",
-                "lst1_magic_event_coincidence.py",
+                "lst1_magic_event_coincidence",
                 "-l",
                 f"{lst_dir_name}/dl1_LST-1.Run{lst_RunID}.{lst_run}.h5",
                 "-m",
diff --git a/magicctapipe/scripts/lst1_magic/lst1_magic_event_coincidence.py b/magicctapipe/scripts/lst1_magic/lst1_magic_event_coincidence.py
index c3102a2c..abd1d50d 100644
--- a/magicctapipe/scripts/lst1_magic/lst1_magic_event_coincidence.py
+++ b/magicctapipe/scripts/lst1_magic/lst1_magic_event_coincidence.py
@@ -272,20 +272,16 @@ def event_coincidence(
 
         if input_dir_toff is not None:
             files_toff = glob.glob(f"{input_dir_toff}/*M{str(tel_id - 1)}*_detail.npy")
+            files_npy = []
             for i, file_toff in enumerate(files_toff):  # files_toff:
                 file_npy = np.load(file_toff)
-                if i == 0:
-                    file_npy_ = file_npy
-                else:
-                    file_npy_ = np.hstack([file_npy, file_npy_])
-
+                files_npy.append(file_npy)
+            file_npy_ = np.hstack(files_npy)
             df_toff = pd.DataFrame(file_npy_.T, columns=["timestamp", "toff1", "n"])
-
             timestamps_lst_min, timestamps_lst_max = (  # noqa: F841
                 float(min(timestamps_lst_org.value)) / 1e9,
                 float(max(timestamps_lst_org.value)) / 1e9,
             )
-
             df_toff = df_toff.query(
                 f"{timestamps_lst_min} < timestamp < {timestamps_lst_max}"
             )