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

[1]:
%matplotlib inline
[2]:
import numpy as np
import astropy.units as u
from astropy.coordinates import Angle
from gammapy.irf import load_cta_irfs
from gammapy.maps import WcsGeom, MapAxis, WcsNDMap
from gammapy.modeling.models import PowerLawSpectralModel
from gammapy.modeling.models import GaussianSpatialModel
from gammapy.modeling.models import SkyModel, BackgroundModel
from gammapy.cube import MapDataset, PSFKernel
from gammapy.cube import make_map_exposure_true_energy, make_map_background_irf
from gammapy.modeling import Fit
from gammapy.data import FixedPointingInfo
[3]:
!gammapy info --no-envvar --no-dependencies --no-system

Gammapy package:

        version                : 0.15.dev640+g0153d954a
        path                   : /Users/deil/work/code/gammapy-docs/build/dev/gammapy/gammapy

Simulate

[4]:
irfs = load_cta_irfs(
    "$GAMMAPY_DATA/cta-1dc/caldb/data/cta/1dc/bcf/South_z20_50h/irf_file.fits"
)
[5]:
# Define sky model to simulate the data
spatial_model = GaussianSpatialModel(
    lon_0="0.2 deg", lat_0="0.1 deg", sigma="0.3 deg", frame="galactic"
)
spectral_model = PowerLawSpectralModel(
    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

   name     value   error      unit         min        max    frozen
--------- --------- ----- -------------- ---------- --------- ------
    lon_0 2.000e-01   nan            deg        nan       nan  False
    lat_0 1.000e-01   nan            deg -9.000e+01 9.000e+01  False
    sigma 3.000e-01   nan            deg  0.000e+00       nan  False
        e 0.000e+00   nan                 0.000e+00 1.000e+00   True
      phi 0.000e+00   nan            deg        nan       nan   True
    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
[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]
)
[7]:
# Define some observation parameters
# We read in the pointing info from one of the 1dc event list files as an example
pointing = FixedPointingInfo.read(
    "$GAMMAPY_DATA/cta-1dc/data/baseline/gps/gps_baseline_110380.fits"
)

livetime = 1 * u.hour
offset_max = 2 * u.deg
offset = Angle("2 deg")
[8]:
exposure = make_map_exposure_true_energy(
    pointing=pointing.radec, livetime=livetime, aeff=irfs["aeff"], geom=geom
)
exposure.slice_by_idx({"energy": 3}).plot(add_cbar=True);
../_images/notebooks_simulate_3d_11_0.png
[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);
WARNING: Tried to get polar motions for times after IERS data is valid. Defaulting to polar motion from the 50-yr mean for those. This may affect precision at the 10s of arcsec level [astropy.coordinates.builtin_frames.utils]
../_images/notebooks_simulate_3d_12_1.png
[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
[11]:
energy = axis.edges
edisp = irfs["edisp"].to_energy_dispersion(
    offset, e_reco=energy, e_true=energy
)
edisp.plot_matrix();
/Users/deil/work/code/gammapy-docs/build/dev/gammapy/gammapy/utils/interpolation.py:159: Warning: Interpolated values reached float32 precision limit
  "Interpolated values reached float32 precision limit", Warning
../_images/notebooks_simulate_3d_14_1.png

Now we have to compute npred maps, i.e. “predicted counts per pixel” given the model and the observation infos: exposure, background, PSF and EDISP. For this we use the gammapy.cube.MapDataset object:

[12]:
background_model = BackgroundModel(background)
dataset = MapDataset(
    model=sky_model,
    exposure=exposure,
    background_model=background_model,
    psf=psf_kernel,
    edisp=edisp,
)
[13]:
npred = dataset.npred()
[14]:
npred.sum_over_axes().plot(add_cbar=True);
../_images/notebooks_simulate_3d_18_0.png
[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.data)
counts_map = WcsNDMap(geom, counts)
[16]:
counts_map.sum_over_axes().plot();
../_images/notebooks_simulate_3d_20_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, …

[17]:
# Define sky model to fit the data
spatial_model = GaussianSpatialModel(
    lon_0="0.1 deg", lat_0="0.1 deg", sigma="0.5 deg", frame="galactic"
)
spectral_model = PowerLawSpectralModel(
    index=2, amplitude="1e-11 cm-2 s-1 TeV-1", reference="1 TeV"
)
model = SkyModel(spatial_model=spatial_model, spectral_model=spectral_model)
print(model)
SkyModel

   name     value   error      unit         min        max    frozen
--------- --------- ----- -------------- ---------- --------- ------
    lon_0 1.000e-01   nan            deg        nan       nan  False
    lat_0 1.000e-01   nan            deg -9.000e+01 9.000e+01  False
    sigma 5.000e-01   nan            deg  0.000e+00       nan  False
        e 0.000e+00   nan                 0.000e+00 1.000e+00   True
      phi 0.000e+00   nan            deg        nan       nan   True
    index 2.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
[18]:
# We do not want to fit the background in this case, so we will freeze the parameters
background_model.parameters["norm"].value = 1.0
background_model.parameters["norm"].frozen = True
background_model.parameters["tilt"].frozen = True

print(background_model)
BackgroundModel

   name     value   error unit    min    max frozen
--------- --------- ----- ---- --------- --- ------
     norm 1.000e+00   nan      0.000e+00 nan   True
     tilt 0.000e+00   nan            nan nan   True
reference 1.000e+00   nan  TeV       nan nan   True
[19]:
dataset = MapDataset(
    model=model,
    exposure=exposure,
    counts=counts_map,
    background_model=background_model,
    psf=psf_kernel,
    edisp=edisp,
)
[20]:
%%time
fit = Fit([dataset])
result = fit.run(optimize_opts={"print_level": 1})
------------------------------------------------------------------
| FCN = 2.34E+05                |     Ncalls=214 (214 total)     |
| EDM = 6.27E-06 (Goal: 1E-05)  |            up = 1.0            |
------------------------------------------------------------------
|  Valid Min.   | Valid Param.  | Above EDM | Reached call limit |
------------------------------------------------------------------
|     True      |     True      |   False   |       False        |
------------------------------------------------------------------
| Hesse failed  |   Has cov.    | Accurate  | Pos. def. | Forced |
------------------------------------------------------------------
|     False     |     True      |   True    |   True    | False  |
------------------------------------------------------------------
CPU times: user 10.3 s, sys: 716 ms, total: 11.1 s
Wall time: 11.1 s

True model:

[21]:
print(sky_model)
SkyModel

   name     value   error      unit         min        max    frozen
--------- --------- ----- -------------- ---------- --------- ------
    lon_0 2.000e-01   nan            deg        nan       nan  False
    lat_0 1.000e-01   nan            deg -9.000e+01 9.000e+01  False
    sigma 3.000e-01   nan            deg  0.000e+00       nan  False
        e 0.000e+00   nan                 0.000e+00 1.000e+00   True
      phi 0.000e+00   nan            deg        nan       nan   True
    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:

[22]:
print(model)
SkyModel

   name     value   error      unit         min        max    frozen
--------- --------- ----- -------------- ---------- --------- ------
    lon_0 1.859e-01   nan            deg        nan       nan  False
    lat_0 8.768e-02   nan            deg -9.000e+01 9.000e+01  False
    sigma 2.963e-01   nan            deg  0.000e+00       nan  False
        e 0.000e+00   nan                 0.000e+00 1.000e+00   True
      phi 0.000e+00   nan            deg        nan       nan   True
    index 3.056e+00   nan                       nan       nan  False
amplitude 9.060e-12   nan cm-2 s-1 TeV-1        nan       nan  False
reference 1.000e+00   nan            TeV        nan       nan   True

To get the errors on the model, we can check the parameter table, or the covariance matrix

[23]:
result.parameters.to_table()
[23]:
Table length=11
namevalueerrorunitminmaxfrozen
str9float64float64str14float64float64bool
lon_01.859e-018.818e-03degnannanFalse
lat_08.768e-028.747e-03deg-9.000e+019.000e+01False
sigma2.963e-016.099e-03deg0.000e+00nanFalse
e0.000e+000.000e+000.000e+001.000e+00True
phi0.000e+000.000e+00degnannanTrue
index3.056e+003.055e-02nannanFalse
amplitude9.060e-124.611e-13cm-2 s-1 TeV-1nannanFalse
reference1.000e+000.000e+00TeVnannanTrue
norm1.000e+000.000e+000.000e+00nanTrue
tilt0.000e+000.000e+00nannanTrue
reference1.000e+000.000e+00TeVnannanTrue
[24]:
result.parameters.covariance[:3, :3]
[24]:
array([[ 7.77564445e-05,  3.64438368e-07, -8.77372121e-08],
       [ 3.64438368e-07,  7.65067462e-05,  1.83707774e-06],
       [-8.77372121e-08,  1.83707774e-06,  3.71937623e-05]])
[25]:
result.parameters.correlation[:3, :3]
/Users/deil/work/code/gammapy-docs/build/dev/gammapy/gammapy/modeling/parameter.py:511: RuntimeWarning: invalid value encountered in true_divide
  return self.covariance / np.outer(err, err)
[25]:
array([[ 1.        ,  0.00472505, -0.00163148],
       [ 0.00472505,  1.        ,  0.03443837],
       [-0.00163148,  0.03443837,  1.        ]])
[26]:
# Or, to see the value of and error on an individual parameter, say index:
print(result.parameters["index"].value, result.parameters.error("index"))
3.0555770657748837 0.030549883113329776
[27]:
# TODO: show e.g. how to make a residual image