diff --git a/pylot/correlation/parameters_adriaarray.yaml b/pylot/correlation/parameters_adriaarray.yaml
index 74d36b20..ed3b36f0 100644
--- a/pylot/correlation/parameters_adriaarray.yaml
+++ b/pylot/correlation/parameters_adriaarray.yaml
@@ -59,6 +59,9 @@ P:
   filter_type: bandpass
   sampfreq: 20.0
 
+  # ignore if autopylot fails to pick master-trace (not recommended if absolute onset times matter)
+  ignore_autopylot_fail_on_master: True
+
 # S-phase
 S:
   min_corr_stacking: 0.7
@@ -93,3 +96,6 @@ S:
   filter_type: bandpass
   sampfreq: 20.0
 
+  # ignore if autopylot fails to pick master-trace (not recommended if absolute onset times matter)
+  ignore_autopylot_fail_on_master: True
+
diff --git a/pylot/correlation/pick_correlation_correction.py b/pylot/correlation/pick_correlation_correction.py
index f23c3e03..70ceb988 100644
--- a/pylot/correlation/pick_correlation_correction.py
+++ b/pylot/correlation/pick_correlation_correction.py
@@ -436,14 +436,15 @@ def correlation_main(database_path_dmt: str, pylot_infile_path: str, params: dic
 
     # iterate over all events in "database_path_dmt"
     for eventindex, eventdir in enumerate(eventdirs):
-        if not istart <= eventindex + 1 < istop:
+        eventindex += 1
+        if not istart <= eventindex < istop:
             continue
 
         if select_events and not os.path.split(eventdir)[-1] in select_events:
             continue
 
         logging.info('\n' + 100 * '#')
-        logging.info('Working on event {} ({}/{})'.format(eventdir, eventindex + 1, len(eventdirs)))
+        logging.info('Working on event {} ({}/{})'.format(eventdir, eventindex, len(eventdirs)))
         if event_blacklist and get_event_id(eventdir) in event_blacklist:
             logging.info('Event on blacklist. Continue')
 
@@ -526,7 +527,7 @@ def cut_stream_to_same_length(wfdata: Stream) -> None:
         st.trim(tstart, tend, pad=True, fill_value=0.)
         # check for stream length
         if len(st) < 3:
-            logging.info('Not enough traces in stream, remove it from dataset:', str(st))
+            logging.info(f'Not enough traces in stream, remove it from dataset: {st}')
             remove(wfdata, st)
             continue
         if st[0].stats.starttime != st[1].stats.starttime or st[1].stats.starttime != st[2].stats.starttime:
@@ -684,19 +685,56 @@ def correlate_event(eventdir: str, pylot_parameter: PylotParameter, params: dict
         if not params[phase_type]['use_stacked_trace']:
             # first stack mastertrace by using stations with high correlation on that station in station list with the
             # highest correlation coefficient
-            stack_result = stack_mastertrace(wfdata_lowf, wfdata_highf, wfdata, picks, params=params[phase_type],
-                                             channels=channels_list, method=method, fig_dir=fig_dir)
-        else:
-            stack_result = load_stacked_trace(eventdir, params[phase_type]['min_corr_stacking'])
+            logging.info('Searching for master trace. ' + 20 * '-')
+            best_result = None
+            for stack_result in stack_mastertrace(wfdata_lowf, wfdata_highf, wfdata, picks, params=params[phase_type],
+                                                  channels=channels_list, method=method, fig_dir=fig_dir):
+
+                if not stack_result:
+                    continue
+
+                # save best result in case it is still needed
+                if not best_result:
+                    best_result = stack_result
+
+                # extract stack result
+                correlations_dict, nwst_id, trace_master, nstack = stack_result
+
+                # now pick stacked trace with PyLoT for a more precise pick (use raw trace, gets filtered by autoPyLoT)
+                pick_stacked = repick_master_trace(wfdata_lowf, trace_master, pylot_parameter, event, event_id, metadata,
+                                                   phase_type, correlation_out_dir)
+
+                if not pick_stacked or not pick_stacked.time_errors.uncertainty:
+                    logging.info(f'Invalid autoPyLoT pick on master trace. Try next one.')
+                    continue
+                else:
+                    break
+
+            else:
+                logging.info('Did not find autoPyLoT pick for any stacked trace.')
+
+                if params[phase_type]['ignore_autopylot_fail_on_master'] is True:
+                    # Fallback in case autopylot failed
+                    logging.warning(f'Could not pick stacked trace. Using reference pick instead. If you need '
+                                    f'absolute onsets you need to account for a time-shift.')
+
+                    correlations_dict, nwst_id, trace_master, nstack = best_result
+                    pick_stacked = get_pick4station(picks, network_code=trace_master.stats.network,
+                                                    station_code=trace_master.stats.station, method='auto')
+                else:
+                    pick_stacked = None
+
+                if not pick_stacked:
+                    logging.info('Could not find reference pick also. Continue with next phase.')
+                    continue
+
+        else:
+            raise NotImplementedError('Loading stacked trace currently not implemented')
+            # stack_result = load_stacked_trace(eventdir, params[phase_type]['min_corr_stacking'])
 
-        if not stack_result:
-            logging.info('No stack result. Continue.')
-            continue
 
         #############################################
         # NOW DO THE FINAL CORRELATION
-        # extract stack result
-        correlations_dict, nwst_id, trace_master, nstack = stack_result
 
         if params[phase_type]['plot']:
             # plot correlations of traces used to generate stacked trace
@@ -713,11 +751,6 @@ def correlate_event(eventdir: str, pylot_parameter: PylotParameter, params: dict
         # write unfiltered trace
         trace_master.write(os.path.join(correlation_out_dir, '{}_stacked.mseed'.format(trace_master.id)))
 
-        # now pick stacked trace with PyLoT for a more precise pick (use raw trace, gets filtered by autoPyLoT)
-        pick_stacked = repick_master_trace(wfdata_lowf, trace_master, pylot_parameter, event, event_id, metadata,
-                                           phase_type, correlation_out_dir)
-        if not pick_stacked:
-            continue
 
         # correlate stations with repicked and stacked master trace
         fig_dir_traces = make_figure_dirs(fig_dir, trace_master.id)
@@ -1178,38 +1211,38 @@ def stack_mastertrace(wfdata_lowf: Stream, wfdata_highf: Stream, wfdata_raw: Str
     A master trace will be created by stacking well correlating traces onto this station.
     """
 
-    def get_best_station4stack(sta_result):
-        """ return station with maximum mean_ccc"""
+    def get_best_stations4stack(sta_result, n_max=4):
+        """ return stations sorted after maximum mean_ccc"""
         ccc_means = {nwst_id: value['mean_ccc'] for nwst_id, value in sta_result.items() if
                      not np.isnan(value['mean_ccc'])}
         if len(ccc_means) < 1:
             logging.warning('No valid station found for stacking! Return.')
             return
-        best_station_id = max(ccc_means, key=ccc_means.get)
-        logging.info(
-            'Found highest mean correlation for station {} ({})'.format(best_station_id, max(ccc_means.values())))
-        return best_station_id
+        best_station_ids = sorted(ccc_means, key=ccc_means.get, reverse=True)[:n_max]
+        logging.info(f'Found mean correlations for stations: {ccc_means}')
+        return best_station_ids, ccc_means
 
     station_results = iterate_correlation(wfdata_lowf, wfdata_highf, channels, picks, method, params, fig_dir=fig_dir)
-    nwst_id_master = get_best_station4stack(station_results)
+    nwst_ids_master, ccc_means = get_best_stations4stack(station_results)
 
     # in case no stream with a valid pick is found
-    if not nwst_id_master:
-        logging.info('No mastertrace found! Will skip this event.')
+    if not nwst_ids_master:
+        logging.info('No mastertrace found! Continue with next.')
         return None
 
-    trace_master = station_results[nwst_id_master]['trace']
-    stations4stack = station_results[nwst_id_master]['stations4stack']
-    correlations_dict = station_results[nwst_id_master]['correlations_dict']
+    for nwst_id_master in nwst_ids_master:
+        trace_master = station_results[nwst_id_master]['trace']
+        stations4stack = station_results[nwst_id_master]['stations4stack']
+        correlations_dict = station_results[nwst_id_master]['correlations_dict']
 
-    wfdata_highf += trace_master
+        wfdata_highf += trace_master
 
-    dt_pre, dt_post = params['dt_stacking']
-    trace_master, nstack = apply_stacking(trace_master, stations4stack, wfdata_raw, picks, method=method,
-                                          do_rms_check=params['check_RMS'], plot=params['plot'], fig_dir=fig_dir,
-                                          dt_pre=dt_pre, dt_post=dt_post)
+        dt_pre, dt_post = params['dt_stacking']
+        trace_master, nstack = apply_stacking(trace_master, stations4stack, wfdata_raw, picks, method=method,
+                                              do_rms_check=params['check_RMS'], plot=params['plot'], fig_dir=fig_dir,
+                                              dt_pre=dt_pre, dt_post=dt_post)
 
-    return correlations_dict, nwst_id_master, trace_master, nstack
+        yield correlations_dict, nwst_id_master, trace_master, nstack
 
 
 def iterate_correlation(wfdata_lowf: Stream, wfdata_highf: Stream, channels: list, picks: list, method: str,
@@ -1998,4 +2031,4 @@ if __name__ == "__main__":
     correlation_main(ARGS.dmt_path, ARGS.pylot_infile, params=CORR_PARAMS, istart=int(ARGS.istart),
                      istop=int(ARGS.istop),
                      channel_config=CHANNELS, update=ARGS.update, event_blacklist=ARGS.blacklist,
-                     select_events=['20250107_010516.a'])
+                     select_events=None) # can add a list of event ids if only specific events shall be picked
diff --git a/pylot/correlation/submit_pick_corr_correction.sh b/pylot/correlation/submit_pick_corr_correction.sh
index 1531f801..40bf3869 100755
--- a/pylot/correlation/submit_pick_corr_correction.sh
+++ b/pylot/correlation/submit_pick_corr_correction.sh
@@ -32,9 +32,10 @@ export PYTHONPATH="$PYTHONPATH:/home/marcel/git/pylot/"
 #python pick_correlation_correction.py '/data/AlpArray_Data/dmt_database_mantle_M5.8-6.0' '/home/marcel/.pylot/pylot_alparray_mantle_corr_stack_0.03-0.5.in' -pd -n ${NSLOTS:=1} -istart 100 -istop 200
 #python pick_correlation_correction.py 'H:\sciebo\dmt_database' 'H:\Sciebo\dmt_database\pylot_alparray_mantle_corr_S_0.01-0.2.in' -pd -n 4 -t
 
-pylot_infile='/home/marcel/.pylot/pylot_alparray_syn_fwi_mk6_it3.in'
-#pylot_infile='/home/marcel/.pylot/pylot_adriaarray_corr_P_and_S.in'
+#pylot_infile='/home/marcel/.pylot/pylot_alparray_syn_fwi_mk6_it3.in'
+pylot_infile='/home/marcel/.pylot/pylot_adriaarray_corr_P_and_S.in'
 
 # THIS SCRIPT SHOLD BE CALLED BY "submit_to_grid_engine.py" using the following line:
-python pick_correlation_correction.py $1 $pylot_infile -pd -n ${NSLOTS:=1} -istart $2 --params 'parameters_fwi_mk6_it3.yaml'
+# use -pd for detailed plots in eventdir/correlation_XX_XX/figures
+python pick_correlation_correction.py $1 $pylot_infile -n ${NSLOTS:=1} -istart $2 --params 'parameters_adriaarray.yaml' # -pd
 #--event_blacklist eventlist.txt
diff --git a/pylot/correlation/submit_to_grid_engine.py b/pylot/correlation/submit_to_grid_engine.py
index a5410d07..fb5ea728 100755
--- a/pylot/correlation/submit_to_grid_engine.py
+++ b/pylot/correlation/submit_to_grid_engine.py
@@ -3,21 +3,26 @@
 import subprocess
 
 fnames = [
-    ('/data/AlpArray_Data/dmt_database_synth_model_mk6_it3_no_rotation', 0),
+    ('/data/AdriaArray_Data/dmt_database_mantle_M5.0-5.4', 0),
+    ('/data/AdriaArray_Data/dmt_database_mantle_M5.4-5.7', 0),
+    ('/data/AdriaArray_Data/dmt_database_mantle_M5.7-6.0', 0),
+    ('/data/AdriaArray_Data/dmt_database_mantle_M6.0-6.3', 0),
+    ('/data/AdriaArray_Data/dmt_database_mantle_M6.3-10.0', 0),
+    # ('/data/AdriaArray_Data/dmt_database_ISC_mantle_M5.0-5.4', 0),
+    # ('/data/AdriaArray_Data/dmt_database_ISC_mantle_M5.4-5.7', 0),
+    # ('/data/AdriaArray_Data/dmt_database_ISC_mantle_M5.7-6.0', 0),
+    # ('/data/AdriaArray_Data/dmt_database_ISC_mantle_M6.0-10.0', 0),
           ]
 
 #fnames = [('/data/AlpArray_Data/dmt_database_mantle_0.01-0.2_SKS-phase', 0),
 #          ('/data/AlpArray_Data/dmt_database_mantle_0.01-0.2_S-phase', 0),]
 
 ####
-script_location = '/home/marcel/VersionCtrl/git/code_base/correlation_picker/submit_pick_corr_correction.sh'
+script_location = '/home/marcel/VersionCtrl/git/pylot/pylot/correlation/submit_pick_corr_correction.sh'
 ####
 
 for fnin, istart in fnames:
     input_cmds = f'qsub -q low.q@minos15,low.q@minos14,low.q@minos13,low.q@minos12,low.q@minos11 {script_location} {fnin} {istart}'
 
     print(input_cmds)
-    print(subprocess.check_output(input_cmds.split()))
-
-
-
+    print(subprocess.check_output(input_cmds.split()))
\ No newline at end of file