Repository: cta-observatory/pyirf

Issue:
angular resolution does not work when events are outside of the defined bins
```
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-15-873b0be41b34> in <module>
----> 1 angular_resolution(events, bins)

~/Uni/CTA/pyirf/pyirf/benchmarks/angular_resolution.py in angular_resolution(events, energy_bins, energy_type)
     56     index = by_bin.groups.keys["bin_index"]
     57     result["angular_resolution"][index] = by_bin["theta"].groups.aggregate(
---> 58         lambda x: np.quantile(x, ONE_SIGMA_QUANTILE)
     59     )
     60     return result

~/.local/anaconda/envs/pyirf/lib/python3.7/site-packages/astropy/table/column.py in __setitem__(self, index, value)
   1114         # Set items using a view of the underlying data, as it gives an
   1115         # order-of-magnitude speed-up. [#2994]
-> 1116         self.data[index] = value
   1117 
   1118     def _make_compare(oper):

IndexError: index 24 is out of bounds for axis 0 with size 24
```

Discussion:

Patch:
diff --git a/pyirf/benchmarks/angular_resolution.py b/pyirf/benchmarks/angular_resolution.py
index 7ae02d8..f98478e 100644
--- a/pyirf/benchmarks/angular_resolution.py
+++ b/pyirf/benchmarks/angular_resolution.py
@@ -6,7 +6,7 @@ import astropy.units as u
 from ..binning import calculate_bin_indices
 
 
-ONE_SIGMA_PERCENTILE = norm.cdf(1) - norm.cdf(-1)
+ONE_SIGMA_QUANTILE = norm.cdf(1) - norm.cdf(-1)
 
 
 def angular_resolution(
@@ -34,7 +34,6 @@ def angular_resolution(
         Table containing the 68% containment of the angular
         distance distribution per each reconstructed energy bin.
     """
-
     # create a table to make use of groupby operations
     table = Table(events[[f"{energy_type}_energy", "theta"]])
 
@@ -42,6 +41,9 @@ def angular_resolution(
         table[f"{energy_type}_energy"].quantity, energy_bins
     )
 
+    n_bins =  len(energy_bins) - 1
+    mask = (table["bin_index"] >= 0) & (table["bin_index"] < n_bins)
+
     result = Table()
     result[f"{energy_type}_energy_low"] = energy_bins[:-1]
     result[f"{energy_type}_energy_high"] = energy_bins[1:]
@@ -50,10 +52,10 @@ def angular_resolution(
     result["angular_resolution"] = np.nan * u.deg
 
     # use groupby operations to calculate the percentile in each bin
-    by_bin = table.group_by("bin_index")
+    by_bin = table[mask].group_by("bin_index")
 
     index = by_bin.groups.keys["bin_index"]
     result["angular_resolution"][index] = by_bin["theta"].groups.aggregate(
-        lambda x: np.percentile(x, 100 * ONE_SIGMA_PERCENTILE)
+        lambda x: np.quantile(x, ONE_SIGMA_QUANTILE)
     )
     return result
diff --git a/pyirf/benchmarks/energy_bias_resolution.py b/pyirf/benchmarks/energy_bias_resolution.py
index 43585b2..bb077ec 100644
--- a/pyirf/benchmarks/energy_bias_resolution.py
+++ b/pyirf/benchmarks/energy_bias_resolution.py
@@ -88,6 +88,8 @@ def energy_bias_resolution(
     table["bin_index"] = calculate_bin_indices(
         table[f"{energy_type}_energy"].quantity, energy_bins
     )
+    n_bins =  len(energy_bins) - 1
+    mask = (table["bin_index"] >= 0) & (table["bin_index"] < n_bins)
 
     result = Table()
     result[f"{energy_type}_energy_low"] = energy_bins[:-1]
@@ -98,7 +100,7 @@ def energy_bias_resolution(
     result["resolution"] = np.nan
 
     # use groupby operations to calculate the percentile in each bin
-    by_bin = table.group_by("bin_index")
+    by_bin = table[mask].group_by("bin_index")
 
     index = by_bin.groups.keys["bin_index"]
     result["bias"][index] = by_bin["rel_error"].groups.aggregate(bias_function)


Test patch:
diff --git a/pyirf/benchmarks/tests/test_angular_resolution.py b/pyirf/benchmarks/tests/test_angular_resolution.py
index ad26656..6623764 100644
--- a/pyirf/benchmarks/tests/test_angular_resolution.py
+++ b/pyirf/benchmarks/tests/test_angular_resolution.py
@@ -13,7 +13,12 @@ def test_angular_resolution():
     true_resolution = np.append(np.full(1000, TRUE_RES_1), np.full(1000, TRUE_RES_2))
 
     events = QTable({
-        'true_energy': np.append(np.full(1000, 5.0), np.full(1000, 50.0)) * u.TeV,
+        'true_energy': np.concatenate([
+            [0.5], # below bin 1 to test with underflow
+            np.full(999, 5.0),
+            np.full(999, 50.0),
+            [500], # above bin 2 to test with overflow
+        ]) * u.TeV,
         'theta': np.abs(np.random.normal(0, true_resolution)) * u.deg
     })
 
diff --git a/pyirf/benchmarks/tests/test_bias_resolution.py b/pyirf/benchmarks/tests/test_bias_resolution.py
index a0c13a2..77c3083 100644
--- a/pyirf/benchmarks/tests/test_bias_resolution.py
+++ b/pyirf/benchmarks/tests/test_bias_resolution.py
@@ -16,7 +16,12 @@ def test_energy_bias_resolution():
     true_bias = np.append(np.full(1000, TRUE_BIAS_1), np.full(1000, TRUE_BIAS_2))
     true_resolution = np.append(np.full(1000, TRUE_RES_1), np.full(1000, TRUE_RES_2))
 
-    true_energy = np.append(np.full(1000, 5.0), np.full(1000, 50.0)) * u.TeV
+    true_energy = np.concatenate([
+        [0.5], # below bin 1 to test with underflow
+        np.full(999, 5.0),
+        np.full(999, 50.0),
+        [500], # above bin 2 to test with overflow
+    ]) * u.TeV
     reco_energy = true_energy * (1 + np.random.normal(true_bias, true_resolution))
 
     events = QTable({