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 SkyCoord
from gammapy.irf import load_cta_irfs
from gammapy.maps import WcsGeom, MapAxis
from gammapy.modeling.models import PowerLawSpectralModel
from gammapy.modeling.models import GaussianSpatialModel
from gammapy.modeling.models import SkyModel
from gammapy.cube import MapDataset, MapDatasetMaker, SafeMaskMaker
from gammapy.modeling import Fit
from gammapy.data import Observation
[3]:
!gammapy info --no-envvar --no-dependencies --no-system

Gammapy package:

        version                : 0.16.dev144+gb770832e4
        path                   : /Users/adonath/github/adonath/gammapy/gammapy

Simulation

We will simulate using the CTA-1DC IRFs shipped with gammapy. Note that for dedictaed CTA simulations, you can simply use `Observation.from_caldb() <>`__ without having to externally load the IRFs

[4]:
# Loading IRFs
irfs = load_cta_irfs(
    "$GAMMAPY_DATA/cta-1dc/caldb/data/cta/1dc/bcf/South_z20_50h/irf_file.fits"
)
[5]:
# Define the observation parameters (typically the observation duration and the pointing position):
livetime = 2.0 * u.hr
pointing = SkyCoord(0, 0, unit="deg", frame="galactic")
[6]:
# Define map geometry for binned simulation
energy_reco = 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=(6, 6),
    frame="galactic",
    axes=[energy_reco],
)
# It is usually useful to have a separate binning for the true energy axis
energy_true = MapAxis.from_edges(
    np.logspace(-1.5, 1.5, 30), unit="TeV", name="energy", interp="log"
)
[7]:
# Define sky model to used simulate the data.
# Here we use a Gaussian spatial model and a Power Law spectral model.
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"
)
model_simu = SkyModel(
    spatial_model=spatial_model,
    spectral_model=spectral_model,
    name="model_simu",
)
print(model_simu)
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

Now, comes the main part of dataset simulation. We create an in-memory observation and an empty dataset. We then predict the number of counts for the given model, and Poission fluctuate it using fake() to make a simulated counts maps. Keep in mind that it is important to specify the selection of the maps that you want to produce

[8]:
# Create an in-memory observation
obs = Observation.create(pointing=pointing, livetime=livetime, irfs=irfs)
print(obs)
# Make the MapDataset
empty = MapDataset.create(geom)
maker = MapDatasetMaker(selection=["exposure", "background", "psf", "edisp"])
maker_safe_mask = SafeMaskMaker(methods=["offset-max"], offset_max=4.0 * u.deg)
dataset = maker.run(empty, obs)
dataset = maker_safe_mask.run(dataset, obs)
print(dataset)
Info for OBS_ID = 0
- Pointing pos: RA 266.40 deg / Dec -28.94 deg
- Livetime duration: 7200.0 s

WARNING: AstropyDeprecationWarning: The truth value of a Quantity is ambiguous. In the future this will raise a ValueError. [astropy.units.quantity]
MapDataset

    Name                            : obs_0

    Total counts                    : nan
    Total predicted counts          : 161422.07
    Total background counts         : 161422.07

    Exposure min                    : 6.41e+07 m2 s
    Exposure max                    : 2.53e+10 m2 s

    Number of total bins            : 0
    Number of fit bins              : 804492

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

    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


[9]:
# Add the model on the dataset and Poission fluctuate
dataset.models = model_simu
dataset.fake()
# Do a print on the dataset - there is now a counts maps
print(dataset)
MapDataset

    Name                            : obs_0

    Total counts                    : 169811
    Total predicted counts          : 169652.54
    Total background counts         : 161422.07

    Exposure min                    : 6.41e+07 m2 s
    Exposure max                    : 2.53e+10 m2 s

    Number of total bins            : 810000
    Number of fit bins              : 804492

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

    Number of models                : 2
    Number of parameters            : 11
    Number of free parameters       : 6

    Component 0:
        Name                        : model_simu
        Type                        : SkyModel
        Spatial  model type         : GaussianSpatialModel
        Spectral model type         : PowerLawSpectralModel
        Parameters:
            lon_0                   : 0.200  deg
            lat_0                   : 0.100  deg
            sigma                   : 0.300  deg
            e            (frozen)   : 0.000
            phi          (frozen)   : 0.000  deg
            index                   : 3.000
            amplitude               : 1.00e-11  1 / (cm2 s TeV)
            reference    (frozen)   : 1.000  TeV

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


Now use this dataset as you would in all standard analysis. You can plot the maps, or proceed with your custom analysis. In the next section, we show the standard 3D fitting as in analysis_3d.

[10]:
# To plot, eg, counts:
dataset.counts.smooth(0.05 * u.deg).plot_interactive(
    add_cbar=True, stretch="linear"
)
../_images/notebooks_simulate_3d_16_0.png

Fit

In this section, we do a usual 3D fit with the same model used to simulated the data and see the stability of the simulations. Often, it is useful to simulate many such datasets and look at the distribution of the reconstructed parameters.

[11]:
# Make a copy of the dataset
dataset1 = dataset.copy()
[12]:
# Define sky model to fit the data
spatial_model1 = GaussianSpatialModel(
    lon_0="0.1 deg", lat_0="0.1 deg", sigma="0.5 deg", frame="galactic"
)
spectral_model1 = PowerLawSpectralModel(
    index=2, amplitude="1e-11 cm-2 s-1 TeV-1", reference="1 TeV"
)
model_fit = SkyModel(
    spatial_model=spatial_model1,
    spectral_model=spectral_model1,
    name="model_fit",
)

dataset1.models = model_fit
print(model_fit)
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
[13]:
# We do not want to fit the background in this case, so we will freeze the parameters
background_model = dataset1.background_model
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
[14]:
%%time
fit = Fit([dataset1])
result = fit.run(optimize_opts={"print_level": 1})
------------------------------------------------------------------
| FCN = 5.637E+05               |     Ncalls=222 (222 total)     |
| EDM = 6.18E-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 21.1 s, sys: 1.32 s, total: 22.4 s
Wall time: 22.6 s
[15]:
dataset1.plot_residuals(method="diff/sqrt(model)", vmin=-0.5, vmax=0.5)
[15]:
(<matplotlib.axes._subplots.WCSAxesSubplot at 0x1192b54e0>, None)
../_images/notebooks_simulate_3d_22_1.png

Compare the injected and fitted models:

[16]:
print("True model: \n", model_simu, "\n\n Fitted model: \n", model_fit)
True 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

 Fitted model:
 SkyModel

   name     value   error      unit         min        max    frozen
--------- --------- ----- -------------- ---------- --------- ------
    lon_0 2.095e-01   nan            deg        nan       nan  False
    lat_0 9.582e-02   nan            deg -9.000e+01 9.000e+01  False
    sigma 3.041e-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.040e+00   nan                       nan       nan  False
amplitude 9.637e-12   nan cm-2 s-1 TeV-1        nan       nan  False
reference 1.000e+00   nan            TeV        nan       nan   True

Get the errors on the fitted parameters from the parameter table

[17]:
result.parameters.to_table()
[17]:
Table length=11
namevalueerrorunitminmaxfrozen
str9float64float64str14float64float64bool
lon_02.095e-016.166e-03degnannanFalse
lat_09.582e-026.100e-03deg-9.000e+019.000e+01False
sigma3.041e-014.167e-03deg0.000e+00nanFalse
e0.000e+000.000e+000.000e+001.000e+00True
phi0.000e+000.000e+00degnannanTrue
index3.040e+002.103e-02nannanFalse
amplitude9.637e-123.330e-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