Return reconstructed energy in .get_bias_energy()

gammapy/irf/effective_area.py

@@ -247,7 +247,11 @@ def max_area(self):
         return cleaned_data.max()
 
     def find_energy(self, aeff, emin=None, emax=None):
-        """Find energy for given effective area.
+        """Find energy for a given effective area.
+
+        In case the solution is not unique, provide the `emin` or `emax` arguments
+        to limit the solution to the given range. By default the peak energy of the
+        effective area is chosen as `emax`.
 
         Parameters
         ----------
@@ -261,14 +265,16 @@ def find_energy(self, aeff, emin=None, emax=None):
         Returns
         -------
         energy : `~astropy.units.Quantity`
-            Energy corresponding to aeff
+            Energy corresponding to the given aeff.
         """
         from ..spectrum.models import TableModel
+        energy = self.energy.nodes
 
         if emin is None:
-            emin = self.energy.nodes[0]
+            emin = energy[0]
         if emax is None:
-            emax = self.energy.nodes[-1]
+            # use the peak effective area as a default for the energy maximum
+            emax = energy[np.argmax(self.data.data)]
 
         aeff_spectrum = TableModel(self.energy.nodes, self.data.data, values_scale="lin")
         return aeff_spectrum.inverse(aeff, emin=emin, emax=emax)

gammapy/irf/energy_dispersion.py

@@ -442,7 +442,11 @@ def get_bias(self, e_true):
         return bias
 
     def get_bias_energy(self, bias, emin=None, emax=None):
-        """Find energy correspoding to a given bias
+        """Find energy corresponding to a given bias.
+
+        In case the solution is not unique, provide the `emin` or `emax` arguments
+        to limit the solution to the given range.  By default the peak energy of the
+        bias is chosen as `emin`.
 
         Parameters
         ----------
@@ -456,19 +460,23 @@ def get_bias_energy(self, bias, emin=None, emax=None):
         Returns
         -------
         bias_energy : `~astropy.unit.Quantity`
-            Bias energy
+            Reconstructed energy corresponding to the given bias.
         """
         from ..spectrum.models import TableModel
-        energy = self.e_true.nodes
-        values = self.get_bias(energy)
+        e_true = self.e_true.nodes
+        values = self.get_bias(e_true)
 
         if emin is None:
-            emin = energy[np.nanargmax(values)]
+            # use the peak bias energy as default minimum
+            emin = e_true[np.nanargmax(values)]
         if emax is None:
-            emax = energy[-1]
+            emax = e_true[-1]
+
+        bias_spectrum = TableModel(e_true, values)
+        e_true = bias_spectrum.inverse(Quantity(bias), emin=emin, emax=emax)
 
-        bias_spectrum = TableModel(energy, values)
-        return bias_spectrum.inverse(Quantity(bias), emin=emin, emax=emax)
+        # return reconstructed energy
+        return e_true * (1 + bias)
 
     def get_mean(self, e_true):
         """Get mean reconstructed energy for a given true energy."""

gammapy/irf/tests/test_effective_area.py

@@ -101,7 +101,7 @@ def test_EffectiveAreaTable(tmpdir, aeff):
         assert elo_threshold.unit == "TeV"
         assert_allclose(elo_threshold.value, 0.554086, rtol=1e-3)
 
-        ehi_threshold = arf.find_energy(0.9 * arf.max_area, emin=30 * u.TeV)
+        ehi_threshold = arf.find_energy(0.9 * arf.max_area, emin=30 * u.TeV, emax=100 * u.TeV)
         assert ehi_threshold.unit == "TeV"
         assert_allclose(ehi_threshold.value, 53.347217, rtol=1e-3)
 

gammapy/spectrum/observation.py

@@ -201,12 +201,11 @@ def compute_energy_threshold(
         # vector, otherwise Sherpa will not understand the files
 
         # Low threshold
-        e_center = np.sqrt(self.e_true[0] * self.e_true[-1])
         if method_lo == "area_max":
             aeff_thres = kwargs["area_percent_lo"] / 100 * self.aeff.max_area
-            thres_lo = self.aeff.find_energy(aeff_thres, emax=e_center)
+            thres_lo = self.aeff.find_energy(aeff_thres)
         elif method_lo == "energy_bias":
-            thres_lo = self.edisp.get_bias_energy(kwargs["bias_percent_lo"] / 100, emax=e_center)
+            thres_lo = self.edisp.get_bias_energy(kwargs["bias_percent_lo"] / 100)
         elif method_lo == "none":
             thres_lo = self.e_true[0]
         else:

gammapy/spectrum/tests/test_observation.py

@@ -28,7 +28,7 @@ def test_spectrum_observation_1():
         npred=109.014552,
         excess=169.916667,
         excess_safe_range=116.33,
-        lo_threshold=8.912509e+08,
+        lo_threshold=1.e+09,
         hi_threshold=1e11,
     )
     tester = SpectrumObservationTester(obs, pars)
@@ -58,7 +58,7 @@ def test_spectrum_observation_2():
         npred=292.00223031875987,
         excess=824,
         excess_safe_range=824,
-        lo_threshold=0.012045,
+        lo_threshold=0.013219,
         hi_threshold=100,
     )
     tester = SpectrumObservationTester(obs, pars)
@@ -193,8 +193,8 @@ def test_energy_thresholds(self):
                 bias_percent_hi=10
             )
 
-            assert_allclose(self.obs.lo_threshold.value, self.vals["lo_threshold"], rtol=1e-5)
-            assert_allclose(self.obs.hi_threshold.value, self.vals["hi_threshold"], rtol=1e-5)
+            assert_allclose(self.obs.lo_threshold.value, self.vals["lo_threshold"], rtol=1e-4)
+            assert_allclose(self.obs.hi_threshold.value, self.vals["hi_threshold"], rtol=1e-4)
             self.obs.reset_thresholds()