diff --git a/pyirf/cut_optimization.py b/pyirf/cut_optimization.py
index 0c8e1073..d8c1a514 100644
--- a/pyirf/cut_optimization.py
+++ b/pyirf/cut_optimization.py
@@ -75,7 +75,7 @@ def optimize_cuts(
     fill_value = signal['gh_score'].max()
 
     sensitivities = []
-    cut_indicies = []
+    cut_indices = []
     n_theta_cuts = len(theta_cut_efficiencies)
     n_gh_cuts = len(gh_cut_efficiencies)
     n_cuts = len(multiplicity_cuts) * n_theta_cuts * n_gh_cuts
@@ -156,7 +156,7 @@ def optimize_cuts(
                     signal_hist, background_hist, alpha=alpha,
                     **kwargs,
                 )
-                cut_indicies.append((multiplicity_index, theta_index, gh_index))
+                cut_indices.append((multiplicity_index, theta_index, gh_index))
                 sensitivities.append(sensitivity)
                 bar.update(1)
 
@@ -192,7 +192,7 @@ def optimize_cuts(
             # if all are invalid, just use the first one
             best = 0
 
-        multiplicity_index, theta_index, gh_index = cut_indicies[best]
+        multiplicity_index, theta_index, gh_index = cut_indices[best]
 
         best_sensitivity[bin_id] = sensitivities[best][bin_id]
 
diff --git a/pyirf/tests/test_cuts.py b/pyirf/tests/test_cuts.py
index e5dfca6b..8648a043 100644
--- a/pyirf/tests/test_cuts.py
+++ b/pyirf/tests/test_cuts.py
@@ -178,3 +178,37 @@ def test_calculate_percentile_cuts_table():
         [dist1.ppf(0.68), dist2.ppf(0.68)],
         rtol=0.1,
     )
+
+
+
+
+def test_calculate_percentile_cuts_multiple():
+    from pyirf.cuts import calculate_percentile_cut
+
+    np.random.seed(0)
+
+    dist1 = norm(0, 1)
+    dist2 = norm(10, 1)
+    N = int(1e4)
+
+    values = np.append(dist1.rvs(size=N), dist2.rvs(size=N)) * u.deg
+    bin_values = np.append(np.zeros(N), np.ones(N)) * u.m
+    # add some values outside of binning to test that under/overflow are ignored
+    bin_values[10] = 5 * u.m
+    bin_values[30] = -1 * u.m
+
+    bins = [-0.5, 0.5, 1.5] * u.m
+
+    cuts = calculate_percentile_cut(values, bin_values, bins, fill_value=np.nan * u.deg, percentile=[50, 68, 95])
+    assert np.all(cuts["low"] == bins[:-1])
+    assert np.all(cuts["high"] == bins[1:])
+
+    np.testing.assert_allclose(
+        cuts["cut"].to_value(u.deg),
+        [
+            [dist1.ppf(0.5), dist1.ppf(0.68), dist1.ppf(0.95)],
+            [dist2.ppf(0.5), dist2.ppf(0.68), dist2.ppf(0.95)],
+        ],
+        rtol=0.1,
+        atol=0.1, # dist1.ppf(0.5) == 0, so we also need a non-zero atol 
+    )