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

3D simulation and fitting

This tutorial shows how to do a 3D map-based simulation and fit.

For a tutorial on how to do a 3D map analyse of existing data, see the analysis_3d tutorial.

This can be useful to do a performance / sensitivity study, or to evaluate the capabilities of Gammapy or a given analysis method. Note that is is a binned simulation as is e.g. done also in Sherpa for Chandra, not an event sampling and anbinned analysis as is done e.g. in the Fermi ST or ctools.

Imports and versions

In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
In [2]:
import numpy as np
import astropy.units as u
from astropy.coordinates import SkyCoord, Angle
from gammapy.irf import load_cta_irfs
from gammapy.maps import WcsGeom, MapAxis, WcsNDMap, Map
from gammapy.spectrum.models import PowerLaw
from gammapy.image.models import SkyGaussian
from gammapy.cube.models import SkyModel, SkyModels
from gammapy.cube import MapFit, MapEvaluator, PSFKernel
from gammapy.cube import make_map_exposure_true_energy, make_map_background_irf
In [3]:
!gammapy info --no-envvar --no-dependencies --no-system

Gammapy package:

        path                   : /Users/deil/work/code/gammapy-docs/build/dev/gammapy/gammapy
        version                : 0.9.dev7911
        githash                : ec8bb9fc2056c3cb2c6ab01a353cffa13d6d2430

Simulate

In [4]:
filename = (
    "$GAMMAPY_DATA/cta-1dc/caldb/data/cta/1dc/bcf/South_z20_50h/irf_file.fits"
)
irfs = load_cta_irfs(filename)
In [5]:
# Define sky model to simulate the data
spatial_model = SkyGaussian(lon_0="0.2 deg", lat_0="0.1 deg", sigma="0.3 deg")
spectral_model = PowerLaw(
    index=3, amplitude="1e-11 cm-2 s-1 TeV-1", reference="1 TeV"
)
sky_model = SkyModel(
    spatial_model=spatial_model, spectral_model=spectral_model
)
print(sky_model)
SkyModel

Parameters:

           name     value   error      unit      min max frozen
        --------- --------- ----- -------------- --- --- ------
            lon_0 2.000e-01   nan            deg nan nan  False
            lat_0 1.000e-01   nan            deg nan nan  False
            sigma 3.000e-01   nan            deg nan nan  False
            index 3.000e+00   nan                nan nan  False
        amplitude 1.000e-11   nan cm-2 s-1 TeV-1 nan nan  False
        reference 1.000e+00   nan            TeV nan nan   True
In [6]:
# Define map geometry
axis = MapAxis.from_edges(
    np.logspace(-1.0, 1.0, 10), unit="TeV", name="energy", interp="log"
)
geom = WcsGeom.create(
    skydir=(0, 0), binsz=0.02, width=(5, 4), coordsys="GAL", axes=[axis]
)
In [7]:
# Define some observation parameters
# Here we just have a single observation,
# we are not simulating many pointings / observations
pointing = SkyCoord(1, 0.5, unit="deg", frame="galactic")
livetime = 1 * u.hour
offset_max = 2 * u.deg
offset = Angle("2 deg")
In [8]:
exposure = make_map_exposure_true_energy(
    pointing=pointing, livetime=livetime, aeff=irfs["aeff"], geom=geom
)
exposure.slice_by_idx({"energy": 3}).plot(add_cbar=True);
../_images/notebooks_simulate_3d_11_0.png
In [9]:
background = make_map_background_irf(
    pointing=pointing, ontime=livetime, bkg=irfs["bkg"], geom=geom
)
background.slice_by_idx({"energy": 3}).plot(add_cbar=True);
../_images/notebooks_simulate_3d_12_0.png
In [10]:
psf = irfs["psf"].to_energy_dependent_table_psf(theta=offset)
psf_kernel = PSFKernel.from_table_psf(psf, geom, max_radius=0.3 * u.deg)
psf_kernel.psf_kernel_map.sum_over_axes().plot(stretch="log");
../_images/notebooks_simulate_3d_13_0.png
In [11]:
energy = axis.edges * axis.unit
edisp = irfs["edisp"].to_energy_dispersion(
    offset, e_reco=energy, e_true=energy
)
edisp.plot_matrix();
../_images/notebooks_simulate_3d_14_0.png
In [12]:
%%time
# The idea is that we have this class that can compute `npred`
# maps, i.e. "predicted counts per pixel" given the model and
# the observation infos: exposure, background, PSF and EDISP
evaluator = MapEvaluator(
    model=sky_model,
    exposure=exposure,
    background=background,
    psf=psf_kernel,
    edisp=edisp,
)
CPU times: user 20 µs, sys: 4 µs, total: 24 µs
Wall time: 11 µs
In [13]:
# Accessing and saving a lot of the following maps is for debugging.
# Just for a simulation one doesn't need to store all these things.
# dnde = evaluator.compute_dnde()
# flux = evaluator.compute_flux()
npred = evaluator.compute_npred()
npred_map = WcsNDMap(geom, npred)
In [14]:
npred_map.sum_over_axes().plot(add_cbar=True);
../_images/notebooks_simulate_3d_17_0.png
In [15]:
# This one line is the core of how to simulate data when
# using binned simulation / analysis: you Poisson fluctuate
# npred to obtain simulated observed counts.
# Compute counts as a Poisson fluctuation
rng = np.random.RandomState(seed=42)
counts = rng.poisson(npred)
counts_map = WcsNDMap(geom, counts)
In [16]:
counts_map.sum_over_axes().plot();
../_images/notebooks_simulate_3d_19_0.png

Fit

Now let’s analyse the simulated data. Here we just fit it again with the same model we had before, but you could do any analysis you like here, e.g. fit a different model, or do a region-based analysis, …

In [17]:
# Define sky model to fit the data
spatial_model = SkyGaussian(lon_0="0 deg", lat_0="0 deg", sigma="1 deg")
spectral_model = PowerLaw(
    index=2, amplitude="1e-11 cm-2 s-1 TeV-1", reference="1 TeV"
)
model = SkyModel(spatial_model=spatial_model, spectral_model=spectral_model)

# Impose that specrtal index remains within limits
spectral_model.parameters["index"].min = 0.0
spectral_model.parameters["index"].max = 10.0

print(model)
SkyModel

Parameters:

           name     value   error      unit         min       max    frozen
        --------- --------- ----- -------------- --------- --------- ------
            lon_0 0.000e+00   nan            deg       nan       nan  False
            lat_0 0.000e+00   nan            deg       nan       nan  False
            sigma 1.000e+00   nan            deg       nan       nan  False
            index 2.000e+00   nan                0.000e+00 1.000e+01  False
        amplitude 1.000e-11   nan cm-2 s-1 TeV-1       nan       nan  False
        reference 1.000e+00   nan            TeV       nan       nan   True
In [18]:
%%time
fit = MapFit(
    model=model,
    counts=counts_map,
    exposure=exposure,
    background=background,
    psf=psf_kernel,
)

result = fit.run(optimize_opts={"print_level": 1})

FCN = 237085.25425881322 TOTAL NCALL = 228 NCALLS = 228
EDM = 3.4943207004158642e-06 GOAL EDM = 1e-05 UP = 1.0
Valid Valid Param Accurate Covar PosDef Made PosDef
True True True True False
Hesse Fail HasCov Above EDM Reach calllim
False True False False
+ Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed?
0 par_000_lon_0 0.198526 0.00848876 No
1 par_001_lat_0 0.0940424 0.00845611 No
2 par_002_sigma -0.297599 0.00572867 No
3 par_003_index 3.00415 0.0281002 0 10 No
4 par_004_amplitude 0.978582 0.0456543 No
5 par_005_reference 1 1 Yes

CPU times: user 10.1 s, sys: 484 ms, total: 10.6 s
Wall time: 10.6 s

True model:

In [19]:
print(sky_model)
SkyModel

Parameters:

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

Best-fit model:

In [20]:
print(result.model)
SkyModel

Parameters:

           name     value      error        unit         min       max    frozen
        --------- ---------- --------- -------------- --------- --------- ------
            lon_0  1.985e-01 8.489e-03            deg       nan       nan  False
            lat_0  9.404e-02 8.456e-03            deg       nan       nan  False
            sigma -2.976e-01 5.729e-03            deg       nan       nan  False
            index  3.004e+00 2.810e-02                0.000e+00 1.000e+01  False
        amplitude  9.786e-12 4.565e-13 cm-2 s-1 TeV-1       nan       nan  False
        reference  1.000e+00 0.000e+00            TeV       nan       nan   True

Covariance:

           name     lon_0      lat_0      sigma      index    amplitude  reference
        --------- ---------- ---------- ---------- ---------- ---------- ---------
            lon_0  7.206e-05 -3.427e-06 -3.679e-07  2.135e-06 -4.423e-17 0.000e+00
            lat_0 -3.427e-06  7.151e-05  1.172e-06  3.088e-06 -8.719e-17 0.000e+00
            sigma -3.679e-07  1.172e-06  3.282e-05 -1.719e-06 -7.109e-16 0.000e+00
            index  2.135e-06  3.088e-06 -1.719e-06  7.896e-04 -1.070e-14 0.000e+00
        amplitude -4.423e-17 -8.719e-17 -7.109e-16 -1.070e-14  2.084e-25 0.000e+00
        reference  0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00 0.000e+00
In [21]:
# TODO: show e.g. how to make a residual image