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

%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

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.

# Loading IRFs
irfs = load_cta_irfs(
    "$GAMMAPY_DATA/cta-1dc/caldb/data/cta/1dc/bcf/South_z20_50h/irf_file.fits"
)

# 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", interp="log"
)

on_region_radius = Angle("0.11 deg")
on_region = CircleSkyRegion(center=center, radius=on_region_radius)

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

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

Table length=5

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

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

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

datasets = Datasets()

empty = SpectrumDataset.create(
    e_reco=energy_axis.edges,
    e_true=energy_axis_true.edges,
    region=on_region,
    name="empty",
)

for i in range(n_obs):
    obs = Observation.create(
        pointing=pointing,
        livetime=lvtm[i],
        tstart=tstart[i],
        irfs=irfs,
        reference_time=gti_t0,
    )
    empty_i = empty.copy(name=f"dataset_{i}")
    maker = SpectrumDatasetMaker(selection=["aeff", "background", "edisp"])
    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	livetime	n_on	background	excess	significance	background_rate	gamma_rate
	s					1 / s	1 / s
str9	float64	int64	float64	float64	float64	float64	float64
stacked	3299.999999999999	796	20.30377224384548	775.6962277561545	65.49294001748362	0.006152658255710753	0.2350594629564105
dataset_1	1500.0	343	9.228987383566126	333.77101261643384	42.57475095546916	0.006152658255710751	0.2225140084109559
dataset_2	1560.0	287	9.598146878908771	277.4018531210912	37.35770173570651	0.006152658255710751	0.17782170071864822
dataset_3	2400.0	324	14.766379813705804	309.23362018629416	37.186151161016724	0.006152658255710751	0.12884734174428925
dataset_4	2399.9999999999964	191	14.766379813705779	176.23362018629425	25.00844612269254	0.00615265825571075	0.07343067507762271
dataset_5	3000.0	182	18.457974767132253	163.54202523286776	22.492968959365548	0.006152658255710751	0.05451400841095592
dataset_6	2400.0	38	14.766379813705804	23.233620186294196	5.036912012105179	0.006152658255710751	0.009680675077622582
dataset_7	3120.0	51	19.196293757817543	31.803706242182457	6.004801820565906	0.006152658255710751	0.010193495590443096
dataset_8	2580.0	26	15.873858299733739	10.126141700266261	2.3250178763136047	0.006152658255710751	0.00392486112413421
dataset_9	2820.0	29	17.35049628110432	11.649503718895682	2.548357572372354	0.006152658255710752	0.00413102968755166

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.

# 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 each dataset
for dataset in datasets:
    dataset.models = model_fit

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

CPU times: user 1.16 s, sys: 16.6 ms, total: 1.17 s
Wall time: 1.2 s

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

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

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

model.parameters.to_table()

Table length=5

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

for dataset in datasets:
    dataset.models = model

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

CPU times: user 8.36 s, sys: 111 ms, total: 8.47 s
Wall time: 9.03 s

result.parameters.to_table()

Table length=5

name	value	unit	min	max	frozen	error
str9	float64	str14	float64	float64	bool	int64
index	2.983e+00		nan	nan	False	0.000e+00
amplitude	9.907e-12	cm-2 s-1 TeV-1	nan	nan	False	0.000e+00
reference	1.000e+00	TeV	nan	nan	True	0.000e+00
t0	2.497e-01	d	nan	nan	False	0.000e+00
t_ref	5.891e+04	d	nan	nan	True	0.000e+00

We see that the fitted parameters 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.