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

3D analysis

This tutorial shows how to run a 3D map-based analysis (two spatial and one energy axis).

The example data is three observations of the Galactic center region with CTA.

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
from gammapy.extern.pathlib import Path
from gammapy.data import DataStore
from gammapy.irf import EnergyDispersion
from gammapy.maps import WcsGeom, MapAxis, Map
from gammapy.cube import MapMaker, PSFKernel, MapFit
from gammapy.cube.models import SkyModel
from gammapy.spectrum.models import PowerLaw
from gammapy.image.models import SkyGaussian, SkyPointSource
from regions import CircleSkyRegion
In [3]:
!gammapy info --no-envvar --no-dependencies --no-system

Gammapy package:

        path                   : /Users/deil/work/code/gammapy/gammapy
        version                : 0.8.dev7377
        githash                : 148b1c4221a90bff80e0601c91260bcc80682152

Make maps

In [4]:
# Define which data to use
data_store = DataStore.from_dir("$GAMMAPY_DATA/cta-1dc/index/gps/")
obs_ids = [110380, 111140, 111159]
obs_list = data_store.obs_list(obs_ids)
In [5]:
# Define map geometry (spatial and energy binning)
axis = MapAxis.from_edges(
    np.logspace(-1., 1., 10), unit="TeV", name="energy", interp="log"
)
geom = WcsGeom.create(
    skydir=(0, 0),
    binsz=0.02,
    width=(10, 8),
    coordsys="GAL",
    proj="CAR",
    axes=[axis],
)
In [6]:
%%time
maker = MapMaker(geom, 4. * u.deg)
maps = maker.run(obs_list)
CPU times: user 11.4 s, sys: 2.38 s, total: 13.8 s
Wall time: 13.8 s
In [7]:
%%time
# Make 2D images (for plotting, analysis will be 3D)
images = maker.make_images()
CPU times: user 16.1 ms, sys: 8.88 ms, total: 25 ms
Wall time: 23 ms
In [8]:
images["counts"].smooth(width=0.1 * u.deg).plot(stretch="sqrt");
../_images/notebooks_analysis_3d_11_0.png
In [9]:
images["background"].plot(stretch="sqrt");
../_images/notebooks_analysis_3d_12_0.png
In [10]:
residual = images["counts"].copy()
residual.data -= images["background"].data
residual.smooth(5).plot(stretch="sqrt");
../_images/notebooks_analysis_3d_13_0.png

Compute PSF kernel

For the moment we rely on the ObservationList.make_mean_psf() method.

In [11]:
%%time
obs_list = data_store.obs_list(obs_ids)
src_pos = SkyCoord(0, 0, unit="deg", frame="galactic")

table_psf = obs_list.make_mean_psf(src_pos)
psf_kernel = PSFKernel.from_table_psf(
    table_psf, maps["exposure"].geom, max_radius="0.3 deg"
)
CPU times: user 292 ms, sys: 5.62 ms, total: 297 ms
Wall time: 295 ms

Compute energy dispersion

In [12]:
%%time
energy_axis = geom.get_axis_by_name("energy")
energy = energy_axis.edges * energy_axis.unit
edisp = obs_list.make_mean_edisp(
    position=src_pos, e_true=energy, e_reco=energy
)
CPU times: user 219 ms, sys: 14.3 ms, total: 233 ms
Wall time: 231 ms

Save maps

It’s common to run the “precompute” step and the “likelihood fit” step separately, because often the “precompute” of maps, PSF and EDISP is slow if you have a lot of data.

Here it woudn’t really be necessary, because the precompute step (everything above this cell) takes less than a minute.

But usually you would do it like this: write precomputed things to FITS files, and then read them from your script that does the likelihood fitting without having to run the precomputations again.

In [13]:
# Write
path = Path("analysis_3d")
path.mkdir(exist_ok=True)
maps["counts"].write(str(path / "counts.fits"), overwrite=True)
maps["background"].write(str(path / "background.fits"), overwrite=True)
maps["exposure"].write(str(path / "exposure.fits"), overwrite=True)
psf_kernel.write(str(path / "psf.fits"), overwrite=True)
edisp.write(str(path / "edisp.fits"), overwrite=True)
In [14]:
# Read
maps = {
    "counts": Map.read(str(path / "counts.fits")),
    "background": Map.read(str(path / "background.fits")),
    "exposure": Map.read(str(path / "exposure.fits")),
}
psf_kernel = PSFKernel.read(str(path / "psf.fits"))
edisp = EnergyDispersion.read(str(path / "edisp.fits"))

Cutout

Let’s cut out only part of the map, so that we can have a faster fit

In [15]:
cmaps = {
    name: m.cutout(SkyCoord(0, 0, unit="deg", frame="galactic"), 1.5 * u.deg)
    for name, m in maps.items()
}
cmaps["counts"].sum_over_axes().plot(stretch="sqrt");
../_images/notebooks_analysis_3d_22_0.png

Fit mask

To select a certain region and/or energy range for the fit we can create a fit mask.

In [16]:
mask = Map.from_geom(cmaps["counts"].geom)

region = CircleSkyRegion(center=src_pos, radius=0.6 * u.deg)
mask.data = mask.geom.region_mask([region])

mask.get_image_by_idx((0,)).plot();
../_images/notebooks_analysis_3d_24_0.png

In addition we also exclude the range below 0.3 TeV for the fit:

In [17]:
coords = mask.geom.get_coord()
mask.data &= coords["energy"] > 0.3

Model fit

  • TODO: Add diffuse emission model? (it’s 800 MB, maybe prepare a cutout)
  • TODO: compare against true model known for DC1
In [18]:
spatial_model = SkyPointSource(lon_0="0.01 deg", lat_0="0.01 deg")
spectral_model = PowerLaw(
    index=2.2, amplitude="3e-12 cm-2 s-1 TeV-1", reference="1 TeV"
)
model = SkyModel(spatial_model=spatial_model, spectral_model=spectral_model)
In [19]:
%%time
fit = MapFit(
    model=model,
    counts=cmaps["counts"],
    exposure=cmaps["exposure"],
    background=cmaps["background"],
    mask=mask,
    psf=psf_kernel,
    edisp=edisp,
)

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

FCN = 14281.494005819377 TOTAL NCALL = 162 NCALLS = 162
EDM = 7.274632570778809e-05 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 -4.82862 0.219073 No
1 par_001_lat_0 -5.2709 0.217932 No
2 par_002_index 2.37661 0.0602355 No
3 par_003_amplitude 0.27121 0.0149826 No
4 par_004_reference 1 1 0 Yes

CPU times: user 2.58 s, sys: 71 ms, total: 2.65 s
Wall time: 2.59 s

Check model fit

  • plot counts spectrum for some on region (e.g. the one used in 1D spec analysis, 0.2 deg)
  • plot residual image for some energy band (e.g. the total one used here)
In [20]:
# Parameter error are not synched back to
# sub model components automatically yet
spec = result.model.spectral_model
print(spec)
PowerLaw

Parameters:

           name     value   error       unit         min    max
        --------- --------- ----- --------------- --------- ---
            index 2.377e+00   nan                       nan nan
        amplitude 2.712e-12   nan 1 / (cm2 s TeV)       nan nan
        reference 1.000e+00   nan             TeV 0.000e+00 nan
In [21]:
# For now, we can copy the parameter error manually
spec.parameters.set_parameter_errors(
    {
        "index": result.model.parameters.error("index"),
        "amplitude": result.model.parameters.error("amplitude"),
    }
)
print(spec)
PowerLaw

Parameters:

           name     value     error         unit         min    max
        --------- --------- --------- --------------- --------- ---
            index 2.377e+00 6.024e-02                       nan nan
        amplitude 2.712e-12 1.498e-13 1 / (cm2 s TeV)       nan nan
        reference 1.000e+00 0.000e+00             TeV 0.000e+00 nan

Covariance:

           name           index               amplitude       reference
        --------- --------------------- --------------------- ---------
            index 0.0036283176318476075                   0.0       0.0
        amplitude                   0.0 2.244783432817743e-26       0.0
        reference                   0.0                   0.0       0.0
In [22]:
energy_range = [0.1, 10] * u.TeV
spec.plot(energy_range, energy_power=2)
spec.plot_error(energy_range, energy_power=2);
../_images/notebooks_analysis_3d_33_0.png
In [23]:
result.model.parameters.covariance = 10

Exercises

  • Analyse the second source in the field of view: G0.9+0.1
  • Run the model fit with energy dispersion (pass edisp to MapFit)