.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "tutorials/analysis-3d/analysis_mwl.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-3d_analysis_mwl.py: Multi instrument joint 3D and 1D analysis ========================================= Joint 3D analysis using 3D Fermi datasets, a H.E.S.S. reduced spectrum and HAWC flux points. Prerequisites ------------- - Handling of Fermi-LAT data with Gammapy see the :doc:`/tutorials/data/fermi_lat` tutorial. - Knowledge of spectral analysis to produce 1D On-Off datasets, see the following :doc:`/tutorials/analysis-1d/spectral_analysis` tutorial. - Using flux points to directly fit a model (without forward-folding) from the :doc:`/tutorials/analysis-1d/sed_fitting` tutorial. Context ------- Some science studies require to combine heterogeneous data from various instruments to extract physical information. In particular, it is often useful to add flux measurements of a source at different energies to an analysis to better constrain the wide-band spectral parameters. This can be done using a joint fit of heterogeneous datasets. **Objectives: Constrain the spectral parameters of the gamma-ray emission from the Crab nebula between 10 GeV and 100 TeV, using a 3D Fermi dataset, a H.E.S.S. reduced spectrum and HAWC flux points.** Proposed approach ----------------- This tutorial illustrates how to perform a joint modeling and fitting of the Crab Nebula spectrum using different datasets. The spectral parameters are optimized by combining a 3D analysis of Fermi-LAT data, a ON/OFF spectral analysis of HESS data, and flux points from HAWC. In this tutorial we are going to use pre-made datasets. We prepared maps of the Crab region as seen by Fermi-LAT using the same event selection than the `3FHL catalog `__ (7 years of data with energy from 10 GeV to 2 TeV). For the HESS ON/OFF analysis we used two observations from the `first public data release `__ with a significant signal from energy of about 600 GeV to 10 TeV. These observations have an offset of 0.5° and a zenith angle of 45-48°. The HAWC flux points data are taken from a `recent analysis `__ based on 2.5 years of data with energy between 300 Gev and 300 TeV. The setup --------- .. GENERATED FROM PYTHON SOURCE LINES 53-63 .. code-block:: Python from pathlib import Path from astropy import units as u import matplotlib.pyplot as plt from gammapy.datasets import Datasets, FluxPointsDataset, SpectrumDatasetOnOff from gammapy.estimators import FluxPoints, FluxPointsEstimator from gammapy.maps import MapAxis from gammapy.modeling import Fit from gammapy.modeling.models import Models, create_crab_spectral_model .. GENERATED FROM PYTHON SOURCE LINES 64-66 Check setup ----------- .. GENERATED FROM PYTHON SOURCE LINES 66-72 .. code-block:: Python from gammapy.utils.check import check_tutorials_setup from gammapy.utils.scripts import make_path check_tutorials_setup() .. rst-class:: sphx-glr-script-out .. code-block:: none System: python_executable : /home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/bin/python python_version : 3.9.19 machine : x86_64 system : Linux Gammapy package: version : 1.3.dev241+g0271bebfc path : /home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/lib/python3.9/site-packages/gammapy Other packages: numpy : 1.26.4 scipy : 1.13.0 astropy : 5.2.2 regions : 0.8 click : 8.1.7 yaml : 6.0.1 IPython : 8.18.1 jupyterlab : not installed matplotlib : 3.8.4 pandas : not installed healpy : 1.16.6 iminuit : 2.25.2 sherpa : 4.16.0 naima : 0.10.0 emcee : 3.1.6 corner : 2.2.2 ray : 2.20.0 Gammapy environment variables: GAMMAPY_DATA : /home/runner/work/gammapy-docs/gammapy-docs/gammapy-datasets/dev .. GENERATED FROM PYTHON SOURCE LINES 73-82 Data and models files --------------------- The datasets serialization produce YAML files listing the datasets and models. In the following cells we show an example containing only the Fermi-LAT dataset and the Crab model. Fermi-LAT-3FHL_datasets.yaml: .. GENERATED FROM PYTHON SOURCE LINES 82-89 .. code-block:: Python path = make_path("$GAMMAPY_DATA/fermi-3fhl-crab/Fermi-LAT-3FHL_datasets.yaml") with path.open("r") as f: print(f.read()) .. rst-class:: sphx-glr-script-out .. code-block:: none datasets: - name: Fermi-LAT type: MapDataset filename: Fermi-LAT-3FHL_data_Fermi-LAT.fits .. GENERATED FROM PYTHON SOURCE LINES 90-98 We used as model a point source with a log-parabola spectrum. The initial parameters were taken from the latest Fermi-LAT catalog `4FGL `__, then we have re-optimized the spectral parameters for our dataset in the 10 GeV - 2 TeV energy range (fixing the source position). Fermi-LAT-3FHL_models.yaml: .. GENERATED FROM PYTHON SOURCE LINES 98-105 .. code-block:: Python path = make_path("$GAMMAPY_DATA/fermi-3fhl-crab/Fermi-LAT-3FHL_models.yaml") with path.open("r") as f: print(f.read()) .. rst-class:: sphx-glr-script-out .. code-block:: none components: - name: Crab Nebula type: SkyModel spectral: type: LogParabolaSpectralModel parameters: - name: amplitude value: 0.018182745349064267 unit: cm-2 s-1 TeV-1 min: .nan max: .nan frozen: false error: 0.003026327991562108 - name: reference value: 5.054833602905273e-05 unit: TeV min: .nan max: .nan frozen: true error: 0.0 - name: alpha value: 1.652368617859867 unit: '' min: .nan max: .nan frozen: false error: 0.05762513693893088 - name: beta value: 0.03921700077803329 unit: '' min: .nan max: .nan frozen: false error: 0.00521472221220211 spatial: type: PointSpatialModel frame: icrs parameters: - name: lon_0 value: 83.63310241699219 unit: deg min: .nan max: .nan frozen: true error: 0.0 - name: lat_0 value: 22.019899368286133 unit: deg min: -90.0 max: 90.0 frozen: true error: 0.0 - type: FoVBackgroundModel datasets_names: - Fermi-LAT spectral: type: PowerLawNormSpectralModel parameters: - name: norm value: 1.3004625872247901 unit: '' min: 0.0 max: .nan frozen: false error: 0.07512322002655547 - name: tilt value: 0.0 unit: '' min: .nan max: .nan frozen: true error: 0.0 - name: reference value: 1.0 unit: TeV min: .nan max: .nan frozen: true error: 0.0 .. GENERATED FROM PYTHON SOURCE LINES 106-115 Reading different datasets -------------------------- Fermi-LAT 3FHL: map dataset for 3D analysis ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ For now we let’s use the datasets serialization only to read the 3D `MapDataset` associated to Fermi-LAT 3FHL data and models. .. GENERATED FROM PYTHON SOURCE LINES 115-125 .. code-block:: Python path = Path("$GAMMAPY_DATA/fermi-3fhl-crab") filename = path / "Fermi-LAT-3FHL_datasets.yaml" datasets = Datasets.read(filename=filename) models = Models.read(path / "Fermi-LAT-3FHL_models.yaml") print(models) .. rst-class:: sphx-glr-script-out .. code-block:: none /home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/lib/python3.9/site-packages/gammapy/utils/scripts.py:65: UserWarning: Checksum verification failed for /home/runner/work/gammapy-docs/gammapy-docs/gammapy-datasets/dev/fermi-3fhl-crab/Fermi-LAT-3FHL_datasets.yaml. warnings.warn(f"Checksum verification failed for {filename}.", UserWarning) Models Component 0: SkyModel Name : Crab Nebula Datasets names : None Spectral model type : LogParabolaSpectralModel Spatial model type : PointSpatialModel Temporal model type : Parameters: amplitude : 1.82e-02 +/- 3.0e-03 1 / (cm2 s TeV) reference (frozen): 0.000 TeV alpha : 1.652 +/- 0.06 beta : 0.039 +/- 0.01 lon_0 (frozen): 83.633 deg lat_0 (frozen): 22.020 deg Component 1: FoVBackgroundModel Name : Fermi-LAT-bkg Datasets names : ['Fermi-LAT'] Spectral model type : PowerLawNormSpectralModel Parameters: norm : 1.300 +/- 0.08 tilt (frozen): 0.000 reference (frozen): 1.000 TeV .. GENERATED FROM PYTHON SOURCE LINES 126-128 We get the Crab model in order to share it with the other datasets .. GENERATED FROM PYTHON SOURCE LINES 128-132 .. code-block:: Python print(models["Crab Nebula"]) .. rst-class:: sphx-glr-script-out .. code-block:: none SkyModel Name : Crab Nebula Datasets names : None Spectral model type : LogParabolaSpectralModel Spatial model type : PointSpatialModel Temporal model type : Parameters: amplitude : 1.82e-02 +/- 3.0e-03 1 / (cm2 s TeV) reference (frozen): 0.000 TeV alpha : 1.652 +/- 0.06 beta : 0.039 +/- 0.01 lon_0 (frozen): 83.633 deg lat_0 (frozen): 22.020 deg .. GENERATED FROM PYTHON SOURCE LINES 133-142 HESS-DL3: 1D ON/OFF dataset for spectral fitting ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The ON/OFF datasets can be read from PHA files following the `OGIP standards `__. We read the PHA files from each observation, and compute a stacked dataset for simplicity. Then the Crab spectral model previously defined is added to the dataset. .. GENERATED FROM PYTHON SOURCE LINES 142-158 .. code-block:: Python datasets_hess = Datasets() for obs_id in [23523, 23526]: dataset = SpectrumDatasetOnOff.read( f"$GAMMAPY_DATA/joint-crab/spectra/hess/pha_obs{obs_id}.fits" ) datasets_hess.append(dataset) dataset_hess = datasets_hess.stack_reduce(name="HESS") datasets.append(dataset_hess) print(datasets) .. rst-class:: sphx-glr-script-out .. code-block:: none Datasets -------- Dataset 0: Type : MapDataset Name : Fermi-LAT Instrument : Models : Dataset 1: Type : SpectrumDatasetOnOff Name : HESS Instrument : Models : .. GENERATED FROM PYTHON SOURCE LINES 159-166 HAWC: 1D dataset for flux point fitting ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The HAWC flux point are taken from https://arxiv.org/pdf/1905.12518.pdf Then these flux points are read from a pre-made FITS file and passed to a `FluxPointsDataset` together with the source spectral model. .. GENERATED FROM PYTHON SOURCE LINES 166-180 .. code-block:: Python # read flux points from https://arxiv.org/pdf/1905.12518.pdf filename = "$GAMMAPY_DATA/hawc_crab/HAWC19_flux_points.fits" flux_points_hawc = FluxPoints.read( filename, reference_model=create_crab_spectral_model("meyer") ) dataset_hawc = FluxPointsDataset(data=flux_points_hawc, name="HAWC") datasets.append(dataset_hawc) print(datasets) .. rst-class:: sphx-glr-script-out .. code-block:: none Datasets -------- Dataset 0: Type : MapDataset Name : Fermi-LAT Instrument : Models : Dataset 1: Type : SpectrumDatasetOnOff Name : HESS Instrument : Models : Dataset 2: Type : FluxPointsDataset Name : HAWC Instrument : Models : .. GENERATED FROM PYTHON SOURCE LINES 181-189 Datasets serialization ---------------------- The `datasets` object contains each dataset previously defined. It can be saved on disk as datasets.yaml, models.yaml, and several data files specific to each dataset. Then the `datasets` can be rebuild later from these files. .. GENERATED FROM PYTHON SOURCE LINES 189-203 .. code-block:: Python path = Path("crab-3datasets") path.mkdir(exist_ok=True) filename = path / "crab_10GeV_100TeV_datasets.yaml" datasets.write(filename, overwrite=True) datasets = Datasets.read(filename) datasets.models = models print(datasets) .. rst-class:: sphx-glr-script-out .. code-block:: none /home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/lib/python3.9/site-packages/gammapy/utils/scripts.py:65: UserWarning: Checksum verification failed for crab-3datasets/crab_10GeV_100TeV_datasets.yaml. warnings.warn(f"Checksum verification failed for {filename}.", UserWarning) Datasets -------- Dataset 0: Type : MapDataset Name : Fermi-LAT Instrument : Models : ['Crab Nebula', 'Fermi-LAT-bkg'] Dataset 1: Type : SpectrumDatasetOnOff Name : HESS Instrument : Models : ['Crab Nebula'] Dataset 2: Type : FluxPointsDataset Name : HAWC Instrument : Models : ['Crab Nebula'] .. GENERATED FROM PYTHON SOURCE LINES 204-210 Joint analysis -------------- We run the fit on the `Datasets` object that include a dataset for each instrument .. GENERATED FROM PYTHON SOURCE LINES 212-217 .. code-block:: Python fit_joint = Fit() results_joint = fit_joint.run(datasets=datasets) print(results_joint) .. rst-class:: sphx-glr-script-out .. code-block:: none OptimizeResult backend : minuit method : migrad success : True message : Optimization terminated successfully. nfev : 302 total stat : -12697.22 CovarianceResult backend : minuit method : hesse success : True message : Hesse terminated successfully. .. GENERATED FROM PYTHON SOURCE LINES 218-220 Let’s display only the parameters of the Crab spectral model .. GENERATED FROM PYTHON SOURCE LINES 220-225 .. code-block:: Python crab_spec = datasets[0].models["Crab Nebula"].spectral_model print(crab_spec) .. rst-class:: sphx-glr-script-out .. code-block:: none LogParabolaSpectralModel type name value unit error ... frozen is_norm link prior ---- --------- ---------- -------------- --------- ... ------ ------- ---- ----- amplitude 3.9741e-03 cm-2 s-1 TeV-1 3.125e-04 ... False True reference 5.0548e-05 TeV 0.000e+00 ... True False alpha 1.2634e+00 1.707e-02 ... False False beta 6.1321e-02 9.448e-04 ... False False .. GENERATED FROM PYTHON SOURCE LINES 226-229 We can compute flux points for Fermi-LAT and HESS datasets in order plot them together with the HAWC flux point. .. GENERATED FROM PYTHON SOURCE LINES 229-246 .. code-block:: Python # compute Fermi-LAT and HESS flux points energy_edges = MapAxis.from_energy_bounds("10 GeV", "2 TeV", nbin=5).edges flux_points_fermi = FluxPointsEstimator( energy_edges=energy_edges, source="Crab Nebula", ).run([datasets["Fermi-LAT"]]) energy_edges = MapAxis.from_bounds(1, 15, nbin=6, interp="log", unit="TeV").edges flux_points_hess = FluxPointsEstimator( energy_edges=energy_edges, source="Crab Nebula", selection_optional=["ul"] ).run([datasets["HESS"]]) .. GENERATED FROM PYTHON SOURCE LINES 247-250 Now, let’s plot the Crab spectrum fitted and the flux points of each instrument. .. GENERATED FROM PYTHON SOURCE LINES 250-267 .. code-block:: Python # display spectrum and flux points fig, ax = plt.subplots(figsize=(8, 6)) energy_bounds = [0.01, 300] * u.TeV sed_type = "e2dnde" crab_spec.plot(ax=ax, energy_bounds=energy_bounds, sed_type=sed_type, label="Model") crab_spec.plot_error(ax=ax, energy_bounds=energy_bounds, sed_type=sed_type) flux_points_fermi.plot(ax=ax, sed_type=sed_type, label="Fermi-LAT") flux_points_hess.plot(ax=ax, sed_type=sed_type, label="HESS") flux_points_hawc.plot(ax=ax, sed_type=sed_type, label="HAWC") ax.set_xlim(energy_bounds) ax.legend() plt.show() .. image-sg:: /tutorials/analysis-3d/images/sphx_glr_analysis_mwl_001.png :alt: analysis mwl :srcset: /tutorials/analysis-3d/images/sphx_glr_analysis_mwl_001.png :class: sphx-glr-single-img .. _sphx_glr_download_tutorials_analysis-3d_analysis_mwl.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-3d/analysis_mwl.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: analysis_mwl.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: analysis_mwl.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_