From c33d6ff2d90b36f785a8107540632afdce196ffc Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Maximilian=20N=C3=B6the?= <maximilian.noethe@tu-dortmund.de>
Date: Mon, 5 Oct 2020 17:43:03 +0200
Subject: [PATCH] Fix column order and make sure spectrum is the same

---
 .../comparison_with_EventDisplay.ipynb        | 10 ++----
 examples/calculate_eventdisplay_irfs.py       | 10 ++++--
 pyirf/sensitivity.py                          | 35 ++++++++++---------
 3 files changed, 27 insertions(+), 28 deletions(-)

diff --git a/docs/notebooks/comparison_with_EventDisplay.ipynb b/docs/notebooks/comparison_with_EventDisplay.ipynb
index d8906fe00..4288f82c9 100644
--- a/docs/notebooks/comparison_with_EventDisplay.ipynb
+++ b/docs/notebooks/comparison_with_EventDisplay.ipynb
@@ -268,6 +268,7 @@
     "sensitivity['reco_energy_low'].info.format = '.3g'\n",
     "sensitivity['reco_energy_high'].info.format = '.3g'\n",
     "sensitivity['reco_energy_center'].info.format = '.3g'\n",
+    "sensitivity['n_signal'].info.format = '.1f'\n",
     "sensitivity['n_signal_weighted'].info.format = '.1f'\n",
     "sensitivity['n_background_weighted'].info.format = '.1f'\n",
     "sensitivity['n_background'].info.format = '.1f'\n",
@@ -679,13 +680,6 @@
     "\n",
     "None # to remove clutter by mpl objects"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {
@@ -704,7 +698,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.9"
+   "version": "3.7.7"
   }
  },
  "nbformat": 4,
diff --git a/examples/calculate_eventdisplay_irfs.py b/examples/calculate_eventdisplay_irfs.py
index c88e4c5b1..a58130d82 100644
--- a/examples/calculate_eventdisplay_irfs.py
+++ b/examples/calculate_eventdisplay_irfs.py
@@ -179,6 +179,7 @@ def main():
     )
     gammas["selected"] = gammas["selected_theta"] & gammas["selected_gh"]
 
+    # calculate sensitivity
     signal_hist = create_histogram_table(
         gammas[gammas["selected"]], bins=sensitivity_bins
     )
@@ -189,12 +190,15 @@ def main():
         alpha=ALPHA,
         background_radius=MAX_BG_RADIUS,
     )
-    sensitivity = calculate_sensitivity(signal_hist, background_hist, alpha=ALPHA)
+    sensitivity = calculate_sensitivity(
+        signal_hist, background_hist, alpha=ALPHA
+    )
 
     # scale relative sensitivity by Crab flux to get the flux sensitivity
+    spectrum = particles['gamma']['target_spectrum']
     for s in (sensitivity_step_2, sensitivity):
-        s["flux_sensitivity"] = s["relative_sensitivity"] * CRAB_HEGRA(
-            s["reco_energy_center"]
+        s["flux_sensitivity"] = (
+            s["relative_sensitivity"] * spectrum(s['reco_energy_center'])
         )
 
     log.info('Calculating IRFs')
diff --git a/pyirf/sensitivity.py b/pyirf/sensitivity.py
index e35f71d30..58bf3beda 100644
--- a/pyirf/sensitivity.py
+++ b/pyirf/sensitivity.py
@@ -155,17 +155,8 @@ def calculate_sensitivity(
     """
     check_histograms(signal_hist, background_hist)
 
-    s = QTable()
-    for key in ("low", "high", "center"):
-        k = "reco_energy_" + key
-        s[k] = signal_hist[k]
-
-    key = 'relative_sensitivity'
-    # add event number information
-    s["n_background"] = background_hist["n"]
-    s["n_background_weighted"] = background_hist["n_weighted"]
-
-    s[key] = [
+    # calculate sensitivity in each bin
+    rel_sens = np.array([
         relative_sensitivity(
             n_on=n_signal + alpha * n_background,
             n_off=n_background,
@@ -174,14 +165,22 @@ def calculate_sensitivity(
         for n_signal, n_background in zip(
             signal_hist["n_weighted"], background_hist["n_weighted"]
         )
-    ]
+    ])
+
+    # fill output table
+    s = QTable()
+    for key in ("low", "high", "center"):
+        k = "reco_energy_" + key
+        s[k] = signal_hist[k]
+
+    s["n_signal"] = signal_hist["n"] * rel_sens
+    s["n_background"] = background_hist["n"]
+    s["n_signal_weighted"] = signal_hist["n_weighted"] * rel_sens
+    s["n_background_weighted"] = background_hist["n_weighted"]
 
     # perform safety checks
     # we use the number of signal events at the flux level that yields
     # the target significance
-    s["n_signal"] = signal_hist["n"] * s[key]
-    s["n_signal_weighted"] = signal_hist["n_weighted"] * s[key]
-
     # safety checks according to the IRF document
     # at least ten signal events and the number of signal events
     # must be larger then five percent of the remaining background
@@ -190,10 +189,12 @@ def calculate_sensitivity(
     s['failed_checks'] = (
         (s['n_signal_weighted'] < 10)
         | ((s['n_signal_weighted'] < min_excess) << 1)
-        | ((s[key] > 1) << 2)
+        | ((rel_sens > 1) << 2)
     )
 
-    s[key][s['failed_checks'] > 0] = np.nan
+    rel_sens[s['failed_checks'] > 0] = np.nan
+
+    s["relative_sensitivity"] = rel_sens
 
     return s