diff --git a/gammapy/modeling/models/spectral.py b/gammapy/modeling/models/spectral.py index c2340f04d2..e1200de187 100644 --- a/gammapy/modeling/models/spectral.py +++ b/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): diff --git a/gammapy/modeling/models/tests/test_spectral.py b/gammapy/modeling/models/tests/test_spectral.py index 06ede9b903..20a2a92d5f 100644 --- a/gammapy/modeling/models/tests/test_spectral.py +++ b/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") diff --git a/gammapy/modeling/tests/test_serialize_yaml.py b/gammapy/modeling/tests/test_serialize_yaml.py index 7f5d6f2dd4..b7e8285913 100644 --- a/gammapy/modeling/tests/test_serialize_yaml.py +++ b/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"), diff --git a/gammapy/spectrum/flux_point.py b/gammapy/spectrum/flux_point.py index b7345f4559..4e1a8d0370 100644 --- a/gammapy/spectrum/flux_point.py +++ b/gammapy/spectrum/flux_point.py @@ -1289,13 +1289,15 @@ 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 @@ -1303,7 +1305,7 @@ def peek(self, method="diff/model", **kwargs): 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 diff --git a/gammapy/spectrum/tests/test_flux_point.py b/gammapy/spectrum/tests/test_flux_point.py index ed77c5016c..610fde83f5 100644 --- a/gammapy/spectrum/tests/test_flux_point.py +++ b/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)