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

Modeling and fitting 2D images using Gammapy

Prerequisites:

  • To understand how a generel modelling and fiiting works in gammapy, please refer to the analysis_3d tutorial

Context:

We often want the determine the position and morphology of an object. To do so, we don’t necessarily have to resort to a full 3D fitting but can perform a simple image fitting, in particular, in an energy range where the PSF does not vary strongly, or if we want to explore a possible energy dependance of the morphology.

Objective:

To localize a source and/or constrain its morphology.

Proposed approach:

The first step here, as in most analysis with DL3 data, is to create reduced datasets. For this, we will use the Analysis class to create a single set of stacked maps with a single bin in energy (thus, an image which behaves as a cube). This, we will then model with a spatial model of our choice, while keeping the spectral model fixed to an integrated power law.

Setup

As usual, we’ll start with some general imports…

[1]:
%matplotlib inline
import astropy.units as u
import numpy as np
from astropy.coordinates import SkyCoord
from astropy.time import Time
from regions import CircleSkyRegion
from astropy.coordinates import Angle

import logging

log = logging.getLogger(__name__)

Now let’s import gammapy specific classes and functions

[2]:
from gammapy.analysis import Analysis, AnalysisConfig

Creating the config file

Now, we create a config file for out analysis. You may load this from disc if you have a pre-defined config file.

Here, we use 3 simulated CTA runs of the galactic center.

[3]:
config = AnalysisConfig()
# Selecting the observations
config.observations.datastore = "$GAMMAPY_DATA/cta-1dc/index/gps/"
config.observations.obs_ids = [110380, 111140, 111159]

Technically, gammapy implements 2D analysis as a special case of 3D analysis (one one bin in energy). So, we must specify the type of analysis as 3D, and define the geometry of the analysis.

[4]:
config.datasets.type = "3d"
config.datasets.geom.wcs.skydir = {
    "lon": "0 deg",
    "lat": "0 deg",
    "frame": "galactic",
}  # The WCS geometry - centered on the galactic center
config.datasets.geom.wcs.fov = {"width": "10 deg", "height": "8 deg"}
config.datasets.geom.wcs.binsize = "0.02 deg"

# The FoV offset cut
config.datasets.geom.selection.offset_max = 3.5 * u.deg

# We now fix the energy axis for the counts map - (the reconstructed energy binning)
config.datasets.geom.axes.energy.min = "0.1 TeV"
config.datasets.geom.axes.energy.max = "10 TeV"
config.datasets.geom.axes.energy.nbins = 1

config.datasets.geom.wcs.binsize_irf = 0.2 * u.deg
[5]:
print(config)
AnalysisConfig

    general:
        log: {level: info, filename: null, filemode: null, format: null, datefmt: null}
        outdir: .
    observations:
        datastore: $GAMMAPY_DATA/cta-1dc/index/gps
        obs_ids: [110380, 111140, 111159]
        obs_file: null
        obs_cone: {frame: null, lon: null, lat: null, radius: null}
        obs_time: {start: null, stop: null}
    datasets:
        type: 3d
        stack: true
        geom:
            wcs:
                skydir: {frame: galactic, lon: 0.0 deg, lat: 0.0 deg}
                binsize: 0.02 deg
                fov: {width: 10.0 deg, height: 8.0 deg}
                binsize_irf: 0.2 deg
            selection: {offset_max: 3.5 deg}
            axes:
                energy: {min: 0.1 TeV, max: 10.0 TeV, nbins: 1}
                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}

Getting the reduced dataset

We now use the config file and create a single MapDataset containing counts, background, exposure, psf and edisp maps.

[6]:
%%time
analysis = Analysis(config)
analysis.get_observations()
analysis.get_datasets()
Setting logging config: {'level': 'INFO', 'filename': None, 'filemode': None, 'format': None, 'datefmt': None}
Fetching observations.
Number of selected observations: 3
Creating geometry.
Creating datasets.
Processing observation 110380
Processing observation 111140
Processing observation 111159
CPU times: user 6.84 s, sys: 1.64 s, total: 8.47 s
Wall time: 8.68 s
[7]:
print(analysis.datasets["stacked"])
MapDataset

    Name                            : stacked

    Total counts                    : 137137
    Total predicted counts          : 122006.90
    Total background counts         : 122006.90

    Exposure min                    : 1.59e+07 m2 s
    Exposure max                    : 1.87e+10 m2 s

    Number of total bins            : 200000
    Number of fit bins              : 164958

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

    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


The counts and background maps have only one bin in reconstructed energy. The exposure and IRF maps are in true energy, and hence, have a different binning based upon the binning of the IRFs. We need not bother about them presently.

[8]:
print(analysis.datasets["stacked"].counts)
WcsNDMap

        geom  : WcsGeom
        axes  : ['lon', 'lat', 'energy']
        shape : (500, 400, 1)
        ndim  : 3
        unit  :
        dtype : float32

[9]:
print(analysis.datasets["stacked"].background_model.map)
WcsNDMap

        geom  : WcsGeom
        axes  : ['lon', 'lat', 'energy']
        shape : (500, 400, 1)
        ndim  : 3
        unit  :
        dtype : float64

[10]:
print(analysis.datasets["stacked"].exposure)
WcsNDMap

        geom  : WcsGeom
        axes  : ['lon', 'lat', 'energy']
        shape : (500, 400, 30)
        ndim  : 3
        unit  : m2 s
        dtype : float32

Modelling

Now, we define a model to be fitted to the dataset. The important thing to note here is the dummy spectral model - an integrated powerlaw with only free normalisation. Here, we use its YAML definition to load it:

[11]:
model_config = """
components:
- name: GC-1
  type: SkyModel
  spatial:
    type: PointSpatialModel
    frame: galactic
    parameters:
    - name: lon_0
      value: 0.02
      unit: deg
    - name: lat_0
      value: 0.01
      unit: deg
  spectral:
    type: PowerLaw2SpectralModel
    parameters:
    - name: amplitude
      value: 1.0e-12
      unit: cm-2 s-1
    - name: index
      value: 2.0
      unit: ''
      frozen: true
    - name: emin
      value: 0.1
      unit: TeV
      frozen: true
    - name: emax
      value: 10.0
      unit: TeV
      frozen: true
"""
[12]:
analysis.set_models(model_config)
Reading model.
SkyModels

Component 0: SkyModel

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




We will freeze the parameters of the background

[13]:
analysis.datasets["stacked"].background_model.parameters["norm"].frozen = True
analysis.datasets["stacked"].background_model.parameters["tilt"].frozen = True
[14]:
# To run the fit
analysis.run_fit()
Fitting datasets.
OptimizeResult

        backend    : minuit
        method     : minuit
        success    : True
        message    : Optimization terminated successfully.
        nfev       : 146
        total stat : 275451.82

[15]:
# To see the best fit values along with the errors
analysis.fit_result.parameters.to_table()
[15]:
Table length=9
namevalueerrorunitminmaxfrozen
str9float64float64str8float64float64bool
lon_0-5.413e-022.341e-03degnannanFalse
lat_0-4.787e-022.328e-03deg-9.000e+019.000e+01False
amplitude3.240e-111.522e-12cm-2 s-1nannanFalse
index2.000e+000.000e+00nannanTrue
emin1.000e-010.000e+00TeVnannanTrue
emax1.000e+010.000e+00TeVnannanTrue
norm1.000e+000.000e+000.000e+00nanTrue
tilt0.000e+000.000e+00nannanTrue
reference1.000e+000.000e+00TeVnannanTrue
[ ]: