This is a fixed-text formatted version of a Jupyter notebook

First analysis

This notebook shows a simple example of a Crab analysis using the H.E.S.S. DL3 data release 1. It reduces the data to cube datasets and performs a simple 3D model fitting of the Crab nebula.

It uses the high level Analysis class to orchestrate data reduction. In its current state, Analysis supports the standard analysis cases of joint or stacked 3D and 1D analyses. It is instantiated with an AnalysisConfig object that gives access to analysis parameters either directly or via a YAML config file.

To see what is happening under-the-hood and to get an idea of the internal API, a second notebook performs the same analysis without using the Analysis class.

We will first show how to configure and run a stacked 3D analysis. The structure of the tutorial follows a typical analysis:

  • Analysis configuration

  • Observation selection

  • Data reduction

  • Model fitting

  • Estimating flux points

Finally we will compare the results against a reference model.

Setup

[1]:
%matplotlib inline
import matplotlib.pyplot as plt
[2]:
from pathlib import Path
from astropy import units as u
from gammapy.analysis import Analysis, AnalysisConfig
from gammapy.modeling.models import create_crab_spectral_model

Analysis configuration

For configuration of the analysis we use the YAML data format. YAML is a machine readable serialisation format, that is also friendly for humans to read. In this tutorial we will write the configuration file just using Python strings, but of course the file can be created and modified with any text editor of your choice.

Here is what the configuration for our analysis looks like:

[3]:
config = AnalysisConfig()
# the AnalysisConfig gives access to the various parameters used from logging to reduced dataset geometries
print(config)
AnalysisConfig

    general:
        log: {level: info, filename: null, filemode: null, format: null, datefmt: null}
        outdir: .
    observations:
        datastore: $GAMMAPY_DATA/hess-dl3-dr1
        obs_ids: []
        obs_file: null
        obs_cone: {frame: null, lon: null, lat: null, radius: null}
        obs_time: {start: null, stop: null}
    datasets:
        type: 1d
        stack: true
        geom:
            wcs:
                skydir: {frame: null, lon: null, lat: null}
                binsize: 0.02 deg
                fov: {width: 5.0 deg, height: 5.0 deg}
                binsize_irf: 0.2 deg
            selection: {offset_max: 2.5 deg}
            axes:
                energy: {min: 0.1 TeV, max: 10.0 TeV, nbins: 30}
                energy_true: {min: 0.1 TeV, max: 10.0 TeV, nbins: 30}
        map_selection: [counts, exposure, background, psf, edisp]
        background: {method: reflected, exclusion: null}
        on_region: {frame: null, lon: null, lat: null, radius: null}
        containment_correction: true
    fit:
        fit_range: {min: 0.1 TeV, max: 10.0 TeV}
    flux_points:
        energy: {min: 0.1 TeV, max: 10.0 TeV, nbins: 30}

Setting the data to use

We want to use Crab runs from the H.E.S.S. DL3-DR1. We define here the datastore and a cone search of observations pointing with 5 degrees of the Crab nebula. Parameters can be set directly or as a python dict.

[4]:
# We define the datastore containing the data
config.observations.datastore = "$GAMMAPY_DATA/hess-dl3-dr1"

# We define the cone search parameters
config.observations.obs_cone.frame = "icrs"
config.observations.obs_cone.lon = "83.633 deg"
config.observations.obs_cone.lat = "22.014 deg"
config.observations.obs_cone.radius = "5 deg"

# Equivalently we could have set parameters with a python dict
# config.observations.obs_cone = {"frame": "icrs", "lon": "83.633 deg", "lat": "22.014 deg", "radius": "5 deg"}

Setting the reduced datasets geometry

[5]:
# We want to perform a 3D analysis
config.datasets.type = "3d"
# We want to stack the data into a single reduced dataset
config.datasets.stack = True

# We fiw the WCS geometry of the datasets
config.datasets.geom.wcs.skydir = {
    "lon": "83.633 deg",
    "lat": "22.014 deg",
    "frame": "icrs",
}
config.datasets.geom.wcs.fov = {"width": "2 deg", "height": "2 deg"}
config.datasets.geom.wcs.binsize = "0.02 deg"

# We now fix the energy axis for the counts map
config.datasets.geom.axes.energy.min = "1 TeV"
config.datasets.geom.axes.energy.max = "10 TeV"
config.datasets.geom.axes.energy.nbins = 4

# We now fix the energy axis for the IRF maps (exposure, etc)
config.datasets.geom.axes.energy_true.min = "0.5 TeV"
config.datasets.geom.axes.energy_true.max = "20 TeV"
config.datasets.geom.axes.energy.nbins = 10

Setting modeling and fitting parameters

Analysis can perform a few modeling and fitting tasks besides data reduction. Parameters have then to be passed to the configuration object.

[6]:
config.fit.fit_range.min = 1 * u.TeV
config.fit.fit_range.max = 10 * u.TeV
config.flux_points.energy = {"min": "1 TeV", "max": "10 TeV", "nbins": 3}

We’re all set. But before we go on let’s see how to save or import AnalysisConfig objects though YAML files.

Using YAML configuration files

One can export/import the AnalysisConfig to/from a YAML file.

[7]:
config.write("config.yaml", overwrite=True)
[8]:
config = AnalysisConfig.read("config.yaml")
print(config)
AnalysisConfig

    general:
        log: {level: info, filename: null, filemode: null, format: null, datefmt: null}
        outdir: .
    observations:
        datastore: $GAMMAPY_DATA/hess-dl3-dr1
        obs_ids: []
        obs_file: null
        obs_cone: {frame: icrs, lon: 83.633 deg, lat: 22.014 deg, radius: 5.0 deg}
        obs_time: {start: null, stop: null}
    datasets:
        type: 3d
        stack: true
        geom:
            wcs:
                skydir: {frame: icrs, lon: 83.633 deg, lat: 22.014 deg}
                binsize: 0.02 deg
                fov: {width: 2.0 deg, height: 2.0 deg}
                binsize_irf: 0.2 deg
            selection: {offset_max: 2.5 deg}
            axes:
                energy: {min: 1.0 TeV, max: 10.0 TeV, nbins: 10}
                energy_true: {min: 0.5 TeV, max: 20.0 TeV, nbins: 30}
        map_selection: [counts, exposure, background, psf, edisp]
        background: {method: reflected, exclusion: null}
        on_region: {frame: null, lon: null, lat: null, radius: null}
        containment_correction: true
    fit:
        fit_range: {min: 1.0 TeV, max: 10.0 TeV}
    flux_points:
        energy: {min: 1.0 TeV, max: 10.0 TeV, nbins: 3}

Running the analysis

We first create an gammapy.analysis.Analysis object from our configuration.

[9]:
analysis = Analysis(config)
Setting logging config: {'level': 'INFO', 'filename': None, 'filemode': None, 'format': None, 'datefmt': None}

Observation selection

We can directly select and load the observations from disk using gammapy.analysis.Analysis.get_observations():

[10]:
analysis.get_observations()
Fetching observations.
Number of selected observations: 4

The observations are now available on the Analysis object. The selection corresponds to the following ids:

[11]:
analysis.observations.ids
[11]:
['23592', '23523', '23526', '23559']

To see how to explore observations, please refer to the following notebook: CTA with Gammapy or HESS with Gammapy

Data reduction

Now we proceed to the data reduction. In the config file we have chosen a WCS map geometry, energy axis and decided to stack the maps. We can run the reduction using .get_datasets():

[12]:
%%time
analysis.get_datasets()
Creating geometry.
Creating datasets.
Processing observation 23592
Processing observation 23523
Processing observation 23526
Processing observation 23559
CPU times: user 2.29 s, sys: 200 ms, total: 2.49 s
Wall time: 2.49 s

As we have chosen to stack the data, there is finally one dataset contained which we can print:

[13]:
print(analysis.datasets["stacked"])
MapDataset

    Name                            : stacked

    Total counts                    : 2479
    Total predicted counts          : 2010.78
    Total background counts         : 2010.78

    Exposure min                    : 2.38e+08 m2 s
    Exposure max                    : 3.53e+09 m2 s

    Number of total bins            : 100000
    Number of fit bins              : 100000

    Fit statistic type              : cash
    Fit statistic value (-2 log(L)) : 22303.58

    Number of models                : 1
    Number of parameters            : 3
    Number of free parameters       : 1

    Component 0:
        Name                        : background
        Type                        : BackgroundModel
        Parameters:
            norm                    : 1.000
            tilt         (frozen)   : 0.000
            reference    (frozen)   : 1.000  TeV


As you can see the dataset comes with a predefined background model out of the data reduction, but no source model has been set yet.

The counts, exposure and background model maps are directly available on the dataset and can be printed and plotted:

[14]:
counts = analysis.datasets["stacked"].counts
counts.smooth("0.05 deg").plot_interactive()
../_images/notebooks_analysis_1_30_0.png

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:

[15]:
path = Path("analysis_1")
path.mkdir(exist_ok=True)

And then write the maps and IRFs to disk by calling the dedicated write() method:

[16]:
filename = path / "crab-stacked-dataset.fits.gz"
analysis.datasets[0].write(filename, overwrite=True)

Model fitting

Now we define a model to be fitted to the dataset. Here we use its YAML definition to load it:

[17]:
model_config = """
components:
- name: crab
  type: SkyModel
  spatial:
    type: PointSpatialModel
    frame: icrs
    parameters:
    - name: lon_0
      value: 83.63
      unit: deg
    - name: lat_0
      value: 22.14
      unit: deg
  spectral:
    type: PowerLawSpectralModel
    parameters:
    - name: amplitude
      value: 1.0e-12
      unit: cm-2 s-1 TeV-1
    - name: index
      value: 2.0
      unit: ''
    - name: reference
      value: 1.0
      unit: TeV
      frozen: true
"""

Now we set the model on the analysis object:

[18]:
analysis.set_models(model_config)
Reading model.
SkyModels

Component 0: SkyModel

   name     value   error      unit         min        max    frozen
--------- --------- ----- -------------- ---------- --------- ------
    lon_0 8.363e+01   nan            deg        nan       nan  False
    lat_0 2.214e+01   nan            deg -9.000e+01 9.000e+01  False
    index 2.000e+00   nan                       nan       nan  False
amplitude 1.000e-12   nan cm-2 s-1 TeV-1        nan       nan  False
reference 1.000e+00   nan            TeV        nan       nan   True




Finally we run the fit:

[19]:
analysis.run_fit()
Fitting datasets.
OptimizeResult

        backend    : minuit
        method     : minuit
        success    : True
        message    : Optimization terminated successfully.
        nfev       : 331
        total stat : 20722.27

[20]:
print(analysis.fit_result)
OptimizeResult

        backend    : minuit
        method     : minuit
        success    : True
        message    : Optimization terminated successfully.
        nfev       : 331
        total stat : 20722.27

This is how we can write the model back to file again:

[21]:
filename = path / "model-best-fit.yaml"
analysis.models.write(filename, overwrite=True)
[22]:
!cat analysis_1/model-best-fit.yaml
components:
-   name: crab
    type: SkyModel
    spectral:
        type: PowerLawSpectralModel
        parameters:
        - {name: index, value: 2.5949081491767294, unit: '', min: .nan, max: .nan,
            frozen: false}
        - {name: amplitude, value: 4.736520791604545e-11, unit: cm-2 s-1 TeV-1, min: .nan,
            max: .nan, frozen: false}
        - {name: reference, value: 1.0, unit: TeV, min: .nan, max: .nan, frozen: true}
    spatial:
        type: PointSpatialModel
        frame: icrs
        parameters:
        - {name: lon_0, value: 83.61906728143654, unit: deg, min: .nan, max: .nan,
            frozen: false}
        - {name: lat_0, value: 22.02364538552178, unit: deg, min: -90.0, max: 90.0,
            frozen: false}

Flux points

[23]:
analysis.get_flux_points(source="crab")
Calculating flux points.

      e_ref               ref_flux        ...        dnde_err        is_ul
       TeV              1 / (cm2 s)       ...    1 / (cm2 s TeV)
------------------ ---------------------- ... ---------------------- -----
1.4125375446227544 1.9829261209401868e-11 ...   1.28723694449008e-12 False
3.1622776601683795  7.597133159858373e-12 ... 2.1412700426918383e-13 False
  7.07945784384138  1.516599223581321e-12 ... 4.8784109538905096e-14 False
[24]:
plt.figure(figsize=(8, 5))
ax_sed, ax_residuals = analysis.flux_points.peek()
../_images/notebooks_analysis_1_47_0.png

The flux points can be exported to a fits table following the format defined here

[25]:
filename = path / "flux-points.fits"
analysis.flux_points.write(filename, overwrite=True)

What’s next

You can look at the same analysis without the high level interface in analysis_2

You can see how to perform a 1D spectral analysis of the same data in spectrum analysis

[ ]: