.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "tutorials/analysis-1d/sed_fitting.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. or to run this example in your browser via Binder .. rst-class:: sphx-glr-example-title .. _sphx_glr_tutorials_analysis-1d_sed_fitting.py: Flux point fitting ================== Fit spectral models to combined Fermi-LAT and IACT flux points tables. Prerequisites ------------- - Some knowledge about retrieving information from catalogs, see :doc:`/tutorials/api/catalog` tutorial. Context ------- Some 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 approach ----------------- Here 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` .. GENERATED FROM PYTHON SOURCE LINES 53-58 Setup ----- Let us start with the usual IPython notebook and Python imports: .. GENERATED FROM PYTHON SOURCE LINES 58-73 .. code-block:: Python from astropy import units as u # %matplotlib inline import matplotlib.pyplot as plt from gammapy.catalog import CATALOG_REGISTRY from gammapy.datasets import Datasets, FluxPointsDataset from gammapy.modeling import Fit from gammapy.modeling.models import ( ExpCutoffPowerLawSpectralModel, LogParabolaSpectralModel, PowerLawSpectralModel, SkyModel, ) .. GENERATED FROM PYTHON SOURCE LINES 74-82 Load spectral points -------------------- For 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: .. GENERATED FROM PYTHON SOURCE LINES 82-92 .. code-block:: Python 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"] .. GENERATED FROM PYTHON SOURCE LINES 93-96 The corresponding flux points data can be accessed with ``.flux_points`` attribute: .. GENERATED FROM PYTHON SOURCE LINES 96-107 .. code-block:: Python 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) .. raw:: html
Table length=5
e_refe_mine_maxdndednde_errpdnde_errndnde_ulsqrt_tsis_ul
GeVGeVGeV1 / (cm2 GeV s)1 / (cm2 GeV s)1 / (cm2 GeV s)1 / (cm2 GeV s)
float64float64float64float64float64float64float64float32bool
14.14210.00020.0009.288e-122.343e-122.128e-12nan5.660False
31.62320.00050.0002.777e-126.572e-135.818e-13nan6.940False
86.60350.000150.0002.335e-131.055e-138.554e-14nan3.835False
273.861150.000500.0006.411e-142.697e-142.133e-14nan5.697False
1000.000500.0002000.0009.188e-214.034e-15nan8.068e-150.000True


.. GENERATED FROM PYTHON SOURCE LINES 108-114 Power Law Fit ------------- First we start with fitting a simple `~gammapy.modeling.models.PowerLawSpectralModel`. .. GENERATED FROM PYTHON SOURCE LINES 114-121 .. code-block:: Python pwl = PowerLawSpectralModel( index=2, amplitude="1e-12 cm-2 s-1 TeV-1", reference="1 TeV" ) model = SkyModel(spectral_model=pwl, name="j1507-pl") .. GENERATED FROM PYTHON SOURCE LINES 122-125 After creating the model we run the fit by passing the ``flux_points`` and ``model`` objects: .. GENERATED FROM PYTHON SOURCE LINES 125-134 .. code-block:: Python datasets = Datasets([dataset_gammacat, dataset_3fgl, dataset_3fhl]) datasets.models = model print(datasets) fitter = Fit() result_pwl = fitter.run(datasets=datasets) .. rst-class:: sphx-glr-script-out .. code-block:: none Datasets -------- Dataset 0: Type : FluxPointsDataset Name : gammacat Instrument : Models : ['j1507-pl'] Dataset 1: Type : FluxPointsDataset Name : 3fgl Instrument : Models : ['j1507-pl'] Dataset 2: Type : FluxPointsDataset Name : 3fhl Instrument : Models : ['j1507-pl'] .. GENERATED FROM PYTHON SOURCE LINES 135-137 And print the result: .. GENERATED FROM PYTHON SOURCE LINES 137-143 .. code-block:: Python print(result_pwl) print(model) .. rst-class:: sphx-glr-script-out .. code-block:: none OptimizeResult backend : minuit method : migrad success : True message : Optimization terminated successfully. nfev : 40 total stat : 28.29 CovarianceResult backend : minuit method : hesse success : True message : Hesse terminated successfully. SkyModel Name : j1507-pl Datasets names : None Spectral model type : PowerLawSpectralModel Spatial model type : Temporal model type : Parameters: index : 1.985 +/- 0.03 amplitude : 1.28e-12 +/- 1.6e-13 1 / (cm2 s TeV) reference (frozen): 1.000 TeV .. GENERATED FROM PYTHON SOURCE LINES 144-146 Finally we plot the data points and the best fit model: .. GENERATED FROM PYTHON SOURCE LINES 146-164 .. code-block:: Python ax = plt.subplot() ax.yaxis.set_units(u.Unit("TeV cm-2 s-1")) kwargs = {"ax": ax, "sed_type": "e2dnde"} for d in datasets: d.data.plot(label=d.name, **kwargs) energy_bounds = [1e-4, 1e2] * u.TeV pwl.plot(energy_bounds=energy_bounds, color="k", **kwargs) pwl.plot_error(energy_bounds=energy_bounds, **kwargs) ax.set_ylim(1e-13, 1e-11) ax.set_xlim(energy_bounds) ax.legend() plt.show() .. image-sg:: /tutorials/analysis-1d/images/sphx_glr_sed_fitting_001.png :alt: sed fitting :srcset: /tutorials/analysis-1d/images/sphx_glr_sed_fitting_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 165-172 Exponential Cut-Off Powerlaw Fit -------------------------------- Next we fit an `~gammapy.modeling.models.ExpCutoffPowerLawSpectralModel` law to the data. .. GENERATED FROM PYTHON SOURCE LINES 172-182 .. code-block:: Python 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") .. GENERATED FROM PYTHON SOURCE LINES 183-186 We run the fitter again by passing the flux points and the model instance: .. GENERATED FROM PYTHON SOURCE LINES 186-192 .. code-block:: Python datasets.models = model result_ecpl = fitter.run(datasets=datasets) print(model) .. rst-class:: sphx-glr-script-out .. code-block:: none SkyModel Name : j1507-ecpl Datasets names : None Spectral model type : ExpCutoffPowerLawSpectralModel Spatial model type : Temporal model type : Parameters: index : 1.894 +/- 0.05 amplitude : 1.96e-12 +/- 3.9e-13 1 / (cm2 s TeV) reference (frozen): 1.000 TeV lambda_ : 0.078 +/- 0.05 1 / TeV alpha (frozen): 1.000 .. GENERATED FROM PYTHON SOURCE LINES 193-195 We plot the data and best fit model: .. GENERATED FROM PYTHON SOURCE LINES 195-212 .. code-block:: Python ax = plt.subplot() kwargs = {"ax": ax, "sed_type": "e2dnde"} ax.yaxis.set_units(u.Unit("TeV cm-2 s-1")) for d in datasets: d.data.plot(label=d.name, **kwargs) ecpl.plot(energy_bounds=energy_bounds, color="k", **kwargs) ecpl.plot_error(energy_bounds=energy_bounds, **kwargs) ax.set_ylim(1e-13, 1e-11) ax.set_xlim(energy_bounds) ax.legend() plt.show() .. image-sg:: /tutorials/analysis-1d/images/sphx_glr_sed_fitting_002.png :alt: sed fitting :srcset: /tutorials/analysis-1d/images/sphx_glr_sed_fitting_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 213-219 Log-Parabola Fit ---------------- Finally we try to fit a `~gammapy.modeling.models.LogParabolaSpectralModel` model: .. GENERATED FROM PYTHON SOURCE LINES 219-245 .. code-block:: Python 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"} ax.yaxis.set_units(u.Unit("TeV cm-2 s-1")) for d in datasets: d.data.plot(label=d.name, **kwargs) log_parabola.plot(energy_bounds=energy_bounds, color="k", **kwargs) log_parabola.plot_error(energy_bounds=energy_bounds, **kwargs) ax.set_ylim(1e-13, 1e-11) ax.set_xlim(energy_bounds) ax.legend() plt.show() .. image-sg:: /tutorials/analysis-1d/images/sphx_glr_sed_fitting_003.png :alt: sed fitting :srcset: /tutorials/analysis-1d/images/sphx_glr_sed_fitting_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none SkyModel Name : j1507-lp Datasets names : None Spectral model type : LogParabolaSpectralModel Spatial model type : Temporal model type : Parameters: amplitude : 1.88e-12 +/- 2.8e-13 1 / (cm2 s TeV) reference (frozen): 1.000 TeV alpha : 2.144 +/- 0.07 beta : 0.049 +/- 0.02 .. GENERATED FROM PYTHON SOURCE LINES 246-256 Exercises --------- - Fit a `~gammapy.modeling.models.PowerLaw2SpectralModel` and `~gammapy.modeling.models.ExpCutoffPowerLaw3FGLSpectralModel` to the same data. - Fit a `~gammapy.modeling.models.ExpCutoffPowerLawSpectralModel` model to Vela X (‘HESS J0835-455’) only and check if the best fit values correspond to the values given in the Gammacat catalog .. GENERATED FROM PYTHON SOURCE LINES 259-271 What next? ---------- This was an introduction to SED fitting in Gammapy. - If you would like to learn how to perform a full Poisson maximum likelihood spectral fit, please check out the :doc:`/tutorials/analysis-1d/spectral_analysis` tutorial. - To learn how to combine heterogeneous datasets to perform a multi-instrument forward-folding fit see the :doc:`/tutorials/analysis-3d/analysis_mwl` tutorial. .. _sphx_glr_download_tutorials_analysis-1d_sed_fitting.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: binder-badge .. image:: images/binder_badge_logo.svg :target: https://mybinder.org/v2/gh/gammapy/gammapy-webpage/main?urlpath=lab/tree/notebooks/dev/tutorials/analysis-1d/sed_fitting.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: sed_fitting.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: sed_fitting.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_