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

Try online
You may download all the notebooks in the documentation as a tar file.
Source files: spectrum_simulation.ipynb | spectrum_simulation.py

1D spectrum simulation¶

Prerequisites¶

Knowledge of spectral extraction and datasets used in gammapy, see for instance the spectral analysis tutorial

Context¶

To simulate a specific observation, it is not always necessary to simulate the full photon list. For many uses cases, simulating directly a reduced binned dataset is enough: the IRFs reduced in the correct geometry are combined with a source model to predict an actual number of counts per bin. The latter is then used to simulate a reduced dataset using Poisson probability distribution.

This can be done to check the feasibility of a measurement, to test whether fitted parameters really provide a good fit to the data etc.

Here we will see how to perform a 1D spectral simulation of a CTA observation, in particular, we will generate OFF observations following the template background stored in the CTA IRFs.

Objective: simulate a number of spectral ON-OFF observations of a source with a power-law spectral model with CTA using the CTA 1DC response, fit them with the assumed spectral model and check that the distribution of fitted parameters is consistent with the input values.

Proposed approach:¶

We will use the following classes:

Setup¶

[1]:

%matplotlib inline
import matplotlib.pyplot as plt

[2]:

import numpy as np
import astropy.units as u
from astropy.coordinates import SkyCoord, Angle
from regions import CircleSkyRegion
from gammapy.datasets import SpectrumDatasetOnOff, SpectrumDataset, Datasets
from gammapy.makers import SpectrumDatasetMaker
from gammapy.modeling import Fit
from gammapy.modeling.models import (
    PowerLawSpectralModel,
    SkyModel,
)
from gammapy.irf import load_cta_irfs
from gammapy.data import Observation
from gammapy.maps import MapAxis, RegionGeom

Simulation of a single spectrum¶

To do a simulation, we need to define the observational parameters like the livetime, the offset, the assumed integration radius, the energy range to perform the simulation for and the choice of spectral model. We then use an in-memory observation which is convolved with the IRFs to get the predicted number of counts. This is Poission fluctuated using the fake() to get the simulated counts for each observation.

[3]:

# Define simulation parameters parameters
livetime = 1 * u.h

pointing = SkyCoord(0, 0, unit="deg", frame="galactic")
offset = 0.5 * u.deg

# Reconstructed and true energy axis
energy_axis = MapAxis.from_edges(
    np.logspace(-0.5, 1.0, 10), unit="TeV", name="energy", interp="log"
)
energy_axis_true = MapAxis.from_edges(
    np.logspace(-1.2, 2.0, 31), unit="TeV", name="energy_true", interp="log"
)

on_region_radius = Angle("0.11 deg")

center = pointing.directional_offset_by(
    position_angle=0 * u.deg, separation=offset
)
on_region = CircleSkyRegion(center=center, radius=on_region_radius)

[4]:

# Define spectral model - a simple Power Law in this case
model_simu = PowerLawSpectralModel(
    index=3.0,
    amplitude=2.5e-12 * u.Unit("cm-2 s-1 TeV-1"),
    reference=1 * u.TeV,
)
print(model_simu)
# we set the sky model used in the dataset
model = SkyModel(spectral_model=model_simu, name="source")

PowerLawSpectralModel

  type      name     value         unit        error   min max frozen link
-------- --------- ---------- -------------- --------- --- --- ------ ----
spectral     index 3.0000e+00                0.000e+00 nan nan  False
spectral amplitude 2.5000e-12 cm-2 s-1 TeV-1 0.000e+00 nan nan  False
spectral reference 1.0000e+00            TeV 0.000e+00 nan nan   True

[5]:

# Load the IRFs
# In this simulation, we use the CTA-1DC irfs shipped with gammapy.
irfs = load_cta_irfs(
    "$GAMMAPY_DATA/cta-1dc/caldb/data/cta/1dc/bcf/South_z20_50h/irf_file.fits"
)

Invalid unit found in background table! Assuming (s-1 MeV-1 sr-1)

[6]:

obs = Observation.create(pointing=pointing, livetime=livetime, irfs=irfs)
print(obs)

Observation

        obs id            : 0
        tstart            : 51544.00
        tstop             : 51544.04
        duration          : 3600.00 s
        pointing (icrs)   : 266.4 deg, -28.9 deg

        deadtime fraction : 0.0%

[7]:

# Make the SpectrumDataset
geom = RegionGeom.create(region=on_region, axes=[energy_axis])

dataset_empty = SpectrumDataset.create(
    geom=geom, energy_axis_true=energy_axis_true, name="obs-0"
)
maker = SpectrumDatasetMaker(selection=["exposure", "edisp", "background"])

dataset = maker.run(dataset_empty, obs)

[8]:

# Set the model on the dataset, and fake
dataset.models = model
dataset.fake(random_state=42)
print(dataset)

SpectrumDataset
---------------

  Name                            : obs-0

  Total counts                    : 298
  Total background counts         : 22.29
  Total excess counts             : 275.71

  Predicted counts                : 303.66
  Predicted background counts     : 22.29
  Predicted excess counts         : 281.37

  Exposure min                    : 2.53e+08 m2 s
  Exposure max                    : 1.77e+10 m2 s

  Number of total bins            : 9
  Number of fit bins              : 9

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

  Number of models                : 1
  Number of parameters            : 3
  Number of free parameters       : 2

  Component 0: SkyModel

    Name                      : source
    Datasets names            : None
    Spectral model type       : PowerLawSpectralModel
    Spatial  model type       :
    Temporal model type       :
    Parameters:
      index                   :      3.000   +/-    0.00
      amplitude               :   2.50e-12   +/- 0.0e+00 1 / (cm2 s TeV)
      reference    (frozen)   :      1.000       TeV

You can see that background counts are now simulated

On-Off analysis¶

To do an on off spectral analysis, which is the usual science case, the standard would be to use SpectrumDatasetOnOff, which uses the acceptance to fake off-counts

[9]:

dataset_on_off = SpectrumDatasetOnOff.from_spectrum_dataset(
    dataset=dataset, acceptance=1, acceptance_off=5
)
dataset_on_off.fake(npred_background=dataset.npred_background())
print(dataset_on_off)

SpectrumDatasetOnOff
--------------------

  Name                            : zDADetr1

  Total counts                    : 298
  Total background counts         : 19.00
  Total excess counts             : 279.00

  Predicted counts                : 300.39
  Predicted background counts     : 19.02
  Predicted excess counts         : 281.37

  Exposure min                    : 2.53e+08 m2 s
  Exposure max                    : 1.77e+10 m2 s

  Number of total bins            : 9
  Number of fit bins              : 9

  Fit statistic type              : wstat
  Fit statistic value (-2 log(L)) : 8.03

  Number of models                : 1
  Number of parameters            : 3
  Number of free parameters       : 2

  Component 0: SkyModel

    Name                      : source
    Datasets names            : None
    Spectral model type       : PowerLawSpectralModel
    Spatial  model type       :
    Temporal model type       :
    Parameters:
      index                   :      3.000   +/-    0.00
      amplitude               :   2.50e-12   +/- 0.0e+00 1 / (cm2 s TeV)
      reference    (frozen)   :      1.000       TeV

    Total counts_off                : 95
  Acceptance                      : 9
  Acceptance off                  : 45

You can see that off counts are now simulated as well. We now simulate several spectra using the same set of observation conditions.

[10]:

%%time

n_obs = 100
datasets = Datasets()

for idx in range(n_obs):
    dataset_on_off.fake(
        random_state=idx, npred_background=dataset.npred_background()
    )
    dataset_fake = dataset_on_off.copy(name=f"obs-{idx}")
    dataset_fake.meta_table["OBS_ID"] = [idx]
    datasets.append(dataset_fake)

CPU times: user 362 ms, sys: 11.4 ms, total: 374 ms
Wall time: 367 ms

[11]:

table = datasets.info_table()
table

[11]:

Table length=100

name	counts	background	excess	sqrt_ts	npred	npred_background	npred_signal	exposure_min	exposure_max	livetime	ontime	counts_rate	background_rate	excess_rate	n_bins	n_fit_bins	stat_type	stat_sum	counts_off	acceptance	acceptance_off	alpha
								m2 s	m2 s	s	s	1 / s	1 / s	1 / s
str7	int64	float64	float64	float64	float64	float64	float64	float64	float64	float64	float64	float64	float64	float64	int64	int64	str5	float64	int64	float64	float64	float64
stacked	317	18.399999618530273	298.6000061035156	27.082402135529488	68.16666751313541	68.16666751313541	nan	252718176.0	17719697408.0	3600.0	3600.000000000002	0.08805555555555555	0.005111111005147298	0.08294444613986546	9	9	wstat	738.7245834271496	92	9.0	45.000000932942285	0.20000000298023224
obs-1	275	22.0	253.0	23.76785365487285	64.16666666666669	64.16666666666669	nan	252718170.97287682	17719697919.599277	3600.0	3600.0	0.0763888888888889	0.006111111111111111	0.07027777777777777	9	9	wstat	575.779512784738	110	9.0	45.0	0.2
obs-2	293	20.6	272.4	25.17110555404655	66.0	66.0	nan	252718170.97287682	17719697919.599277	3600.0	3600.0	0.08138888888888889	0.005722222222222222	0.07566666666666666	9	9	wstat	645.4993824075303	103	9.0	45.0	0.2
obs-3	280	22.4	257.6	23.982951737405376	65.33333333333334	65.33333333333334	nan	252718170.97287682	17719697919.599277	3600.0	3600.0	0.07777777777777778	0.006222222222222222	0.07155555555555557	9	9	wstat	585.9241546985872	112	9.0	45.0	0.2
obs-4	337	20.6	316.4	27.682709945184747	73.33333333333334	73.33333333333334	nan	252718170.97287682	17719697919.599277	3600.0	3600.0	0.09361111111111112	0.005722222222222222	0.08788888888888888	9	9	wstat	787.314723949448	103	9.0	45.0	0.2
obs-5	283	24.400000000000002	258.6	23.727154782347895	67.50000000000001	67.50000000000001	nan	252718170.97287682	17719697919.599277	3600.0	3600.0	0.07861111111111112	0.006777777777777778	0.07183333333333335	9	9	wstat	574.8525737727499	122	9.0	44.99999999999999	0.2
obs-6	330	22.400000000000006	307.6	26.889184475727866	73.66666666666667	73.66666666666667	nan	252718170.97287682	17719697919.599277	3600.0	3600.0	0.09166666666666666	0.006222222222222224	0.08544444444444445	9	9	wstat	734.2414755285937	112	9.0	44.999999999999986	0.2
obs-7	283	25.8	257.2	23.43087974795853	68.66666666666667	68.66666666666667	nan	252718170.97287682	17719697919.599277	3600.0	3600.0	0.07861111111111112	0.007166666666666667	0.07144444444444444	9	9	wstat	552.3698779448044	129	9.0	45.0	0.2
obs-8	308	23.400000000000002	284.6	25.42049273328333	70.83333333333334	70.83333333333334	nan	252718170.97287682	17719697919.599277	3600.0	3600.0	0.08555555555555555	0.006500000000000001	0.07905555555555556	9	9	wstat	652.2621572567633	117	9.0	44.99999999999999	0.2
obs-9	299	20.4	278.6	25.57085071486863	66.83333333333334	66.83333333333334	nan	252718170.97287682	17719697919.599277	3600.0	3600.0	0.08305555555555555	0.005666666666666666	0.0773888888888889	9	9	wstat	666.9062670260786	102	9.0	45.0	0.2
obs-10	310	15.2	294.8	27.488972161356774	64.33333333333334	64.33333333333334	nan	252718170.97287682	17719697919.599277	3600.0	3600.0	0.08611111111111111	0.004222222222222222	0.08188888888888889	9	9	wstat	768.5404492529331	76	9.0	45.0	0.2
obs-11	285	24.000000000000004	261.0	23.933833745454685	67.50000000000001	67.50000000000001	nan	252718170.97287682	17719697919.599277	3600.0	3600.0	0.07916666666666666	0.006666666666666668	0.0725	9	9	wstat	583.448753292733	120	9.0	44.99999999999999	0.2
obs-12	299	21.000000000000004	278.0	25.4325263895117	67.33333333333336	67.33333333333336	nan	252718170.97287682	17719697919.599277	3600.0	3600.0	0.08305555555555555	0.0058333333333333345	0.07722222222222222	9	9	wstat	655.2768246480357	105	9.0	44.99999999999999	0.2
obs-13	309	21.6	287.4	25.877406235012618	69.50000000000001	69.50000000000001	nan	252718170.97287682	17719697919.599277	3600.0	3600.0	0.08583333333333333	0.006	0.07983333333333333	9	9	wstat	672.355748274336	108	9.0	45.0	0.2
obs-14	320	22.599999999999998	297.4	26.282826028192567	72.16666666666667	72.16666666666667	nan	252718170.97287682	17719697919.599277	3600.0	3600.0	0.08888888888888889	0.006277777777777777	0.08261111111111111	9	9	wstat	697.5949745415169	113	9.0	45.00000000000001	0.2
obs-15	283	22.0	261.0	24.25408979304137	65.50000000000001	65.50000000000001	nan	252718170.97287682	17719697919.599277	3600.0	3600.0	0.07861111111111112	0.006111111111111111	0.0725	9	9	wstat	592.8754563739889	110	9.0	45.0	0.2
obs-16	298	24.200000000000003	273.8	24.664218201697352	69.83333333333336	69.83333333333336	nan	252718170.97287682	17719697919.599277	3600.0	3600.0	0.08277777777777778	0.006722222222222223	0.07605555555555556	9	9	wstat	619.2584188386805	121	9.0	44.99999999999999	0.2
obs-17	301	28.200000000000003	272.8	24.01180151925227	73.66666666666667	73.66666666666667	nan	252718170.97287682	17719697919.599277	3600.0	3600.0	0.08361111111111111	0.007833333333333335	0.07577777777777778	9	9	wstat	588.4976379504056	141	9.0	44.99999999999999	0.2
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
obs-81	301	19.400000000000002	281.6	25.922347414834622	66.33333333333334	66.33333333333334	nan	252718170.97287682	17719697919.599277	3600.0	3600.0	0.08361111111111111	0.005388888888888889	0.07822222222222223	9	9	wstat	678.9813616402096	97	9.0	44.99999999999999	0.2
obs-82	280	22.000000000000004	258.0	24.072587044084237	65.0	65.0	nan	252718170.97287682	17719697919.599277	3600.0	3600.0	0.07777777777777778	0.006111111111111112	0.07166666666666667	9	9	wstat	581.8610083160877	110	9.0	44.99999999999999	0.2
obs-83	303	22.200000000000003	280.8	25.394928282299105	69.00000000000001	69.00000000000001	nan	252718170.97287682	17719697919.599277	3600.0	3600.0	0.08416666666666667	0.0061666666666666675	0.078	9	9	wstat	654.4650896470343	111	9.0	44.99999999999999	0.2
obs-84	269	21.200000000000003	247.8	23.580247140026987	62.500000000000014	62.500000000000014	nan	252718170.97287682	17719697919.599277	3600.0	3600.0	0.07472222222222222	0.00588888888888889	0.06883333333333333	9	9	wstat	577.2362891766269	106	9.0	44.99999999999999	0.2
obs-85	297	22.4	274.6	24.998935845735573	68.16666666666667	68.16666666666667	nan	252718170.97287682	17719697919.599277	3600.0	3600.0	0.0825	0.006222222222222222	0.07627777777777778	9	9	wstat	631.6298361198053	112	9.0	45.0	0.2
obs-86	286	17.800000000000004	268.2	25.423408969515453	62.50000000000001	62.50000000000001	nan	252718170.97287682	17719697919.599277	3600.0	3600.0	0.07944444444444444	0.004944444444444446	0.0745	9	9	wstat	655.2082535886425	89	9.0	44.999999999999986	0.2
obs-87	333	19.8	313.2	27.646851451147217	72.00000000000001	72.00000000000001	nan	252718170.97287682	17719697919.599277	3600.0	3600.0	0.0925	0.0055000000000000005	0.087	9	9	wstat	770.1313412321031	99	9.0	45.0	0.2
obs-88	315	18.200000000000003	296.8	27.0172059251097	67.66666666666669	67.66666666666669	nan	252718170.97287682	17719697919.599277	3600.0	3600.0	0.0875	0.005055555555555556	0.08244444444444445	9	9	wstat	734.6476098262428	91	9.0	44.99999999999999	0.2
obs-89	287	22.000000000000004	265.0	24.49456558680515	66.16666666666667	66.16666666666667	nan	252718170.97287682	17719697919.599277	3600.0	3600.0	0.07972222222222222	0.006111111111111112	0.07361111111111111	9	9	wstat	611.2334063685089	110	9.0	44.99999999999999	0.2
obs-90	286	19.000000000000004	267.0	25.131221887043438	63.50000000000001	63.50000000000001	nan	252718170.97287682	17719697919.599277	3600.0	3600.0	0.07944444444444444	0.005277777777777779	0.07416666666666667	9	9	wstat	635.8304116110406	95	9.0	44.99999999999999	0.2
obs-91	285	25.200000000000003	259.8	23.67754591931069	68.50000000000001	68.50000000000001	nan	252718170.97287682	17719697919.599277	3600.0	3600.0	0.07916666666666666	0.007000000000000001	0.07216666666666667	9	9	wstat	573.3564742963381	126	9.0	44.99999999999999	0.2
obs-92	313	23.6	289.4	25.664935420194176	71.83333333333336	71.83333333333336	nan	252718170.97287682	17719697919.599277	3600.0	3600.0	0.08694444444444445	0.006555555555555556	0.08038888888888888	9	9	wstat	669.4496142825196	118	9.0	45.0	0.2
obs-93	302	18.8	283.2	26.123867522605497	66.00000000000001	66.00000000000001	nan	252718170.97287682	17719697919.599277	3600.0	3600.0	0.08388888888888889	0.005222222222222223	0.07866666666666666	9	9	wstat	700.5066738339824	94	9.0	45.0	0.2
obs-94	322	21.8	300.2	26.574029757473202	71.83333333333334	71.83333333333334	nan	252718170.97287682	17719697919.599277	3600.0	3600.0	0.08944444444444444	0.006055555555555555	0.08338888888888889	9	9	wstat	715.2570498037147	109	9.0	45.0	0.2
obs-95	305	24.400000000000006	280.6	25.030731260493152	71.16666666666667	71.16666666666667	nan	252718170.97287682	17719697919.599277	3600.0	3600.0	0.08472222222222223	0.006777777777777779	0.07794444444444446	9	9	wstat	643.2846879404793	122	9.0	44.99999999999999	0.2
obs-96	301	23.6	277.4	24.969845969421428	69.83333333333334	69.83333333333334	nan	252718170.97287682	17719697919.599277	3600.0	3600.0	0.08361111111111111	0.006555555555555556	0.07705555555555554	9	9	wstat	626.7396842428624	118	9.0	45.0	0.2
obs-97	290	18.8	271.2	25.417982194454826	64.00000000000001	64.00000000000001	nan	252718170.97287682	17719697919.599277	3600.0	3600.0	0.08055555555555556	0.005222222222222223	0.07533333333333334	9	9	wstat	650.9762484493856	94	9.0	45.0	0.2
obs-98	301	20.400000000000002	280.6	25.687832964675007	67.16666666666667	67.16666666666667	nan	252718170.97287682	17719697919.599277	3600.0	3600.0	0.08361111111111111	0.005666666666666667	0.07794444444444446	9	9	wstat	664.5620099404952	102	9.0	44.99999999999999	0.2
obs-99	323	20.8	302.2	26.85707842376871	71.16666666666669	71.16666666666669	nan	252718170.97287682	17719697919.599277	3600.0	3600.0	0.08972222222222222	0.005777777777777778	0.08394444444444445	9	9	wstat	733.1576598768017	104	9.0	45.0	0.2

Before moving on to the fit let’s have a look at the simulated observations.

[12]:

fix, axes = plt.subplots(1, 3, figsize=(12, 4))
axes[0].hist(table["counts"])
axes[0].set_xlabel("Counts")
axes[1].hist(table["counts_off"])
axes[1].set_xlabel("Counts Off")
axes[2].hist(table["excess"])
axes[2].set_xlabel("excess");

../../../_images/tutorials_analysis_1D_spectrum_simulation_19_0.png

Now, we fit each simulated spectrum individually

[13]:

%%time
results = []

fit = Fit()

for dataset in datasets:
    dataset.models = model.copy()
    result = fit.optimize(dataset)
    results.append(
        {
            "index": result.parameters["index"].value,
            "amplitude": result.parameters["amplitude"].value,
        }
    )

CPU times: user 6.79 s, sys: 78.8 ms, total: 6.87 s
Wall time: 6.81 s

We take a look at the distribution of the fitted indices. This matches very well with the spectrum that we initially injected.

[14]:

index = np.array([_["index"] for _ in results])
plt.hist(index, bins=10, alpha=0.5)
plt.axvline(x=model_simu.parameters["index"].value, color="red")
print(f"index: {index.mean()} += {index.std()}")

index: 3.0036925550380533 += 0.08081469527527077

../../../_images/tutorials_analysis_1D_spectrum_simulation_23_1.png

Exercises¶

Change the observation time to something longer or shorter. Does the observation and spectrum results change as you expected?
Change the spectral model, e.g. add a cutoff at 5 TeV, or put a steep-spectrum source with spectral index of 4.0
Simulate spectra with the spectral model we just defined. How much observation duration do you need to get back the injected parameters?

[ ]: