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

• Try online

• You may download all the notebooks as a tar file.

• Source files: light_curve_simulation.ipynb | light_curve_simulation.py

Simulating and fitting a time varying source

Prerequisites

• To understand how a single binned simulation works, please refer to spectrum_simulation simulate_3d for 1D and 3D simulations respectively.

• For details of light curve extraction using gammapy, refer to the two tutorials light_curve and light_curve_flare

Context

Frequently, studies of variable sources (eg: decaying GRB light curves, AGN flares, etc) require time variable simulations. For most use cases, generating an event list is an overkill, and it suffices to use binned simulations using a temporal model.

Objective: Simulate and fit a time decaying light curve of a source with CTA using the CTA 1DC response

Proposed approach

We will simulate 10 spectral datasets within given time intervals (Good Time Intervals) following a given spectral (a power law) and temporal profile (an exponential decay, with a decay time of 6 hr ). These are then analysed using the light curve estimator to obtain flux points.

Modelling and fitting of lightcurves can be done either - directly on the output of the LighCurveEstimator (at the DL5 level) - fit the simulated datasets (at the DL4 level)

In summary, necessary steps are:

• Choose observation parameters including a list of gammapy.data.GTI

• Define temporal and spectral models from :ref:model-gallery as per science case

• Perform the simulation (in 1D or 3D)

• Extract the light curve from the reduced dataset as shown in light curve notebook

• Optionally, we show here how to fit the simulated datasets using a source model

Setup

As usual, we’ll start with some general imports…

Setup

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import astropy.units as u
from astropy.coordinates import SkyCoord
from astropy.time import Time

import logging

log = logging.getLogger(__name__)

And some gammapy specific imports

from gammapy.data import Observation
from gammapy.irf import load_cta_irfs
from gammapy.datasets import SpectrumDataset, Datasets, FluxPointsDataset
from gammapy.modeling.models import (
    PowerLawSpectralModel,
    ExpDecayTemporalModel,
    SkyModel,
)
from gammapy.maps import MapAxis, RegionGeom, TimeMapAxis
from gammapy.estimators import LightCurveEstimator
from gammapy.makers import SpectrumDatasetMaker
from gammapy.modeling import Fit
from gammapy.data import observatory_locations

We first define our preferred time format:

TimeMapAxis.time_format = "iso"

Simulating a light curve

We will simulate 10 spectra between 300 GeV and 10 TeV using an PowerLawSpectralModel and a ExpDecayTemporalModel. The important thing to note here is how to attach a different GTI to each dataset. Since we use spectrum datasets here, we will use a RegionGeom.

# Loading IRFs
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)

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

geom = RegionGeom.create("galactic;circle(0, 0, 0.11)", axes=[energy_axis])

# Pointing position
pointing = SkyCoord(0.5, 0.5, unit="deg", frame="galactic")

Note that observations are usually conducted in Wobble mode, in which the source is not in the center of the camera. This allows to have a symmetrical sky position from which background can be estimated.

# Define the source model: A combination of spectral and temporal model

gti_t0 = Time("2020-03-01")
spectral_model = PowerLawSpectralModel(
    index=3, amplitude="1e-11 cm-2 s-1 TeV-1", reference="1 TeV"
)
temporal_model = ExpDecayTemporalModel(t0="6 h", t_ref=gti_t0.mjd * u.d)

model_simu = SkyModel(
    spectral_model=spectral_model,
    temporal_model=temporal_model,
    name="model-simu",
)

/Users/adonath/software/mambaforge/envs/gammapy-dev/lib/python3.9/site-packages/astropy/units/quantity.py:613: RuntimeWarning: overflow encountered in exp
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)
/Users/adonath/software/mambaforge/envs/gammapy-dev/lib/python3.9/site-packages/astropy/units/quantity.py:613: RuntimeWarning: invalid value encountered in subtract
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)

# Look at the model
model_simu.parameters.to_table()

Table length=5

type	name	value	unit	error	min	max	frozen	is_norm	link
str8	str9	float64	str14	int64	float64	float64	bool	bool	str1
spectral	index	3.0000e+00		0.000e+00	nan	nan	False	False
spectral	amplitude	1.0000e-11	cm-2 s-1 TeV-1	0.000e+00	nan	nan	False	True
spectral	reference	1.0000e+00	TeV	0.000e+00	nan	nan	True	False
temporal	t0	6.0000e+00	h	0.000e+00	nan	nan	False	False
temporal	t_ref	5.8909e+04	d	0.000e+00	nan	nan	True	False

Now, define the start and observation livetime wrt to the reference time, gti_t0

n_obs = 10

tstart = gti_t0 + [1, 2, 3, 5, 8, 10, 20, 22, 23, 24] * u.h
lvtm = [55, 25, 26, 40, 40, 50, 40, 52, 43, 47] * u.min

Now perform the simulations

datasets = Datasets()

empty = SpectrumDataset.create(
    geom=geom, energy_axis_true=energy_axis_true, name="empty"
)

maker = SpectrumDatasetMaker(selection=["exposure", "background", "edisp"])


for idx in range(n_obs):
    obs = Observation.create(
        pointing=pointing,
        livetime=lvtm[idx],
        tstart=tstart[idx],
        irfs=irfs,
        reference_time=gti_t0,
        obs_id=idx,
        location=observatory_locations["cta_south"],
    )
    empty_i = empty.copy(name=f"dataset-{idx}")
    dataset = maker.run(empty_i, obs)
    dataset.models = model_simu
    dataset.fake()
    datasets.append(dataset)

The reduced datasets have been successfully simulated. Let’s take a quick look into our datasets.

datasets.info_table()

Table length=10

name	counts	excess	sqrt_ts	background	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
								m2 s	m2 s	s	s	1 / s	1 / s	1 / s
str9	int64	float64	float64	float64	float64	float64	float64	float64	float64	float64	float64	float64	float64	float64	int64	int64	str4	float64
dataset-0	788	767.6812133789062	65.03515595380219	20.318769454956055	20.318769738078117	20.318769454956055	nan	216137904.0	16025275392.0	3299.9999999999973	3300.0	0.238787878787879	0.006157202865138204	0.23263067072088087	9	9	cash	nan
dataset-1	303	293.764195792777	39.08711484570209	9.235804207223024	332.14961280980896	9.235804207223024	322.9138086025859	98244500.93556362	7284216297.312286	1500.0000000000036	1500.0000000000036	0.2019999999999995	0.006157202804815335	0.1958427971951842	9	9	cash	-1804.171041139897
dataset-2	272	262.39476362448806	35.97319869287688	9.605236375511943	293.48960258443094	9.605236375511943	283.88436620891895	102174280.97298616	7575584949.204778	1560.0000000000036	1560.0000000000036	0.17435897435897396	0.006157202804815334	0.16820177155415864	9	9	cash	-1574.3589928828487
dataset-3	335	320.2227132684432	38.08741057470171	14.777286731556805	321.78385941877633	14.777286731556805	307.0065726872195	157191201.49690142	11654746075.69963	2400.0	2400.0	0.13958333333333334	0.0061572028048153355	0.133426130528518	9	9	cash	-2150.529986042023
dataset-4	182	167.2227132684432	24.073382946346577	14.777286731556805	200.9861857996506	14.777286731556805	186.2088990680938	157191201.49690142	11654746075.69963	2400.0	2400.0	0.07583333333333334	0.0061572028048153355	0.069676130528518	9	9	cash	-875.6042880128907
dataset-5	154	135.528391585554	19.547999322228286	18.471608414446006	182.99944272932325	18.471608414446006	164.52783431487723	196489001.87112677	14568432594.624538	3000.0	3000.0	0.051333333333333335	0.0061572028048153355	0.045176130528518	9	9	cash	-708.2822693043212
dataset-6	58	43.222713268443194	8.495137984520849	14.777286731556805	39.97792082811522	14.777286731556805	25.20063409655841	157191201.49690142	11654746075.69963	2400.0	2400.0	0.024166666666666666	0.0061572028048153355	0.01800946386185133	9	9	cash	-147.1795589597059
dataset-7	45	25.789527248976157	5.002955809473021	19.210472751023843	42.30482881204804	19.210472751023843	23.094356061024197	204348561.94597185	15151169898.40952	3120.0	3120.0	0.014423076923076924	0.006157202804815335	0.00826587411826159	9	9	cash	-126.82594551244152
dataset-8	23	7.114416763576436	1.6718101499279763	15.885583236423564	32.24990050264766	15.885583236423564	16.36431726622409	168980541.609169	12528852031.377102	2580.0	2580.0	0.008914728682170543	0.006157202804815335	0.002757525877355208	9	9	cash	-24.280591740673064
dataset-9	28	10.636688090420755	2.3422096643225636	17.363311909579245	32.42183519708895	17.363311909579245	15.05852328750971	184699661.75885916	13694326638.947065	2820.0	2820.0	0.009929078014184398	0.0061572028048153355	0.003771875209369062	9	9	cash	-28.06690147791762

Extract the lightcurve

This section uses standard light curve estimation tools for a 1D extraction. Only a spectral model needs to be defined in this case. Since the estimator returns the integrated flux separately for each time bin, the temporal model need not be accounted for at this stage. We extract the lightcurve in 3 energy binsç

# Define the model:
spectral_model = PowerLawSpectralModel(
    index=3, amplitude="1e-11 cm-2 s-1 TeV-1", reference="1 TeV"
)
model_fit = SkyModel(spectral_model=spectral_model, name="model-fit")

# Attach model to all datasets
datasets.models = model_fit

%%time
lc_maker_1d = LightCurveEstimator(
    energy_edges=[0.3, 0.6, 1.0, 10] * u.TeV,
    source="model-fit",
    selection_optional=["ul"],
)
lc_1d = lc_maker_1d.run(datasets)

CPU times: user 10.7 s, sys: 203 ms, total: 10.9 s
Wall time: 14.1 s

ax = lc_1d.plot(marker="o", axis_name="time", sed_type="flux")

Fitting temporal models

We have the reconstructed lightcurve at this point. Now we want to fit a profile to the obtained light curves, using a joint fitting across the different datasets, while simultaneously minimising across the temporal model parameters as well. The temporal models can be applied

• directly on the obtained lightcurve

• on the simulated datasets

Fitting the obtained light curve

The fitting will proceed through a joint fit of the flux points. First, we need to obtain a set of FluxPointDatasets, one for each time bin

## Create the datasets by iterating over the returned lightcurve
datasets = Datasets()

for idx, fp in enumerate(lc_1d.iter_by_axis(axis_name="time")):
    dataset = FluxPointsDataset(data=fp, name=f"time-bin-{idx}")
    datasets.append(dataset)

We will fit the amplitude, spectral index and the decay time scale. Note that t_ref should be fixed by default for the ExpDecayTemporalModel.

# Define the model:
spectral_model1 = PowerLawSpectralModel(
    index=2.0, amplitude="1e-12 cm-2 s-1 TeV-1", reference="1 TeV"
)
temporal_model1 = ExpDecayTemporalModel(t0="10 h", t_ref=gti_t0.mjd * u.d)

model = SkyModel(
    spectral_model=spectral_model1,
    temporal_model=temporal_model1,
    name="model-test",
)

/Users/adonath/software/mambaforge/envs/gammapy-dev/lib/python3.9/site-packages/astropy/units/quantity.py:613: RuntimeWarning: overflow encountered in exp
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)
/Users/adonath/software/mambaforge/envs/gammapy-dev/lib/python3.9/site-packages/astropy/units/quantity.py:613: RuntimeWarning: invalid value encountered in subtract
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)

datasets.models = model

%%time
# Do a joint fit
fit = Fit()
result = fit.run(datasets=datasets)

CPU times: user 6.8 s, sys: 96.2 ms, total: 6.9 s
Wall time: 8.54 s

/Users/adonath/software/mambaforge/envs/gammapy-dev/lib/python3.9/site-packages/astropy/units/quantity.py:613: RuntimeWarning: overflow encountered in exp
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)
/Users/adonath/software/mambaforge/envs/gammapy-dev/lib/python3.9/site-packages/astropy/units/quantity.py:613: RuntimeWarning: invalid value encountered in subtract
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)

Now let’s plot model and data. We plot only the normalisation of the temporal model in relative units for one particular energy range

lc_1TeV_10TeV = lc_1d.slice_by_idx({"energy": slice(2, 3)})
ax = lc_1TeV_10TeV.plot(sed_type="norm", axis_name="time")

time_range = lc_1TeV_10TeV.geom.axes["time"].time_bounds
temporal_model1.plot(ax=ax, time_range=time_range, label="Best fit model")

ax.set_yscale("linear")
plt.legend()

<matplotlib.legend.Legend at 0x16c5c0f70>

Fit the datasets

Here, we apply the models directly to the simulated datasets.

For modelling and fitting more complex flares, you should attach the relevant model to each group of datasets. The parameters of a model in a given group of dataset will be tied. For more details on joint fitting in gammapy, see here.

# Define the model:
spectral_model2 = PowerLawSpectralModel(
    index=2.0, amplitude="1e-12 cm-2 s-1 TeV-1", reference="1 TeV"
)
temporal_model2 = ExpDecayTemporalModel(t0="10 h", t_ref=gti_t0.mjd * u.d)

model2 = SkyModel(
    spectral_model=spectral_model2,
    temporal_model=temporal_model2,
    name="model-test2",
)

/Users/adonath/software/mambaforge/envs/gammapy-dev/lib/python3.9/site-packages/astropy/units/quantity.py:613: RuntimeWarning: overflow encountered in exp
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)
/Users/adonath/software/mambaforge/envs/gammapy-dev/lib/python3.9/site-packages/astropy/units/quantity.py:613: RuntimeWarning: invalid value encountered in subtract
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)

model2.parameters.to_table()

Table length=5

type	name	value	unit	error	min	max	frozen	is_norm	link
str8	str9	float64	str14	int64	float64	float64	bool	bool	str1
spectral	index	2.0000e+00		0.000e+00	nan	nan	False	False
spectral	amplitude	1.0000e-12	cm-2 s-1 TeV-1	0.000e+00	nan	nan	False	True
spectral	reference	1.0000e+00	TeV	0.000e+00	nan	nan	True	False
temporal	t0	1.0000e+01	h	0.000e+00	nan	nan	False	False
temporal	t_ref	5.8909e+04	d	0.000e+00	nan	nan	True	False

datasets.models = model2

%%time
# Do a joint fit
fit = Fit()
result = fit.run(datasets=datasets)

CPU times: user 5.75 s, sys: 71.6 ms, total: 5.83 s
Wall time: 7.26 s

/Users/adonath/software/mambaforge/envs/gammapy-dev/lib/python3.9/site-packages/astropy/units/quantity.py:613: RuntimeWarning: overflow encountered in exp
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)
/Users/adonath/software/mambaforge/envs/gammapy-dev/lib/python3.9/site-packages/astropy/units/quantity.py:613: RuntimeWarning: invalid value encountered in subtract
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)

result.parameters.to_table()

Table length=5

type	name	value	unit	error	min	max	frozen	is_norm	link
str8	str9	float64	str14	float64	float64	float64	bool	bool	str1
spectral	index	2.9926e+00		2.868e-02	nan	nan	False	False
spectral	amplitude	9.3341e-12	cm-2 s-1 TeV-1	3.377e-13	nan	nan	False	True
spectral	reference	1.0000e+00	TeV	0.000e+00	nan	nan	True	False
temporal	t0	6.0565e+00	h	2.562e-01	nan	nan	False	False
temporal	t_ref	5.8909e+04	d	0.000e+00	nan	nan	True	False

We see that the fitted parameters are consistent between fitting flux points and datasets, and match well with the simulated ones

Exercises

1. Re-do the analysis with MapDataset instead of SpectralDataset

2. Model the flare of PKS 2155-304 which you obtained using the light curve flare tutorial. Use a combination of a Gaussian and Exponential flare profiles, and fit using scipy.optimize.curve_fit

3. Do a joint fitting of the datasets.