.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "tutorials/starting/analysis_2.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_starting_analysis_2.py: Low level API ============= Introduction to Gammapy analysis using the low level API. Prerequisites ------------- - Understanding the gammapy data workflow, in particular what are DL3 events and instrument response functions (IRF). - Understanding of the data reduction and modeling fitting process as shown in the analysis with the high level interface tutorial :doc:`/tutorials/starting/analysis_1` Context ------- This notebook is an introduction to gammapy analysis this time using the lower level classes and functions the library. This allows to understand what happens during two main gammapy analysis steps, data reduction and modeling/fitting. **Objective: Create a 3D dataset of the Crab using the H.E.S.S. DL3 data release 1 and perform a simple model fitting of the Crab nebula using the lower level gammapy API.** Proposed approach ----------------- Here, we have to interact with the data archive (with the `~gammapy.data.DataStore`) to retrieve a list of selected observations (`~gammapy.data.Observations`). Then, we define the geometry of the `~gammapy.datasets.MapDataset` object we want to produce and the maker object that reduce an observation to a dataset. We can then proceed with data reduction with a loop over all selected observations to produce datasets in the relevant geometry and stack them together (i.e.sum them all). In practice, we have to: - Create a `~gammapy.data.DataStore` pointing to the relevant data - Apply an observation selection to produce a list of observations, a `~gammapy.data.Observations` object. - Define a geometry of the Map we want to produce, with a sky projection and an energy range. - Create a `~gammapy.maps.MapAxis` for the energy - Create a `~gammapy.maps.WcsGeom` for the geometry - Create the necessary makers: - the map dataset maker `~gammapy.makers.MapDatasetMaker` - the background normalization maker, here a `~gammapy.makers.FoVBackgroundMaker` - and usually the safe range maker : `~gammapy.makers.SafeMaskMaker` - Perform the data reduction loop. And for every observation: - Apply the makers sequentially to produce the current `~gammapy.datasets.MapDataset` - Stack it on the target one. - Define the`~gammapy.modeling.models.SkyModel` to apply to the dataset. - Create a `~gammapy.modeling.Fit` object and run it to fit the model parameters - Apply a `~gammapy.estimators.FluxPointsEstimator` to compute flux points for the spectral part of the fit. Setup ----- First, we setup the analysis by performing required imports. .. GENERATED FROM PYTHON SOURCE LINES 73-96 .. code-block:: Python from pathlib import Path from astropy import units as u from astropy.coordinates import SkyCoord from regions import CircleSkyRegion # %matplotlib inline import matplotlib.pyplot as plt from gammapy.data import DataStore from gammapy.datasets import MapDataset from gammapy.estimators import FluxPointsEstimator from gammapy.makers import FoVBackgroundMaker, MapDatasetMaker, SafeMaskMaker from gammapy.maps import MapAxis, WcsGeom from gammapy.modeling import Fit from gammapy.modeling.models import ( FoVBackgroundModel, PointSpatialModel, PowerLawSpectralModel, SkyModel, ) from gammapy.utils.check import check_tutorials_setup from gammapy.visualization import plot_npred_signal .. GENERATED FROM PYTHON SOURCE LINES 97-99 Check setup ----------- .. GENERATED FROM PYTHON SOURCE LINES 99-103 .. code-block:: Python 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.20 machine : x86_64 system : Linux Gammapy package: version : 1.3 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.1 astropy : 5.2.2 regions : 0.8 click : 8.1.7 yaml : 6.0.2 IPython : 8.18.1 jupyterlab : not installed matplotlib : 3.9.2 pandas : not installed healpy : 1.17.3 iminuit : 2.30.1 sherpa : 4.16.1 naima : 0.10.0 emcee : 3.1.6 corner : 2.2.3 ray : 2.39.0 Gammapy environment variables: GAMMAPY_DATA : /home/runner/work/gammapy-docs/gammapy-docs/gammapy-datasets/1.3 .. GENERATED FROM PYTHON SOURCE LINES 104-110 Defining the datastore and selecting observations ------------------------------------------------- We first use the `~gammapy.data.DataStore` object to access the observations we want to analyse. Here the H.E.S.S. DL3 DR1. .. GENERATED FROM PYTHON SOURCE LINES 110-114 .. code-block:: Python data_store = DataStore.from_dir("$GAMMAPY_DATA/hess-dl3-dr1") .. GENERATED FROM PYTHON SOURCE LINES 115-122 We can now define an observation filter to select only the relevant observations. Here we use a cone search which we define with a python dict. We then filter the `ObservationTable` with `~gammapy.data.ObservationTable.select_observations`. .. GENERATED FROM PYTHON SOURCE LINES 122-133 .. code-block:: Python selection = dict( type="sky_circle", frame="icrs", lon="83.633 deg", lat="22.014 deg", radius="5 deg", ) selected_obs_table = data_store.obs_table.select_observations(selection) .. GENERATED FROM PYTHON SOURCE LINES 134-138 We can now retrieve the relevant observations by passing their ``obs_id`` to the `~gammapy.data.DataStore.get_observations` method. .. GENERATED FROM PYTHON SOURCE LINES 138-142 .. code-block:: Python observations = data_store.get_observations(selected_obs_table["OBS_ID"]) .. GENERATED FROM PYTHON SOURCE LINES 143-150 Preparing reduced datasets geometry ----------------------------------- Now we define a reference geometry for our analysis, We choose a WCS based geometry with a binsize of 0.02 deg and also define an energy axis: .. GENERATED FROM PYTHON SOURCE LINES 150-168 .. code-block:: Python energy_axis = MapAxis.from_energy_bounds(1.0, 10.0, 4, unit="TeV") geom = WcsGeom.create( skydir=(83.633, 22.014), binsz=0.02, width=(2, 2), frame="icrs", proj="CAR", axes=[energy_axis], ) # Reduced IRFs are defined in true energy (i.e. not measured energy). energy_axis_true = MapAxis.from_energy_bounds( 0.5, 20, 10, unit="TeV", name="energy_true" ) .. GENERATED FROM PYTHON SOURCE LINES 169-171 Now we can define the target dataset with this geometry. .. GENERATED FROM PYTHON SOURCE LINES 171-177 .. code-block:: Python stacked = MapDataset.create( geom=geom, energy_axis_true=energy_axis_true, name="crab-stacked" ) .. GENERATED FROM PYTHON SOURCE LINES 178-190 Data reduction -------------- Create the maker classes to be used ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The `~gammapy.makers.MapDatasetMaker` object is initialized as well as the `~gammapy.makers.SafeMaskMaker` that carries here a maximum offset selection. The `~gammapy.makers.FoVBackgroundMaker` utilised here has the default ``spectral_model`` but it is possible to set your own. For further details see the :doc:`FoV background `. .. GENERATED FROM PYTHON SOURCE LINES 190-202 .. code-block:: Python offset_max = 2.5 * u.deg maker = MapDatasetMaker() maker_safe_mask = SafeMaskMaker( methods=["offset-max", "aeff-max"], offset_max=offset_max ) circle = CircleSkyRegion(center=SkyCoord("83.63 deg", "22.14 deg"), radius=0.2 * u.deg) exclusion_mask = ~geom.region_mask(regions=[circle]) maker_fov = FoVBackgroundMaker(method="fit", exclusion_mask=exclusion_mask) .. GENERATED FROM PYTHON SOURCE LINES 203-206 Perform the data reduction loop ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 208-229 .. code-block:: Python for obs in observations: # First a cutout of the target map is produced cutout = stacked.cutout( obs.get_pointing_icrs(obs.tmid), width=2 * offset_max, name=f"obs-{obs.obs_id}" ) # A MapDataset is filled in this cutout geometry dataset = maker.run(cutout, obs) # The data quality cut is applied dataset = maker_safe_mask.run(dataset, obs) # fit background model dataset = maker_fov.run(dataset) print( f"Background norm obs {obs.obs_id}: {dataset.background_model.spectral_model.norm.value:.2f}" ) # The resulting dataset cutout is stacked onto the final one stacked.stack(dataset) print(stacked) .. rst-class:: sphx-glr-script-out .. code-block:: none Background norm obs 23523: 0.99 Background norm obs 23526: 1.08 Background norm obs 23559: 0.99 Background norm obs 23592: 1.10 MapDataset ---------- Name : crab-stacked Total counts : 2479 Total background counts : 2112.97 Total excess counts : 366.03 Predicted counts : 2112.97 Predicted background counts : 2112.97 Predicted excess counts : nan Exposure min : 3.75e+08 m2 s Exposure max : 3.48e+09 m2 s Number of total bins : 40000 Number of fit bins : 40000 Fit statistic type : cash Fit statistic value (-2 log(L)) : nan Number of models : 0 Number of parameters : 0 Number of free parameters : 0 .. GENERATED FROM PYTHON SOURCE LINES 230-233 Inspect the reduced dataset ~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 233-238 .. code-block:: Python stacked.counts.sum_over_axes().smooth(0.05 * u.deg).plot(stretch="sqrt", add_cbar=True) plt.show() .. image-sg:: /tutorials/starting/images/sphx_glr_analysis_2_001.png :alt: analysis 2 :srcset: /tutorials/starting/images/sphx_glr_analysis_2_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 239-246 Save dataset to disk -------------------- It is common to run the preparation step independent of the likelihood fit, because often the preparation of maps, PSF and energy dispersion is slow if you have a lot of data. We first create a folder: .. GENERATED FROM PYTHON SOURCE LINES 246-251 .. code-block:: Python path = Path("analysis_2") path.mkdir(exist_ok=True) .. GENERATED FROM PYTHON SOURCE LINES 252-255 And then write the maps and IRFs to disk by calling the dedicated `~gammapy.datasets.MapDataset.write` method: .. GENERATED FROM PYTHON SOURCE LINES 255-260 .. code-block:: Python filename = path / "crab-stacked-dataset.fits.gz" stacked.write(filename, overwrite=True) .. GENERATED FROM PYTHON SOURCE LINES 261-267 Define the model ---------------- We first define the model, a `~gammapy.modeling.models.SkyModel`, as the combination of a point source `~gammapy.modeling.models.SpatialModel` with a powerlaw `~gammapy.modeling.models.SpectralModel`: .. GENERATED FROM PYTHON SOURCE LINES 267-286 .. code-block:: Python target_position = SkyCoord(ra=83.63308, dec=22.01450, unit="deg") spatial_model = PointSpatialModel( lon_0=target_position.ra, lat_0=target_position.dec, frame="icrs" ) spectral_model = PowerLawSpectralModel( index=2.702, amplitude=4.712e-11 * u.Unit("1 / (cm2 s TeV)"), reference=1 * u.TeV, ) sky_model = SkyModel( spatial_model=spatial_model, spectral_model=spectral_model, name="crab" ) bkg_model = FoVBackgroundModel(dataset_name="crab-stacked") .. GENERATED FROM PYTHON SOURCE LINES 287-289 Now we assign this model to our reduced dataset: .. GENERATED FROM PYTHON SOURCE LINES 289-293 .. code-block:: Python stacked.models = [sky_model, bkg_model] .. GENERATED FROM PYTHON SOURCE LINES 294-303 Fit the model ------------- The `~gammapy.modeling.Fit` class is orchestrating the fit, connecting the ``stats`` method of the dataset to the minimizer. By default, it uses ``iminuit``. Its constructor takes a list of dataset as argument. .. GENERATED FROM PYTHON SOURCE LINES 305-309 .. code-block:: Python fit = Fit(optimize_opts={"print_level": 1}) result = fit.run([stacked]) .. GENERATED FROM PYTHON SOURCE LINES 310-313 The `~gammapy.modeling.FitResult` contains information about the optimization and parameter error calculation. .. GENERATED FROM PYTHON SOURCE LINES 313-317 .. code-block:: Python print(result) .. rst-class:: sphx-glr-script-out .. code-block:: none OptimizeResult backend : minuit method : migrad success : True message : Optimization terminated successfully. nfev : 129 total stat : 16240.82 CovarianceResult backend : minuit method : hesse success : True message : Hesse terminated successfully. .. GENERATED FROM PYTHON SOURCE LINES 318-321 The fitted parameters are visible from the `~astropy.modeling.models.Models` object. .. GENERATED FROM PYTHON SOURCE LINES 321-325 .. code-block:: Python print(stacked.models.to_parameters_table()) .. rst-class:: sphx-glr-script-out .. code-block:: none model type name value ... max frozen link prior ---------------- ---- --------- ---------- ... --------- ------ ---- ----- crab index 2.6024e+00 ... nan False crab amplitude 4.5924e-11 ... nan False crab reference 1.0000e+00 ... nan True crab lon_0 8.3619e+01 ... nan False crab lat_0 2.2024e+01 ... 9.000e+01 False crab-stacked-bkg norm 9.3514e-01 ... nan False crab-stacked-bkg tilt 0.0000e+00 ... nan True crab-stacked-bkg reference 1.0000e+00 ... nan True .. GENERATED FROM PYTHON SOURCE LINES 326-330 Here we can plot the number of predicted counts for each model and for the background in our dataset. In order to do this, we can use the `~gammapy.visualization.plot_npred_signal` function. .. GENERATED FROM PYTHON SOURCE LINES 330-335 .. code-block:: Python plot_npred_signal(stacked) plt.show() .. image-sg:: /tutorials/starting/images/sphx_glr_analysis_2_002.png :alt: analysis 2 :srcset: /tutorials/starting/images/sphx_glr_analysis_2_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 336-344 Inspecting residuals ~~~~~~~~~~~~~~~~~~~~ For any fit it is useful to inspect the residual images. We have a few options on the dataset object to handle this. First we can use `~gammapy.datasets.MapDataset.plot_residuals_spatial` to plot a residual image, summed over all energies: .. GENERATED FROM PYTHON SOURCE LINES 344-349 .. code-block:: Python stacked.plot_residuals_spatial(method="diff/sqrt(model)", vmin=-0.5, vmax=0.5) plt.show() .. image-sg:: /tutorials/starting/images/sphx_glr_analysis_2_003.png :alt: analysis 2 :srcset: /tutorials/starting/images/sphx_glr_analysis_2_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 350-353 In addition, we can also specify a region in the map to show the spectral residuals: .. GENERATED FROM PYTHON SOURCE LINES 353-363 .. code-block:: Python region = CircleSkyRegion(center=SkyCoord("83.63 deg", "22.14 deg"), radius=0.5 * u.deg) stacked.plot_residuals( kwargs_spatial=dict(method="diff/sqrt(model)", vmin=-0.5, vmax=0.5), kwargs_spectral=dict(region=region), ) plt.show() .. image-sg:: /tutorials/starting/images/sphx_glr_analysis_2_004.png :alt: analysis 2 :srcset: /tutorials/starting/images/sphx_glr_analysis_2_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 364-367 We can also directly access the ``.residuals()`` to get a map, that we can plot interactively: .. GENERATED FROM PYTHON SOURCE LINES 367-374 .. code-block:: Python residuals = stacked.residuals(method="diff") residuals.smooth("0.08 deg").plot_interactive( cmap="coolwarm", vmin=-0.2, vmax=0.2, stretch="linear", add_cbar=True ) plt.show() .. image-sg:: /tutorials/starting/images/sphx_glr_analysis_2_005.png :alt: analysis 2 :srcset: /tutorials/starting/images/sphx_glr_analysis_2_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none interactive(children=(SelectionSlider(continuous_update=False, description='Select energy:', layout=Layout(width='50%'), options=('1.00 TeV - 1.78 TeV', '1.78 TeV - 3.16 TeV', '3.16 TeV - 5.62 TeV', '5.62 TeV - 10.0 TeV'), style=SliderStyle(description_width='initial'), value='1.00 TeV - 1.78 TeV'), RadioButtons(description='Select stretch:', options=('linear', 'sqrt', 'log'), style=DescriptionStyle(description_width='initial'), value='linear'), Output()), _dom_classes=('widget-interact',)) .. GENERATED FROM PYTHON SOURCE LINES 375-378 Plot the fitted spectrum ------------------------ .. GENERATED FROM PYTHON SOURCE LINES 381-388 Making a butterfly plot ~~~~~~~~~~~~~~~~~~~~~~~ The `~gammapy.modeling.models.SpectralModel` component can be used to produce a, so-called, butterfly plot showing the envelope of the model taking into account parameter uncertainties: .. GENERATED FROM PYTHON SOURCE LINES 388-392 .. code-block:: Python spec = sky_model.spectral_model .. GENERATED FROM PYTHON SOURCE LINES 393-395 Now we can actually do the plot using the ``plot_error`` method: .. GENERATED FROM PYTHON SOURCE LINES 395-404 .. code-block:: Python energy_bounds = [1, 10] * u.TeV fig, ax = plt.subplots(figsize=(8, 6)) spec.plot(ax=ax, energy_bounds=energy_bounds, sed_type="e2dnde") spec.plot_error(ax=ax, energy_bounds=energy_bounds, sed_type="e2dnde") plt.show() .. image-sg:: /tutorials/starting/images/sphx_glr_analysis_2_006.png :alt: analysis 2 :srcset: /tutorials/starting/images/sphx_glr_analysis_2_006.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 405-415 Computing flux points ~~~~~~~~~~~~~~~~~~~~~ We can now compute some flux points using the `~gammapy.estimators.FluxPointsEstimator`. Besides the list of datasets to use, we must provide it the energy intervals on which to compute flux points as well as the model component name. .. GENERATED FROM PYTHON SOURCE LINES 415-419 .. code-block:: Python energy_edges = [1, 2, 4, 10] * u.TeV fpe = FluxPointsEstimator(energy_edges=energy_edges, source="crab") .. GENERATED FROM PYTHON SOURCE LINES 420-427 .. code-block:: Python flux_points = fpe.run(datasets=[stacked]) fig, ax = plt.subplots(figsize=(8, 6)) spec.plot(ax=ax, energy_bounds=energy_bounds, sed_type="e2dnde") spec.plot_error(ax=ax, energy_bounds=energy_bounds, sed_type="e2dnde") flux_points.plot(ax=ax, sed_type="e2dnde") plt.show() .. image-sg:: /tutorials/starting/images/sphx_glr_analysis_2_007.png :alt: analysis 2 :srcset: /tutorials/starting/images/sphx_glr_analysis_2_007.png :class: sphx-glr-single-img .. _sphx_glr_download_tutorials_starting_analysis_2.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/v1.3?urlpath=lab/tree/notebooks/1.3/tutorials/starting/analysis_2.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: analysis_2.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: analysis_2.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: analysis_2.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_