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");

In [9]:

images["background"].plot(stretch="sqrt");

In [10]:

residual = images["counts"].copy()
residual.data -= images["background"].data
residual.smooth(5).plot(stretch="sqrt");


## 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(
)

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 = {
}


## 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");


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)



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

In [17]:

coords = mask.geom.get_coord()


## 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"],
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);

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)