Extracting flux points with a compound or a naima model in gammapy v0.19

[registerrier]

As you may have experienced, a bug in gammapy v0.19 (#3690) make it impossible to run FluxPointsEstimator on a model that contains a CompoundSpectralModel or a NaimaSpectralModel. While this will be corrected in the next gammapy release, this can be solved with the following trick:

If your compound model or NaimaSpectralModel is called spectrum simply add the amplitude attribute to it doing:

spectrum.amplitude = spectrum.parameters["amplitude"]

Here is a working (stupid) example:

import numpy as np
import astropy.units as u
from gammapy.modeling import Fit
from gammapy.modeling.models import SkyModel, PowerLawSpectralModel, EBLAbsorptionNormSpectralModel
from gammapy.datasets import SpectrumDatasetOnOff
from gammapy.estimators import FluxPointsEstimator

# read dataset from disk
dataset=SpectrumDatasetOnOff.read("$GAMMAPY_DATA/joint-crab/spectra/hess/pha_obs23523.fits")

# create spectral model
pl = PowerLawSpectralModel()
ebl = EBLAbsorptionNormSpectralModel.read_builtin(reference='dominguez', redshift=0.01)

compound_model = ebl*pl
compound_model.amplitude = compound_model.parameters["amplitude"] 

model = SkyModel(spectral_model = compound_model, name="compound")

dataset.models = [model]

fit_result = Fit().run([dataset])

fpe = FluxPointsEstimator(energy_edges=np.geomspace(0.2,10,6)*u.TeV, source="compound")
fpe.run([dataset])