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

# Simulating and fitting a time varying source#

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

[1]:

%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

[2]:

from gammapy.data import Observation
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:

[3]:

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.

[4]:

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

[5]:

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

[6]:

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

[7]:

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

[8]:

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

[8]:

Table length=5
str8str9float64str14int64float64float64boolboolstr1
spectralindex3.0000e+000.000e+00nannanFalseFalse
spectralamplitude1.0000e-11cm-2 s-1 TeV-10.000e+00nannanFalseTrue
spectralreference1.0000e+00TeV0.000e+00nannanTrueFalse
temporalt06.0000e+00h0.000e+00nannanFalseFalse
temporalt_ref5.8909e+04d0.000e+00nannanTrueFalse

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

[9]:

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

[10]:

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.

[11]:

datasets.info_table()

[11]:

Table length=10
namecountsexcesssqrt_tsbackgroundnprednpred_backgroundnpred_signalexposure_minexposure_maxlivetimeontimecounts_ratebackground_rateexcess_raten_binsn_fit_binsstat_typestat_sum
m2 sm2 sss1 / s1 / s1 / s
str9int64float64float64float64float64float64float64float64float64float64float64float64float64float64int64int64str4float64
dataset-0788767.681213378906265.0351559538021920.31876945495605520.31876973807811720.318769454956055nan216137904.016025275392.03299.99999999999733300.00.2387878787878790.0061572028651382040.2326306707208808799cashnan
dataset-1303293.76419579277739.087114845702099.235804207223024332.149612809808969.235804207223024322.913808602585998244500.935563627284216297.3122861500.00000000000361500.00000000000360.20199999999999950.0061572028048153350.195842797195184299cash-1804.171041139897
dataset-2272262.3947636244880635.973198692876889.605236375511943293.489602584430949.605236375511943283.88436620891895102174280.972986167575584949.2047781560.00000000000361560.00000000000360.174358974358973960.0061572028048153340.1682017715541586499cash-1574.3589928828487
dataset-3335320.222713268443238.0874105747017114.777286731556805321.7838594187763314.777286731556805307.0065726872195157191201.4969014211654746075.699632400.02400.00.139583333333333340.00615720280481533550.13342613052851899cash-2150.529986042023
dataset-4182167.222713268443224.07338294634657714.777286731556805200.986185799650614.777286731556805186.2088990680938157191201.4969014211654746075.699632400.02400.00.075833333333333340.00615720280481533550.06967613052851899cash-875.6042880128907
dataset-5154135.52839158555419.54799932222828618.471608414446006182.9994427293232518.471608414446006164.52783431487723196489001.8711267714568432594.6245383000.03000.00.0513333333333333350.00615720280481533550.04517613052851899cash-708.2822693043212
dataset-65843.2227132684431948.49513798452084914.77728673155680539.9779208281152214.77728673155680525.20063409655841157191201.4969014211654746075.699632400.02400.00.0241666666666666660.00615720280481533550.0180094638618513399cash-147.1795589597059
dataset-74525.7895272489761575.00295580947302119.21047275102384342.3048288120480419.21047275102384323.094356061024197204348561.9459718515151169898.409523120.03120.00.0144230769230769240.0061572028048153350.0082658741182615999cash-126.82594551244152
dataset-8237.1144167635764361.671810149927976315.88558323642356432.2499005026476615.88558323642356416.36431726622409168980541.60916912528852031.3771022580.02580.00.0089147286821705430.0061572028048153350.00275752587735520899cash-24.280591740673064
dataset-92810.6366880904207552.342209664322563617.36331190957924532.4218351970889517.36331190957924515.05852328750971184699661.7588591613694326638.9470652820.02820.00.0099290780141843980.00615720280481533550.00377187520936906299cash-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ç

[12]:

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

[13]:

# Attach model to all datasets
datasets.models = model_fit

[14]:

%%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

[15]:

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

[16]:

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

[17]:

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

[18]:

datasets.models = model

[19]:

%%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

[20]:

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

[20]:

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

[21]:

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

[22]:

model2.parameters.to_table()

[22]:

Table length=5
str8str9float64str14int64float64float64boolboolstr1
spectralindex2.0000e+000.000e+00nannanFalseFalse
spectralamplitude1.0000e-12cm-2 s-1 TeV-10.000e+00nannanFalseTrue
spectralreference1.0000e+00TeV0.000e+00nannanTrueFalse
temporalt01.0000e+01h0.000e+00nannanFalseFalse
temporalt_ref5.8909e+04d0.000e+00nannanTrueFalse
[23]:

datasets.models = model2

[24]:

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

[25]:

result.parameters.to_table()

[25]:

Table length=5
str8str9float64str14float64float64float64boolboolstr1
spectralindex2.9926e+002.868e-02nannanFalseFalse
spectralamplitude9.3341e-12cm-2 s-1 TeV-13.377e-13nannanFalseTrue
spectralreference1.0000e+00TeV0.000e+00nannanTrueFalse
temporalt06.0565e+00h2.562e-01nannanFalseFalse
temporalt_ref5.8909e+04d0.000e+00nannanTrueFalse

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.