Merge pull request #2454 from JouvinLea/change_for_ebl

Add a norm to the EBL model

gammapy/modeling/models/spectral.py

@@ -1275,6 +1275,8 @@ def evaluate(self, energy, norm):
 class Absorption:
     r"""Gamma-ray absorption models.
 
+    Usually used as part of `AbsorbedSpectralModel`.
+
     Parameters
     ----------
     energy : `~astropy.units.Quantity`
@@ -1410,25 +1412,43 @@ def evaluate(self, energy, parameter):
 
 
 class AbsorbedSpectralModel(SpectralModel):
-    """Spectral model with EBL absorption.
+    r"""Spectral model with EBL absorption.
+
+    The spectral model is evaluated, and then multiplied with an EBL
+    absorption factor given by
+
+    .. math::
+        \exp{ \left ( -\alpha \times \tau(E, z) \right )}
+
+    where :math:`\tau(E, z)` is the optical depth predicted by the model
+    (`Absorption`), which depends on the energy
+    of the gamma-rays and the redshift z of the source, and :math:`\alpha`
+    is a scale factor (default: 1) for the optical depth.
 
     Parameters
     ----------
-    spectral_model : `~gammapy.modeling.models.SpectralModel`
+    spectral_model : `SpectralModel`
         Spectral model.
-    absorption : `~gammapy.modeling.models.Absorption`
+    absorption : `Absorption`
         Absorption model.
     parameter : float
         parameter value for absorption model
     parameter_name : str, optional
         parameter name
+    alpha_norm: float
+        Norm of the EBL model
     """
 
     __slots__ = ["spectral_model", "absorption", "parameter", "parameter_name"]
     tag = "AbsorbedSpectralModel"
 
     def __init__(
-        self, spectral_model, absorption, parameter, parameter_name="redshift"
+        self,
+        spectral_model,
+        absorption,
+        parameter,
+        parameter_name="redshift",
+        alpha_norm=1.0,
     ):
         self.spectral_model = spectral_model
         self.absorption = absorption
@@ -1437,22 +1457,25 @@ def __init__(
         min_ = self.absorption.param.min()
         max_ = self.absorption.param.max()
         par = Parameter(parameter_name, parameter, min=min_, max=max_, frozen=True)
+        self.alpha_norm = Parameter("alpha_norm", alpha_norm, frozen=True)
 
-        parameters = spectral_model.parameters.parameters.copy()
-        parameters.append(par)
-
-        super().__init__(parameters)
+        super().__init__(
+            spectral_model.parameters.parameters.copy() + [par, self.alpha_norm]
+        )
 
     def evaluate(self, energy, **kwargs):
         """Evaluate the model at a given energy."""
-        # assign redshift value and remove it from dictionnary
+        # assign redshift value and remove it from dictionary
         # since it does not belong to the spectral model
         parameter = kwargs[self.parameter_name]
         del kwargs[self.parameter_name]
+        del kwargs["alpha_norm"]
 
-        flux = self.spectral_model.evaluate(energy=energy, **kwargs)
+        dnde = self.spectral_model.evaluate(energy=energy, **kwargs)
         absorption = self.absorption.evaluate(energy=energy, parameter=parameter)
-        return flux * absorption
+        # Power rule: (e ^ a) ^ b = e ^ (a * b)
+        absorption = np.power(absorption, self.alpha_norm.value)
+        return dnde * absorption
 
     def to_dict(self):
         return {

gammapy/modeling/models/tests/test_spectral.py

@@ -393,6 +393,21 @@ def test_absorption():
     )
     desired = u.Quantity(5.140765e-13, "TeV-1 s-1 cm-2")
     assert_quantity_allclose(model(1 * u.TeV), desired, rtol=1e-3)
+    assert model.alpha_norm.value == 1.0
+
+    # EBL + PWL model: test if norm of EBL=0: it mean model =pwl
+    model = AbsorbedSpectralModel(
+        spectral_model=pwl, absorption=absorption, alpha_norm=0, parameter=redshift
+    )
+    assert_quantity_allclose(model(1 * u.TeV), pwl(1 * u.TeV), rtol=1e-3)
+
+    # EBL + PWL model: Test with a norm different of 1
+    model = AbsorbedSpectralModel(
+        spectral_model=pwl, absorption=absorption, alpha_norm=1.5, parameter=redshift
+    )
+    desired = u.Quantity(2.739695e-13, "TeV-1 s-1 cm-2")
+    assert model.alpha_norm.value == 1.5
+    assert_quantity_allclose(model(1 * u.TeV), desired, rtol=1e-3)
 
 
 @requires_dependency("uncertainties")

gammapy/modeling/tests/test_serialize_yaml.py

@@ -163,10 +163,12 @@ def test_absorption_io(tmpdir):
     new_model = AbsorbedSpectralModel.from_dict(model_dict)
     assert new_model.parameter == 0.5
     assert new_model.parameter_name == "redshift"
+    assert new_model.alpha_norm.name == "alpha_norm"
+    assert new_model.alpha_norm.value == 1
     assert new_model.spectral_model.tag == "PowerLawSpectralModel"
     assert_allclose(new_model.absorption.energy, dominguez.energy)
     assert_allclose(new_model.absorption.param, dominguez.param)
-    assert len(new_model.parameters.parameters) == 4
+    assert len(new_model.parameters.parameters) == 5
 
     test_absorption = Absorption(
         u.Quantity(range(3), "keV"),

gammapy/spectrum/flux_point.py

@@ -1425,9 +1425,13 @@ def plot_spectrum(self, ax=None, fp_kwargs=None, model_kwargs=None):
 
         plot_kwargs.setdefault("color", ax.lines[-1].get_color())
         del plot_kwargs["label"]
-
+        
         if self.model.parameters.covariance is not None:
-            self.model.plot_error(ax=ax, **plot_kwargs)
+            try:
+                self.model.plot_error(ax=ax, **plot_kwargs)
+            except AttributeError:
+                log.debug("Model does not support evaluation of errors")
+                pass
 
         # format axes
         ax.set_xlim(self._e_range.to_value(self._e_unit))