diff --git a/dhi/plots/exclusion.py b/dhi/plots/exclusion.py
index 0df034c61949d7007af354f640d1d3d086a317e8..1a1d02f420017286ed4b3cc747ac50df305c807d 100644
--- a/dhi/plots/exclusion.py
+++ b/dhi/plots/exclusion.py
@@ -7,12 +7,10 @@ Exclusion result plots using ROOT.
 import math
 
 import numpy as np
-import scipy.interpolate
 
 from dhi.config import poi_data, campaign_labels, colors, br_hh_names
-from dhi.util import (
-    import_ROOT, DotDict, to_root_latex, create_tgraph, minimize_1d, try_int,
-)
+from dhi.util import import_ROOT, to_root_latex, create_tgraph, try_int
+from dhi.plots.limits import evaluate_limit_scan_1d
 from dhi.plots.likelihoods import evaluate_likelihood_scan_1d, evaluate_likelihood_scan_2d
 from dhi.plots.util import (
     use_style, draw_model_parameters, invert_graph, get_graph_points, get_contours, get_text_extent,
@@ -533,63 +531,6 @@ def plot_exclusion_and_bestfit_2d(
     canvas.SaveAs(path)
 
 
-def evaluate_limit_scan_1d(scan_values, limit_values):
-    """
-    Takes the results of an upper limit scan given by the *scan_values* and the corresponding *limit*
-    values, performs an interpolation and returns certain results of the scan in a dict.
-
-    The returned fields are:
-
-    - ``interp``: The generated interpolation function.
-    - ``excluded_ranges``: A list of 2-tuples denoting ranges in units of the poi where limits are
-      below one.
-    """
-    scan_values = np.array(scan_values)
-    limit_values = np.array(limit_values)
-
-    # first, obtain an interpolation function
-    mask = ~np.isnan(limit_values)
-    scan_values = scan_values[mask]
-    limit_values = limit_values[mask]
-    # interp = scipy.interpolate.interp1d(scan_values, limit_values, kind="cubic")
-    interp = scipy.interpolate.interp1d(scan_values, limit_values, kind="linear")
-
-    # interpolation bounds
-    bounds = (scan_values.min() + 1e-4, scan_values.max() - 1e-4)
-
-    # helper to find intersections with one given a starting point
-    def get_intersection(start):
-        objective = lambda x: (interp(x) - 1) ** 2.0
-        res = minimize_1d(objective, bounds, start=start)
-        return res.x[0] if res.status == 0 and (bounds[0] < res.x[0] < bounds[1]) else None
-
-    # get exclusion range edges from intersections and remove duplicates
-    rnd = lambda v: round(v, 3)
-    edges = [rnd(scan_values.min()), rnd(scan_values.max())]
-    for start in np.linspace(bounds[0], bounds[1], 100):
-        edge = get_intersection(start)
-        if edge is None:
-            continue
-        edge = rnd(edge)
-        if edge not in edges:
-            edges.append(edge)
-    edges.sort()
-
-    # create ranges consisting of two adjacent edges
-    ranges = [(edges[i - 1], edges[i]) for i in range(1, len(edges))]
-
-    # select those ranges whose central value is below 1
-    excluded_ranges = [
-        r for r in ranges
-        if interp(0.5 * (r[1] + r[0])) < 1
-    ]
-
-    return DotDict(
-        interp=interp,
-        excluded_ranges=excluded_ranges,
-    )
-
-
 def get_auto_contour_levels(values, steps=(1,)):
     min_value = min(values)
     max_value = max(values)
diff --git a/dhi/plots/likelihoods.py b/dhi/plots/likelihoods.py
index 2f838567257ef9794450e94a058ed8cbc26a5bdf..ed141eadc1f35c4d1249045d685f5f31b0702ae6 100644
--- a/dhi/plots/likelihoods.py
+++ b/dhi/plots/likelihoods.py
@@ -285,26 +285,9 @@ def plot_likelihood_scans_1d(
             r.setup_line(line, props={"LineColor": colors.black, "LineStyle": 2, "NDC": False})
             draw_objs.append(line)
 
-    # theory prediction with uncertainties
-    if theory_value:
-        has_thy_err = len(theory_value) == 3
-        if has_thy_err:
-            # theory graph
-            g_thy = create_tgraph(1, theory_value[0], y_min, theory_value[2], theory_value[1],
-                0, y_max_line)
-            r.setup_graph(g_thy, props={"LineColor": colors.red, "FillStyle": 1001,
-                "FillColor": colors.red_trans_50})
-            draw_objs.append((g_thy, "SAME,02"))
-            legend_entries.append((g_thy, "Theory prediction", "LF"))
-        # theory line
-        line_thy = ROOT.TLine(theory_value[0], y_min, theory_value[0], y_max_line)
-        r.setup_line(line_thy, props={"NDC": False}, color=colors.red)
-        draw_objs.append(line_thy)
-        if not has_thy_err:
-            legend_entries.append((line_thy, "Theory prediction", "L"))
-
     # perform scans and draw nll curves
-    for d, col, ms in zip(data, color_sequence[:len(data)], marker_sequence[:len(data)]):
+    for d, col, ms in zip(data[::-1], color_sequence[:len(data)][::-1],
+            marker_sequence[:len(data)][::-1]):
         # evaluate the scan, run interpolation and error estimation
         scan = evaluate_likelihood_scan_1d(d["values"][poi], d["values"]["dnll2"],
             poi_min=d["poi_min"])
@@ -315,7 +298,7 @@ def plot_likelihood_scans_1d(
         r.setup_graph(g_nll, props={"LineWidth": 2, "MarkerStyle": ms, "MarkerSize": 1.2},
             color=colors[col])
         draw_objs.append((g_nll, "SAME,CP" if show_points else "SAME,C"))
-        legend_entries.append((g_nll, to_root_latex(br_hh_names.get(d["name"], d["name"])),
+        legend_entries.insert(-1, (g_nll, to_root_latex(br_hh_names.get(d["name"], d["name"])),
             "LP" if show_points else "L"))
 
         # line for best fit value
@@ -323,6 +306,24 @@ def plot_likelihood_scans_1d(
         r.setup_line(line_fit, props={"LineWidth": 2, "NDC": False}, color=colors[col])
         draw_objs.append(line_fit)
 
+    # theory prediction with uncertainties
+    if theory_value:
+        has_thy_err = len(theory_value) == 3
+        if has_thy_err:
+            # theory graph
+            g_thy = create_tgraph(1, theory_value[0], y_min, theory_value[2], theory_value[1],
+                0, y_max_line)
+            r.setup_graph(g_thy, props={"LineColor": colors.red, "FillStyle": 1001,
+                "FillColor": colors.red_trans_50})
+            draw_objs.insert(-len(data), (g_thy, "SAME,02"))
+            legend_entries.append((g_thy, "Theory prediction", "LF"))
+        # theory line
+        line_thy = ROOT.TLine(theory_value[0], y_min, theory_value[0], y_max_line)
+        r.setup_line(line_thy, props={"NDC": False}, color=colors.red)
+        draw_objs.append(line_thy)
+        if not has_thy_err:
+            legend_entries.insert(-len(data), (line_thy, "Theory prediction", "L"))
+
     # legend
     legend_cols = min(int(math.ceil(len(legend_entries) / 4.)), 3)
     legend_rows = int(math.ceil(len(legend_entries) / float(legend_cols)))
@@ -590,7 +591,7 @@ def plot_likelihood_scans_2d(
     draw_objs.append((h_dummy, "HIST"))
 
     # loop through data entries
-    for d, (cont1, cont2), col in zip(data, contours, color_sequence[:len(data)]):
+    for d, (cont1, cont2), col in zip(data[::-1], contours[::-1], color_sequence[:len(data)][::-1]):
         # evaluate the scan
         scan = evaluate_likelihood_scan_2d(
             d["values"][poi1], d["values"][poi2], d["values"]["dnll2"],
@@ -605,7 +606,7 @@ def plot_likelihood_scans_2d(
             r.setup_graph(g2, props={"LineWidth": 2, "LineStyle": 2, "LineColor": colors[col]})
             draw_objs.append((g2, "SAME,C"))
         name = to_root_latex(br_hh_names.get(d["name"], d["name"]))
-        legend_entries.append((g1, name, "L"))
+        legend_entries.insert(-1, (g1, name, "L"))
 
         # best fit point
         g_fit = create_tgraph(1, scan.num1_min(), scan.num2_min())
diff --git a/dhi/plots/limits.py b/dhi/plots/limits.py
index abfe18b489cd25ba386bebc9d795d3a2faefc979..3c25f8156bc0e3d9cbf06989b6a2ce454e27c9c1 100644
--- a/dhi/plots/limits.py
+++ b/dhi/plots/limits.py
@@ -7,11 +7,14 @@ Limit plots using ROOT.
 import math
 
 import numpy as np
+import scipy.interpolate
 
 from dhi.config import (
     poi_data, br_hh_names, campaign_labels, colors, color_sequence, marker_sequence,
 )
-from dhi.util import import_ROOT, to_root_latex, create_tgraph, try_int, colored
+from dhi.util import (
+    import_ROOT, DotDict, to_root_latex, create_tgraph, try_int, colored, minimize_1d,
+)
 from dhi.plots.util import (
     use_style, draw_model_parameters, create_hh_process_label, determine_limit_digits,
     get_graph_points,
@@ -342,7 +345,7 @@ def plot_limit_scans(
             color=colors[col])
         draw_objs.append((g_exp, "SAME,CP" if show_points else "SAME,C"))
         name = names[n_graphs - i - 1]
-        legend_entries.append((g_exp, to_root_latex(br_hh_names.get(name, name)),
+        legend_entries.insert(0, (g_exp, to_root_latex(br_hh_names.get(name, name)),
             "LP" if show_points else "L"))
         y_max_value = max(y_max_value, max(ev["limit"]))
         y_min_value = min(y_min_value, min(ev["limit"]))
@@ -863,3 +866,60 @@ def plot_benchmark_limits(
     # save
     r.update_canvas(canvas)
     canvas.SaveAs(path)
+
+
+def evaluate_limit_scan_1d(scan_values, limit_values):
+    """
+    Takes the results of an upper limit scan given by the *scan_values* and the corresponding *limit*
+    values, performs an interpolation and returns certain results of the scan in a dict.
+
+    The returned fields are:
+
+    - ``interp``: The generated interpolation function.
+    - ``excluded_ranges``: A list of 2-tuples denoting ranges in units of the poi where limits are
+      below one.
+    """
+    scan_values = np.array(scan_values)
+    limit_values = np.array(limit_values)
+
+    # first, obtain an interpolation function
+    mask = ~np.isnan(limit_values)
+    scan_values = scan_values[mask]
+    limit_values = limit_values[mask]
+    # interp = scipy.interpolate.interp1d(scan_values, limit_values, kind="cubic")
+    interp = scipy.interpolate.interp1d(scan_values, limit_values, kind="linear")
+
+    # interpolation bounds
+    bounds = (scan_values.min() + 1e-4, scan_values.max() - 1e-4)
+
+    # helper to find intersections with one given a starting point
+    def get_intersection(start):
+        objective = lambda x: (interp(x) - 1) ** 2.0
+        res = minimize_1d(objective, bounds, start=start)
+        return res.x[0] if res.status == 0 and (bounds[0] < res.x[0] < bounds[1]) else None
+
+    # get exclusion range edges from intersections and remove duplicates
+    rnd = lambda v: round(v, 3)
+    edges = [rnd(scan_values.min()), rnd(scan_values.max())]
+    for start in np.linspace(bounds[0], bounds[1], 100):
+        edge = get_intersection(start)
+        if edge is None:
+            continue
+        edge = rnd(edge)
+        if edge not in edges:
+            edges.append(edge)
+    edges.sort()
+
+    # create ranges consisting of two adjacent edges
+    ranges = [(edges[i - 1], edges[i]) for i in range(1, len(edges))]
+
+    # select those ranges whose central value is below 1
+    excluded_ranges = [
+        r for r in ranges
+        if interp(0.5 * (r[1] + r[0])) < 1
+    ]
+
+    return DotDict(
+        interp=interp,
+        excluded_ranges=excluded_ranges,
+    )
diff --git a/dhi/plots/significances.py b/dhi/plots/significances.py
index fd0e79d63e535aa1b5779e48bb730db4272d63f9..b6a3a3918c91f5f40d3dbd7cb5c56d09fd5146b6 100644
--- a/dhi/plots/significances.py
+++ b/dhi/plots/significances.py
@@ -246,7 +246,7 @@ def plot_significance_scans(
             "MarkerSize": 1.2}, color=colors[col])
         draw_objs.append((g_exp, "SAME,CP" if show_points else "SAME,C"))
         name = names[n_graphs - i - 1]
-        legend_entries.append((g_exp, to_root_latex(br_hh_names.get(name, name)),
+        legend_entries.insert(0, (g_exp, to_root_latex(br_hh_names.get(name, name)),
             "LP" if show_points else "L"))
         y_max_value = max(y_max_value, max(ev["significance"]))
         y_min_value = min(y_min_value, min(ev["significance"]))
diff --git a/dhi/tasks/limits.py b/dhi/tasks/limits.py
index 5903b1064bb0db78a8a7ebcd3cc1db69583dd39f..589447f761d3cc664412caa77bad08a8dfd49247 100644
--- a/dhi/tasks/limits.py
+++ b/dhi/tasks/limits.py
@@ -415,6 +415,8 @@ class PlotMultipleUpperLimitsByModel(PlotUpperLimits, MultiHHModelTask):
 
     unblinded = None
 
+    allow_empty_hh_model = True
+
     def requires(self):
         return [
             [