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

Try online You can also contribute with your own notebooks in this GitHub repository. Source files: cube_analysis_part2.ipynb | cube_analysis_part2.py

# Cube analysis with Gammapy (part 2)¶

## Introduction¶

In this tutorial we will learn how to compute a morphological and spectral fit simultanously.

This is part 2 of 2 for IACT cube analysis. If you haven’t prepared your cubes, PSF and EDISP yet, do part 1 first.

The fitting is done using Sherpa.

We will use the following classes:

We will use the cubes built on the 4 Crab observations of gammapy-extra.

You could use the cubes we just created with the notebook cube_analysis_part1.ipynb

## Setup¶

In [1]:

%matplotlib inline
import matplotlib.pyplot as plt

import logging
logging.getLogger("sherpa.fit").setLevel(logging.ERROR)

In [2]:

import numpy as np
from astropy import units as u
from astropy.io import fits
from astropy.coordinates import SkyCoord, Angle
from regions import CircleSkyRegion

from gammapy.extern.pathlib import Path
from gammapy.irf import EnergyDispersion
from gammapy.cube import SkyCube
from gammapy.cube.sherpa_ import (
CombinedModel3DInt,
CombinedModel3DIntConvolveEdisp,
NormGauss2DInt,
)

from sherpa.models import PowLaw1D, TableModel
from sherpa.estmethods import Covariance
from sherpa.stats import Cash
from sherpa.fit import Fit

WARNING: imaging routines will not be available,
failed to import sherpa.image.ds9_backend due to
'RuntimeErr: DS9Win unusable: Could not find xpaget on your PATH'

WARNING:sherpa.image:imaging routines will not be available,
failed to import sherpa.image.ds9_backend due to
'RuntimeErr: DS9Win unusable: Could not find xpaget on your PATH'

WARNING: failed to import sherpa.astro.xspec; XSPEC models will not be available

WARNING:sherpa.astro.all:failed to import sherpa.astro.xspec; XSPEC models will not be available


## 3D analysis assuming that there is no energy dispersion (perfect energy resolution)¶

### Load the different cubes needed for the analysis¶

We will use the Cubes build on the 4 Crab observations of gammapy-extra. You could use the Cubes we just created with the notebook cube_analysis.ipynb by changing the cube_directory by your local path.

• Counts cube
In [3]:

cube_dir = Path('\$GAMMAPY_EXTRA/test_datasets/cube')
# Transformation to a sherpa object
cube = counts_3d.to_sherpa_data3d(dstype='Data3DInt')

• Background Cube
In [4]:

bkg_3d = SkyCube.read(cube_dir / 'bkg_cube.fits')
bkg = TableModel('bkg')
bkg.ampl = 1
bkg.ampl.freeze()

• Exposure Cube
In [5]:

exposure_3d = SkyCube.read(cube_dir / 'exposure_cube.fits')
i_nan = np.where(np.isnan(exposure_3d.data))
exposure_3d.data[i_nan] = 0
# In order to have the exposure in cm2 s
exposure_3d.data = exposure_3d.data * u.Unit('m2 / cm2').to('')

• PSF Cube
In [6]:

psf_3d = SkyCube.read(cube_dir / 'psf_cube.fits')


### Setup Combined spatial and spectral model¶

In [7]:

# Define a 2D gaussian for the spatial model
spatial_model = NormGauss2DInt('spatial-model')

# Define a power law for the spectral model
spectral_model = PowLaw1D('spectral-model')

# Combine spectral and spatial model
coord = counts_3d.sky_image_ref.coordinates(mode="edges")
energies = counts_3d.energies(mode='edges').to("TeV")
# Here the source model will be convolve by the psf:
# PSF * (source_model * exposure)
source_model = CombinedModel3DInt(
coord=coord,
energies=energies,
use_psf=True,
exposure=exposure_3d,
psf=psf_3d,
spatial_model=spatial_model,
spectral_model=spectral_model,
)


### Set starting value¶

In [8]:

center = SkyCoord(83.633083, 22.0145, unit="deg").galactic
source_model.gamma = 2.2
source_model.xpos = center.l.value
source_model.ypos = center.b.value
source_model.fwhm = 0.12
source_model.ampl = 1.0


### Fit¶

In [9]:

# Define the model
flux_factor = 1e-11
model = bkg + flux_factor * source_model

In [10]:

# Fit to the counts Cube sherpa object
fit = Fit(
data=cube,
model=model,
stat=Cash(),
estmethod=Covariance(),
)
fit_results = fit.fit()
print(fit_results.format())

Method                = neldermead
Statistic             = cash
Initial fit statistic = 8614.61
Final fit statistic   = 8028.8 at function evaluation 619
Data points           = 12500
Degrees of freedom    = 12495
Change in statistic   = 585.806
spatial-model.xpos   184.195
spatial-model.ypos   -6.16917
spatial-model.ampl   6.16628
spatial-model.fwhm   0.0765675
spectral-model.gamma   2.30648

In [11]:

err_est_results = fit.est_errors()
print(err_est_results.format())

Confidence Method     = covariance
Iterative Fit Method  = None
Statistic             = cash
covariance 1-sigma (68.2689%) bounds:
Param            Best-Fit  Lower Bound  Upper Bound
-----            --------  -----------  -----------
spatial-model.xpos      184.195   -0.0128654    0.0128654
spatial-model.ypos     -6.16917   -0.0124289    0.0124289
spatial-model.ampl      6.16628    -0.308859     0.308859
spatial-model.fwhm    0.0765675     -0.01881      0.01881
spectral-model.gamma      2.30648   -0.0599405    0.0599405


For example if you want to exclude a region in the FoV of your cube, you just have to provide a SkyCube with the same dimension than the Counts cube and with 0 in the region you want to exlude and 1 outside. With this SkyCube mask you can select only some energy band for the fit or just some spatial region whatever the energy band or both.

Here this is a test case, there is no source to exlude in our FOV but we will create a mask that remove some events from the source. You just see at the end that the amplitude fitted in the 3D analysis is lower than the one when you use all the events from the Crab. The principle works even if this is not usefull here, just a testcase.

In [12]:

exclusion_region = CircleSkyRegion(
center=SkyCoord(83.60, 21.88, unit='deg'),
)



In [13]:

# Set the counts
cube = counts_3d.to_sherpa_data3d(dstype='Data3DInt')


Select only the background pixels of the cube of the selected region to create the background model

In [14]:

# Set the bkg and select only the data points of the selected region
bkg = TableModel('bkg')
bkg.ampl = 1
bkg.ampl.freeze()


Add the indices of the selected region of the Cube to combine the model

In [15]:

# The model is evaluated on all the points then it is compared with the data only on the selected_region
source_model = CombinedModel3DInt(
coord=coord,
energies=energies,
use_psf=True,
exposure=exposure_3d,
psf=psf_3d,
spatial_model=spatial_model,
spectral_model=spectral_model,
select_region=True,
index_selected_region=index_region_selected_3d,
)

# Set starting values
source_model.gamma = 2.2
source_model.xpos = center.l.value
source_model.ypos = center.b.value
source_model.fwhm = 0.12
source_model.ampl = 1.0

# Define the model
model = bkg + flux_factor * (source_model)

In [16]:

fit = Fit(
data=cube,
model=model,
stat=Cash(),
estmethod=Covariance(),
)
fit_results = fit.fit()
print(fit_results.format())

Method                = neldermead
Statistic             = cash
Initial fit statistic = 670.276
Final fit statistic   = 448.967 at function evaluation 553
Data points           = 390
Degrees of freedom    = 385
Change in statistic   = 221.309
spatial-model.xpos   184.607
spatial-model.ypos   -5.81938
spatial-model.ampl   10.4384
spatial-model.fwhm   0.00276986
spectral-model.gamma   2.41052

In [17]:

err_est_results = fit.est_errors()

print(err_est_results.format())

Confidence Method     = covariance
Iterative Fit Method  = None
Statistic             = cash
covariance 1-sigma (68.2689%) bounds:
Param            Best-Fit  Lower Bound  Upper Bound
-----            --------  -----------  -----------
spatial-model.xpos      184.607        -----        -----
spatial-model.ypos     -5.81938        -----        -----
spatial-model.ampl      10.4384        -----        -----
spatial-model.fwhm   0.00276986        -----        -----
spectral-model.gamma      2.41052        -----        -----


The fitted flux is less than in the previous example since here the mask remove some events from the Crab

## 3D analysis taking into account the energy dispersion¶

In [18]:

# Set the counts
cube = counts_3d.to_sherpa_data3d(dstype='Data3DInt')

# Set the bkg
bkg = TableModel('bkg')
bkg.ampl = 1
bkg.ampl.freeze()

# Set the exposure
i_nan = np.where(np.isnan(exposure_3d.data))
exposure_3d.data[i_nan] = 0
exposure_3d.data = exposure_3d.data * 1e4

# Set the mean psf model

# Load the mean rmf calculated for the 4 Crab runs

In [19]:

# Setup combined spatial and spectral model
spatial_model = NormGauss2DInt('spatial-model')
spectral_model = PowLaw1D('spectral-model')


### Add the mean RMF to the Combine3DInt object¶

The model is evaluated on the true energy bin then it is convolved by the energy dispersion to compare to the counts data in reconstructed energy

In [20]:

coord = counts_3d.sky_image_ref.coordinates(mode="edges")
energies = counts_3d.energies(mode='edges').to("TeV")
source_model = CombinedModel3DIntConvolveEdisp(
coord=coord,
energies=energies,
use_psf=True,
exposure=exposure_3d,
psf=psf_3d,
spatial_model=spatial_model,
spectral_model=spectral_model,
edisp=rmf.data.data,
)

# Set starting values
center = SkyCoord(83.633083, 22.0145, unit="deg").galactic
source_model.gamma = 2.2
source_model.xpos = center.l.value
source_model.ypos = center.b.value
source_model.fwhm = 0.12
source_model.ampl = 1.0

# Define the model
model = bkg + flux_factor * source_model

In [21]:

fit = Fit(
data=cube,
model=model,
stat=Cash(),
estmethod=Covariance(),
)
fit_results = fit.fit()
print(fit_results.format())

Method                = neldermead
Statistic             = cash
Initial fit statistic = 6415.62
Final fit statistic   = 5848.02 at function evaluation 704
Data points           = 12500
Degrees of freedom    = 12495
Change in statistic   = 567.598
spatial-model.xpos   184.192
spatial-model.ypos   -6.17565
spatial-model.ampl   6.22776
spatial-model.fwhm   0.0712976
spectral-model.gamma   2.26903

In [22]:

err_est_results = fit.est_errors()
print(err_est_results.format())

Confidence Method     = covariance
Iterative Fit Method  = None
Statistic             = cash
covariance 1-sigma (68.2689%) bounds:
Param            Best-Fit  Lower Bound  Upper Bound
-----            --------  -----------  -----------
spatial-model.xpos      184.192   -0.0117762    0.0117762
spatial-model.ypos     -6.17565   -0.0136609    0.0136609
spatial-model.ampl      6.22776    -0.312187     0.312187
spatial-model.fwhm    0.0712976   -0.0205746    0.0205746
spectral-model.gamma      2.26903   -0.0515044    0.0515044


## Discussion¶

Here the Cubes are constructed from dummy data. We don’t expect to find the good spectral shape or amplitude since there is a lot of problem with the PFS, rmf etc.. in these dummy data. On real data, we find the good Crab position and spectral shape for a power law or a power law with an exponential cutoff….

In [23]:

# TODO: there should be some visualisation of results.
# E.g. a plotted spectral model + butterfly
# Or a counts, model and residual image