{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\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.17?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/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": [], "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\", 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": [], "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", "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", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
namevalueunitminmaxfrozenerror
str9float64str14float64float64boolint64
index3.000e+00nannanFalse0.000e+00
amplitude1.000e-11cm-2 s-1 TeV-1nannanFalse0.000e+00
reference1.000e+00TeVnannanTrue0.000e+00
t02.500e-01dnannanFalse0.000e+00
t_ref5.891e+04dnannanTrue0.000e+00
" ], "text/plain": [ "\n", " name value unit min max frozen error \n", " str9 float64 str14 float64 float64 bool int64 \n", "--------- --------- -------------- ------- ------- ------ ---------\n", " index 3.000e+00 nan nan False 0.000e+00\n", "amplitude 1.000e-11 cm-2 s-1 TeV-1 nan nan False 0.000e+00\n", "reference 1.000e+00 TeV nan nan True 0.000e+00\n", " t0 2.500e-01 d nan nan False 0.000e+00\n", " t_ref 5.891e+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.edges,\n", " e_true=energy_axis_true.edges,\n", " region=on_region,\n", " name=\"empty\",\n", ")\n", "\n", "for i in range(n_obs):\n", " obs = Observation.create(\n", " pointing=pointing,\n", " livetime=lvtm[i],\n", " tstart=tstart[i],\n", " irfs=irfs,\n", " reference_time=gti_t0,\n", " )\n", " empty_i = empty.copy(name=f\"dataset_{i}\")\n", " maker = SpectrumDatasetMaker(selection=[\"aeff\", \"background\", \"edisp\"])\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", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
namelivetimen_onbackgroundexcesssignificancebackground_rategamma_rate
s1 / s1 / s
str9float64int64float64float64float64float64float64
stacked3299.99999999999979620.30377224384548775.696227756154565.492940017483620.0061526582557107530.2350594629564105
dataset_11500.03439.228987383566126333.7710126164338442.574750955469160.0061526582557107510.2225140084109559
dataset_21560.02879.598146878908771277.401853121091237.357701735706510.0061526582557107510.17782170071864822
dataset_32400.032414.766379813705804309.2336201862941637.1861511610167240.0061526582557107510.12884734174428925
dataset_42399.999999999996419114.766379813705779176.2336201862942525.008446122692540.006152658255710750.07343067507762271
dataset_53000.018218.457974767132253163.5420252328677622.4929689593655480.0061526582557107510.05451400841095592
dataset_62400.03814.76637981370580423.2336201862941965.0369120121051790.0061526582557107510.009680675077622582
dataset_73120.05119.19629375781754331.8037062421824576.0048018205659060.0061526582557107510.010193495590443096
dataset_82580.02615.87385829973373910.1261417002662612.32501787631360470.0061526582557107510.00392486112413421
dataset_92820.02917.3504962811043211.6495037188956822.5483575723723540.0061526582557107520.00413102968755166
" ], "text/plain": [ "\n", " name livetime n_on ... background_rate gamma_rate \n", " s ... 1 / s 1 / s \n", " str9 float64 int64 ... float64 float64 \n", "--------- ------------------ ----- ... -------------------- --------------------\n", " stacked 3299.999999999999 796 ... 0.006152658255710753 0.2350594629564105\n", "dataset_1 1500.0 343 ... 0.006152658255710751 0.2225140084109559\n", "dataset_2 1560.0 287 ... 0.006152658255710751 0.17782170071864822\n", "dataset_3 2400.0 324 ... 0.006152658255710751 0.12884734174428925\n", "dataset_4 2399.9999999999964 191 ... 0.00615265825571075 0.07343067507762271\n", "dataset_5 3000.0 182 ... 0.006152658255710751 0.05451400841095592\n", "dataset_6 2400.0 38 ... 0.006152658255710751 0.009680675077622582\n", "dataset_7 3120.0 51 ... 0.006152658255710751 0.010193495590443096\n", "dataset_8 2580.0 26 ... 0.006152658255710751 0.00392486112413421\n", "dataset_9 2820.0 29 ... 0.006152658255710752 0.00413102968755166" ] }, "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 1.16 s, sys: 16.6 ms, total: 1.17 s\n", "Wall time: 1.2 s\n" ] } ], "source": [ "%%time\n", "lc_maker_1d = LightCurveEstimator(\n", " energy_range=[energy_axis.edges[0], energy_axis.edges[-1]],\n", " source=\"model-fit\",\n", " reoptimize=False,\n", ")\n", "lc_1d = lc_maker_1d.run(datasets)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAERCAYAAAB2CKBkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAdzklEQVR4nO3dcZQdZZnn8e+P2Eo7qI2mFegQw2CMogjhtEDM7B4ExkBkITo6wFEU9UzARWUcJ1ni7uqwzqzZwXEGRAhRUVkdRGOmh8EAMgYExhDo0KEDBjCKE9KJEsEOZO0xIXn2j6rW2ze3763uvnVv367f55x7uuqtt249lU7f5771Vr2vIgIzMyuug5odgJmZNZcTgZlZwTkRmJkVnBOBmVnBORGYmRWcE4GZWcG1ZCKQdL2kpyQ9XKf3u03SoKRbyso/ImmLpJA0vR7HMjObbFoyEQBfA86o4/tdAVxQofzfgNOBf6/jsczMJpWWTAQRcTfwTGmZpKPTb/YbJN0j6XVjeL8fAM9VKO+LiJ9POGAzs0nsBc0OoI5WAhdHxE8knQRcA5za5JjMzCa9KZEIJB0CvAX4jqTh4hel294J/K8Kuw1ExILGRGhmNnlNiURAcolrMCKOL98QEauB1Y0PycysNbRkH0G5iHgWeELSuwGUOK7JYZmZtYSWTASSbgTWAXMkbZP0IeA9wIckPQQ8Apwzhve7B/gOcFr6fgvS8o9J2gbMAPolfbne52Jm1mzyMNRmZsXWki0CMzOrn5brLJ4+fXrMmjWr2WGYmbWUDRs2/CoiOitta7lEMGvWLHp7e5sdhplZS5E06ggJvjRkZlZwuScCSdMk9ZUP6JZuk6Sr0oHd+iWdkHc8ZmY2UiNaBJcCm0fZdiYwO30tBq5tQDxmZlYi10QgaQbwdmC0++/PAW6IxH1Ah6TD84zJzMxGyrtF8A/AUmD/KNu7gCdL1relZSNIWiypV1Lvzp076x+lmVmB5ZYIJJ0FPBURG6pVq1B2wBNuEbEyIrojoruzs+LdT2ZmNk55tgjmA2dL+jnwLeBUSd8oq7MNOLJkfQawPceYzMysTG6JICKWRcSMiJgFnAesjYj3llW7GXhfevfQycCuiNiRV0xmZnaghj9QJuligIhYAawBFgJbgN8AH2h0POdetw6Amy6a1+hDm5lNCg1JBBFxF3BXuryipDyASxoRg5mZVeYni83MCs6JwMys4AqdCHr6BujbOsj6J55h/vK19PQNNDskM7OGK2wi6OkbYNnqTezZlzzrNjA4xLLVm5wMzKxwWm4Y6okYvkMIoG/r4O+SwLChvftYuqqfG+/fCvhOIjMrhsK2CMqTQK3yc69bNyKRmJlNFYVqEZR+w5+/fC0Dg0MH1OnqaHdLwMwKpbAtgiUL5tDeNm1EWXvbNJYsmNOkiMzMmqNQLYJSi+Ymg5wuXdXPnn376epoZ8mCOb8rNzMrisImAkiSgTuGzazoCntpyMzMEk4EZmYFV+hLQ+BLQmZmbhFk4KEozGwqcyKowUNRmNlU50RQw9JV/Qzt3TeibGjvPq64/bEmRWRmVl95Tl5/sKT7JT0k6RFJl1eoc4qkXZI2pq9P5RXPeI025MT2Ck8lm5m1ojw7i38LnBoRuyW1AfdKujUi7iurd09EnJVjHBPS1dFecSiKIzramxCNmVn95Tl5fUTE7nS1LX1FXsfLi4eiMLOpLtc+AknTJG0EngLuiIj1FarNSy8f3SrpDaO8z2JJvZJ6d+7cmWfIB1g0t4vPvvNYXjgt+afq6mjns+881kNRmNmUoWT++JwPInUA/wR8NCIeLil/KbA/vXy0ELgyImZXe6/u7u7o7e3NN+AKhoeg9nMHZtaKJG2IiO5K2xpy11BEDAJ3AWeUlT87fPkoItYAbZKmNyImMzNL5HnXUGfaEkBSO3A68GhZncMkKV0+MY3n6bxiMjOzA+V519DhwNclTSP5gP92RNwi6WKAiFgBvAv4sKTngSHgvGjEtSozM/ud3BJBRPQDcyuUryhZvhq4Oq8YzMysNj9ZbGZWcIUffTQr3y1kZlOVWwRmZgXnRGBmVnBOBGZmBedEYGZWcE4EZmYF50RgZlZwTgRmZgXnRGBmVnBOBGZmBedEYGZWcE4EZmYF50RgZlZwTgRmZgXnRGBmVnBOBGZmBZfnnMUHS7pf0kOSHpF0eYU6knSVpC2S+iWdkFc8ZmZWWZ4T0/wWODUidktqA+6VdGtE3FdS50xgdvo6Cbg2/TklnHvdOsCT2pjZ5JZbiyASu9PVtvRVPjH9OcANad37gA5Jh+cVk5mZHSjXPgJJ0yRtBJ4C7oiI9WVVuoAnS9a3pWXl77NYUq+k3p07d+YXsJlZAeWaCCJiX0QcD8wATpT0xrIqqrRbhfdZGRHdEdHd2dmZR6hmZoXVkLuGImIQuAs4o2zTNuDIkvUZwPZGxGRmZok87xrqlNSRLrcDpwOPllW7GXhfevfQycCuiNiRV0xmZnagPFsEhwN3SuoHHiDpI7hF0sWSLk7rrAF+BmwBvgT81xzjaaievgH6tg6y/olnmL98LT19A80OycysotxuH42IfmBuhfIVJcsBXJJXDM3S0zfAstWb2LNvPwADg0MsW70JgEVzD+gLNzNrKj9ZnIOlq/oZ2rtvRNnQ3n1ccftjTYrIzGx0TgQ5GG4JlNs+ONTgSMzManMiyEFXR3vF8iNGKTczayYnghwsWTCH9rZpI8ra26axZMGcJkVkZja6PMcaKqzhDuGlq/rZs28/XR3tLFkwxx3FZjYpORHkZNHcLm68fyvgQefMbHLzpSEzs4JzIjAzKzgnAjOzgnMiMDMrOHcW58idxGbWCjIlAkmvBOYDRwBDwMNAb0RUfoTWzMxaRtVEIOmtwGXAy4E+kpnGDgYWAUdLWgX8XUQ8m3egZmaWj1otgoXAn0XE1vINkl4AnAX8MfDdHGIzM7MGqJoIImJJlW3PAz11j8jMzBpq3HcNSfpAPQMxM7PmmMjto5dX2yjpSEl3Stos6RFJl1aoc4qkXZI2pq9PTSAeMzMbh1qdxf2jbQJeVeO9nwc+EREPSnoJsEHSHRHx47J690TEWdnCNTOzeqvVWfwqYAHw67JyAT+qtmM6Cf2OdPk5SZuBLqA8EZiZWRPVSgS3AIdExMbyDZLuynoQSbNI5i9eX2HzPEkPAduBv4yIRyrsvxhYDDBz5syshzUzswyUzB+f4wGkQ4AfAn8TEavLtr0U2B8RuyUtBK6MiNnV3q+7uzt6e3vzC9jMbAqStCEiuittG3NncfrtPGvdNpJnDL5ZngQAIuLZiNidLq8B2iRNH2tMZmY2fuO5a+jiLJUkCfgKsDkiPj9KncPSekg6MY3n6XHEZGZm4zSeQeeUsd584AJgk6ThPoZPAjMBImIF8C7gw5KeJxnD6LzI+1qVmZmNMJ5E8F+yVIqIe6mRNCLiauDqccRgZmZ1kunSkKRLJb00vYxzuaQHJb0t59jMzKwBsvYRfDAdYfRtQCfwAWB5blGZmVnDZE0Ew5d4FgJfjYiHyN5XYGZmk1jWRLBB0vdJEsHt6ZARnpTGzGwKyNpZ/CHgeOBnEfEbSa8guTxkZmYtLlMiSKekfLBk/Wl8v7+Z2ZQwkWGozcxsCnAiMDMrOCcCM7OCq5oIJB0r6T5JT0paKenQkm335x+emZnlrVaL4Frgr4BjgceBeyUdnW5ryzEuG6Nzr1vHudeta3YYZtaCat01dEhE3JYuf07SBuA2SRcAHhzOzGwKqJUIJOllEbELICLulPQnJHMMvDz36MzMLHe1Lg39H+D1pQUR0Q+cBhww0Yw1R0/fAH1bB1n/xDPMX76Wnr6BZodkZi2kaosgIv6xvEzSYRGxFfiz3KKyzHr6Bli2ehN79iUjfgwMDrFs9SYAFs3tamZoZtYixnP76Jq6R2Hjcu5161i6qp+hvftGlA/t3cfSVf1NisrMWs14EoFHHZ1EhlsCWcvNzMqNJxF8KUslSUdKulPSZkmPSLq0Qh1JukrSFkn9kk4YRzyFddNF8+jqaK+4bbRyM7NymROBpEMlvQm4T9IJGT60nwc+ERGvB04GLpF0TFmdM4HZ6WsxyXMLNgZLFsyhvW3aiLL2tmksWTCnSRGZWavJNPqopM8AFwI/5ffPDwRw6mj7RMQOYEe6/JykzUAX8OOSaucAN6QT1t8nqUPS4em+lsFwh/DSVf3s2befro52liyY445iM8ss63wEfwocHRF7xnMQSbOAucD6sk1dwJMl69vSshGJQNJikhYDM2fOHE8IU9qiuV3ceP9WILlcZGY2FlkvDT0MdIznAJIOIXkA7c/TeY9HbK6wywFPLEfEyojojojuzs7O8YRhZmajyNoi+CzQJ+lh4LfDhRFxdrWdJLWRJIFvRkSlB9C2AUeWrM8AtmeMyUq4JWBm45U1EXyd5CnjTWScq1iSgK8AmyPi86NUuxn4iKRvAScBu9w/YGbWWFkTwa8i4qoxvvd84AJgk6SNadkngZkAEbGC5OG0hcAW4Dd4HmQzs4bLmgg2SPosyTf40ktDD462Q0TcS42Hz9K7hS7JGIOZmeUgayKYm/48uaSs6u2jZmbWGjIlgoh4a96BmJlZc2S6fVTS/5bUUbJ+qKS/zi8sMzNrlKzPEZwZEYPDKxHxa5JOXjMza3FZE8E0SS8aXpHUDryoSn0zM2sRWTuLvwH8QNJXSTqJP0jybIGZmbW4rJ3FfyupHzid5JbQz0TE7blGZmZmDVE1EUhSeq8/EXEbcFu1OmZm1npq9RHcKemjkkYM+SnphZJOlfR14P35hWdmZnmrdWnoDJL+gBslHQUMAu0kCeT7wN9HxMYq+5uZ2SRXNRFExH8A1wDXpCOJTgeGSm8lNTOz1pb1riEiYi9lE8aYmVnrG8/k9WZmNoU4EZiZFVzWsYaOqVB2St2jMTOzhsvaIvi2pP+mRLukL5BMX2lmZi0uayI4iWRu4R8BD5DMKzy/2g6Srpf0VDrPcaXtp0jaJWlj+vrUWAI3M7P6yHrX0F5giOQZgoOBJyKi1tzFXwOuBm6oUueeiDgrYwxmZpaDrC2CB0gSwZuBPwLOl7Sq2g4RcTfwzMTCMzOzvGVNBB+KiE9FxN6I+EVEnAP8cx2OP0/SQ5JulfSG0SpJWiypV1Lvzp0763BYMzMblvXS0FPl4w0BP5zgsR8EXh0RuyUtBHqA2ZUqRsRKYCVAd3e3B7gzM6ujrIngeyTzEIikj+Ao4DFg1G/xtUTEsyXLayRdI2l6RPxqvO9pZmZjl3U+gmNL1yWdAFw0kQNLOgz4ZUSEpBNJLlM9PZH3NDOzscs81lCpiHhQ0pur1ZF0I3AKMF3SNuDTQFu6/wrgXcCHJT1P0hF9nuc1MDNrvEyJQNJflKweBJwAVO21jYjza2y/muT2UjMza6KsLYKXlCw/T9Jn8N36h2NmZo2WtY/g8rwDMTOz5qg1Z/G/kNwtVFFEnF33iMzMrKFqtQg+15AozMysaWolgiciYmtDIjEzs6aoNcREz/CCJHcOm5lNQbUSgUqW/zDPQGxqOve6dZx73bpmh2FmVdRKBDHKspmZTRG1EsFxkp6V9BzwpnT5WUnPSXq2xr5WcD19A/RtHWT9E88wf/laevoGmh2SmVVQtbM4IqY1KhCbWnr6Bli2ehN79iXzFw0MDrFs9SYAFs3tamZoZlYm63wEZpmde906lq7qZ2jvvhHlQ3v3sXRVv/sMzCYZJwLLxXBLIGu5mTWPE4HV3U0XzaOro73itq6Odm66aF6DIzKzapwILBdLFsyhvW1kF1N72zSWLJjTpIjMbDTjmo/ArJbhDuGlq/rZs28/XR3tLFkwxx3FZpOQE4HlZtHcLn/wm7WA3C4NSbpe0lOSHh5luyRdJWmLpP50+kszM2uwPPsIvgacUWX7mcDs9LUYuDbHWMzMbBS5JYKIuBt4pkqVc4AbInEf0CHp8LziMTOzypp511AX8GTJ+ra0zMzMGqiZiUAVyioObCdpsaReSb07d+7MOSwzs2JpZiLYBhxZsj4D2F6pYkSsjIjuiOju7OxsSHBmZkXRzERwM/C+9O6hk4FdEbGjifFYi/PcB2bjk9tzBJJuBE4BpkvaBnwaaAOIiBXAGmAhsAX4DfCBvGKxqW94yOs9+/Yzf/laP7xmNga5JYKIOL/G9gAuyev4Vhwe8tpsYvxksbWs4ctAwy2BUsNDXt94/1YPcmdWgweds5bnIa/NJsYtAmtZw9/05y9fy8Dg0AHbPeS1WTZuEVjL85DXZhPjFoG1vOEO4Stuf4ztg0Mc4SGvzcbEicCmBA95bTZ+vjRkZlZwTgRmZgXnRGBmVnBOBGZmBedEYGZWcE4EZmYF50RgZlZwTgRmZgXnRGBmVnBOBGZmBedEYGZWcLkmAklnSHpM0hZJl1XYfoqkXZI2pq9P5RmPmZkdKM85i6cBXwT+GNgGPCDp5oj4cVnVeyLirLziMDOz6vJsEZwIbImIn0XEHuBbwDk5Hs/MzMYhz0TQBTxZsr4tLSs3T9JDkm6V9IZKbyRpsaReSb07d+7MI1Yzs8LKMxGoQlmUrT8IvDoijgO+APRUeqOIWBkR3RHR3dnZWecwzcyKLc9EsA04smR9BrC9tEJEPBsRu9PlNUCbpOk5xmRmZmXyTAQPALMlHSXphcB5wM2lFSQdJknp8olpPE/nGJOZmZXJ7a6hiHhe0keA24FpwPUR8Yiki9PtK4B3AR+W9DwwBJwXEeWXj8zMLEdqtc/d7u7u6O3tbXYYZmYtRdKGiOiutM1PFpuZFZwTgZlZwTkRmJkVnBOBmVnBORGYmRWcE4GZWcE5EZiZFZwTgZlZweX2ZLGZWSP09A1wxe2PsX1wiCM62lmyYA6L5lYa6Li1jtVITgRm1rJ6+gZYtnoTQ3v3ATAwOMSy1ZsA6v4B3chjNZqHmDCzlnTudevo2zrInn37D9j2wmkH8fjfnFmXYwyrdqy5Mzu46aJ5Ez5enjzEhJlNSZU+mKuVt8qxGs2XhsysJd100TzmL1/LwODQAdu6Otrrdoxh1Y412VsDtbhFYGYta8mCObS3TRtR1t42jSUL5kyKY/X0DTB/+VqOuux7zF++lp6+gZrHGc8+E+UWgZm1rOFO2kbcyTPWY42nc7lZHdLuLDYzq6PhDuaxdGRn2WeiHdJN6yyWdIakxyRtkXRZhe2SdFW6vV/SCXnGY2bWKOPpXG5Wh3RuiUDSNOCLwJnAMcD5ko4pq3YmMDt9LQauzSseM7NGuOmiedx00bxRO6wrlWfZJ88O6TxbBCcCWyLiZxGxB/gWcE5ZnXOAGyJxH9Ah6fAcYzIza4jxdC43svO7VJ6dxV3AkyXr24CTMtTpAnaUVpK0mKTFwMyZM+seqJlZvY2nI7uRnd+l8kwEqlBW3jOdpQ4RsRJYCUln8cRDMzPL36K5XWP+EB/PPhOV56WhbcCRJeszgO3jqGNmZjnKMxE8AMyWdJSkFwLnATeX1bkZeF9699DJwK6I2FH+RmZmlp/cLg1FxPOSPgLcDkwDro+IRyRdnG5fAawBFgJbgN8AH8grHjMzqyzXJ4sjYg3Jh31p2YqS5QAuyTMGMzOrzmMNmZkVnBOBmVnBORGYmRVcyw06J2kn8O/AdOBXTQ6nGYp43kU8ZyjmeRfxnKEx5/3qiOistKHlEsEwSb2jjaQ3lRXxvIt4zlDM8y7iOUPzz9uXhszMCs6JwMys4Fo5EaxsdgBNUsTzLuI5QzHPu4jnDE0+75btIzAzs/po5RaBmZnVgROBmVnBtUwikPRuSY9I2i9p1Nusas2T3GokvVzSHZJ+kv48dJR6H0//fR6WdKOkgxsda72M4Zw7JK2S9KikzZLym8uvAbKed1p3mqQ+Sbc0MsZ6y3LOko6UdGf6O35E0qXNiHWiJvMc7i2TCICHgXcCd49WIeM8ya3mMuAHETEb+EG6PoKkLuBjQHdEvJFktNfzGhplfdU859SVwG0R8TrgOGBzg+LLS9bzBriU1j9fyHbOzwOfiIjXAycDl7Ta3/Vkn8O9ZRJBRGyOiMdqVMsyT3KrOQf4err8dWDRKPVeALRLegHwYlp7gp+a5yzppcB/Br4CEBF7ImKwYRHmI9PvWtIM4O3AlxsUV55qnnNE7IiIB9Pl50gSYGOn8Jq4ST2He8skgoxGmwO5lb1qeLKe9OcryytExADwOWAryXzPuyLi+w2Nsr5qnjPwh8BO4KvpJZIvS/qDRgaZgyznDfAPwFJgf6MCy1HWcwZA0ixgLrA+98jqK8tnU9M+v3Kdj2CsJP0rcFiFTf89Iv45y1tUKJv098dWO++M+x9K8m3iKGAQ+I6k90bEN+oXZX1N9JxJ/u+eAHw0ItZLupLkssL/rFOIuajD7/os4KmI2CDplHrGlpc6/K6H3+cQ4LvAn0fEs/WIrYHqNod7HiZVIoiI0yf4Fi05B3K185b0S0mHR8SOtJn4VIVqpwNPRMTOdJ/VwFuASZsI6nDO24BtETH8zXAV1a+pTwp1OO/5wNmSFgIHAy+V9I2IeG9OIU9YHc4ZSW0kSeCbEbE6p1DzNKnncJ9ql4ayzJPcam4G3p8uvx+o1DLaCpws6cWSBJxGa3ck1jzniPgF8KSkOWnRacCPGxNebrKc97KImBERs0j+f6+dzEkgg5rnnP6f/gqwOSI+38DY6mlyz+EeES3xAt5BkjF/C/wSuD0tPwJYU1JvIfA48FOSS0pNj32C5/0KkrspfpL+fPko53058CjJ3VX/F3hRs2NvwDkfD/QC/UAPcGizY2/EeZfUPwW4pdlx533OwB+RXCLpBzamr4XNjn0c53rAZxNwMXBxuiySO4t+CmwiuQuwIbF5iAkzs4KbapeGzMxsjJwIzMwKzonAzKzgnAjMzArOicDMbAIk/ZWkAUkb09fCUepVHBhS0nGS1knaJOlf0uFTkPSKdLC93ZKuHmNMX5C0O2t9JwKb0tI/puE/0F+U/cH+KKdjzpX05XT5Qkkh6bSS7e9Iy96Vrt+ldERdST9PPxA2SfqxpL+W9KJ0W6ek2/KI2bKRdIqkr1XY9PcRcXz6WlNhv2oDQ34ZuCwijgX+CViSlv8HyZPyfznGGLuBjrHs40RgU1pEPD38BwqsYOQf7FtyOuwngS+UrG8Czi9ZPw94qMr+b00/FE4kGVNpJUAkT47vkDS/vuFag4w2MOQcfj+q8h3AnwBExP+LiHtJEsIIkt6WtiIelPSddPiN4VFOryAZiyozJwIrrOGmc/ot74eSvi3pcUnLJb1H0v3pN/Oj03qdkr4r6YH0dcAHsqSXAG+KiNIP+nuAEyW1pX+wryF5KKqqiNhN8sDRIkkvT4t7gPdM6MQtDx9RMofA9aowp0JUHxjyYeDsdPndjBxm4gCSpgP/Azg9Ik4geajyL4bjAG6OMT6R7ERgljiOZIz/Y4ELgNdGxIkkzfaPpnWuJGlRvJnkW1ulYaC7Sf6wSwXwr8ACksEBMw97Esngak+QjFEPyR/9f8q6v9WHpPWSNpL8zs8uuby4gGTegKNJnnTfAfxdhf1LB4Y8AvgDScNDg3yQZI6FDcBLgD01wjmZZE6Df0tjej/waklHkCSSL1TbuZJJNeicWRM9MPwtStJPgeFva5uAt6bLpwPHJEPfAMmAby+JZIz8YYeTDI9d7lsk14hfBnyC5PJRVqWjUj5F8kFiDRQRJ0HSegQujIgLK9WT9CWg0qxxow4MGRGPAm9Ly19LMtdENQLuiIjzRxRKbydpbW5J/4++WNKWiHhNrfNzIjBL/LZkeX/J+n5+/3dyEDAvIoaqvM8QyaigI0TE/ZLeCAxFxOMlyaSq9FLTLJIxakjfu9rxrcGGR09NV9/BgS1CKBkYkuT3dxpJ6w5Jr4yIpyQdRHLJZ0WNQ94HfFHSayJiS/qeMyLie5QM9y1pd5YkAL40ZDYW3ye5BguApOMr1NlM8q2skmWMoSWQ9idcA/RExK/T4tdS+YPGmudv076kfpLW48cBJB0haQ1AJMOlrwIeJGllHkR6EwDJtJWPkwwauR346vAbS/o58HngQknbJB2TtiouBG5Mj3kf8LqJnIBbBGbZfYzkm1g/yd/O3SSdub8TEY9KelmFS0ZExK2jvO8LGNkiuVNJk+EgktsJP1Oy7a3A9yZ2GjZeEXEXcFdZ2QWj1N1OMuLo8PqngU9XqHclSf9TpfeYNUr5WuDNNWI9pNr2Uh591KzOJH0ceC4ias4pnD4jsAV4Y0TsylD/buCckhaC2YT50pBZ/V3LyG/4FaUP/mwErsmYBDqBzzsJWL25RWBmVnBuEZiZFZwTgZlZwTkRmJkVnBOBmVnBORGYmRXc/wfNi/VcMNlMXAAAAABJRU5ErkJggg==\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": 15, "metadata": {}, "outputs": [], "source": [ "# Define the model:\n", "\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", "model = SkyModel(\n", " spectral_model=spectral_model1,\n", " temporal_model=temporal_model1,\n", " name=\"model-test\",\n", ")" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/html": [ "Table length=5\n", "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
namevalueunitminmaxfrozenerror
str9float64str14float64float64boolint64
index2.000e+00nannanFalse0.000e+00
amplitude1.000e-12cm-2 s-1 TeV-1nannanFalse0.000e+00
reference1.000e+00TeVnannanTrue0.000e+00
t04.167e-01dnannanFalse0.000e+00
t_ref5.891e+04dnannanTrue0.000e+00
" ], "text/plain": [ "\n", " name value unit min max frozen error \n", " str9 float64 str14 float64 float64 bool int64 \n", "--------- --------- -------------- ------- ------- ------ ---------\n", " index 2.000e+00 nan nan False 0.000e+00\n", "amplitude 1.000e-12 cm-2 s-1 TeV-1 nan nan False 0.000e+00\n", "reference 1.000e+00 TeV nan nan True 0.000e+00\n", " t0 4.167e-01 d nan nan False 0.000e+00\n", " t_ref 5.891e+04 d nan nan True 0.000e+00" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model.parameters.to_table()" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "for dataset in datasets:\n", " dataset.models = model" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 8.36 s, sys: 111 ms, total: 8.47 s\n", "Wall time: 9.03 s\n" ] } ], "source": [ "%%time\n", "# Do a joint fit\n", "fit = Fit(datasets)\n", "result = fit.optimize()" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/html": [ "Table length=5\n", "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
namevalueunitminmaxfrozenerror
str9float64str14float64float64boolint64
index2.983e+00nannanFalse0.000e+00
amplitude9.907e-12cm-2 s-1 TeV-1nannanFalse0.000e+00
reference1.000e+00TeVnannanTrue0.000e+00
t02.497e-01dnannanFalse0.000e+00
t_ref5.891e+04dnannanTrue0.000e+00
" ], "text/plain": [ "\n", " name value unit min max frozen error \n", " str9 float64 str14 float64 float64 bool int64 \n", "--------- --------- -------------- ------- ------- ------ ---------\n", " index 2.983e+00 nan nan False 0.000e+00\n", "amplitude 9.907e-12 cm-2 s-1 TeV-1 nan nan False 0.000e+00\n", "reference 1.000e+00 TeV nan nan True 0.000e+00\n", " t0 2.497e-01 d nan nan False 0.000e+00\n", " t_ref 5.891e+04 d nan nan True 0.000e+00" ] }, "execution_count": 19, "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 }