Add a norm parameter to the EBL model

gammapy/modeling/models/spectral.py

@@ -1410,7 +1410,22 @@ def evaluate(self, energy, parameter):
 
 
 class AbsorbedSpectralModel(SpectralModel):
-    """Spectral model with EBL absorption.
+    r"""Spectral model with EBL absorption.
+       The absorption of the EBL is described as:
+
+    .. math::
+        \exp{ \left ( -\alpha \times \tau(E,z) \right )}
+
+    where tau(E,z) is the optical depth predicted by the model (`~gammapy.modeling.models.Absorption`), which depends
+    on the energy of the gamma-rays and the redshift z of the source. The class `~gammapy.modeling.models.Absorption`
+    computes this correction
+
+    .. math:: \exp{ \left ( - \tau(E,z) \right )}
+
+    With the optical depth scaled by a factor alpha, the observed spectrum is formed as
+
+    .. math::
+        \mathrm{spectral \, model} \times \exp{ \left ( -\alpha \times \tau(E,z) \right )}
 
     Parameters
     ----------
@@ -1422,13 +1437,20 @@ class AbsorbedSpectralModel(SpectralModel):
         parameter value for absorption model
     parameter_name : str, optional
         parameter name
+    alpha_norm: float
+        Norm of the EBL model, by default 1.
     """
 
     __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,9 +1459,9 @@ 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)
+        parameters.extend([par, self.alpha_norm])
 
         super().__init__(parameters)
 
@@ -1449,9 +1471,11 @@ def evaluate(self, energy, **kwargs):
         # 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)
         absorption = self.absorption.evaluate(energy=energy, parameter=parameter)
+        absorption = np.power(absorption, self.alpha_norm.value)
         return flux * absorption
 
     def to_dict(self):

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

@@ -1289,21 +1289,23 @@ def residuals(self, method="diff"):
         residuals[fp.is_ul] = np.nan
         return residuals
 
-    def peek(self, method="diff/model", **kwargs):
+    def peek(self, method="diff/model", plot_error=True, **kwargs):
         """Plot flux points, best fit model and residuals.
 
         Parameters
         ----------
         method : {"diff", "diff/model", "diff/sqrt(model)"}
             Method used to compute the residuals, see `MapDataset.residuals()`
+        plot_error : bool
+            False to not plot the error band of the model
         """
         from matplotlib.gridspec import GridSpec
         import matplotlib.pyplot as plt
 
         gs = GridSpec(7, 1)
 
         ax_spectrum = plt.subplot(gs[:5, :])
-        self.plot_spectrum(ax=ax_spectrum, **kwargs)
+        self.plot_spectrum(ax=ax_spectrum, plot_error=plot_error, **kwargs)
 
         ax_spectrum.set_xticks([])
 
@@ -1381,7 +1383,9 @@ def plot_residuals(self, ax=None, method="diff", **kwargs):
         ax.set_ylim(-y_max, y_max)
         return ax
 
-    def plot_spectrum(self, ax=None, fp_kwargs=None, model_kwargs=None):
+    def plot_spectrum(
+        self, ax=None, fp_kwargs=None, plot_error=True, model_kwargs=None
+    ):
         """
         Plot spectrum including flux points and model.
 
@@ -1393,6 +1397,8 @@ def plot_spectrum(self, ax=None, fp_kwargs=None, model_kwargs=None):
             Keyword arguments passed to `FluxPoints.plot`.
         model_kwargs : dict
             Keywords passed to `SpectralModel.plot` and `SpectralModel.plot_error`
+        plot_error : bool
+            False to not plot the error band of the model
 
         Returns
         -------
@@ -1426,7 +1432,8 @@ 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:
+        # TODO: remove this plot error option when there is a better manipulation of uncertainties
+        if self.model.parameters.covariance is not None and plot_error:
             self.model.plot_error(ax=ax, **plot_kwargs)
 
         # format axes

gammapy/spectrum/tests/test_flux_point.py

@@ -307,3 +307,4 @@ def test_fp_dataset_peek(fit):
 
         with mpl_plot_check():
             fp_dataset.peek(method="diff/model")
+            fp_dataset.peek(method="diff/model", plot_error=False)