diff --git a/mlair/plotting/postprocessing_plotting.py b/mlair/plotting/postprocessing_plotting.py
index 282dc76bb0f0e09b14c67b7123c311081cd8a72a..a0c701ac5371dd70e56f759e4fd615cc6d9de1c2 100644
--- a/mlair/plotting/postprocessing_plotting.py
+++ b/mlair/plotting/postprocessing_plotting.py
@@ -610,7 +610,7 @@ class PlotCompetitiveSkillScore(AbstractPlotClass):  # pragma: no cover
 class PlotSectorialSkillScore(AbstractPlotClass):  # pragma: no cover
 
     def __init__(self, data: xr.DataArray, plot_folder: str = ".", model_setup: str = "NN", sampling: str = "daily",
-                 model_name_for_plots: Union[str, None] = None, ahead_dim: str = "ahead"):
+                 model_name_for_plots: Union[str, None] = None, ahead_dim: str = "ahead", ):
         """Initialise."""
         super().__init__(plot_folder, f"skill_score_sectorial_{model_setup}")
         self._model_setup = model_setup
@@ -618,7 +618,7 @@ class PlotSectorialSkillScore(AbstractPlotClass):  # pragma: no cover
         self._ahead_dim = ahead_dim
         self._labels = None
         self._model_name_for_plots = model_name_for_plots
-        self._data = self._prepare_data(data)
+        self._data, self._reference_model = self._prepare_data(data)
         self._plot()
         self._save()
         self.plot_name = self.plot_name + "_vertical"
@@ -628,24 +628,21 @@ class PlotSectorialSkillScore(AbstractPlotClass):  # pragma: no cover
     @TimeTrackingWrapper
     def _prepare_data(self, data: xr.DataArray):
         self._labels = [str(i) + self._sampling for i in data.coords[self._ahead_dim].values]
+        reference_model = data.attrs["reference_model"]
         logging.info(f"PlotSectorialSkillScore._prepare_data: shape of data (xarray) is {data.shape}\n dims: {data.dims}")
         data = data.to_dataframe("data")[["data"]].stack(level=0).reset_index(level=3, drop=True).reset_index(name="data")
-        return data
-
-#    def _prepare_data(self, data: xr.DataArray):
-#        self._labels = [str(i) + self._sampling for i in data.coords[self._ahead_dim].values]
-#        data = data.to_dataframe("data")[["data"]].stack(level=0).reset_index(level=3, drop=True).reset_index(name="data")
-#        return data
+        return data, reference_model
 
     def _plot(self):
         size = max([len(np.unique(self._data.sector)), 6])
         fig, ax = plt.subplots(figsize=(size, size * 0.8))
         data = self._data
-        sns.boxplot(x="sector", y="data", hue="ahead", data=data, whis=1.5, ax=ax, palette="Blues_d",
+        sns.boxplot(x="sector", y="data", hue="ahead", data=data, whis=1, ax=ax, palette="Blues_d",
                     showmeans=True, meanprops={"markersize": 3, "markeredgecolor": "k"}, flierprops={"marker": "."},
                     )
         ax.axhline(y=0, color="grey", linewidth=.5)
-        ax.set(ylabel="skill score", xlabel="sector", title="summary of all stations")
+        ax.set(ylabel=f"skill score ({self._model_setup} vs. {self._reference_model})", xlabel="sector",
+               title="summary of all stations")
         handles, _ = ax.get_legend_handles_labels()
         plt.xticks(rotation=45, horizontalalignment="right")
         ax.legend(handles, self._labels)
@@ -659,7 +656,8 @@ class PlotSectorialSkillScore(AbstractPlotClass):  # pragma: no cover
                     showmeans=True, meanprops={"markersize": 3, "markeredgecolor": "k"}, flierprops={"marker": "."},
                     )
         ax.axvline(x=0, color="grey", linewidth=.5)
-        ax.set(xlabel="skill score", ylabel="sector", title="summary of all stations")
+        ax.set(xlabel=f"skill score ({self._model_setup} vs. {self._reference_model})", ylabel="sector",
+               title="summary of all stations")
         handles, _ = ax.get_legend_handles_labels()
         ax.legend(handles, self._labels)
         plt.tight_layout()
diff --git a/mlair/run_modules/post_processing.py b/mlair/run_modules/post_processing.py
index 297f378dae13167d2ef9ccac79a86c3038a35718..7bc60fb063492bb6cb72c785302ed42324b8b488 100644
--- a/mlair/run_modules/post_processing.py
+++ b/mlair/run_modules/post_processing.py
@@ -953,7 +953,14 @@ class PostProcessing(RunEnvironment):
             external_data_expd.to_netcdf(os.path.join(path, f"forecasts_ds_{str(station)}_test.nc"))
     
     @TimeTrackingWrapper
-    def calculate_error_metrics_based_on_upstream_wind_dir(self):
+    def calculate_error_metrics_based_on_upstream_wind_dir(self, ref_name: str = "ols") -> xr.DataArray:
+        """
+        Calculates the error metrics (mse)/skill scores based on the wind sector at time t0.
+
+        The reference model's name is stored as xr attribute
+        :return:
+        :rtype:
+        """
         self.load_forecast_array_and_store_as_dataset()
         path = self.data_store.get("forecast_path")
         files = os.path.join(path, "forecasts_ds_*_test.nc")
@@ -965,15 +972,10 @@ class PostProcessing(RunEnvironment):
             h_sector_skill_scores.append(statistics.skill_score_based_on_mse(
                 ds.where(self.upstream_wind_sector.squeeze() == sec),
                 obs_name=self.observation_indicator, pred_name=self.model_display_name,
-                ref_name="ols").assign_coords({"sector": sec})
+                ref_name=ref_name).assign_coords({"sector": sec})
                                          )
-
-            #sec_coords = da.argwhere(self.upstream_wind_sector.squeeze() == sec)
-            #sector_collector[sec] = dict()
-            #for i, dim in enumerate(self.upstream_wind_sector.squeeze().dims):
-            #    sector_collector[sec][dim] = sec_coords[:,i]
         sector_skill_scores = xr.concat(h_sector_skill_scores, dim="sector")
-        #.to_dataframe("data")[["data"]].stack(level=0).reset_index(level=2, drop=True).reset_index(name="data")
+        sector_skill_scores = sector_skill_scores.assign_attrs({f"reference_model": ref_name})
         return sector_skill_scores