add a norm to the EBL model

gammapy/modeling/models/spectral.py

@@ -1422,13 +1422,23 @@ 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.
+    alpha_norm_frozen: bool
+        False if you want to have it as a free parameter
     """
 
     __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,
+        alpha_norm_frozen=True,
     ):
         self.spectral_model = spectral_model
         self.absorption = absorption
@@ -1440,6 +1450,8 @@ def __init__(
 
         parameters = spectral_model.parameters.parameters.copy()
         parameters.append(par)
+        self.alpha_norm = Parameter("alpha_norm", alpha_norm, frozen=alpha_norm_frozen)
+        parameters.append(self.alpha_norm)
 
         super().__init__(parameters)
 
@@ -1449,9 +1461,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

@@ -130,8 +130,8 @@ def test_datasets_to_io(tmpdir):
     assert dataset0.model.skymodels[1].name == "gll_iem_v06_cutout"
 
     assert (
-        dataset0.model.skymodels[0].parameters["reference"]
-        is dataset1.model.skymodels[1].parameters["reference"]
+            dataset0.model.skymodels[0].parameters["reference"]
+            is dataset1.model.skymodels[1].parameters["reference"]
     )
 
     assert_allclose(
@@ -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 errors of the model since for some of them it is not handle yet properly in Gammapy
         """
         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 errors of the model since for some of them it is not handle yet properly in Gammapy
 
         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)