<jupyter_start><jupyter_text>Flux point fitting Prerequisites- Some knowledge about retrieving information from catalogs, see [the catalogs tutorial](../../api/catalog.ipynb)  ContextSome high level studies do not rely on reduced datasets with their IRFs but directly on higher level products such as flux points. This is not ideal because flux points already contain some hypothesis for the underlying spectral shape and the uncertainties they carry are usually simplified (e.g. symmetric gaussian errors). Yet, this is an efficient way to combine heterogeneous data. **Objective: fit spectral models to combined Fermi-LAT and IACT flux points.** Proposed approachHere we will load, the spectral points from Fermi-LAT and TeV catalogs and fit them with various spectral models to find the best representation of the wide band spectrum. The central class we're going to use for this example analysis is:  - `~gammapy.datasets.FluxPointsDataset`In addition we will work with the following data classes:- `~gammapy.estimators.FluxPoints`- `~gammapy.catalog.SourceCatalogGammaCat`- `~gammapy.catalog.SourceCatalog3FHL`- `~gammapy.catalog.SourceCatalog3FGL`And the following spectral model classes:- `~gammapy.modeling.models.PowerLawSpectralModel`- `~gammapy.modeling.models.ExpCutoffPowerLawSpectralModel`- `~gammapy.modeling.models.LogParabolaSpectralModel`  SetupLet us start with the usual IPython notebook and Python imports:<jupyter_code>%matplotlib inline
import matplotlib.pyplot as plt
from astropy import units as u
from gammapy.modeling.models import (
    PowerLawSpectralModel,
    ExpCutoffPowerLawSpectralModel,
    LogParabolaSpectralModel,
    SkyModel,
)
from gammapy.estimators import FluxPoints
from gammapy.datasets import FluxPointsDataset, Datasets
from gammapy.catalog import CATALOG_REGISTRY
from gammapy.modeling import Fit<jupyter_output><empty_output><jupyter_text>Load spectral pointsFor this analysis we choose to work with the source 'HESS J1507-622' and the associated Fermi-LAT sources '3FGL J1506.6-6219' and '3FHL J1507.9-6228e'. We load the source catalogs, and then access source of interest by name:<jupyter_code>catalog_3fgl = CATALOG_REGISTRY.get_cls("3fgl")()
catalog_3fhl = CATALOG_REGISTRY.get_cls("3fhl")()
catalog_gammacat = CATALOG_REGISTRY.get_cls("gamma-cat")()
source_fermi_3fgl = catalog_3fgl["3FGL J1506.6-6219"]
source_fermi_3fhl = catalog_3fhl["3FHL J1507.9-6228e"]
source_gammacat = catalog_gammacat["HESS J1507-622"]<jupyter_output><empty_output><jupyter_text>The corresponding flux points data can be accessed with `.flux_points` attribute:<jupyter_code>dataset_gammacat = FluxPointsDataset(data=source_gammacat.flux_points, name="gammacat")
dataset_gammacat.data.to_table(sed_type="dnde", formatted=True)
dataset_3fgl = FluxPointsDataset(data=source_fermi_3fgl.flux_points, name="3fgl")
dataset_3fgl.data.to_table(sed_type="dnde", formatted=True)
dataset_3fhl = FluxPointsDataset(data=source_fermi_3fhl.flux_points, name="3fhl")
dataset_3fhl.data.to_table(sed_type="dnde", formatted=True)<jupyter_output><empty_output><jupyter_text>Power Law FitFirst we start with fitting a simple `~gammapy.modeling.models.PowerLawSpectralModel`.<jupyter_code>pwl = PowerLawSpectralModel(
    index=2, amplitude="1e-12 cm-2 s-1 TeV-1", reference="1 TeV"
)
model = SkyModel(spectral_model=pwl, name="j1507-pl")<jupyter_output><empty_output><jupyter_text>After creating the model we run the fit by passing the `'flux_points'` and `'model'` objects:<jupyter_code>datasets = Datasets([dataset_gammacat, dataset_3fgl, dataset_3fhl])
datasets.models = model
print(datasets)
fitter = Fit()
result_pwl = fitter.run(datasets=datasets)<jupyter_output><empty_output><jupyter_text>And print the result:<jupyter_code>print(result_pwl["optimize_result"])
print(model)<jupyter_output><empty_output><jupyter_text>Finally we plot the data points and the best fit model:<jupyter_code>ax = plt.subplot()

kwargs = {
    "ax": ax,
    "sed_type": "e2dnde",
    "yunits": u.Unit("TeV cm-2 s-1")
}

for d in datasets:
    d.data.plot(label=d.name, **kwargs)

pwl.plot(energy_bounds=[1e-4, 1e2] * u.TeV, color="k", **kwargs)
pwl.plot_error(energy_bounds=[1e-4, 1e2] * u.TeV, **kwargs)
ax.set_ylim(1e-13, 1e-11);
ax.legend()<jupyter_output><empty_output><jupyter_text>Exponential Cut-Off Powerlaw FitNext we fit an `~gammapy.modeling.models.ExpCutoffPowerLawSpectralModel` law to the data.<jupyter_code>ecpl = ExpCutoffPowerLawSpectralModel(
    index=1.8,
    amplitude="2e-12 cm-2 s-1 TeV-1",
    reference="1 TeV",
    lambda_="0.1 TeV-1",
)
model = SkyModel(spectral_model=ecpl, name="j1507-ecpl")<jupyter_output><empty_output><jupyter_text>We run the fitter again by passing the flux points and the model instance:<jupyter_code>datasets.models = model
result_ecpl = fitter.run(datasets=datasets)
print(model)<jupyter_output><empty_output><jupyter_text>We plot the data and best fit model:<jupyter_code>ax = plt.subplot()

kwargs = {
    "ax": ax,
    "sed_type": "e2dnde",
    "yunits": u.Unit("TeV cm-2 s-1")
}

for d in datasets:
    d.data.plot(label=d.name,**kwargs)

ecpl.plot(energy_bounds=[1e-4, 1e2] * u.TeV, color="k", **kwargs)
ecpl.plot_error(energy_bounds=[1e-4, 1e2] * u.TeV, **kwargs)
ax.set_ylim(1e-13, 1e-11);
ax.legend()<jupyter_output><empty_output><jupyter_text>Log-Parabola FitFinally we try to fit a `~gammapy.modeling.models.LogParabolaSpectralModel` model:<jupyter_code>log_parabola = LogParabolaSpectralModel(
    alpha=2, amplitude="1e-12 cm-2 s-1 TeV-1", reference="1 TeV", beta=0.1
)
model = SkyModel(spectral_model=log_parabola, name="j1507-lp")
datasets.models = model
result_log_parabola = fitter.run(datasets=datasets)
print(model)
ax = plt.subplot()

kwargs = {
    "ax": ax,
    "sed_type": "e2dnde",
    "yunits": u.Unit("TeV cm-2 s-1")
}

for d in datasets:
    d.data.plot(label=d.name, **kwargs)

log_parabola.plot(energy_bounds=[1e-4, 1e2] * u.TeV, color="k", **kwargs)
log_parabola.plot_error(energy_bounds=[1e-4, 1e2] * u.TeV, **kwargs)
ax.set_ylim(1e-13, 1e-11);
ax.legend()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1169002473206925, 'fraction_numerical': 0.027535037098103874, 'mean_word_length': 6.5825, 'pattern_counts': {'":': 9, '<': 45, '<?xml version=': 0, '>': 45, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '82782', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2328, 'n_tokens_neox': 2230, 'n_words': 486, 'path': 'docs/tutorials/analysis/1D/sed_fitting.ipynb', 'repo_name': 'LauraOlivera/gammapy', 'size': '6065'}