Add function to determine pivot energy for any spectral model

gammapy/modeling/models/spectral.py

@@ -213,6 +213,45 @@ def evaluate_error(self, energy, epsilon=1e-4):
         """
         return self._propagate_error(epsilon=epsilon, fct=self, energy=energy)
 
+    @property
+    def pivot_energy(self):
+        """The pivot or decorrelation energy, for a given spectral model calculated numerically.
+        It is defined as the energy at which the correlation between the spectral parameters is minimized.
+
+        Returns
+        -------
+        pivot energy : `~astropy.units.Quantity`
+            The energy at which the statistical error in the computed flux is smallest.
+            If no minimum is found, NaN will be returned.
+        """
+
+        x_unit = self.reference.unit
+
+        def min_func(x):
+            """Function to minimize"""
+            x = np.exp(x)
+            dnde, dnde_error = self.evaluate_error(x * x_unit)
+            return dnde_error / dnde
+
+        bounds = [np.log(self.reference.value) - 3, np.log(self.reference.value) + 3]
+
+        std = np.std(min_func(x=np.linspace(bounds[0], bounds[1], 100)))
+        if std < 1e-5:
+            log.warning(
+                "The relative error on the flux does not depend on energy. No pivot energy found."
+            )
+            return np.nan * x_unit
+
+        minimizer = scipy.optimize.minimize_scalar(min_func, bounds=bounds)
+
+        if not minimizer.success:
+            log.warning(
+                "No minima found in the relative error on the flux. Pivot energy computation failed."
+            )
+            return np.nan * x_unit
+        else:
+            return np.exp(minimizer.x) * x_unit
+
     def integral(self, energy_min, energy_max, **kwargs):
         r"""Integrate spectral model numerically if no analytical solution defined.
 
@@ -815,13 +854,18 @@ def inverse(self, value, *args):
 
     @property
     def pivot_energy(self):
-        r"""The decorrelation energy is defined as:
+        r"""The pivot or decorrelation energy is defined as:
 
         .. math::
 
             E_D = E_0 * \exp{cov(\phi_0, \Gamma) / (\phi_0 \Delta \Gamma^2)}
 
         Formula (1) in https://arxiv.org/pdf/0910.4881.pdf
+
+        Returns
+        -------
+        pivot energy : `~astropy.units.Quantity`
+            If no minimum is found, NaN will be returned.
         """
         index_err = self.index.error
         reference = self.reference.quantity
@@ -908,13 +952,18 @@ def inverse(self, value, *args):
 
     @property
     def pivot_energy(self):
-        r"""The decorrelation energy is defined as:
+        r"""The pivot or decorrelation energy is defined as:
 
         .. math::
 
             E_D = E_0 * \exp{cov(\phi_0, \Gamma) / (\phi_0 \Delta \Gamma^2)}
 
         Formula (1) in https://arxiv.org/pdf/0910.4881.pdf
+
+        Returns
+        -------
+        pivot energy : `~astropy.units.Quantity`
+            If no minimum is found, NaN will be returned.
         """
         tilt_err = self.tilt.error
         reference = self.reference.quantity

gammapy/modeling/models/tests/test_spectral.py

@@ -614,14 +614,29 @@ def test_ecpl_integrate():
 
 def test_pwl_pivot_energy():
     pwl = PowerLawSpectralModel(amplitude="5.35510540e-11 cm-2 s-1 TeV-1")
+    assert_quantity_allclose(pwl.pivot_energy, np.nan * u.TeV, rtol=1e-5)
 
     pwl.covariance = [
         [0.0318377**2, 6.56889442e-14, 0],
         [6.56889442e-14, 0, 0],
         [0, 0, 0],
     ]
+    assert_quantity_allclose(pwl.pivot_energy, 3.3540034 * u.TeV, rtol=1e-5)
 
-    assert_quantity_allclose(pwl.pivot_energy, 3.3540034240210987 * u.TeV)
+
+def test_num_pivot_energy():
+    lp = LogParabolaSpectralModel(
+        amplitude="5.82442e-11 cm-2 s-1 GeV-1",
+        reference="17.337 GeV",
+        alpha="1.9134",
+        beta="0.2217",
+    )
+    lp.amplitude.error = "2.8804e-12 cm-2 s-1 GeV-1"
+    assert_quantity_allclose(lp.pivot_energy, np.nan * u.GeV, rtol=1e-5)
+
+    lp.alpha.error = "0.1126"
+    lp.beta.error = "0.0670"
+    assert_quantity_allclose(lp.pivot_energy, 17.337042 * u.GeV, rtol=1e-5)
 
 
 def test_template_spectral_model_evaluate_tiny():