This is a fixed-text formatted version of a Jupyter notebook
Try online
You can contribute with your own notebooks in this GitHub repository.
Source files: simulate_3d.ipynb | simulate_3d.py
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.15
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), coordsys="GAL", 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 = 1
- 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_1
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_1
Total counts : 169096
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)) : 562017.57
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"
)
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.62E+05 | Ncalls=231 (231 total) |
| EDM = 3.7E-05 (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 11.4 s, sys: 2.58 s, total: 14 s
Wall time: 14.1 s
[15]:
dataset1.plot_residuals(method="diff/sqrt(model)", vmin=-0.5, vmax=0.5)
[15]:
(<matplotlib.axes._subplots.WCSAxesSubplot at 0x11a4bdef0>, None)
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 1.951e-01 nan deg nan nan False
lat_0 1.014e-01 nan deg -9.000e+01 9.000e+01 False
sigma 2.957e-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.998e+00 nan nan nan False
amplitude 1.012e-11 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
| name | value | error | unit | min | max | frozen |
|---|---|---|---|---|---|---|
| str9 | float64 | float64 | str14 | float64 | float64 | bool |
| lon_0 | 1.951e-01 | 6.024e-03 | deg | nan | nan | False |
| lat_0 | 1.014e-01 | 5.848e-03 | deg | -9.000e+01 | 9.000e+01 | False |
| sigma | 2.957e-01 | 4.150e-03 | deg | 0.000e+00 | nan | False |
| e | 0.000e+00 | 0.000e+00 | 0.000e+00 | 1.000e+00 | True | |
| phi | 0.000e+00 | 0.000e+00 | deg | nan | nan | True |
| index | 2.998e+00 | 2.001e-02 | nan | nan | False | |
| amplitude | 1.012e-11 | 3.315e-13 | cm-2 s-1 TeV-1 | nan | nan | False |
| reference | 1.000e+00 | 0.000e+00 | TeV | nan | nan | True |
| norm | 1.000e+00 | 0.000e+00 | 0.000e+00 | nan | True | |
| tilt | 0.000e+00 | 0.000e+00 | nan | nan | True | |
| reference | 1.000e+00 | 0.000e+00 | TeV | nan | nan | True |