Repository: cta-observatory/pyirf

Issue:
Doc update (or functionality fix) is required in  pyirf.sensitivity.relative_sensitivity
The docstring of the method says:
```
    ----------
    n_on: int or array-like
        Number of signal-like events for the on observations
    n_off: int or array-like
        Number of signal-like events for the off observations
```
However the method itself can't work with arrays, only with scalars. The biggest issue with that is in scipy.optimize.brentq root finding algo, so fixing might be not straightforward. In case of leaving as is a doc update is needed. 

Discussion:

Patch:
diff --git a/pyirf/sensitivity.py b/pyirf/sensitivity.py
index 58bf3be..bda24c4 100644
--- a/pyirf/sensitivity.py
+++ b/pyirf/sensitivity.py
@@ -19,7 +19,7 @@ __all__ = ["relative_sensitivity", "calculate_sensitivity"]
 log = logging.getLogger(__name__)
 
 
-def relative_sensitivity(
+def _relative_sensitivity(
     n_on,
     n_off,
     alpha,
@@ -42,16 +42,19 @@ def relative_sensitivity(
     The reference time should be incorporated by appropriately weighting the events
     before calculating ``n_on`` and ``n_off``.
 
+    ``n_on``, ``n_off``, ``alpha`` and ``target_significance`` must be broadcastable
+    to a common shape.
+
     Parameters
     ----------
     n_on: int or array-like
         Number of signal-like events for the on observations
     n_off: int or array-like
         Number of signal-like events for the off observations
-    alpha: float
+    alpha: float or array-like
         Scaling factor between on and off observations.
         1 / number of off regions for wobble observations.
-    significance: float
+    significance: float or array-like
         Significance necessary for a detection
     significance_function: function
         A function f(n_on, n_off, alpha) -> significance in sigma
@@ -104,6 +107,12 @@ def relative_sensitivity(
     return relative_flux
 
 
+relative_sensitivity = np.vectorize(
+    _relative_sensitivity,
+    excluded=['significance_function']
+)
+
+
 def calculate_sensitivity(
     signal_hist,
     background_hist,
@@ -155,17 +164,15 @@ def calculate_sensitivity(
     """
     check_histograms(signal_hist, background_hist)
 
+    n_on = signal_hist["n_weighted"] + alpha * background_hist["n_weighted"]
+
+    # convert any quantities to arrays,
+    # since quantitites don't work with vectorize
+    n_on = u.Quantity(n_on, copy=False).to_value(u.one)
+    n_off = u.Quantity(background_hist["n_weighted"], copy=False).to_value(u.one)
+
     # calculate sensitivity in each bin
-    rel_sens = np.array([
-        relative_sensitivity(
-            n_on=n_signal + alpha * n_background,
-            n_off=n_background,
-            alpha=alpha,
-        )
-        for n_signal, n_background in zip(
-            signal_hist["n_weighted"], background_hist["n_weighted"]
-        )
-    ])
+    rel_sens = relative_sensitivity(n_on=n_on, n_off=n_off, alpha=alpha)
 
     # fill output table
     s = QTable()
diff --git a/pyirf/version.py b/pyirf/version.py
index 4b57329..260c070 100644
--- a/pyirf/version.py
+++ b/pyirf/version.py
@@ -1,1 +1,1 @@
-__version__ = "0.3.0.post1"
+__version__ = "0.3.1"


Test patch:
diff --git a/pyirf/tests/test_sensitivity.py b/pyirf/tests/test_sensitivity.py
index b8ed8ba..ae1287a 100644
--- a/pyirf/tests/test_sensitivity.py
+++ b/pyirf/tests/test_sensitivity.py
@@ -26,6 +26,13 @@ def test_relative_sensitivity():
     # no background, should work
     assert relative_sensitivity(10, 0, 0.2) > 0
 
+    # test vectorization
+
+    n_on = np.array([10, 20, 30])
+    n_off = np.array([10, 100])
+    s = relative_sensitivity(n_on, n_off[:, np.newaxis], alpha=0.2)
+    assert s.shape == (2, 3)
+
 
 def test_calculate_sensitivity():
     from pyirf.sensitivity import calculate_sensitivity