\n",
"\n",
"**This is a fixed-text formatted version of a Jupyter notebook**\n",
"\n",
"- Try online [![Binder](https://static.mybinder.org/badge.svg)](https://mybinder.org/v2/gh/gammapy/gammapy-webpage/v0.18.1?urlpath=lab/tree/light_curve_simulation.ipynb)\n",
"- You can contribute with your own notebooks in this\n",
"[GitHub repository](https://github.com/gammapy/gammapy/tree/master/docs/tutorials).\n",
"- **Source files:**\n",
"[light_curve_simulation.ipynb](../_static/notebooks/light_curve_simulation.ipynb) |\n",
"[light_curve_simulation.py](../_static/notebooks/light_curve_simulation.py)\n",
"
\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Binned light curve simulation and fitting\n",
"\n",
"## Prerequisites:\n",
"\n",
"- To understand how a single binned simulation works, please refer to [spectrum_simulation](spectrum_simulation.ipynb) [simulate_3d](simulate_3d.ipynb) for 1D and 3D simulations respectively.\n",
"- For details of light curve extraction using gammapy, refer to the two tutorials [light_curve](light_curve.ipynb) and [light_curve_flare](light_curve_flare.ipynb) \n",
"\n",
"## Context\n",
"\n",
"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.\n",
"\n",
"**Objective: Simulate and fit a time decaying light curve of a source with CTA using the CTA 1DC response**\n",
"\n",
"## Proposed approach:\n",
"\n",
"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.\n",
"\n",
"In summary, necessary steps are:\n",
"\n",
"- Choose observation parameters including a list of `gammapy.data.GTI`\n",
"- Define temporal and spectral models from :ref:model-gallery as per science case\n",
"- Perform the simulation (in 1D or 3D)\n",
"- Extract the light curve from the reduced dataset as shown in [light curve notebook](light_curve.ipynb)\n",
"- Optionaly, we show here how to fit the simulated datasets using a source model \n",
"\n",
"\n",
"## Setup \n",
"\n",
"As usual, we'll start with some general imports..."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setup"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import astropy.units as u\n",
"from astropy.coordinates import SkyCoord, Angle\n",
"from astropy.time import Time\n",
"from regions import CircleSkyRegion\n",
"\n",
"import logging\n",
"\n",
"log = logging.getLogger(__name__)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And some gammapy specific imports"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"from gammapy.data import Observation\n",
"from gammapy.irf import load_cta_irfs\n",
"from gammapy.datasets import SpectrumDataset, Datasets\n",
"from gammapy.modeling.models import (\n",
" PowerLawSpectralModel,\n",
" ExpDecayTemporalModel,\n",
" SkyModel,\n",
")\n",
"from gammapy.maps import MapAxis\n",
"from gammapy.estimators import LightCurveEstimator\n",
"from gammapy.makers import SpectrumDatasetMaker\n",
"from gammapy.modeling import Fit"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Simulating a light curve\n",
"\n",
"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."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Invalid unit found in background table! Assuming (s-1 MeV-1 sr-1)\n"
]
}
],
"source": [
"# Loading IRFs\n",
"irfs = load_cta_irfs(\n",
" \"$GAMMAPY_DATA/cta-1dc/caldb/data/cta/1dc/bcf/South_z20_50h/irf_file.fits\"\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"# Reconstructed and true energy axis\n",
"center = SkyCoord(0.0, 0.0, unit=\"deg\", frame=\"galactic\")\n",
"energy_axis = MapAxis.from_edges(\n",
" np.logspace(-0.5, 1.0, 10), unit=\"TeV\", name=\"energy\", interp=\"log\"\n",
")\n",
"energy_axis_true = MapAxis.from_edges(\n",
" np.logspace(-1.2, 2.0, 31), unit=\"TeV\", name=\"energy_true\", interp=\"log\"\n",
")\n",
"\n",
"on_region_radius = Angle(\"0.11 deg\")\n",
"on_region = CircleSkyRegion(center=center, radius=on_region_radius)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"# Pointing position\n",
"pointing = SkyCoord(0.5, 0.5, unit=\"deg\", frame=\"galactic\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/adonath/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/astropy/units/quantity.py:477: RuntimeWarning: overflow encountered in exp\n",
" result = super().__array_ufunc__(function, method, *arrays, **kwargs)\n",
"/Users/adonath/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/astropy/units/quantity.py:477: RuntimeWarning: invalid value encountered in subtract\n",
" result = super().__array_ufunc__(function, method, *arrays, **kwargs)\n"
]
}
],
"source": [
"# Define the source model: A combination of spectral and temporal model\n",
"\n",
"gti_t0 = Time(\"2020-03-01\")\n",
"spectral_model = PowerLawSpectralModel(\n",
" index=3, amplitude=\"1e-11 cm-2 s-1 TeV-1\", reference=\"1 TeV\"\n",
")\n",
"temporal_model = ExpDecayTemporalModel(t0=\"6 h\", t_ref=gti_t0.mjd * u.d)\n",
"\n",
"model_simu = SkyModel(\n",
" spectral_model=spectral_model,\n",
" temporal_model=temporal_model,\n",
" name=\"model-simu\",\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"Table length=5\n",
"
\n",
"
name
value
unit
min
max
frozen
error
\n",
"
str9
float64
str14
float64
float64
bool
int64
\n",
"
index
3.0000e+00
nan
nan
False
0.000e+00
\n",
"
amplitude
1.0000e-11
cm-2 s-1 TeV-1
nan
nan
False
0.000e+00
\n",
"
reference
1.0000e+00
TeV
nan
nan
True
0.000e+00
\n",
"
t0
2.5000e-01
d
nan
nan
False
0.000e+00
\n",
"
t_ref
5.8909e+04
d
nan
nan
True
0.000e+00
\n",
"
"
],
"text/plain": [
"
\n",
" name value unit min max frozen error \n",
" str9 float64 str14 float64 float64 bool int64 \n",
"--------- ---------- -------------- ------- ------- ------ ---------\n",
" index 3.0000e+00 nan nan False 0.000e+00\n",
"amplitude 1.0000e-11 cm-2 s-1 TeV-1 nan nan False 0.000e+00\n",
"reference 1.0000e+00 TeV nan nan True 0.000e+00\n",
" t0 2.5000e-01 d nan nan False 0.000e+00\n",
" t_ref 5.8909e+04 d nan nan True 0.000e+00"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Look at the model\n",
"model_simu.parameters.to_table()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, define the start and observation livetime wrt to the reference time, `gti_t0`"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"n_obs = 10\n",
"tstart = [1, 2, 3, 5, 8, 10, 20, 22, 23, 24] * u.h\n",
"lvtm = [55, 25, 26, 40, 40, 50, 40, 52, 43, 47] * u.min"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now perform the simulations"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"datasets = Datasets()\n",
"\n",
"empty = SpectrumDataset.create(\n",
" e_reco=energy_axis, e_true=energy_axis_true, region=on_region, name=\"empty\"\n",
")\n",
"\n",
"maker = SpectrumDatasetMaker(selection=[\"exposure\", \"background\", \"edisp\"])\n",
"\n",
"for idx in range(n_obs):\n",
" obs = Observation.create(\n",
" pointing=pointing,\n",
" livetime=lvtm[idx],\n",
" tstart=tstart[idx],\n",
" irfs=irfs,\n",
" reference_time=gti_t0,\n",
" obs_id=idx,\n",
" )\n",
" empty_i = empty.copy(name=f\"dataset-{idx}\")\n",
" dataset = maker.run(empty_i, obs)\n",
" dataset.models = model_simu\n",
" dataset.fake()\n",
" datasets.append(dataset)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The reduced datasets have been successfully simulated. Let's take a quick look into our datasets."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"Table length=10\n",
"
\n",
"
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
\n",
"
m2 s
m2 s
s
s
1 / s
1 / s
1 / s
\n",
"
str9
int64
float64
float64
float64
float64
float64
float64
float64
float64
float64
float64
float64
float64
float64
int64
int64
str4
float64
\n",
"
dataset-0
799
20.30377174963533
778.6962282503647
65.66086564575978
825.925411153998
20.30377174963533
805.6216394043627
216137902.05823946
16025275854.086987
3299.999999999999
3299.999999999999
0.24212121212121218
0.006152658105950102
0.2359685540152621
9
9
cash
-6410.135496836574
\n",
"
dataset-1
349
9.228987158925152
339.7710128410748
43.082459985513815
332.14278827322767
9.228987158925152
322.91380111430254
98244500.93556342
7284216297.312269
1500.0
1500.0
0.23266666666666666
0.006152658105950101
0.22651400856071655
9
9
cash
-2211.4834010173
\n",
"
dataset-2
275
9.598146645282158
265.40185335471784
36.256818460507986
293.48250627100055
9.598146645282158
283.8843596257184
102174280.97298595
7575584949.20476
1560.0
1560.0
0.1762820512820513
0.006152658105950101
0.17012939317610118
9
9
cash
-1636.4473425792842
\n",
"
dataset-3
306
14.766379454280242
291.23362054571976
35.674176138470614
321.772945022102
14.766379454280242
307.0065655678218
157191201.49690145
11654746075.69963
2400.0
2400.0
0.1275
0.006152658105950101
0.1213473418940499
9
9
cash
-1820.364859399368
\n",
"
dataset-4
189
14.76637945428022
174.2336205457198
24.803300586459546
200.97527420424075
14.76637945428022
186.2088947499605
157191201.4969012
11654746075.699614
2399.9999999999964
2399.9999999999964
0.07875000000000013
0.0061526581059501
0.07259734189405002
9
9
cash
-974.660229500435
\n",
"
dataset-5
210
18.457974317850304
191.5420256821497
25.26247386864822
182.9858048173724
18.457974317850304
164.52783049952214
196489001.87112683
14568432594.624538
3000.0
3000.0
0.07
0.006152658105950101
0.0638473418940499
9
9
cash
-1114.2590725169475
\n",
"
dataset-6
47
14.766379454280242
32.23362054571976
6.660739995200845
39.96701296644291
14.766379454280242
25.20063351216266
157191201.49690145
11654746075.69963
2400.0
2400.0
0.019583333333333335
0.006152658105950101
0.013430675227383233
9
9
cash
-125.64973733510145
\n",
"
dataset-7
33
19.196293290564316
13.803706709435684
2.854951884968713
42.290648816036764
19.196293290564316
23.094355525472444
204348561.9459719
15151169898.40952
3120.0
3120.0
0.010576923076923078
0.006152658105950101
0.0044242649709729754
9
9
cash
-38.16647740667102
\n",
"
dataset-8
31
15.873857913351259
15.126142086648741
3.353380282915971
32.23817480009135
15.873857913351259
16.364316886740088
168980541.60916907
12528852031.377102
2580.0
2580.0
0.012015503875968992
0.0061526581059501
0.005862845770018892
9
9
cash
-55.061167541004416
\n",
"
dataset-9
37
17.350495858779286
19.649504141220714
4.0915714473025995
32.409018797086
17.350495858779286
15.058522938306712
184699661.75885922
13694326638.947065
2820.0
2820.0
0.013120567375886525
0.006152658105950101
0.0069679092699364235
9
9
cash
-66.46581488391693
\n",
"
"
],
"text/plain": [
"
\n",
" name counts background ... n_fit_bins stat_type stat_sum \n",
" ... \n",
" str9 int64 float64 ... int64 str4 float64 \n",
"--------- ------ ------------------ ... ---------- --------- -------------------\n",
"dataset-0 799 20.30377174963533 ... 9 cash -6410.135496836574\n",
"dataset-1 349 9.228987158925152 ... 9 cash -2211.4834010173\n",
"dataset-2 275 9.598146645282158 ... 9 cash -1636.4473425792842\n",
"dataset-3 306 14.766379454280242 ... 9 cash -1820.364859399368\n",
"dataset-4 189 14.76637945428022 ... 9 cash -974.660229500435\n",
"dataset-5 210 18.457974317850304 ... 9 cash -1114.2590725169475\n",
"dataset-6 47 14.766379454280242 ... 9 cash -125.64973733510145\n",
"dataset-7 33 19.196293290564316 ... 9 cash -38.16647740667102\n",
"dataset-8 31 15.873857913351259 ... 9 cash -55.061167541004416\n",
"dataset-9 37 17.350495858779286 ... 9 cash -66.46581488391693"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"datasets.info_table()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Extract the lightcurve\n",
"\n",
"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."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"# Define the model:\n",
"spectral_model = PowerLawSpectralModel(\n",
" index=3, amplitude=\"1e-11 cm-2 s-1 TeV-1\", reference=\"1 TeV\"\n",
")\n",
"model_fit = SkyModel(spectral_model=spectral_model, name=\"model-fit\")"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"# Attach model to each dataset\n",
"for dataset in datasets:\n",
" dataset.models = model_fit"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"CPU times: user 4 s, sys: 43.7 ms, total: 4.04 s\n",
"Wall time: 4.29 s\n"
]
}
],
"source": [
"%%time\n",
"lc_maker_1d = LightCurveEstimator(\n",
" energy_edges=[energy_axis.edges[0], energy_axis.edges[-1]],\n",
" source=\"model-fit\",\n",
")\n",
"lc_1d = lc_maker_1d.run(datasets)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"lc_1d.table[\"is_ul\"] = lc_1d.table[\"ts\"] < 1"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAERCAYAAAB2CKBkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dfZxdVX3v8c/XcSBjAx0gQWBCDGJIfeAh6UhI0+sr4cFApCZSegP3VkV9NUGlpWgTSdtbi32AFqvyoCbxCbjWiKZppBiCaFDRJoRJJpmAPBihDTOkEqED5jqQZPK7f+x98MzJmZkzM2efM2fO9/16zWv2XnudfX4bMud31tprr6WIwMzM6terqh2AmZlVlxOBmVmdcyIwM6tzTgRmZnXOicDMrM45EZiZ1bmaTASSvizpWUkPl+l8GyR1S7q7oPwqSbskhaQJ5XgvM7PRpiYTAXAbcGEZz3cj8O4i5T8Gzgf+s4zvZWY2qtRkIoiIHwLP55dJOjX9Zr9V0gOSfmsI5/se8Msi5e0R8R8jDtjMbBR7dbUDKKNVwJUR8VNJM4HPAedWOSYzs1FvTCQCSeOB3wG+KSlXfGR67BLgE0Ve1hUR8yoToZnZ6DUmEgFJF1d3RJxVeCAi1gJrKx+SmVltqMl7BIUi4kXgKUl/AKDEmVUOy8ysJtRkIpC0GtgETJPUKekDwP8GPiBpB/AIsGAI53sA+CZwXnq+eWn5n0jqBCYBHZK+WO5rMTOrNnkaajOz+laTLQIzMyufmrtZPGHChJgyZUq1wzAzqylbt279RURMLHas5hLBlClTaGtrq3YYZmY1RVK/MyS4a8jMrM45EZiZ1bnME4GkBknthTN7psck6eZ0hs8OSTOyjsfMzPqqxD2Cq4FHgaOLHLsImJr+zAQ+n/42M6u4AwcO0NnZyUsvvVTtUIZt3LhxTJo0icbGxpJfk2kikDQJeAfwd8BHilRZANwRycMMmyU1SzoxIvZkGZeZWTGdnZ0cddRRTJkyhbx5y2pGRPDcc8/R2dnJKaecUvLrsu4a+gywDDjUz/EW4Om8/c60rA9JiyW1SWrbu3dv+aM0MwNeeukljjvuuJpMAgCSOO6444bcosksEUi6GHg2IrYOVK1I2WGPOkfEqohojYjWiROLDoM1MyuLWk0COcOJP8uuodnAOyXNB8YBR0v6akT8YV6dTuDkvP1JwDMZxmRmVhafvu8JbvreTw8rv/q8qVxzwWlViGj4MksEEbEcWA4gaQ7wZwVJAOAu4CpJXye5SfxCpe8PLFq5CYA7l8yq5NuaWY275oLTuOaC08r+GfLSSy/xtre9jZdffpmDBw9y6aWXct1113HFFVfwgx/8gKOPPpqenh7OOeccrr/+elpaDutNH7KKP0cg6UpJV6a764EngV3AF4APVToeM7Oh+PR9T7yyva69i/bd3Tz41PPMvmEj69q7itYbiiOPPJKNGzeyY8cOtm/fzoYNG9i8eTMAN954Izt27ODxxx9n+vTpzJ07l/3794/sgqhQIoiI70fExen2iohYkW5HRHw4Ik6NiNMjwnNHmNmolusOWtfexfK1O9nfm4yF6eruYfnana8kg2LdRqWQxPjx44FkOOuBAwcO6/eXxDXXXMMJJ5zAPffcM9xLeYWfLDYzG6JFKzexbE0HPQd6+5T3HOhl2ZqOV7qLhqu3t5ezzjqL448/ngsuuICZM4s/XjVjxgwee+yxEb0X1HkiGKhZZ2Y2kFxLoNTyoWhoaGD79u10dnayZcsWHn744aL1yrWeTN0mgsGadYUWrdw04ixvZmPDnUtm0dLcVPRYS3NT2W4cNzc3M2fOHDZs2FD0eHt7O2984xtH/D41Nw31SOR/kLfv7j4sc+eadau37AY8ksjM+rd03jSWr93Zp3uoqbGBpfOmjei8e/fupbGxkebmZnp6evjud7/Lxz72sT6tgojglltuYc+ePVx44YUjej+o4xZBls06Mxu7rj5vKgALp7dw/SWnc0RD8jHa0tzE9ZeczsLpLX3qDdWePXuYO3cuZ5xxBm9961u54IILuPjiiwFYunQpZ555JqeddhoPPfQQ999/P0ccccSIr6muWgT53/Bn37CRru6ew+qUs1lnZmNP/sNiC6e39NuDMNyHys444wza29sPK7/tttuGdb5S1G2LYOm8aTQ1NvQpK0ezzszqw6fve4Ip136bB596ngefep4p136bKdd+e9jPD1RTXbUI8uWab8vWdLC/9xAtzU0snTftlXIzs4HkniweC+o2EcDAzTozq08RUdMTzw1nSGnddg0NhZ83MKsP48aN47nnnivb+PxKy61HMG7cuCG9rq5bBKXo73kDwN1IZmPMpEmT6OzspJbXPcmtUDYUqrXM19raGm1tlZuS6LS/uKfokNKW5iZ+fO25FYvDzGwkJG2NiNZix9w1NIj+nit4psjQUzOzWuREMIj+HiM/qZ9yM7Na40QwCD9vYGZjnW8WD8LPG5jZWJdZIpA0DvghcGT6Pmsi4uMFdeYA3wKeSovWRsQnsoppuPy8gZmNZVm2CF4Gzo2IfZIagR9JuiciNhfUeyC3epmZmVVelovXB7Av3W1Mf2prrKqZWR3I9GaxpAZJ24Fngfsi4sEi1WZJ2iHpHklv7uc8iyW1SWqr5Qc9zMxGo0wTQUT0RsRZwCTgbElvKaiyDXhdRJwJ3AKs6+c8qyKiNSJaJ06cmGXIZmZ1pyLDRyOiG/g+cGFB+YsRsS/dXg80SppQiZjMzCyR5aihicCBiOiW1AScD/xDQZ0TgJ9HREg6myQxPZdVTCPh0UJmNlZlOWroROB2SQ0kH/DfiIi7JV0JEBErgEuBD0o6CPQAl0WtTX5kZlbjshw11AFML1K+Im/7VuDWrGIwM7PBeYoJM7M650RgZlbnnAjMzOqcE4GZWZ1zIjAzq3NOBGZmdc6JwMyszjkRmJnVOScCM7M650RgZlbnnAjMzOqcE4GZWZ1zIjAzq3NOBGZmdc6JwMyszjkRmJnVucwSgaRxkrZI2iHpEUnXFakjSTdL2iWpQ9KMrOIxM7Pislyq8mXg3IjYJ6kR+JGkeyJic16di4Cp6c9M4PPpbzMzq5DMWgSR2JfuNqY/hesRLwDuSOtuBpolnZhVTGZmdrhM7xFIapC0HXgWuC8iHiyo0gI8nbffmZYVnmexpDZJbXv37s0uYDOzOpRpIoiI3og4C5gEnC3pLQVVVOxlRc6zKiJaI6J14sSJWYRqZla3KjJqKCK6ge8DFxYc6gROztufBDxTiZjMzCyR5aihiZKa0+0m4HzgsYJqdwHvSUcPnQO8EBF7sorJzMwOl+WooROB2yU1kCScb0TE3ZKuBIiIFcB6YD6wC/gV8L4M4zEzsyIySwQR0QFML1K+Im87gA9nFUO1LVq5CYA7l8yqciRmZv3zk8VmZnXOicDMrM45EZiZ1TknAjOzOudEkJF17V207+7mwaeeZ/YNG1nX3lXtkMzMinIiyMC69i6Wr93J/t5DAHR197B87U4nAzMblZwIMrBsTQc9B3r7lPUc6OXGex+vUkRmZv1zIshAriVQ6JnungpHYmY2OCeCDLQ0NxUtP6mfcjOzanIiyMDSedNoamzoU9bU2MDSedOqFJGZWf+ynGuobi2cniypsGxNB/t7D9HS3MTSedNeKTczG01KSgSSjgdmAycBPcDDQFtEFO8MNxZOb2H1lt2A5xoys9FtwEQgaS5wLXAs0E6y0tg4YCFwqqQ1wD9FxItZB2pmZtkYrEUwH/ijiNhdeEDSq4GLgQuAf8kgNjMzq4ABE0FELB3g2EFgXdkjMjOzihr2qCFJXkTGzGwMGMmooeuAr/R3UNLJwB3ACcAhYFVE3FRQZw7wLeCptGhtRHxiBDGNKr5JbGa1YLCbxR39HQJeO8i5DwIfjYhtko4Ctkq6LyJ+UlDvgYi4uLRwzcys3AZrEbwWmAf8d0G5gH8f6IXpIvR70u1fSnoUaAEKE4GZmVXRYIngbmB8RGwvPCDp+6W+iaQpJOsXP1jk8CxJO4BngD+LiEeKvH4xsBhg8uTJpb6tmZmVQMn68Rm+gTQe+AHwdxGxtuDY0cChiNgnaT5wU0RMHeh8ra2t0dbWll3AZmZjkKStEdFa7NiQRw2l385LrdtI8ozBPxcmAYCIeDEi9qXb64FGSROGGpOZmQ3fcIaPXllKJUkCvgQ8GhGf6qfOCWk9JJ2dxvPcMGIyM7NhGs7wUZVYbzbwbmCnpNw9hj8HJgNExArgUuCDkg6SzGF0WWTdV2VmZn0MJxH8XimVIuJHDJI0IuJW4NZhxGBmZmVSUteQpKslHZ1241wnaZukt2ccm5mZVUCp9wjen84w+nZgIvA+4IbMojIzs4opNRHkunjmA1+JiB2Ufq/AzMxGsVITwVZJ3yFJBPemU0Z4URozszGg1JvFHwDOAp6MiF9JOo6ke8jMzGpcSYkgXZJyW97+c3i8v5nZmDDs9QjMzGxscCIwM6tzTgRmZnVuwEQg6XRJmyU9LWmVpGPyjm3JPjwzM8vaYC2CzwN/DZwOPAH8SNKp6bHGDOMyM7MKGWzU0PiI2JBuf1LSVmCDpHcDnhzOzGwMGCwRSNJvRsQLABFxv6TfJ1lj4NjMozMzs8wN1jX0D8Ab8wsiogM4DzhsoRkzM6s9A7YIIuJrhWWSToiI3cAfZRaVmZlVzHCGj64vexRmZlY1w0kEJc06KulkSfdLelTSI5KuLlJHkm6WtEtSh6QZw4jHgEUrN7Fo5aZqh2FmNWg4K5R9ocR6B4GPRsS2dLbSrZLui4if5NW5CJia/swkGa46cxgxmZnZMJXcIpB0jKQzgM2SZgz27T0i9kTEtnT7l8CjQEtBtQXAHZHYDDRLOnFol2BmZiNRUotA0t8AVwA/49fPDwRwbomvnwJMBx4sONQCPJ2335mW7Sl4/WJgMcDkyZNLecu6sq69i/bd3ezvPcTsGzaydN40Fk4vzLlmZsWV2jX0P4FTI2L/UN9A0niS5w7+NF3uss/hIi857EG1iFgFrAJobW31g2x51rV3sXztTvb3JusEdXX3sHztTgAnAzMrSaldQw8DzUM9uaRGkiTwzxFR7LmDTuDkvP1JwDNDfZ96tWjlJpat6aDnQG+f8p4DvSxb01GlqMys1pTaIrgeaJf0MPByrjAi3tnfCyQJ+BLwaER8qp9qdwFXSfo6yU3iFyJiTz91rYhcS6DUcjOzQqUmgttJnjLeSelrFc8G3g3slLQ9LftzYDJARKwgeSZhPrAL+BVe/nJI7lwyi9k3bKSru+ewYy3NTVWIyMxqUamJ4BcRcfNQThwRP2KQZw4iIoAPD+W81tfSedNYvnZnn+6hpsYGls6bVsWozKyWlJoItkq6nqQrJ79raFv/L7FKyN0QXramg/29h2hpbvKoITMbklITwfT09zl5ZSUPH7VsLZzewuotu4Gku8jMbChKSgQRMTfrQGxknADMbLhKGj4q6e8lNeftHyPpb7MLy8zMKqXU5wguioju3E5E/DfJaB8zM6txpSaCBklH5nYkNQFHDlDfzMxqRKk3i78KfE/SV0huEr+f5NkCMzOrcaXeLP5HSR3A+STPBvxNRNybaWRmZlYRAyYCSUof+iIiNgAbBqpjZma1Z7B7BPdL+mNJfeZ+lnSEpHMl3Q68N7vwzMwsa4N1DV1Icj9gtaRTgG6giSSBfAf4dERsH+D1ZmY2yg2YCCLiJeBzwOfSKaUnAD35Q0nNzKy2lbxmcUQcoGDlMDMzq30lr1lsZmZjkxOBmVmdK3WuoTcVKZtT9mjMzKziSm0RfEPSx5RoknQLyfKV/ZL0ZUnPpstbFjs+R9ILkranP3811ODNzGzkSk0EM0kWmf934CGSBeZnD/Ka20iGnw7kgYg4K/35RImxmJlZGZWaCA4APSTPEIwDnoqIAdcujogfAs+PLDwzM8taqYngIZJE8Fbgd4HLJa0pw/vPkrRD0j2S3txfJUmLJbVJatu7d28Z3tbMzHJKfY7gAxHRlm7/F7BA0rtH+N7bgNdFxD5J84F1wNRiFSNiFbAKoLW11fMamZmVUaktgmclTc7/AX4wkjeOiBcjYl+6vR5olDRhJOc0M7OhK7VF8G2SdQhEco/gFOBxoN/unMFIOgH4eUSEpLNJktJzwz2fmZkNT6nrEZyevy9pBrBkoNdIWg3MASZI6gQ+DjSm51sBXAp8UNJBkvsPl3k6azOzyit5rqF8EbFN0lsHqXP5IMdvBW4dzvubmVn5lJQIJH0kb/dVwAzAw3fMzMaAUlsER+VtHyS5Z/Av5Q/HzMwqrdR7BNdlHYiZmVXHYGsW/xvJaKGiIuKdZY/IzMwqarAWwScrEoWZmVXNYIngqYjYXZFIzMysKgZ7snhdbkOSbw7bkC1auYlFKzdVOwwzG8BgiUB526/PMhAzM6uOwRJB9LNtNqh17V207+7mwaeeZ/YNG1nX3lXtkMysiMHuEZwp6UWSlkFTuk26HxFxdKbRWc1a197F8rU72d+bLFvR1d3D8rU7AVg4vaWaoZlZgQETQUQ0VCoQGzsWrdxE++7uV5JATs+BXpat6WD1lt3cuWRWlaIzs0KlTkNtNiSFSWCwcjOrHicCK7s7l8yipbmp6LGW5ia3BsxGGScCy8TSedNoauzbs9jU2MDSedOqFJGZ9WdY01CbDSZ3Q3jZmg729x6ipbmJpfOmZXqjOPe8glscZkPjRGCZWTi9hdVbkgfTs/5wzg1V3d97iNk3bMw86ZiNJZklAklfBi4Gno2ItxQ5LuAmYD7wK+CKiNiWVTxWHZX4du6hqmYjk2WL4DaSFcju6Of4RcDU9Gcm8Pn0t1lJcl1BHqpqNjKZ3SyOiB8Czw9QZQFwRyQ2A82STswqHhu7PFTVbGSqOWqoBXg6b78zLTuMpMWS2iS17d3rFTItceeSWR6qalYG1UwEKlJWdD6jiFgVEa0R0Tpx4sSMw7Ja46GqZiNTzVFDncDJefuTgGeqFIvVsGoMVTUbS6qZCO4CrpL0dZKbxC9ExJ4qxmM1bOH0Fn/wmw1TlsNHVwNzgAmSOoGPA40AEbECWE8ydHQXyfDR92UVi5mZ9S+zRBARlw9yPIAPZ/X+ZmZWGs81ZGZW55wIzMzqnBOBmVmdcyIwM6tzTgRmZnXOicDMrM45EZiZ1TknAjOzOudEYGZW55wIzMzqnBOBmVmdcyIwM6tzTgRmZnXOicDMrM5Vc2EaMzMrsK69ixvvfZxnuns4qUKr7TkRmJmNEuvau1i+dic9B3oB6OruYfnanQCZJoNME4GkC4GbgAbgixFxQ8HxOcC3gKfSorUR8YksYzIzG20WrdwEQPvubvb3HupzrOdAL8vWdLB6y27uXDIrk/fPcqnKBuCzwAUkC9U/JOmuiPhJQdUHIuLirOIwM6uWoXbzFCaBwcrLJcsWwdnAroh4EiBdpH4BUJgIzMzGnKF08+S+6c++YSNd3T2HnauluSmz1gBkO2qoBXg6b78zLSs0S9IOSfdIenOG8ZiZZW7Ryk0sWrmJZWs6XkkCOblunv4snTeNpsaGPmVNjQ0snTctk1hzsmwRqEhZFOxvA14XEfskzQfWAVMPO5G0GFgMMHny5HLHaWZWdsPp5sm1FJat6WB/7yFaKjRqKMsWQSdwct7+JOCZ/AoR8WJE7Eu31wONkiYUnigiVkVEa0S0Tpw4McOQzcxG5s4ls7hzySxampuKHu+vPGfh9BamT25m5inH8uNrz808CUC2ieAhYKqkUyQdAVwG3JVfQdIJkpRun53G81yGMZmZVUS1unmGI7OuoYg4KOkq4F6S4aNfjohHJF2ZHl8BXAp8UNJBoAe4LCIKu4/MzGpOtbp5hiPT5wjS7p71BWUr8rZvBW7NMgYzs3IZ6nDQhdNbWL1lN0Cmo35Gyk8Wm5mVoFpP/VaCE4GZWT9yT/xC9Z76rQTPPmpmVoJqPfVbCW4RmJn1I/9bfrWe+q0EtwjMzEpQS8NBh8otAjOzEtTScNChciIws5q2rr2rYh/OtTIcdKjcNWRmNSs3pDN3wzY3pHNde1eVI6stbhGYWU1atHLTgEM6x0KXTaW4RWBmNWssD+msJLcIzKwm3blk1oBDOmvRJ/7tESLgnNcfB8Cn73uCF3sOIIm/+r03Zfa+TgRmVrOWzpvWZ9oHqO0hneOPfDVfeOBJeg78ukXT1NjA4re9PtP3ddeQmdWshdNbuP6S02lpbkIkLYHrLzm9Zu8PfGjuGxg/rrFP2VHjXs0H55ya6fu6RWBmNW3h9JZR+cE/nG6ecY0N/OOlZ/Chr26j50AvTY0N/MOlZzCu4EG2cnMiMDPLwHC7eeZOO57fft0x/PvPfkHrlGOYO+34rEN115CZWRZG0s1z/SWnc3rLb/L37zo9q/D6cCIwM8tArpsnNz/RULp5Tj72NXzrqt/l5GNfk3WYQMaJQNKFkh6XtEvStUWOS9LN6fEOSTOyjMfMrJJy3TyvEhXr5hmOzBKBpAbgs8BFwJuAyyUV3iG5CJia/iwGPp9VPGZm1VDpbp7hyLJFcDawKyKejIj9wNeBBQV1FgB3RGIz0CzpxAxjMjOrqEp38wxHlqOGWoCn8/Y7gZkl1GkB9uRXkrSYpMXA5MmTyx6omdlgqvXUbyVkmQhUpCyGUYeIWAWsAmhtbT3suJlZ1qr11G8lZNk11AmcnLc/CXhmGHXMzKquWk/9VkKWieAhYKqkUyQdAVwG3FVQ5y7gPenooXOAFyJiT+GJzMyqbSTDQUe7zBJBRBwErgLuBR4FvhERj0i6UtKVabX1wJPALuALwIeyisfMbKRqZTjoUCmitrrcW1tbo62trdphmFmdevr5X3HV17Zx6/+aMapHAhWStDUiWosd81xDZmZDkBsOOpZ4igkzszrnRGBmVuecCMzM6lzN3SyWtBf4T2AC8Isqh1MN9Xjd9XjNUJ/XXY/XDJW57tdFxMRiB2ouEeRIauvvDvhYVo/XXY/XDPV53fV4zVD963bXkJlZnXMiMDOrc7WcCFZVO4Aqqcfrrsdrhvq87nq8ZqjyddfsPQIzMyuPWm4RmJlZGTgRmJnVuZpJBJL+QNIjkg5J6neYlaQLJT0uaZekaysZYxYkHSvpPkk/TX8f00+9a9L/Pg9LWi1pXKVjLZchXHOzpDWSHpP0qKRZlY61nEq97rRug6R2SXdXMsZyK+WaJZ0s6f70//Ejkq6uRqwjNdhnUzod/83p8Q5JMyoVW80kAuBh4BLgh/1VkNQAfBa4CHgTcLmk2l5DDq4FvhcRU4Hvpft9SGoB/gRojYi3AA0k6z/UqkGvOXUTsCEifgs4k2S681pW6nUDXE3tXy+Uds0HgY9GxBuBc4AP19rfdYmfTRcBU9OfxcDnKxVfzSSCiHg0Ih4fpNrZwK6IeDIi9gNfBxZkH12mFgC3p9u3Awv7qfdqoEnSq4HXUNsrvQ16zZKOBt4GfAkgIvZHRHfFIsxGSf+vJU0C3gF8sUJxZWnQa46IPRGxLd3+JUkCbKlYhOVRymfTAuCOSGwGmiWdWIngaiYRlKgFeDpvv5Pa+wdT6LW5VdvS34ethBERXcAngd3AHpKV3r5T0SjLa9BrBl4P7AW+knaRfFHSb1QyyAyUct0AnwGWAYf6OV5LSr1mACRNAaYDD2YeWXmV8tlUtc+vUbUegaTvAicUOfQXEfGtUk5RpGzUj48d6LpLfP0xJN8mTgG6gW9K+sOI+Gr5oiyvkV4zyb/dGcAfR8SDkm4i6Vb4P2UKMRNl+H99MfBsRGyVNKecsWWlDP+vc+cZD/wL8KcR8WI5YqugUj6bqvb5NaoSQUScP8JTdAIn5+1Poga6SAa6bkk/l3RiROxJm4nPFql2PvBUROxNX7MW+B1g1CaCMlxzJ9AZEblvhmsYuE99VCjDdc8G3ilpPjAOOFrSVyPiDzMKecTKcM1IaiRJAv8cEWszCjVLpXw2Ve3za6x1DT0ETJV0iqQjSG6Y3lXlmEbqLuC96fZ7gWIto93AOZJeI0nAedT2jcRBrzki/gt4WtK0tOg84CeVCS8zpVz38oiYFBFTSP59bxzNSaAEg15z+m/6S8CjEfGpCsZWTqV8Nt0FvCcdPXQOSRfvnopEFxE18QO8iyRjvgz8HLg3LT8JWJ9Xbz7wBPAzki6lqsc+wus+jmQ0xU/T38f2c93XAY+RjK76v8CR1Y69Atd8FtAGdADrgGOqHXslrjuv/hzg7mrHnfU1A79L0kXSAWxPf+ZXO/ZhXOthn03AlcCV6bZIRhb9DNhJMgqwIrF5igkzszo31rqGzMxsiJwIzMzqnBOBmVmdcyIwM6tzTgRmZiMg6a8ldUnanv7M76de0YkhJZ0paZOknZL+LZ0+BUnHpZPt7ZN06xBjukXSvlLrOxHYmJb+MeX+QP8r7w92n6TPZfSefyrpPen2bZJ+JemovOM3SQpJE9L9fenvKZJ60ikzHpW0RdJ78153saTrsojZSiNpjqTbihz6dESclf6sL/K6gSaG/CJwbUScDvwrsDQtf4nkSfk/G2KMrUDzUF7jRGBjWkQ8l/sDBVbw6z/Y8RHxoXK/Xzrp3/uBr+UV7yKdYEzSq4C5QFc/p/hZREyPZKbNy4BrJL0vPfZtkqeKX1PuuK0i+psYchq/nlX5PuD3ASLi/0XEj0gSQh+S3p62IrZJ+mY6/UZultMbSeaiKpkTgdWl9Jvd3en2X0u6XdJ3JP2HpEsk/WPaVN+QTm+ApN+W9ANJWyXd28/MkOcC2yLiYF7ZamBRuj0H+DHJ1MoDiogngY+QfJMkkod+vg9cPKyLtixdpWQNgS+ryJoKMfDEkA8D70y3/4C+00wcJm1J/iVwfkTMIHmo8iO5OIC7YohPJDsRmCVOJZnaeQHJHE33p031HuAdaTK4Bbg0In4b+DLwd0XOMxvYWlD2U2Bi+gFxOckUxKXaBvxW3n4b8D+G8HorA0kPStpO0o3zzrzuxnkk6wacSvKk+x7gn4q8Pn9iyJOA35CUmxrk/SRrLGwFjgL2DxLOOSRrGvw4jem9wOsknUSSSG4Z6vWNqknnzKronog4IGknSf/thrR8JzCFpPn+FuC+ZOobGhc3O/oAAAIOSURBVEj+6AudSPF5ntaSdPXMBJYMIa7CGSmfJfkgsQqKiJmQtCSBKyLiimL1JH0BKLZqXL8TQ0bEY8Db0/LTSL6QDETAfRFxecF7vwN4A7Ar/Tf6Gkm7IuINg12fE4FZ4mWAiDgk6UD8eu6VQyR/JwIeiYjBlsPsIZkVtNDXSb7d356+R6lxTadvYhmXvoeNErnZU9Pdd5F09RR6ZWJIkv9/55G07pB0fEQ8m94/+kuSe1kD2Qx8VtIbImJXes5JEfFt8qb7lrSvlCQA7hoyK9XjJN07syCZFlnSm4vUe5TkW1kfEbGbZP79kkcqKVmE5ZP0beqfRvEPGque3P2kDpKBANcASDpJ0nqASKZLX0PyZWAnyWfvqvT1l0t6gmTSyGeAr+ROLOk/gE8BV0jqlPSmtFVxBbA6fc/N9O0+HDK3CMxKEBH7JV0K3CzpN0n+dj4DPFJQ9R6S2V+LnWNlYVk6guTlvKJTJbWTfPP/JXBLRHwl7/hcYPmwL8RGJCK+T3LDPr/s3f3UfYZkxtHc/seBjxepdxPJ+tvFzjGln/KNwFsHiXX8QMfzefZRszKT9K/Asoj4aQl1zwS+EBFnl1D3tcDXIuK8MoRp9gp3DZmV37UkN40HJOlKkqGlf1nieScDHx1BXGZFuUVgZlbn3CIwM6tzTgRmZnXOicDMrM45EZiZ1TknAjOzOvf/AUznnHmMrp4eAAAAAElFTkSuQmCC\n",
"text/plain": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"ax = lc_1d.plot(marker=\"o\", label=\"3D\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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`.\n",
"\n",
"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`. \n",
"\n",
"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](modeling.ipynb)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Fit the datasets"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/adonath/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/astropy/units/quantity.py:477: RuntimeWarning: overflow encountered in exp\n",
" result = super().__array_ufunc__(function, method, *arrays, **kwargs)\n",
"/Users/adonath/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/astropy/units/quantity.py:477: RuntimeWarning: invalid value encountered in subtract\n",
" result = super().__array_ufunc__(function, method, *arrays, **kwargs)\n"
]
}
],
"source": [
"# Define the model:\n",
"spectral_model1 = PowerLawSpectralModel(\n",
" index=2.0, amplitude=\"1e-12 cm-2 s-1 TeV-1\", reference=\"1 TeV\"\n",
")\n",
"temporal_model1 = ExpDecayTemporalModel(t0=\"10 h\", t_ref=gti_t0.mjd * u.d)\n",
"\n",
"model = SkyModel(\n",
" spectral_model=spectral_model1,\n",
" temporal_model=temporal_model1,\n",
" name=\"model-test\",\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"Table length=5\n",
"
\n",
"
name
value
unit
min
max
frozen
error
\n",
"
str9
float64
str14
float64
float64
bool
int64
\n",
"
index
2.0000e+00
nan
nan
False
0.000e+00
\n",
"
amplitude
1.0000e-12
cm-2 s-1 TeV-1
nan
nan
False
0.000e+00
\n",
"
reference
1.0000e+00
TeV
nan
nan
True
0.000e+00
\n",
"
t0
4.1667e-01
d
nan
nan
False
0.000e+00
\n",
"
t_ref
5.8909e+04
d
nan
nan
True
0.000e+00
\n",
"
"
],
"text/plain": [
"
\n",
" name value unit min max frozen error \n",
" str9 float64 str14 float64 float64 bool int64 \n",
"--------- ---------- -------------- ------- ------- ------ ---------\n",
" index 2.0000e+00 nan nan False 0.000e+00\n",
"amplitude 1.0000e-12 cm-2 s-1 TeV-1 nan nan False 0.000e+00\n",
"reference 1.0000e+00 TeV nan nan True 0.000e+00\n",
" t0 4.1667e-01 d nan nan False 0.000e+00\n",
" t_ref 5.8909e+04 d nan nan True 0.000e+00"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"model.parameters.to_table()"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"datasets.models = model"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"CPU times: user 6.64 s, sys: 86.5 ms, total: 6.72 s\n",
"Wall time: 7.02 s\n"
]
}
],
"source": [
"%%time\n",
"# Do a joint fit\n",
"fit = Fit(datasets)\n",
"result = fit.run()"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"Table length=5\n",
"
\n",
"
name
value
unit
min
max
frozen
error
\n",
"
str9
float64
str14
float64
float64
bool
float64
\n",
"
index
3.0062e+00
nan
nan
False
3.211e-02
\n",
"
amplitude
9.6409e-12
cm-2 s-1 TeV-1
nan
nan
False
3.339e-13
\n",
"
reference
1.0000e+00
TeV
nan
nan
True
0.000e+00
\n",
"
t0
2.5751e-01
d
nan
nan
False
8.796e-03
\n",
"
t_ref
5.8909e+04
d
nan
nan
True
0.000e+00
\n",
"
"
],
"text/plain": [
"
\n",
" name value unit min max frozen error \n",
" str9 float64 str14 float64 float64 bool float64 \n",
"--------- ---------- -------------- ------- ------- ------ ---------\n",
" index 3.0062e+00 nan nan False 3.211e-02\n",
"amplitude 9.6409e-12 cm-2 s-1 TeV-1 nan nan False 3.339e-13\n",
"reference 1.0000e+00 TeV nan nan True 0.000e+00\n",
" t0 2.5751e-01 d nan nan False 8.796e-03\n",
" t_ref 5.8909e+04 d nan nan True 0.000e+00"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"result.parameters.to_table()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We see that the fitted parameters match well with the simulated ones!"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exercises\n",
"\n",
"1. Re-do the analysis with `MapDataset` instead of `SpectralDataset`\n",
"2. Model the flare of PKS 2155-304 which you obtained using the [light curve flare tutorial](light_curve_flare.ipynb). Use a combination of a Gaussian and Exponential flare profiles, and fit using `scipy.optimize.curve_fit`\n",
"3. Do a joint fitting of the datasets."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.0"
},
"nbsphinx": {
"orphan": true
}
},
"nbformat": 4,
"nbformat_minor": 4
}