{ "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.19?urlpath=lab/tree/tutorials/analysis/time/light_curve_simulation.ipynb)\n", "- You may download all the notebooks in the documentation as a\n", "[tar file](../../../_downloads/notebooks-0.19.tar).\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": [ "# Simulating and fitting a time varying source\n", "\n", "## Prerequisites:\n", "\n", "- To understand how a single binned simulation works, please refer to [spectrum_simulation](../1D/spectrum_simulation.ipynb) [simulate_3d](../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", "- Optionally, 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": { "execution": { "iopub.execute_input": "2021-11-22T21:05:30.314867Z", "iopub.status.busy": "2021-11-22T21:05:30.313354Z", "iopub.status.idle": "2021-11-22T21:05:30.689080Z", "shell.execute_reply": "2021-11-22T21:05:30.689280Z" } }, "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\n", "from astropy.time import Time\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": { "execution": { "iopub.execute_input": "2021-11-22T21:05:30.691851Z", "iopub.status.busy": "2021-11-22T21:05:30.691458Z", "iopub.status.idle": "2021-11-22T21:05:31.009334Z", "shell.execute_reply": "2021-11-22T21:05:31.009526Z" } }, "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, RegionGeom\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": { "execution": { "iopub.execute_input": "2021-11-22T21:05:31.012052Z", "iopub.status.busy": "2021-11-22T21:05:31.011682Z", "iopub.status.idle": "2021-11-22T21:05:31.059776Z", "shell.execute_reply": "2021-11-22T21:05:31.059947Z" } }, "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": { "execution": { "iopub.execute_input": "2021-11-22T21:05:31.070907Z", "iopub.status.busy": "2021-11-22T21:05:31.070610Z", "iopub.status.idle": "2021-11-22T21:05:31.071746Z", "shell.execute_reply": "2021-11-22T21:05:31.071934Z" } }, "outputs": [], "source": [ "# Reconstructed and true energy axis\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", "geom = RegionGeom.create(\"galactic;circle(0, 0, 0.11)\", axes=[energy_axis])" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "execution": { "iopub.execute_input": "2021-11-22T21:05:31.073975Z", "iopub.status.busy": "2021-11-22T21:05:31.073571Z", "iopub.status.idle": "2021-11-22T21:05:31.075126Z", "shell.execute_reply": "2021-11-22T21:05:31.075362Z" } }, "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": { "execution": { "iopub.execute_input": "2021-11-22T21:05:31.079171Z", "iopub.status.busy": "2021-11-22T21:05:31.078875Z", "iopub.status.idle": "2021-11-22T21:05:31.086571Z", "shell.execute_reply": "2021-11-22T21:05:31.086793Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/adonath/software/mambaforge/envs/gammapy-dev/lib/python3.9/site-packages/astropy/units/quantity.py:486: RuntimeWarning: overflow encountered in exp\n", " result = super().__array_ufunc__(function, method, *arrays, **kwargs)\n", "/Users/adonath/software/mambaforge/envs/gammapy-dev/lib/python3.9/site-packages/astropy/units/quantity.py:486: 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": { "execution": { "iopub.execute_input": "2021-11-22T21:05:31.090804Z", "iopub.status.busy": "2021-11-22T21:05:31.090379Z", "iopub.status.idle": "2021-11-22T21:05:31.092402Z", "shell.execute_reply": "2021-11-22T21:05:31.092576Z" } }, "outputs": [ { "data": { "text/html": [ "
Table length=5\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
typenamevalueuniterrorminmaxfrozenlink
str8str9float64str14int64float64float64boolstr1
spectralindex3.0000e+000.000e+00nannanFalse
spectralamplitude1.0000e-11cm-2 s-1 TeV-10.000e+00nannanFalse
spectralreference1.0000e+00TeV0.000e+00nannanTrue
temporalt06.0000e+00h0.000e+00nannanFalse
temporalt_ref5.8909e+04d0.000e+00nannanTrue
" ], "text/plain": [ "\n", " type name value unit error min max frozen link\n", " str8 str9 float64 str14 int64 float64 float64 bool str1\n", "-------- --------- ---------- -------------- --------- ------- ------- ------ ----\n", "spectral index 3.0000e+00 0.000e+00 nan nan False \n", "spectral amplitude 1.0000e-11 cm-2 s-1 TeV-1 0.000e+00 nan nan False \n", "spectral reference 1.0000e+00 TeV 0.000e+00 nan nan True \n", "temporal t0 6.0000e+00 h 0.000e+00 nan nan False \n", "temporal t_ref 5.8909e+04 d 0.000e+00 nan nan True " ] }, "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": { "execution": { "iopub.execute_input": "2021-11-22T21:05:31.094619Z", "iopub.status.busy": "2021-11-22T21:05:31.094331Z", "iopub.status.idle": "2021-11-22T21:05:31.095404Z", "shell.execute_reply": "2021-11-22T21:05:31.095630Z" } }, "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": { "execution": { "iopub.execute_input": "2021-11-22T21:05:31.100461Z", "iopub.status.busy": "2021-11-22T21:05:31.100119Z", "iopub.status.idle": "2021-11-22T21:05:32.062086Z", "shell.execute_reply": "2021-11-22T21:05:32.062266Z" } }, "outputs": [], "source": [ "datasets = Datasets()\n", "\n", "empty = SpectrumDataset.create(\n", " geom=geom, energy_axis_true=energy_axis_true, 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": { "execution": { "iopub.execute_input": "2021-11-22T21:05:32.065323Z", "iopub.status.busy": "2021-11-22T21:05:32.064931Z", "iopub.status.idle": "2021-11-22T21:05:32.089276Z", "shell.execute_reply": "2021-11-22T21:05:32.089460Z" } }, "outputs": [ { "data": { "text/html": [ "
Table length=10\n", "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
namecountsbackgroundexcesssqrt_tsnprednpred_backgroundnpred_signalexposure_minexposure_maxlivetimeontimecounts_ratebackground_rateexcess_raten_binsn_fit_binsstat_typestat_sum
m2 sm2 sss1 / s1 / s1 / s
str9int64float64float64float64float64float64float64float64float64float64float64float64float64float64int64int64str4float64
stacked83820.318769454956055817.681213378906267.8120516609811720.31876973807811720.318769454956055nan216137904.016025275392.03299.9999999999993300.00.2539393939393940.00615720286513820.247782185872395999cashnan
dataset-13359.235804207223003325.76419579277741.88639768268303332.149612809878769.235804207223003322.9138086026557698244500.935563397284216297.3122691500.01500.00.223333333333333330.00615720280481533550.21717613052851899cash-2096.519112853783
dataset-23199.605236375511922309.3947636244880640.20011982055861293.48960258449289.605236375511922283.88436620898085102174280.972985927575584949.204761560.01560.00.204487179487179480.0061572028048153350.1983299766823641699cash-1999.167487932759
dataset-333114.777286731556805316.222713268443237.75884642088334321.783859418843814.777286731556805307.006572687287157191201.4969014211654746075.699632400.02400.00.137916666666666660.00615720280481533550.1317594638618513399cash-2093.961371303704
dataset-419114.777286731556783176.2227132684432125.003242551912294200.9861857996912414.777286731556783186.20889906813449157191201.4969011811654746075.6996142399.99999999999642399.99999999999640.079583333333333450.00615720280481533550.0734261305285181299cash-950.0232173832845
dataset-515618.471608414446006137.52839158555419.764438032194874182.999442729359418.471608414446006164.52783431491338196489001.8711267714568432594.6245383000.03000.00.0520.00615720280481533550.0458427971951846799cash-679.5591099561027
dataset-63714.77728673155680522.2227132684431944.844970711706384539.9779208281207614.77728673155680525.200634096563952157191201.4969014211654746075.699632400.02400.00.0154166666666666670.00615720280481533550.0092594638618513399cash-83.5158051456156
dataset-73419.21047275102384314.7895272489761573.040144950576299642.3048288120531219.21047275102384323.094356061029266204348561.9459718515151169898.409523120.03120.00.0108974358974358970.0061572028048153350.00474023309262056399cash-70.91103544138612
dataset-83115.88558323642356415.1144167635764363.350049281905792832.2499005026512415.88558323642356416.36431726622769168980541.60916912528852031.3771022580.02580.00.0120155038759689920.0061572028048153350.00585830107115365799cash-42.09535699583043
dataset-94117.36331190957924523.6366880904207554.81477484534861332.42183519709226417.36331190957924515.058523287513015184699661.7588591613694326638.9470652820.02820.00.0145390070921985820.00615720280481533550.00838180428738324799cash-96.29813771797772
" ], "text/plain": [ "\n", " name counts background excess sqrt_ts ... excess_rate n_bins n_fit_bins stat_type stat_sum \n", " ... 1 / s \n", " str9 int64 float64 float64 float64 ... float64 int64 int64 str4 float64 \n", "--------- ------ ------------------ ------------------ ------------------ ... -------------------- ------ ---------- --------- ------------------\n", " stacked 838 20.318769454956055 817.6812133789062 67.81205166098117 ... 0.2477821858723959 9 9 cash nan\n", "dataset-1 335 9.235804207223003 325.764195792777 41.88639768268303 ... 0.217176130528518 9 9 cash -2096.519112853783\n", "dataset-2 319 9.605236375511922 309.39476362448806 40.20011982055861 ... 0.19832997668236416 9 9 cash -1999.167487932759\n", "dataset-3 331 14.777286731556805 316.2227132684432 37.75884642088334 ... 0.13175946386185133 9 9 cash -2093.961371303704\n", "dataset-4 191 14.777286731556783 176.22271326844321 25.003242551912294 ... 0.07342613052851812 9 9 cash -950.0232173832845\n", "dataset-5 156 18.471608414446006 137.528391585554 19.764438032194874 ... 0.04584279719518467 9 9 cash -679.5591099561027\n", "dataset-6 37 14.777286731556805 22.222713268443194 4.8449707117063845 ... 0.00925946386185133 9 9 cash -83.5158051456156\n", "dataset-7 34 19.210472751023843 14.789527248976157 3.0401449505762996 ... 0.004740233092620563 9 9 cash -70.91103544138612\n", "dataset-8 31 15.885583236423564 15.114416763576436 3.3500492819057928 ... 0.005858301071153657 9 9 cash -42.09535699583043\n", "dataset-9 41 17.363311909579245 23.636688090420755 4.814774845348613 ... 0.008381804287383247 9 9 cash -96.29813771797772" ] }, "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": { "execution": { "iopub.execute_input": "2021-11-22T21:05:32.094052Z", "iopub.status.busy": "2021-11-22T21:05:32.093618Z", "iopub.status.idle": "2021-11-22T21:05:32.094887Z", "shell.execute_reply": "2021-11-22T21:05:32.095103Z" } }, "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": { "execution": { "iopub.execute_input": "2021-11-22T21:05:32.102311Z", "iopub.status.busy": "2021-11-22T21:05:32.101423Z", "iopub.status.idle": "2021-11-22T21:05:32.103599Z", "shell.execute_reply": "2021-11-22T21:05:32.104070Z" } }, "outputs": [], "source": [ "# Attach model to all datasets\n", "datasets.models = model_fit" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "execution": { "iopub.execute_input": "2021-11-22T21:05:32.114646Z", "iopub.status.busy": "2021-11-22T21:05:32.114333Z", "iopub.status.idle": "2021-11-22T21:05:33.576561Z", "shell.execute_reply": "2021-11-22T21:05:33.576749Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 1.45 s, sys: 18.5 ms, total: 1.47 s\n", "Wall time: 1.47 s\n" ] } ], "source": [ "%%time\n", "lc_maker_1d = LightCurveEstimator(\n", " energy_edges=[0.3, 10] * u.TeV,\n", " source=\"model-fit\",\n", " selection_optional=[\"ul\"],\n", ")\n", "lc_1d = lc_maker_1d.run(datasets)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "execution": { "iopub.execute_input": "2021-11-22T21:05:33.595949Z", "iopub.status.busy": "2021-11-22T21:05:33.595646Z", "iopub.status.idle": "2021-11-22T21:05:33.800689Z", "shell.execute_reply": "2021-11-22T21:05:33.800933Z" }, "nbsphinx-thumbnail": { "tooltip": "Simulate and fit a time decaying light curve of a source using the CTA 1DC response." } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbYAAAE6CAYAAAB6e3fTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAnd0lEQVR4nO3de7xVdZ3/8debm4hHIRFTQQRBM20g9aCRlpJhmpJdLNRpJMcGnamptGxsstEuTv20kKymUEktG7XL5KiJmml0IxVUvDeaph1zfiAFZYiifuaP79qy3R7w3PZel/1+Ph7rcfa67LPfZ+999md/1/qu71JEYGZmVhWD8g5gZmY2kFzYzMysUlzYzMysUlzYzMysUlzYzMysUlzYzMysUobkHcBgm222iQkTJuQdw8ysVJYtW/ZERIxpXO7CVgATJkxg6dKleccwMysVSY90t9y7Is3MrFJc2MzMrFJc2MzMrFJ8jM3MrMLWr19PV1cX69atyztKnw0fPpxx48YxdOjQHm3vwmZmVmFdXV1sueWWTJgwAUl5x+m1iGDVqlV0dXUxceLEHt3HuyLNzCps3bp1jB49upRFDUASo0eP7lWL04XNzKzielvUZi9YwuwFS5qUpvd6m9+FrWIeXbWWmfMWM+kT1zBz3mIeXbU270hm1sbWrVvHPvvsw9SpU9ljjz04/fTTAXjf+97HxIkTmTp1KrvuuivHHnssjz322IA8po+xlVzjt6rlXatZt/55AB5Y8SQHz1/M1HGjXlh/+QnTWxnPzNrcZpttxo033khHRwfr169n//3359BDDwXg7LPP5sgjjyQimD9/PjNmzODuu+9m2LBh/XpMt9gqplbUNjZvZrYpj65ay/Ku1dz88B8HZK+PJDo6OoDUQ3P9+vUv2bUoiZNOOontttuORYsW9evxwC220mtsgc2ct5gHVjwJwCDBpDEdbqWZ2SbV7/nZ1F6fvn6WPPfcc+y99948+OCDfOADH2Dffffl61//+ku222uvvbj//vs54ogj+vQ4NW6xVczCOdPYZdsOBktMGtPBwjnT8o5kZiXSjL0+gwcP5o477qCrq4tbbrmFu+++u9vtIqLfjwVusVXO+NEj+PHJB/Ro29q3NLfozNpb/WdAM/f6jBo1igMPPJBrr7222/W33347Bx10UL8fxy22NjXQ+9HNrBoWzpnG8KGpNAzEXp+VK1eyevVqAJ566iluuOEGdttttxdtExGce+65PP744xxyyCH9ejxwYWtLsxcs4eD5i1+0H/34i2/NOZWZFcH40SOYOm4U+07cmh+ffADjR4/o1+97/PHHmTFjBlOmTGHatGnMnDmTww8/HIBTTjnlhe7+t956KzfddFO/e0SCd0W2rcb95g+t/GtOScysyqZMmcLtt9/+kuUXXXRR0x7Tha0NXX7C9JfsR995zBY9uq+Py5lVX9n/v70rsk31pfekj8uZWRm4xdametN7ElJLrfH8luMvvrVXv8PMrBXcYrMe83E5s3IaqPPD8tLb/C5s1iOXnzCdXbbteGG+N8flzCw/w4cPZ9WqVaUtbrXrsQ0fPrzH9/GuyH6StDPwSWBkRByZLXs7cBiwLfC1iLg+v4QDZ+GcaRx/8a08tPKv7DxmC49qYlYC48aNo6uri5UrV+Ydpc9qV9DuKZW1ig8ESd8EDgdWRMRr6pYfAnwZGAxcEBFf6MHv+n6tsNUtewXwxYg4flP37ezsjKVLl/blTzAza1uSlkVEZ+Pydm+xXQR8FfhWbYGkwcDXgJlAF3CrpCtJRe7zDff/+4hYsYnff1r2u8zMrEXaurBFxM8kTWhYvA/wYEQ8BCDpMuCIiPg8qXX3spSuyfAFYFFE3LaRbeYCcwHGjx/ftz/AzMxewp1HXmos8Pu6+a5sWbckjZb0DWBPSZ/IFv8z8GbgSEkndne/iDgvIjojonPMmDEDFN3MzNq6xbYR6mbZRg9ERsQq4MSGZecC5w5wLjMz6wG32F6qC9ixbn4c8IecspiZWS+5sL3UrcAukiZKGgYcBVyZcyYzM+uhti5ski4FlgCvktQl6fiIeBb4IHAdcB/w3Yi4p0mPP0vSeWvWrGnGrzcza0ttfR5bUfg8NjOz3tvYeWxt3WIzM7PqcWEzM7NKcWEzM7NK8XlsOZI0C5g1efLkvKMUzqOr1nLw/MWsW/88u2ybLoQ6fvSIvGOZWQm480gBuPNIMnvBkhdu11/UFGD40EFMHTeq9JesN7OB484jViqNFzVtnDcz2xjvirTCqG+NzZy3mAdWPAmki5pOGtPh1pqZ9YhbbFZIC+dMY/jQ9PacNKbDFzU1sx5zi80KafzoEdz/2UPzjmFmJeQWW448pJaZ2cBzYctRRFwVEXNHjhyZdxQzs8pwYTMzs0pxYTMzs0pxYTMzs0pxYTMzs0pxYcuRe0WamQ08F7YcuVekmdnAc2EzM7NKcWEzM7NKcWEzM7NKcWEzM7NKcWEzM7NKKf3o/pIGAVOBHYCngHsi4v/nm8rMzPJS2sImaRLwL8CbgQeAlcBwYFdJa4EFwMURUdhLL0uaBcyaPHly3lHMzCpDEZF3hj6RdCnwdeDn0fBHSNoWOAb4U0RcnEe+3ujs7IylS5fmHcPMrFQkLYuIzsblpW2xRcTRm1i3ApjfujRmZlYUpS1sNVnrbD82HGO7G1ha5F2QZmbWPKUtbJJmAKcCWwO3AytIx9jeDkyS9H3gSxHx59xCmplZy5W2sAFvBf4hIh5tXCFpCHA4MBP4QauDWbXNXrAEgMtPmJ5zEjPrTpnPY/tid0UNICKejYgrIsJFzQbUo6vWsrxrNTc//EdmzlvMo6vW5h3JzBqUucW2XNJdwKXADyLC136xpqm10pZ3rWbd+nT49oEVT3Lw/MVMHTfKrTezAilzi20s8EXgDcD/SLpC0mxJm+ecyyqsVtQ2Nm9m+SttYYuI5yLiuog4DtgRuJDUceRhSd/JNZxVzuUnTOfyE6azy7YdLywbJNhl2w631swKprSFrV5EPAPcC9wH/BnYPd9EPeMraJfPwjnTGD40/dtMGtPBwjnTck5kZo1KO/IIgKTxwGzgaGAL4DLgsoi4L9dgveSRR8zMeq9yI49I+hXpONv3gLkR4cpgZmblLWzAJ4CfNY4TaWZm7a20x9giYnFEhKRdJf1E0t0AkqZIOi3vfGZmlo/SFrY655Nab+sBIuJO4KhcE5mZWW6qUNhGRMQtDcuezSWJmZnlrrSFLesRCfBEdtHRyJYfCTyeWzAzM8tVmTuPXAHsBXyQdLXs3SQ9BjwMvDfHXGZmlqMyFzYBRMRvgTdL2gIYFBF/yTeWmZnlqcyFbaykcxsXSgIgIj7U8kRmZpa7Mhe2p4BleYcwM7NiKXNhWxURF+cdoj8kzQJmTZ48Oe8oZmaVUdpekcAzeQfor4i4KiLmjhw5Mu8oZmaVUdrCFhGvyzuDWbPNXrDkhYucmlnPlLawmZmZdceFzaygHl21luVdq7n54T8yc95iHl21Nu9IZqVQ5s4jAGSjjnRFxNOSDgSmAN+KiNV55jLrrcZdjsu7VrNu/fMAPLDiSQ6ev5ip40YB+KrdZptQhRbbD4DnJE0GFgITgf/MN5JZ/9WK2sbmzax7pW+xAc9HxLOS3gHMj4ivSLo971BmvdXYCps5bzEPrHgSgEGCSWM63FIz64EqtNjWSzoamANcnS0bmmMeswGxcM40hg9N/6KTxnSwcM60nBOZlUMVWmzHAScCZ0bEw5ImApfknMms38aPHuFjamZ9oIjIO0Pb6+zsjKVLl+Ydw8ysVCQti4jOxuVV2BVpZmb2Ahc2MzOrFBc2MzOrlNIWNkmDJZ0g6bOS9mtYd1peuczMLF+lLWzAAuAAYBVwrqR5devemU8kMzPLW5kL2z4RcUxEzAf2BTok/ZekzQDlG83MzPJS5sI2rHYjIp6NiLnAHcCNQEdeoczMLF9lLmxLJR1SvyAiPgNcCEzIJVEvSZol6bw1a9bkHcXMrDJ8gnYB+ARtM7Pe29gJ2qUfUkvSYOAwUivthb8nIuZt7D5mZlZdpS9swFXAOuAuwNf1MDNrc1UobOMiYkreIczMrBjK3HmkZpGkg/MOYWZmxVCFFtuvgR9KGgSsJ53DFhGxVb6xzMwsD1UobF8CpgN3hbt4mpm1vSrsinwAuNtFzczMoBottseBn0paBDxdW+ju/mZm7akKhe3hbBpG3TBbZmbWnkpf2CLi03lnMDOz4ij9MTZJP5Y0qm7+FZKuyzGSmZnlqPSFDRgTEatrMxHxJ2Db/OKYmVmeqlDYnpM0vjYjaSfAPSTNzNpU6Y+xAZ8EfiFpcTb/RmBujnnMzCxHpS9sEXGtpL2A15FGHTkpIp7IOZaZmeWktIVN0oSI+B1AVsiublgvYGxEdOUQz8zMclLawgacnY0P+d/AMmAlMByYDMwADgJOB1zYzMzaSGkLW0S8W9LuwN8Cfw9sD6wF7gOuAc6MiHU5RjQzsxyUtrABRMS9pM4jZmZmQDW6+5uZmb3Ahc3MzCqltIVNUql3o5qZWXOUuTj8WlIXcC1wba3rv5mZtbfSFraI6MyGzzoUmC9pLPALYBGwOCKe3uQvMDOzSirtrkiAiHgkIr4REW8HXg9cBbwZ+LmkH7Uig6SdJS2U9P26Za+W9A1J35f0j63IYWZmSakLW72IWB8RN0bExyNiH3owXqSkb0paIenuhuWHSPqNpAclnfoyj/tQRBzfsOy+iDgReA/Q2Yc/x8zM+qgyha1RRDzWg80uAg6pXyBpMPA10i7O3YGjJe0u6W8kXd0wbfTyOJLeRto1+pM+/xFmZtZrpT3GNhAi4meSJjQs3gd4MCIeApB0GXBERHweOLwXv/tK4Mpsl+h/DlBkMzN7GZVtsfXDWOD3dfNd2bJuSRot6RvAnpI+kS07UNK5khaQhvfq7n5zJS2VtHTlypUDGN/MrL1VssUm6byI6Os12dTNso1euDQiVgEnNiz7KfDTTT1IRJwHnAfQ2dnpC6OamQ2Q0hY2SVtvbBXw1n786i5gx7r5ccAf+vH7zMyshUpb2EiXqXmEF7ewIpvfaKeOHrgV2EXSROAx4CjgmH78PjMza6EyF7aHgIMi4tHGFZJ+3832LyHpUuBAYJtsFJPTI2KhpA8C1wGDgW9GxD0DF9vMzJqpzIVtPvAK4CWFDTirJ78gIo7eyPJr2Einj4EkaRYwa/Lkyc1+KDOztqEI91vIW2dnZyxdujTvGGZmpSJpWUS8ZBCM0nb3l7T/y6zfStJrWpXHzMyKocy7It8l6SzS6P7LSJ1JhgOTgRnATsBH84tnZmZ5KG1hi4iTJL0COBJ4N7A98BRwH7AgIn6RZz4zM8tHaQsbQET8CTg/m0rHnUfMzAZeaY+xVUFEXBURc0eOHJl3FDOzynBhMzOzSnFhMzOzSil9YZM0QtKnJJ2fze8iqceXlzEzs2opfWEDLgSeBqZn813A5/KL03OSZkk6b82aNXlHMTOrjCoUtkkRcRawHiAinqL7S88UjjuPmJkNvCoUtmckbU52zTRJk0gtODMza0OlPo8tczpp9JEdJX0H2A94X66JzMz6afaCJQBcfsL0l9nSGpW+sEXEjyXdBryOtAvywxHxRM6xzMwsJ6UtbJL2alj0ePZzvKTxEXFbqzOZmVn+SlvYgC9lP4cDncByUottCnAzsMnR/83MrJpK23kkImZExAzgEWCviOiMiL2BPYEH803XM+7ub2btYvaCJS8cN2y20ha2OrtFxF21mYi4G3htfnF6zt39zcwGXpl3RdbcJ+kC4BJSl//3ki5dY2ZmbagKhe044B+BD2fzPwO+nl8cMzPLU+kLW0SsA87JJjOz0nt01VqWd61m3frnmTlvMQvnTGP86BF5xyqN0hc2SfsBZwA7Uff3RMTOeWUyM+ut+o4VtaIG8MCKJzl4/mKmjhvlk7V7qPSFDVgInAQsA57LOYuZWb/VitrG5m3TqlDY1kTEorxDmJn1R31rbOa8xTyw4kkABgkmjelwa60XqtDd/yZJZ0uaLmmv2pR3qJ7weWxm1p2Fc6YxfGj6eJ40poOFc6blnKh/ascMb374j8yct5hHV61t6uMpIpr6AM0m6aZuFkdEvKnlYfqos7Mzli5dmncMMyuQMg+C3Hgidv0xQ4DhQwcxddwooH9/n6RlEdHZuLz0uyKz0UfMzKygWn3MsLSFTdLJm1ofEfNalcXMzDZobIW1+phhmY+xbZlNnaQTtMdm04nA7jnmMjOzOq0+ZljaFltEfBpA0vWkQZD/ks2fAXwvx2hmZpXVl2N/40ePGJBjaj1V5hZbzXjgmbr5Z4AJ+UQxM7O8lbbFVufbwC2SfkgaBPkdwMX5RjIzs7yUvrBFxJmSFgFvyBYdFxG355nJzMzyU/rCBhARtwG35Z3DzMzyV4VjbKXlkUfMzAaeC1uOfAVtM7OBV4ldkWZmVVPGobSKwi02M7M2NnvBkpeM7Vh2LmxmZlYpLmxmZlYpPsZmZmZN18pjhm6xmZlZj7T6gqF95RabmZltUq1zSf0FQx9Y8SQHz1/c0sGNe8otNjMz65FWXzC0r9xiMzOzTaq1xlp9wdC+covNzMx6pNUXDO0rt9hyJGkWMGvy5Ml5RzEze1mtvmBoX7nFliOPFWlmNvBc2MzMrFJc2MzMrFJc2MzM2lRZTrjuLXceMTNrI/Uj+W/qhGsodgeRTXGLzcysTZXlhOvecovNzKyN1LfCynLCdW+5xWZm1qbKcsJ1b7nFZmbWpspywnVvucVmZmaV4sJmZmaV4sJmZmaV4sJmZmaV4sJmZmaV4sJmZmaV4sJmZmaV4sJmZmaV4hO0c+QraJtZ2ZThRG632HLkK2ibmQ08FzYzM6sUFzYzM6sUFzYzM6sUFzYzM6sUFzYzM6sUFzYzM6sUFzYzM6sUFzYzM6sUFzYzM6sUD6llZtbGyjBEVm+5xWZmZpXiwmZmZpXiwmZmZpXiwmZmZpXiwmZmZpXiwmZmZpXiwmZmZpXiwmZmZpXiwmZmZpWiiMg7Q9uTtBJ4BNgGeCLnON0pai4obrai5oLiZnOu3itqtlbl2ikixjQudGErEElLI6Iz7xyNipoLiputqLmguNmcq/eKmi3vXN4VaWZmleLCZmZmleLCVizn5R1gI4qaC4qbrai5oLjZnKv3ipot11w+xmZmZpXiFpuZmVWKC5uZmVWKC1uBSFLeGTamqNmcq/eKmq2ouaC42Zyrey5sOZI0XdI8SUcCRIEOeBY1m3P1XlGzFTUXFDebc/WMC1tOJB0EnA88Cpwo6WxJ2+QcCyhuNufqvaJmK2ouKG425+o5F7b8TAVuiIj5wPuAscBhkrbIM1SmqNmcq/eKmq2ouaC42Zyrh1zYWkTSoZLeLWl0tui3wPOSRkdEF3AlsC+wi7M5V9WyFTVXkbM5V9+5sDWRks0kXQR8BjgE+IqkTuD3wHDgVdnm3wO2AHbP7tvU16ao2ZyrOtmKmqvI2ZxrYLiwNVF2AFXAEOCQiDgeuAk4PyJuA9YCr5c0PiKeA34J/F123+fbMZtzVSdbUXMVOZtzDQwXtuabCIwC1ksaFBHnZ7dPAL4KjAdOzrbdBrjR2ZyrQtmKmqvI2ZyrvyLC0wBPpC8Mg+rmfw58qG5+L+ABYBiwPfBt0pvgFmBSO2ZzrupkK2quImdzrgHOndcDV2kCJgHzgaOA4XXLX5n9nA6sAIbWrfsOcEB2ezNgQjtlc67qZCtqriJnc67mTt4V2U+SppAOlq4G3gpcJGkHSZsBn5O0e0QsAX5COtg6UtIIYCRwP0BEPB0Rv2uXbM5VnWxFzVXkbM7VfC5s/TcZeCQizgD+HngKmA0EMDci7s22+yfgWWABsAx4DFgtNXXomaJmc67qZCtqriJnc65my7vJWLaJdH7GccDobH4vYCFZ8xuYkc2/sZv7inQy457tlM25qpOtqLmKnM25Wj+5xdZD2XkcZwIXAAcC8yW9A/hf4Dmyczgi4ibgz6Q3DZJGSPqIpF0jWR4Rt7dDNueqTrai5ipyNufKjwtbD0X6mrI1cGxEzAEuAb4MPAE8COwn6dXZ5lcAx2b3Wws8FBH/027ZnKs62Yqaq8jZnCs/Lmw9JGlb0hhoTwJExHXAr4CzgP8AOoAPSxoKjANuyG4TEVe2Yzbnqk62ouYqcjbnylE0cT9nFSZefA7Ht4Gv1s2/gjSi9Y7AVqRvPdcD9wKvb6dswNjsp4qUq6jPV5GzAVsX9bXkxd3Mi5ZNRcxV1OerqX9r3gGKNJHOqp9HGgpmXN3yrUhDyWxL6go7qW7dV4Gj6uYnNinbVsDZwAkNj79lntmyXF8DfgPsULe8I+dcI0nfQN8NbF/A1/JM0rlCOxbstZwPfBPYoiivZd3reS7wJWDfgj1n7wCGZfO14pb3+38UBf0sa8XkXZEZSTuS9icPBfYALpU0KVv9aWC/iFhBevG/LGmHbN1YsnM4ACLi4SZkO460q2AQ6RvVTXWrP0f6RpVLNuD/Ac8Dr4uIP9QtPzOvXEoDs15L6qY8CbhA0oRs9WfI97U8lvRabgW8DlhUtzrP5+x9wD3An0gjS/y1CLmybPsA15F2na0gtTZqcnv/S3obaQDgU0lfoOrl+VoW9rOsZfKurEWZgE7gF3XzXyG9CUYBgxu2nQdcBNwJ/BAYTfZNrQm5RgBHA39Tt+xHwLbZ7a3yyEYqsmOAG+uW7QZsk90enkeu7LGOAC6om/9vUrflzYEhOebaCvgQML1u2Z3A27PbHTlm+xBwT938NrXnKq/3WN3jvRW4sG7+h2StkJyfswOAucB7SMem6vdYbJFjrkJ+lrVyyj1Abn847Ap8AHgl6VjCdqRvglOz9XuQdsnMom7fdLZOwM7AG5qcbduG5eNIJ0TeTNol09HNfZuWrT5X7TnJ/hmOAi4kjRG3KHvOGgtbS3Jl8+8ltSQnZfMfB34KvKXxn7YFr+U4YGTdfC1TbdfVV9nIMYwmP2eNuYaTvuWflb2WlwPfJTvelvNzdlj2vjqZNAbh7cBlpMuiNH5QN/M5m0Q6cfmV9e8j4DWkwwQnbeK+rXj/F+6zLK+pLXdFSvon4GrSuGefBz5I2gXzV2BXSUMi4h7gf4ADI+J5SUMkzZU0DSAiHoqInzc52xckfTRbPgSYQGp5vIF0vaPPZOekDGp2tsZcwEnZql+T/mHuj4g3kS4y+CbSyZ5IOqGFuc6SNBe4hvTt9ExJ55Cer98Cb46IaMVrmb0mZ5IOyB8naVj2WL/Nfj6TbbofsC67j7KfTcvWTa5ab7d1pJEkTgTuAN4PrAdOq7tfq5+zzbJVvyTtctwPuD4i9gR+R2otjc3u28znbIikLwL/BRxM+lJ5XN0m95AK7u61DJIGZz9b+f4v1GdZrvKurHlMpH+Sd2a3p5KGhNmR9E3/HGDvbN1Y4BGylhFwOLB5Dtn2yObre4RtRzqhcrtWZNtIrp2BfYDlwKeydUNI366nZ/NvyynXMNIu3NOyTDuQivBmLXq+diW1FD8C/AB4VTfbvBq4pW6+NgLErGZl21Qu0rf3+t1prwT+wIYBcHN9zkgf3G/JbncAdwG7teA5O5z0hbLWMeQY0u69zeu22Rk4HTglm98h+3lYDu//QnyW5Tm1XYst+xY4kXRGPRGxnPSmXRARlwBPAydIGktqIf2M1AmBiLg6Ip5qcbbzSPvviYj1dZu/ijQY6VPNzraRXBeQnrNbSN/0p2YHqHchFZVnsm2vbHGu84GLIuKZiLg0Ij4XEc8CBwHXRsTT2bZNfS2Bh4BvRcR8YA1wjKTNG7YZBVwraX9Jt5J6vBIRVzUx20ZzRVLfAWgPWvQee7ls2c/BwC7Z7Z1IH+LPZdma+ZzdAHwpsoqQ6YiIp5RdHToiHgKuAg6W9Gfgo9nyH7X4/V+Iz7Lc5V1ZWzmx4RvX6cCihnX3k3Z1DAE+SeqF9QDwrgJkuxc4OLs9lXT8YwlwRAGeswOz23NJ+/XvAo4s0PM1ntRD8hbSrsg83m+7kXqy7kdd5xXg30g9Sm9sxXPWk1ykY22TScfXWvIe60G22vHId5K+RP0SuK2Vz1lDxiNJX55q80Oz6VfZ/8TsFmSoHeMu3GdZEabcAzTpRRdpH3jtpOHam6D+gO9jZLs1svnTgJPr5ncqULZPkh2YJp0z8+GC5DoN+Gjd/CsKkuuTtdeS1HPzxFa+lg3b1PJ+hvRNun7dB9hEh4Mcc725Ge+xgciWLT8g51znAJ/oZv0/NSnXGcDHSN30acyX12dZkafcAzThjfBaUkeBh4Bzullf+2Z6LHArsEs2/03g8AJne1t2+yX/cFV9zvqZa1aer2XddvUfQD8i7Sq9E9iT5pyG0d9c0wr6nN0B7J9nrmzbwaQBCSaQWm7/Dby2SbkOyP7uC0mnYzxC3UU88/wsK/pUxWNsK0gjFLwdGC/pjZB6XAFEOt5CRHyL1Mnh45KWkg7+3lfgbPdk654vWK5mPmf9yXVvd7+wVdlqaq+XpO1Ju4deA3wwIm6P7FOoQLk+EBG3NiHTQGT754j4RZ65MsNIp41cQ9r9/vmIuGOgA2W9KkcC/x4Rx0XEucDPgeOz9cr5s6zY8q6s/Z1I51SdRjr2NDJbNpx0YvNHgG93c5/aLoXBpBNRZ7ZTNudqTba6+w4l9aR7f7vkKnK2fuYaT+o6P+DHqxpyjaibai2z95P1uqy7T0ve/2WbagdqS0nSDNIJrktJPbe2iohj6ta/inRw9fqIuKib+yua9AQUNZtztT5bsxQ1V5Gz9SeXpEGRzgMb8PdaN7k6IuK9Ddt8F7guIhZ2c/+mvf9LKe/K2p+J9I3uS9ntQaSeP0fXrd882+aybH5nGoa6abdszlWdbEXNVeRsZcxF6t04nNRLtXblhYk0jLriacNU9mNs2wArJQ2PtF/+Y8DJdaMpPAVcCgyT9BfSyOVD2zybc1UnW1FzFTlb6XJFOpa2Oem48RRJVwOfAjbb+K9rc3lX1r5MbDjX5VDSOUBb1q27EfiX2PDN5ytAF03oilumbM5VnWxFzVXkbCXOdWp2eybpnMclwAmteC3LPOUeoAcv/FxSt9ra0DmDs5+1g6bXkIbhqb1B3gJ8jw0HXN9Dk3YlFDWbc1UnW1FzFTlbxXJ9n9QxZD9SK21EM17Lqk2F3RUpabqkJaSxBncFviNpVEQ8VzuIm236UdI4ce/J5l8NLI8NXWG/Gy++tlRlszlXdbIVNVeRs1U0150R8RywJCI+GxFrBypXpeVdWbubSCPXf5BsOB/SwdMLefEArTuTRrXenjSg55dJQ+3cSRNGJih6NueqTrai5ipyNufy9KLnPe8AdS/uUNIFBWsjnNePPnABaWSA92Zvgh1JzfYz6rYR8MZ2yuZc1clW1FxFzuZcnjY2FeI8NklvIY1KfQfpHI7LIuIH2bojSFeEvQ44EHhrRLxe0ojImuWSBkdqrrdNNueqTrai5ipyNueyTcq7smaF9VTgqOz2McDFZCN303DNINJQMa/Pbg+iCePtlSGbc1UnW1FzFTmbc3na1DSEYtifdMVXSIOeDgL+VtKLrrEkaQ/SZeGbOW5iWbI5V3WyFTVXkbM5l21Urr0ilV0+nXSNpXdlzfA1pAvi/YHUNRZJkyT9G3AJcFu2TVtmc67qZCtqriJncy7riZYVNkkHSJpSvyxSd1eRTjp8iNR7CGA18Ed4oUW5H7A1cFhEfLFdsjlXdbIVNVeRszmX9Vmz93WSLvT4K9IlF34C/COwXbZus7rtDiA1zXfN5ueTLqMBdVcdbodszlWdbEXNVeRszuWpv1MrWmyTgDsi4g3AvwFjgVMAIuJpSdtI+gipC+z3gHmSfki6iu/SbLtn2yybc1UnW1FzFTmbc1n/NKNaAjsAw7Lb/wDcmN0eSrpy8BXAm0hn4a8EPlV338nAnNr92yWbc1UnW1FzFTmbc3ka0NdtgN8ER5C+mXwPuCZbthmp+T4jm9+S1ISfRzoRcVzd/ZvZDbeQ2ZyrOtmKmqvI2ZzLUzOmfu+K1AZHkM7h+NeIeDewnaRTIuJp4HLS1V+JiL+QegkNIQ002iVpkNSUi/cVMptzVSdbUXMVOZtzWbP1q7BJGhYZ0n7lv4uI67PVnwL2yW7/F9Ah6ZRsfgWwfUQ8Cekcjib8Uxcym3NVJ1tRcxU5m3NZK/T5BG1J/wwcr3TRu19GxCJtOJcD4LXAgwAR8XtJnwYulDSBNNDnPEnK1g/0P3UhszlXdbIVNVeRszmXtUz0Yf8lqTvrzcBewFGkfdH7Z+s2y35+gWxombr77QS8HZjal8ctczbnqk62ouYqcjbn8tTKqTdvgMF1t48EPlM3fyLpekb12/8YGJ9Nn6JuhOsB/yMKms25qpOtqLmKnM25POU1vewxNklDJP078O9KI1cDPAu8sbZNRHwj23Zu9nMCsC3pwnk/zLYZ8LHQiprNuaqTrai5ipzNuSx3m6p6pGb6cuDrwPHArWxopv+GdIC1tu1hwLWkbq+vBp7I7rd1MypyUbM5V3WyFTVXkbM5l6ciTC/3ZnhDwwv+ZeDs7PYs4BFgeDa/F3BOdntHYHJTgxc0m3NVJ1tRcxU5m3N5KsK0yQuNShoBPAc8G2mQz6NJZ9ufGhHPS7oQeAa4AXg38OeIeP9Gf+EAKmo256pOtqLmKnI257Ii2OQxtohYGxFPx4Yrur4F+H1s2Mf8EeBKYDbwm1a+EYqazbmqk62ouYqczbmsEHrSrAMGk4rgImBStuw1QEd2O7ex0Iqazbmqk62ouYqczbk85Tn1dOSR50mDfj4BTJV0FamX0BCAiHimh7+nGYqazbmqk62ouYqczbksPz2tgMDrSG+KXwDH512Ry5DNuaqTrai5ipzNuTzlNW2y80g9SeOAvwPmRRoMtDCKms25eq+o2YqaC4qbzbksLz0ubGZmZmXQiitom5mZtYwLm5mZVYoLm5mZVYoLm5mZVYoLm5mZVYoLm1lJSRot6Y5s+l9Jj2W3n5T0H014vIskPSzpxGz+REnH9uH3zJb0YHbFarMB5+7+ZhUg6QzgyYj4YhMf4yLg6oj4/gD8rgOBj0XE4f39XWaN3GIzqxhJB9ZaQ5LOkHSxpOsl/U7SOyWdJekuSddKGpptt7ekxZKWSbpO0vY9eJwzJH0su/0hSfdKulPSZdmyrSVdkS37taQpzfy7zWpc2MyqbxLp4plHAJcAN0XE3wBPAYdlxe0rwJERsTfwTeDMXj7GqcCeETEFODFb9mng9mzZvwLf6vdfYtYDQ/IOYGZNtygi1ku6izS6/bXZ8ruACcCrSCPc/1gS2TaP9/Ix7gS+I+kK4Ips2f7AuwAi4sbsmODIiFjT9z/F7OW5sJlV39MAkS6ouT42HFh/nvQZIOCeiJjej8c4DHgj8DbgU5L2yH5vIx/Ut6bzrkgz+w0wRtJ0AElDs8LUI5IGATtGxE3Ax4FRQAfwM+Bvs20OBJ6IiD8PaHKzbrjFZtbmIuIZSUcC50oaSfpcmA/c08NfMRi4JLuvgHMiYnXWU/NCSXcCa4E5Ax7erBvu7m9mPeLu/lYW3hVpZj21Bvhs7QTtvpI0G/gP4E8DksqsgVtsZmZWKW6xmZlZpbiwmZlZpbiwmZlZpbiwmZlZpbiwmZlZpfwf5shKpavH2x8AAAAASUVORK5CYII=\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 analysis 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 simultaneously fit the all datasets using a given temporal model. This does a joint fitting across the different datasets, while simultaneously minimising 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 parameters of a model in a given group of dataset will be tied. For more details on joint fitting in gammapy, see [here](../2D/modeling_2D.ipynb)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Fit the datasets" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "execution": { "iopub.execute_input": "2021-11-22T21:05:33.804677Z", "iopub.status.busy": "2021-11-22T21:05:33.804371Z", "iopub.status.idle": "2021-11-22T21:05:33.807636Z", "shell.execute_reply": "2021-11-22T21:05:33.807831Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/adonath/software/mambaforge/envs/gammapy-dev/lib/python3.9/site-packages/astropy/units/quantity.py:486: RuntimeWarning: overflow encountered in exp\n", " result = super().__array_ufunc__(function, method, *arrays, **kwargs)\n", "/Users/adonath/software/mambaforge/envs/gammapy-dev/lib/python3.9/site-packages/astropy/units/quantity.py:486: 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": 16, "metadata": { "execution": { "iopub.execute_input": "2021-11-22T21:05:33.809962Z", "iopub.status.busy": "2021-11-22T21:05:33.808764Z", "iopub.status.idle": "2021-11-22T21:05:33.812022Z", "shell.execute_reply": "2021-11-22T21:05:33.812251Z" } }, "outputs": [ { "data": { "text/html": [ "
Table length=5\n", "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
typenamevalueuniterrorminmaxfrozenlink
str8str9float64str14int64float64float64boolstr1
spectralindex2.0000e+000.000e+00nannanFalse
spectralamplitude1.0000e-12cm-2 s-1 TeV-10.000e+00nannanFalse
spectralreference1.0000e+00TeV0.000e+00nannanTrue
temporalt01.0000e+01h0.000e+00nannanFalse
temporalt_ref5.8909e+04d0.000e+00nannanTrue
" ], "text/plain": [ "\n", " type name value unit error min max frozen link\n", " str8 str9 float64 str14 int64 float64 float64 bool str1\n", "-------- --------- ---------- -------------- --------- ------- ------- ------ ----\n", "spectral index 2.0000e+00 0.000e+00 nan nan False \n", "spectral amplitude 1.0000e-12 cm-2 s-1 TeV-1 0.000e+00 nan nan False \n", "spectral reference 1.0000e+00 TeV 0.000e+00 nan nan True \n", "temporal t0 1.0000e+01 h 0.000e+00 nan nan False \n", "temporal t_ref 5.8909e+04 d 0.000e+00 nan nan True " ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model.parameters.to_table()" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "execution": { "iopub.execute_input": "2021-11-22T21:05:33.814543Z", "iopub.status.busy": "2021-11-22T21:05:33.814268Z", "iopub.status.idle": "2021-11-22T21:05:33.815355Z", "shell.execute_reply": "2021-11-22T21:05:33.815570Z" } }, "outputs": [], "source": [ "datasets.models = model" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "execution": { "iopub.execute_input": "2021-11-22T21:05:35.221573Z", "iopub.status.busy": "2021-11-22T21:05:34.451283Z", "iopub.status.idle": "2021-11-22T21:05:37.053423Z", "shell.execute_reply": "2021-11-22T21:05:37.053590Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 3.22 s, sys: 41.5 ms, total: 3.26 s\n", "Wall time: 3.23 s\n" ] } ], "source": [ "%%time\n", "# Do a joint fit\n", "fit = Fit()\n", "result = fit.run(datasets=datasets)" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "execution": { "iopub.execute_input": "2021-11-22T21:05:37.056902Z", "iopub.status.busy": "2021-11-22T21:05:37.056575Z", "iopub.status.idle": "2021-11-22T21:05:37.057959Z", "shell.execute_reply": "2021-11-22T21:05:37.058136Z" } }, "outputs": [ { "data": { "text/html": [ "
Table length=5\n", "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
typenamevalueuniterrorminmaxfrozenlink
str8str9float64str14float64float64float64boolstr1
spectralindex2.9693e+003.133e-02nannanFalse
spectralamplitude1.0603e-11cm-2 s-1 TeV-13.703e-13nannanFalse
spectralreference1.0000e+00TeV0.000e+00nannanTrue
temporalt05.6466e+00h1.982e-01nannanFalse
temporalt_ref5.8909e+04d0.000e+00nannanTrue
" ], "text/plain": [ "\n", " type name value unit error min max frozen link\n", " str8 str9 float64 str14 float64 float64 float64 bool str1\n", "-------- --------- ---------- -------------- --------- ------- ------- ------ ----\n", "spectral index 2.9693e+00 3.133e-02 nan nan False \n", "spectral amplitude 1.0603e-11 cm-2 s-1 TeV-1 3.703e-13 nan nan False \n", "spectral reference 1.0000e+00 TeV 0.000e+00 nan nan True \n", "temporal t0 5.6466e+00 h 1.982e-01 nan nan False \n", "temporal t_ref 5.8909e+04 d 0.000e+00 nan nan True " ] }, "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.9.0" }, "latex_envs": { "LaTeX_envs_menu_present": true, "autoclose": false, "autocomplete": true, "bibliofile": "biblio.bib", "cite_by": "apalike", "current_citInitial": 1.0, "eqLabelWithNumbers": true, "eqNumInitial": 1.0, "hotkeys": { "equation": "Ctrl-E", "itemize": "Ctrl-I" }, "labels_anchors": false, "latex_user_defs": false, "report_style_numbering": false, "user_envs_cfg": false }, "nbsphinx": { "orphan": true } }, "nbformat": 4, "nbformat_minor": 4 }