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: light_curve_simulation.ipynb | light_curve_simulation.py

Binned light curve simulation and fitting¶

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. Then, we re-fit the simulated datasets to reconstruct back the injected profiles.

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
Optionaly, we show here how to fit the simulated datasets using a source model

Setup¶

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

Setup¶

[1]:

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

import logging

log = logging.getLogger(__name__)

And some gammapy specific imports

[2]:

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

Simulating a light curve¶

We will simulate 10 datasets using an PowerLawSpectralModel and a ExpDecayTemporalModel. The important thing to note here is how to attach a different GTI to each dataset.

[3]:

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

[4]:

# Reconstructed and true energy axis
center = SkyCoord(0.0, 0.0, unit="deg", frame="galactic")
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")
on_region = CircleSkyRegion(center=center, radius=on_region_radius)

[5]:

# 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.

[6]:

# 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/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/astropy/units/quantity.py:477: RuntimeWarning: overflow encountered in exp
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)
/Users/adonath/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/astropy/units/quantity.py:477: RuntimeWarning: invalid value encountered in subtract
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)

[7]:

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

[7]:

Table length=5

name	value	unit	min	max	frozen	error
str9	float64	str14	float64	float64	bool	int64
index	3.0000e+00		nan	nan	False	0.000e+00
amplitude	1.0000e-11	cm-2 s-1 TeV-1	nan	nan	False	0.000e+00
reference	1.0000e+00	TeV	nan	nan	True	0.000e+00
t0	2.5000e-01	d	nan	nan	False	0.000e+00
t_ref	5.8909e+04	d	nan	nan	True	0.000e+00

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

[8]:

n_obs = 10
tstart = [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

[9]:

datasets = Datasets()

empty = SpectrumDataset.create(
    e_reco=energy_axis, e_true=energy_axis_true, region=on_region, 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,
    )
    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.

[10]:

datasets.info_table()

[10]:

Table length=10

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
								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	805	20.30377174963533	784.6962282503647	65.99594588024884	825.925411153998	20.30377174963533	805.6216394043627	216137902.05823946	16025275854.086987	3299.999999999999	3299.999999999999	0.243939393939394	0.006152658105950102	0.23778673583344392	9	9	cash	-6456.367158795519
dataset-1	298	9.228987158925152	288.7710128410748	38.64460511081799	332.14278827322767	9.228987158925152	322.91380111430254	98244500.93556342	7284216297.312269	1500.0	1500.0	0.19866666666666666	0.006152658105950101	0.19251400856071654	9	9	cash	-1797.4710822441302
dataset-2	276	9.598146645282158	266.40185335471784	36.34929048817018	293.48250627100055	9.598146645282158	283.8843596257184	102174280.97298595	7575584949.20476	1560.0	1560.0	0.17692307692307693	0.006152658105950101	0.17077041881712682	9	9	cash	-1559.7025929394501
dataset-3	311	14.766379454280242	296.23362054571976	36.09765120832382	321.772945022102	14.766379454280242	307.0065655678218	157191201.49690145	11654746075.69963	2400.0	2400.0	0.12958333333333333	0.006152658105950101	0.12343067522738323	9	9	cash	-1883.9474995748571
dataset-4	204	14.76637945428022	189.2336205457198	26.32196264783094	200.97527420424075	14.76637945428022	186.2088947499605	157191201.4969012	11654746075.699614	2399.9999999999964	2399.9999999999964	0.08500000000000013	0.0061526581059501	0.07884734189405003	9	9	cash	-1084.5718690277151
dataset-5	212	18.457974317850304	193.5420256821497	25.45462664361616	182.9858048173724	18.457974317850304	164.52783049952214	196489001.87112683	14568432594.624538	3000.0	3000.0	0.07066666666666667	0.006152658105950101	0.06451400856071657	9	9	cash	-1148.7554107527278
dataset-6	40	14.766379454280242	25.23362054571976	5.408778067740319	39.96701296644291	14.766379454280242	25.20063351216266	157191201.49690145	11654746075.69963	2400.0	2400.0	0.016666666666666666	0.006152658105950101	0.010514008560716567	9	9	cash	-79.45903376664398
dataset-7	31	19.196293290564316	11.803706709435684	2.4713006084303153	42.290648816036764	19.196293290564316	23.094355525472444	204348561.9459719	15151169898.40952	3120.0	3120.0	0.009935897435897435	0.006152658105950101	0.0037832393299473345	9	9	cash	-47.659394690074194
dataset-8	29	15.873857913351259	13.126142086648741	2.9495432158022146	32.23817480009135	15.873857913351259	16.364316886740088	168980541.60916907	12528852031.377102	2580.0	2580.0	0.01124031007751938	0.0061526581059501	0.005087651971569279	9	9	cash	-56.277645546824964
dataset-9	38	17.350495858779286	20.649504141220714	4.275784996983856	32.409018797086	17.350495858779286	15.058522938306712	184699661.75885922	13694326638.947065	2820.0	2820.0	0.01347517730496454	0.006152658105950101	0.007322519199014437	9	9	cash	-81.4508661140261

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.

[11]:

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

[12]:

# Attach model to each dataset
for dataset in datasets:
    dataset.models = model_fit

[13]:

%%time
lc_maker_1d = LightCurveEstimator(
    energy_edges=[energy_axis.edges[0], energy_axis.edges[-1]],
    source="model-fit",
)
lc_1d = lc_maker_1d.run(datasets)

CPU times: user 4.93 s, sys: 50.7 ms, total: 4.98 s
Wall time: 5.05 s

[14]:

lc_1d.table["is_ul"] = lc_1d.table["ts"] < 1

[15]:

ax = lc_1d.plot(marker="o", label="3D")

../_images/tutorials_light_curve_simulation_24_0.png

We have the reconstructed lightcurve at this point. Further standard analyis might involve modeling the temporal profiles with an analytical or theoretical model. You may do this using your favourite fitting package, one possible option being curve_fit inside scipy.optimize.

In the next section, we show how to simulatenously fit the all datasets using a given temporal model. This does a joint fitting across the different datasets, while simultaneously miniminsing across the temporal model parameters as well. We will fit the amplitude, spectral index and the decay time scale. Note that t_ref should be fixed by default for the ExpDecayTemporalModel.

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

Fit the datasets¶

[16]:

# 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/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/astropy/units/quantity.py:477: RuntimeWarning: overflow encountered in exp
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)
/Users/adonath/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/astropy/units/quantity.py:477: RuntimeWarning: invalid value encountered in subtract
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)

[17]:

model.parameters.to_table()

[17]:

Table length=5

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

[18]:

datasets.models = model

[19]:

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

CPU times: user 8.73 s, sys: 83.1 ms, total: 8.82 s
Wall time: 9.11 s

[20]:

result.parameters.to_table()

[20]:

Table length=5

name	value	unit	min	max	frozen	error
str9	float64	str14	float64	float64	bool	float64
index	2.9993e+00		nan	nan	False	3.228e-02
amplitude	9.4745e-12	cm-2 s-1 TeV-1	nan	nan	False	3.295e-13
reference	1.0000e+00	TeV	nan	nan	True	0.000e+00
t0	2.5955e-01	d	nan	nan	False	8.852e-03
t_ref	5.8909e+04	d	nan	nan	True	0.000e+00

We see that the fitted parameters match well with the simulated ones!

Exercises¶

Re-do the analysis with MapDataset instead of SpectralDataset
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
Do a joint fitting of the datasets.

[ ]: