\n",
"\n",
"**This is a fixed-text formatted version of a Jupyter notebook**\n",
"\n",
"- Try online [](https://mybinder.org/v2/gh/gammapy/gammapy-webpage/master?urlpath=lab/tree/spectrum_analysis.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",
"[spectrum_analysis.ipynb](../_static/notebooks/spectrum_analysis.ipynb) |\n",
"[spectrum_analysis.py](../_static/notebooks/spectrum_analysis.py)\n",
"
\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Spectral analysis with Gammapy"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Prerequisites \n",
"\n",
"- Understanding how spectral extraction is performed in Cherenkov astronomy, in particular regarding OFF background measurements. \n",
"- Understanding the basics data reduction and modeling/fitting process with the gammapy library API as shown in the [first gammapy analysis with the gammapy library API tutorial](analysis_2.ipynb)\n",
"\n",
"## Context\n",
"\n",
"While 3D analysis allows in principle to deal with complex situations such as overlapping sources, in many cases, it is not required to extract the spectrum of a source. Spectral analysis, where all data inside a ON region are binned into 1D datasets, provides a nice alternative. \n",
"\n",
"In classical Cherenkov astronomy, it is used with a specific background estimation technique that relies on OFF measurements taken in the field-of-view in regions where the background\n",
"rate is assumed to be equal to the one in the ON region. \n",
"\n",
"This allows to use a specific fit statistics for ON-OFF measurements, the wstat (see `~gammapy.stats.fit_statistics`), where no background model is assumed. Background is treated as a set of nuisance parameters. This removes some systematic effects connected\n",
"to the choice or the quality of the background model. But this comes at the expense of larger statistical uncertainties on the fitted model parameters.\n",
"\n",
"**Objective: perform a full region based spectral analysis of 4 Crab observations of H.E.S.S. data release 1 and fit the resulting datasets.**\n",
"\n",
"## Introduction\n",
"\n",
"Here, as usual, we use the `~gammapy.data.DataStore` to retrieve a list of selected observations (`~gammapy.data.Observations`). Then, we define the ON region containing the source and the geometry of the `~gammapy.datasets.SpectrumDataset` object we want to produce. We then create the corresponding dataset Maker. \n",
"\n",
"We have to define the Maker object that will extract the OFF counts from reflected regions in the field-of-view. To ensure we use data in an energy range where the quality of the IRFs is good enough we also create a safe range Maker.\n",
"\n",
"We can then proceed with data reduction with a loop over all selected observations to produce datasets in the relevant geometry.\n",
"\n",
"We can then explore the resulting datasets and look at the cumulative signal and significance of our source. We finally proceed with model fitting. \n",
"\n",
"In practice, we have to:\n",
"- Create a `~gammapy.data.DataStore` poiting to the relevant data \n",
"- Apply an observation selection to produce a list of observations, a `~gammapy.data.Observations` object.\n",
"- Define a geometry of the spectrum we want to produce:\n",
" - Create a `~regions.CircleSkyRegion` for the ON extraction region\n",
" - Create a `~gammapy.maps.MapAxis` for the energy binnings: one for the reconstructed (i.e. measured) energy, the other for the true energy (i.e. the one used by IRFs and models)\n",
"- Create the necessary makers : \n",
" - the spectrum dataset maker : `~gammapy.makers.SpectrumDatasetMaker`\n",
" - the OFF background maker, here a `~gammapy.makers.ReflectedRegionsBackgroundMaker`\n",
" - and the safe range maker : `~gammapy.makers.SafeRangeMaker`\n",
"- Perform the data reduction loop. And for every observation:\n",
" - Apply the makers sequentially to produce a `~gammapy.datasets.SpectrumDatasetOnOff`\n",
" - Append it to list of datasets\n",
"- Define the `~gammapy.modeling.models.SkyModel` to apply to the dataset.\n",
"- Create a `~gammapy.modeling.Fit` object and run it to fit the model parameters\n",
"- Apply a `~gammapy.estimators.FluxPointsEstimator` to compute flux points for the spectral part of the fit.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setup\n",
"\n",
"As usual, we'll start with some setup ..."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"gammapy: 0.18.2\n",
"numpy: 1.18.5\n",
"astropy 4.0.1.post1\n",
"regions 0.4\n"
]
}
],
"source": [
"# Check package versions\n",
"import gammapy\n",
"import numpy as np\n",
"import astropy\n",
"import regions\n",
"\n",
"print(\"gammapy:\", gammapy.__version__)\n",
"print(\"numpy:\", np.__version__)\n",
"print(\"astropy\", astropy.__version__)\n",
"print(\"regions\", regions.__version__)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"from pathlib import Path\n",
"import astropy.units as u\n",
"from astropy.coordinates import SkyCoord, Angle\n",
"from regions import CircleSkyRegion\n",
"from gammapy.maps import Map, MapAxis\n",
"from gammapy.modeling import Fit\n",
"from gammapy.data import DataStore\n",
"from gammapy.datasets import (\n",
" Datasets,\n",
" SpectrumDataset,\n",
" SpectrumDatasetOnOff,\n",
" FluxPointsDataset,\n",
")\n",
"from gammapy.modeling.models import (\n",
" PowerLawSpectralModel,\n",
" create_crab_spectral_model,\n",
" SkyModel,\n",
")\n",
"from gammapy.makers import (\n",
" SafeMaskMaker,\n",
" SpectrumDatasetMaker,\n",
" ReflectedRegionsBackgroundMaker,\n",
")\n",
"from gammapy.estimators import FluxPointsEstimator\n",
"from gammapy.visualization import plot_spectrum_datasets_off_regions"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Load Data\n",
"\n",
"First, we select and load some H.E.S.S. observations of the Crab nebula (simulated events for now).\n",
"\n",
"We will access the events, effective area, energy dispersion, livetime and PSF for containement correction."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"datastore = DataStore.from_dir(\"$GAMMAPY_DATA/hess-dl3-dr1/\")\n",
"obs_ids = [23523, 23526, 23559, 23592]\n",
"observations = datastore.get_observations(obs_ids)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Define Target Region\n",
"\n",
"The next step is to define a signal extraction region, also known as on region. In the simplest case this is just a [CircleSkyRegion](http://astropy-regions.readthedocs.io/en/latest/api/regions.CircleSkyRegion.html), but here we will use the ``Target`` class in gammapy that is useful for book-keeping if you run several analysis in a script."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"target_position = SkyCoord(ra=83.63, dec=22.01, unit=\"deg\", frame=\"icrs\")\n",
"on_region_radius = Angle(\"0.11 deg\")\n",
"on_region = CircleSkyRegion(center=target_position, radius=on_region_radius)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Create exclusion mask\n",
"\n",
"We will use the reflected regions method to place off regions to estimate the background level in the on region.\n",
"To make sure the off regions don't contain gamma-ray emission, we create an exclusion mask.\n",
"\n",
"Using http://gamma-sky.net/ we find that there's only one known gamma-ray source near the Crab nebula: the AGN called [RGB J0521+212](http://gamma-sky.net/#/cat/tev/23) at GLON = 183.604 deg and GLAT = -8.708 deg."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAQ8AAAEHCAYAAACwfMNTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAASsklEQVR4nO3dfbBcdX3H8ffHJEAeQAgBxZB6r7QYIlVJY6HGCeIDQwGRAWSqQkHUtjYiKD4VHavjw6ipKHWqowMiLSiPAR/GB5CmaKxkGkJIiEGp3GCFjMaLkAQjEPn2j/NbsrnZu3vu7+7Zh3s/r5mdPXvO2bPf3bv3s7/z9DuKCMzMxuoZ3S7AzPqTw2MESep2DaPp5dp6Wa9/br1cX7PaHB57eqzbBTSxvdsFjEaSP7d8vVzfqLU5PPbkjUBmJTg8zCyL+nlvy5w5c2JgYKCtyxwaGmJwcLCty2wX15anl2uD3q5vaGgIIIaHh/doaEztfDntMzAwwOrVq7tdhtmEtmjRooYbTb3aYmZZHB5mlsXhYWZZHB5mlsXhYWZZHB5mlsXhYWZZHB5mlsXhYWZZHB5mlsXhYWZZHB5mlsXhYWZZHB5mlsXhYWZZHB5m1lTqEGgPDg8za+qRRx5pON7hYWZZHB5m1tRo/Rw7PMysqYULFzYc7/AwsywODzPL4vAwsywODzPL4vAwsywODzPL4vAwsywODzPL4vAwsywODzPL4vAwsywODzPL4vAwsywODzPL4vAwsywODzPL4vAwsywODzPL4vAwsywODzPL4vAwsywODzPL4vAwsywODzPLUml4SJonaYWkjZI2SLpgxPR3SwpJc+rGLZO0WtKxVdZmZuNTdctjJ3BRRBwBHAMslbQAimABXg38sjazpPlpcAmwtOLazGwcKg2PiNgcEWvS8DZgIzA3Tf4s8F6g/kKYU4Cn0jhVWZuZjU/HtnlIGgCOAlZJOgV4MCLurp8nIjYAM4CVwBc7VZuZlSPp6R/1qR16wVnAjcCFFKsyHwCObzRvRJzfZDkCttcez549u72FmtkehoaGkPRY7bGkWRERlbc8JE2jCI6rI2I5cBgwCNwtaRNwKLBG0rNbLSsKM2u3wcHBKks3M2BwcJD6/7uICKi45ZFaCpcDGyPiEoCIWA8cXDfPJmBRRPy2ylrMrL2qbnksBs4GXiFpbbqdWPFrmlkHVNryiIiVtNhrEhEDVdZgZtXwEaZmlsXhYWZZHB5mlsXhYWZZHB5mlsXhYWZZHB5mlsXhYWZZHB5mlsXhYWZZHB5mlsXhYWZZHB5mlsXhYWZZHB5mlsXhYWZZHB5mlsXhYWZZHB5mlsXhYWZZHB5mlsXhYWZZHB5mlsXhYWZZHB5mlsXhYWZZHB5mlsXhYWZZHB5mlmVqmZkkLQY+DDw3PUdARMTzqivNzHpZqfAALgfeCdwJ/LG6csysX5QNj0cj4ruVVmJmfaVseKyQtAxYDjxeGxkRayqpysx6XtnwODrdL6obF8Ar2luOmfWLUuEREcdVXYiZ9ZdSu2olPVPSJZJWp9tnJD2z6uLMrHeVPc7jK8A24Mx02wpcUVVRZtb7ym7zOCwiTq97/BFJa6soyMz6Q9mWxw5JL6s9SAeN7aimJDPrB2VbHm8DrkzbOQQ8DJxbVVFm1vvK7m1ZC7xI0n7p8dZKqzKzntc0PCSdFRFXSXrXiPEARMQlFdZmZj2sVctjZrrft8G0aHMtZtZHmoZHRHwpDf4gIn5cPy1tNDWzSars3pbPlxxnZpNEq20efwW8FDhoxHaP/YApVRZmZr2t1TaPvYBZab767R5bgTOqKsrMel+rbR63A7dL+mpEPNChmsysD5Q9SOz3qT+PFwD71EZGhE/JN5ukym4wvRq4FxgEPgJsAv6n1ZMkzZO0QtJGSRskXZDGL5N0r6R1km6StH/dc5alM3ePHfO7MbOOKRseB0bE5cCTEXF7RJwHHFPieTuBiyLiiDT/UkkLgFuBIyPihcDPgX8CkDQ/PW8JsHQM78PMOqzsasuT6X6zpJOAh4BDWz0pIjYDm9PwNkkbgbkRcUvdbHewa+PrFOApigPQVLI2M+uCsuHxsXRS3EUUx3fsR9GbemmSBoCjgFUjJp0HXAsQERskzQBWAu8Zy/LNrLMUUf1R5pJmAbcDH4+I5XXjP0DRL+ppUaIQFSfVbK89nj179ozh4eEKKjazmgMPPJCHH37493WjZkVElL3o00HAW4EB6loradtHq+dOA24Erh4RHOcAJwOvLBMc6fWCXefbsGjRIp9fY1axwcFBhoeHZ44cX3a15RvAj4AfMIaLPqWWwuXAxvozcCWdALwPODYifj/a882sd5UNjxkR8b6M5S8GzgbW13VbeDHwr8DewK3p9P47IuIfMpZvZl1SNjy+LenEiPjOWBYeEStpvNdkTMsxs95T9jiPCygCZIekrZK2SXJvYmaTWNluCBt1BmRmk1irU/LnR8S9khY2mu5r1ZpNXq1aHhdR7KL9TINpvlatWRfV+hJupBPHb7U6Jf+t6d7XqjWz3bRabTmt2fT6g77MrHrNWhujzVdVK6TVastrmkwLwOFhVrGygVHm+e0MklarLW9q2yuZ2YRS6jgPSZ8Y0WHPAZI+Vl1ZZgbjb3VUubyyB4n9dUQ8UnsQEb8DTmxbFWa2G0ltD452L7tseEyRtHfdi0+nODfFzCapsue2XAXcJukKig2l5wFXVlaVmfW8soenf1rSOuBVFCe6fTQivl9pZWaTVFWrK41eZzx7X8q2PAA2Ajsj4geSZkjaNyK2Zb+ymfW1sntb3grcANQufD0XuLmqosys95XdYLqUomOfrQARcR9wcFVFmVnvKxsej0fEE7UHkqZSbDg1s0mqbHjcLuliYLqkVwPXA9+qriwz63Vlw+P9wBZgPfD3FN0IfrCqoswmoyoPDKviNcvuqn1K0s3AzRGxJeuVzKyp2m7TTgbIeHbVNm15qPBhSb+luND1zyRtkfSh7Fc0swmh1WrLhRR7WV4SEQdGxGzgaGCxpDFdbtLMJpZW4fG3wOsjYqg2IiLuB85K08xskmoVHtMi4rcjR6btHtOqKcnM+kGr8Hgic5qZTXCt9ra8aJSLOwnYp4J6zCa9iOjIHpfxdknYqhvCKeNauplNWGM5q9bMOqTKYz7a1Qly2SNMzawL2n3ZhHYuz+Fh1uPa8Q8fEW0PIoeHmWVxeJj1gfG0HLp1xTgz6yGduIB1WW55mFkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZak0PCTNk7RC0kZJGyRdkMbPlnSrpPvS/QF1z1kmabWkY6uszczGp+qWx07goog4AjgGWCppAfB+4LaI+DPgtvQYSfPT85YASyuuzczGodLwiIjNEbEmDW8DNgJzgdcCV6bZrgROTcNTgKeAoLgqnZn1qI5t85A0ABwFrAKeFRGboQgY4OA0vAGYAawEvtip2sxs7DrSAbKkWcCNwIURsbXZVbAi4vwmyxGwvfZ49uzZ7SzTzBoYGhpC0mN1o2ZFRFTe8pA0jSI4ro6I5Wn0ryUdkqYfAvymzLKiMLN2GxwcrKZoM3va4OAg9f93kbpwr3pvi4DLgY0RcUndpG8C56Thc4BvVFmHmbVf1asti4GzgfWS1qZxFwOfBK6T9Gbgl8DrKq7DzNqs0vCIiJWMvtfklVW+tplVy0eYmlkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4mFkWRUS3a8g2Z86cGBgYyH7+0NAQg4OD7SuoAv1QI/RHna4xz9DQUAwPD+/R0Ojr8BgvSY9FxMxu19FMP9QI/VGna2wvr7aYWRaHh00m6nYBJfRDjYBXWxQ9/gH0Q439oh8+y36osWZSh4eZ5Zuwqy2SzpV0crfrMJuopna7gHaQtAnYBvwR2BkRi9KkMyWdAPw6Ij5aN/8UYDXwYEScnMa9E3gLEMB64E0R8QdJ5wLHATuAzcA04EjgzIh4ok3171aPpHnAvwPPBp4CvhwRl3ailga17QP8ENib4vtyQ0T8c5q2P3BZqiGA84Dn91iNJwCXAlOAyyLik934HFMtDf+uddNHfg+6UmdpEdH3N2ATMGfEuHOBN6bha0dMexfwNeDb6fFcYAiYnh5fB5xbt5w3pOHb0v3FwFFtrH9kPYcAC9PwvsDPgQWdqKVBbQJmpeFpwCrgmPT4SuAtaXgvYP9eqpEiMH4BPC/Vd3e3Psdmf9cm34Ou1Fn2NmFXW5JH0/3TG3YkHQqcRPGLWW8qMF3SVGAG8FDdtK3pfku6f4LiV27cGtUTEZsjYk0a3gZspAi4SmtpJArb08Np6RaS9gOWAJen+Z6IiEd6qUbgL4H/jYj7o/ilvgZ4bTdqTHWO+ndt8r3seJ1lTZTwCOAWSXdK+rsW834OeC9Fs7F4csSDwL8Av6RoHj4aEbdUVWyreupJGgCOovg17QpJUyStBX4D3BoRqyh+zbcAV0i6S9Jlkrp2cNMoNc4F/q9utl+xK4S7qsHften3oCd1u+nTpubgc9L9wRRN0yWjzHcy8IU0/HJ2NQ8PAP4TOIjiV+tm4KwO1N2wnrrps4A7gdO6/RmnevYHVlCsdy8CdgJHp2mXAh/tsRpfR7GdozbtbODzPVDjbn/XVt+DXr1NiJZHRDyU7n8D3ETRXG1kMXBK2sB6DfAKSVcBrwKGImJLRDwJLAdeWnnho9eDpGnAjcDVEbG8A7W0FMVqyX8BJ1D8iv8qil94gBuAhV0q7WkNapxXN/lQdl8d7bhR/q6jfg96WrfTqw0pPhPYt274v4ETSjzv5exqeRwNbKDY1iGKDYHnd/h91Ncjiq3yn+uBz/cgYP80PB34EXByevwj4Plp+MPAsl6qkWI71v3AILs2mL6gi59ly78rfdTymAi7ap8F3CQJii/L1yLie2NZQESsknQDsIaiKX4X8OV2FzoGiyma2OvTejzAxRHxnS7UcghwZdqN+Azguoj4dpp2PnC1pL0o/knf1IX6mtYo6e3A9yn2vHwlIjZ0qUborb/ruPkIUzPLMiG2eZhZ5zk8zCyLw8PMsjg8zCyLw8PMsjg8zCyLw8PMsjg8JiBJf5S0VtI9kr6V+t1A0nPSwXCtnr99lPGnSlrQ4rl3S/p6XuXtUfZ92vg4PCamHRHx4og4EngYWArFOUARccY4lnsqRX8YDUk6guI7taSbZ9i24X1aCQ6Pie8n7OozYkDSPWl4hqTrJK2TdK2kVZJqPbAh6eOpFXGHpGdJeilwCrAstWoOa/BabwD+A7glzVtb1jsk/TS91jVp3CxJV0han8afnsYfL+knktZIul7SrDR+k6SPpPHrJc1P449N9axNXQPsO+J97lP3OndJOi6NP1fScknfk3SfpE+3+XOf+Lp9co1v7b8B29P9FOB60omCwABwTxp+N/ClNHwkxTk9i9LjAF6Thj8NfDANfxU4o8nr/hx4LnA88M268Q8Be6fh2glsn6LuBDGKbhHmUHQnODONex/woTS8iXSyIvCPpFPtgW8Bi9PwLIrzm+rf50XAFWl4PkWfLftQ9NJ1P/DM9PgBYF63/3b9dHPLY2Kank68GgZmA7c2mOdlFKd/ExH3AOvqpj0B1E5+u5Pin7EpSS8BtkTEA8BtwEJJB6TJ6yhOoDuLIqSg6Abh32rPj4jfUXQduAD4car/HIowqqmdwl5f04+BSyS9gyKYdrK7l1G0hoiIeylC4vA07baIeDQi/gD8dMRrWQsOj4lpR0S8mOKfYS/SNo8Rml1c6MlIP9UUnUqXOfv69cD81CfFL4D9gNPTtJMoguIvgDtTV4+irnvIuppujWJ7zYsjYkFEvLlu+uMja4qIT1J0XD0duKO2OlPyfT5eN1z2fVri8JjAIuJR4B3Au1MnNPVWAmcCpD0of15ikdsoOu7djaRnUPTa9cKIGIiIAYq+Ql+fps2LiBUU3eztT7F6cQvw9rplHADcASyW9Kdp3AxJh9OEpMMiYn1EfIqi5/GR4fFD4I1p3sOBPwF+VuK9WgsOjwkuIu6i6ATnb0ZM+gJwkKR1FNsW1rGrw+jRXAO8J214rN9guoTicgEP1o37IcUqyFzgKknrKfpJ+WwUvX19DDgg7U6+GzguIrZQbIv4eqrrDvYMg5EurFvGDuC7Dd7nlPT611L0iv/4yIXY2Lk/j0kqdZwzLYpr0xxGsZ3i8OiF64FYX/A63uQ1A1iRVmcEvM3BYWPhloeZZfE2DzPL4vAwsywODzPL4vAwsywODzPL8v9RjMCNrrGEkwAAAABJRU5ErkJggg==\n",
"text/plain": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.plot(\n",
" info_table[\"livetime\"].to(\"h\"),\n",
" info_table[\"sqrt_ts\"],\n",
" marker=\"o\",\n",
" ls=\"none\",\n",
")\n",
"plt.xlabel(\"Livetime [h]\")\n",
"plt.ylabel(\"Sqrt(TS)\");"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally you can write the extrated datasets to disk using the OGIP format (PHA, ARF, RMF, BKG, see [here](https://gamma-astro-data-formats.readthedocs.io/en/latest/spectra/ogip/index.html) for details):"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"path = Path(\"spectrum_analysis\")\n",
"path.mkdir(exist_ok=True)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"for dataset in datasets:\n",
" dataset.to_ogip_files(outdir=path, overwrite=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If you want to read back the datasets from disk you can use:"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"datasets = Datasets()\n",
"for obs_id in obs_ids:\n",
" filename = path / f\"pha_obs{obs_id}.fits\"\n",
" datasets.append(SpectrumDatasetOnOff.from_ogip_files(filename))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Fit spectrum\n",
"\n",
"Now we'll fit a global model to the spectrum. First we do a joint likelihood fit to all observations. If you want to stack the observations see below. We will also produce a debug plot in order to show how the global fit matches one of the individual observations."
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"spectral_model = PowerLawSpectralModel(\n",
" index=2, amplitude=2e-11 * u.Unit(\"cm-2 s-1 TeV-1\"), reference=1 * u.TeV\n",
")\n",
"model = SkyModel(spectral_model=spectral_model, name=\"crab\")\n",
"\n",
"for dataset in datasets:\n",
" dataset.models = model\n",
"\n",
"fit_joint = Fit(datasets)\n",
"result_joint = fit_joint.run()\n",
"\n",
"# we make a copy here to compare it later\n",
"model_best_joint = model.copy()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Fit quality and model residuals"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can access the results dictionary to see if the fit converged:"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"OptimizeResult\n",
"\n",
"\tbackend : minuit\n",
"\tmethod : minuit\n",
"\tsuccess : True\n",
"\tmessage : Optimization terminated successfully.\n",
"\tnfev : 47\n",
"\ttotal stat : 124.41\n",
"\n"
]
}
],
"source": [
"print(result_joint)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A simple way to inspect the model residuals is using the function `~SpectrumDataset.plot_fit()`"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(0.1, 40)"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAgIAAAGtCAYAAABk5xxqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deZxU1Zn/8c8jEGgUGRSdiaI2RiESQdDGBTJgRBQF3CAqEiOioGj8KdEkuAUcNyaa4BYwaAiaMa6MRkRiBEXGgCOguLNERcFkIhDTLkBkeX5/dHdZ3V1VfbuWe6vqft+v133RVXXuvU91F+c+dc4955i7IyIiIvG0U9QBiIiISHSUCIiIiMSYEgEREZEYUyIgIiISY0oEREREYqxl1AGEycyGAkPbtWs3pkuXLlGHIyIiEoply5ZtcPc9Ur1mcRw+WFVV5UuXLo06DBERkVCY2TJ3r0r1mroGRAIYNmwYw4YNizoMEYmBsOubWHUNiGRr48aNUYcgIjERdn2jFgEREZEYUyIgIiISY0oEREREYkz3CIgEMGDAgKhDEJGYCLu+0fBBERGRMqfhgyIiIpKSEgGRAE444QROOOGEqMMQkRgIu77RPQIiAWzevDnqEEQkJsKub9QiICIiEmOxSgTMbKiZTa+uro46FBERkaIQq0TA3We7+9j27dtHHYqIiEhR0D0CIgEMGTIk6hBEJCbCrm80j4CIiEiZ0zwCIiIikpISAZEAjj76aI4++uiowxCRGAi7vlEiICIiEmNKBERERGJMiYCIiEiMafigSAFNeXYVt89fnXh86YADGT+wS4QRiYjUp0RAJIDTTz89q/3GD+zC+IFdqJwwhzWTB+c5KhEpR9nWN9lSIiASwEUXXRR1CCISE2HXN7pHQCSATZs2sWnTpqjDEJEYCLu+UYuASAAnnngiAAsWLACa6Puf0h2qP2xwhN/BpAZrXLTfF8a/0aw4dM+BSPlrWN8UmhIBkSxk7Puv/hAmNVjhcsKcxs8lJQZBL/C650BE8k1dAyJFYPzALokL+5rJg/UtX0RCo0RAJM+mbB1W7/GyDz6p92+6ciIiUVAiIJJnt2//6gK/7INPGHnvSwCMvPeleslAcjkRkajoHgGRAEaNGtWs8pUT5jR6bsvWHQybtihPEYlIuWpufZOrWCUCZjYUGHrAAQdEHYqUmOb+x6zr769rEdiydQdtWu3EA+cfyWH7dQBqk4WGIwnyNLpAREqXEoECcvfZwOyqqqoxUccipWXDhg0AdOzYscmyl7aYBdQkAoft14EHzj+SYdMW1UsCoGZkAAObN7pARMpfc+qbfNA9AiIBDB8+nCMH1B+ul+4mwIbqLv7JSQAQeGRA0JsPoWYYooiUtuHDhzN8+PDQzqdEQCSgdZ98NdNXPm4CbHjRTneBD3rzIVBvLgIRkSBi1TUgkqt83gR4+/zViVaBhhf4ht0IuvlQRApFiYBIMwS+CTCgoBf4IOdt7rmbS9Mbi5QndQ2IBNSpQ9vEz3U3AQKNbwJsMSvwMddMHsyayYOZNa4PbVrV/Hds02onZo3rk3gtWabzQu0NiAWi2Q9FypNaBEQCGDduXKPn0t4E2GoWMKPJYyZftDOOLkgahZDpvBD8BkQRKV6p6ptCUiIgEsAZZ5yR92M2vGjnmliISHkoRH2TiboGRAJYu3Yta9eujToMEYmBsOsbJQIiAZx99tmcffbZzd5vyrOrEjfwVU6Yo3H+ItKkbOubbKlrQCSodUuCTwlca/zALrn327ffN9h568rWTkesu/xFJAglAiJBbfsnTPqi/nOppgTOt1TrDKQ7b1JyUJeEVE6Y02j0QTFQoiJSHNQ1ICJpBZ39MJsuDw1HFCkOSgREJK3kb+wZp1XW1MYiJUtdAyIBXH755fC7l6MOo15zeuWEOaE0p2t6Y5FwXX755aGeT4mASABDhw6FZa0Sj6O4IEOebj5spnxPqywimQ0dOjTU8ykREAlg5cqVsGE7XWsf5/uCHFVi0ZTAsx8WcGpjkbhZuXIlAF27dm2iZH4oERAJ4IILLoA1W1hwV2GOH8U3/SACz35YhLGLlKoLLrgAgAULFoRyPiUCIuUk6JwDSfMNiEi8KREQKSdB5xxINRmRiMSShg+KlJHA4/63DgstJhEpbkoERMpI4HH/26NJBAo5QZGIZEddAyIBXHPNNXD/yVGHEUgxj/u/ff7qxI2FDROV5JEIyeVE4uaaa64J9XxKBEQCOPbYY+HF0vjvUuzj/os5UREpBscee2yo51PXgEgAy5cvZ/n/bY86jCalGvcPNB7332JW6LHVWTN5MGsmD2bWuD60aVVTBbVptROzxvVJvCYSZ8uXL2f58uWhna80vuKIROyyyy6rmUcg6kCaEHjcf8dUSyqnkMUww0yrCmqCIpGmXXbZZYDmEQjMzPYHrgbau/vwqOMRKQlBL+5ZDDPMtPyxJigSKT5F2TVgZjPM7GMze7PB84PMbKWZ/dnMJgC4+3vufl40kYqIiJS2okwEgJnAoOQnzKwF8EvgBKAbMMLMuoUfmpSzKc+uonLCnMSmYWxfqfvdAPrdiJSRouwacPeFZlbZ4OnDgT+7+3sAZvYQcDLwdrjRSTnL1Kwdd8W6HkKdTPcmiEh6RZkIpLE3sDbp8TrgCDPbHbgR6GVmV7r7zal2NrOxwFiAfffdt9CxSpm56aab4NcDow5DMlASJ+XipptuCvV8pZQIWIrn3N03Ahc2tbO7TwemA1RVVXmeY5My16dPH/hjKf13EZFS1adPn1DPV0o12zpgn6THnYC/RBSLxMyiRYtg7TbC/e9ZBFKuZpimnFYzFMmLRYtqJtcKKyEopURgCXCgmXUGPgLOBM6KNiSJi6uuugrW/LPo5xHIuwIOMxSR1K666iogvHkEinLUgJk9CCwGuprZOjM7z923AT8AngHeAR5x97eaedyhZja9urq66cIiIiIxUJQtAu4+Is3zTwNP53Dc2cDsqqqqMdkeQ0REpJwUZYuAiIQnH3MnTNlaf1ljLS8sUjqUCIiUqaATAI0f2CUx3G7N5MFZjb2/ffuwmvsEJrVn2U+rGDnteQBGTnueZT+tSryWPM5fRIpDUXYNiBSb2267De7+dtRhNEvYEwBVbvldo+e20JphX16XsnzyBECVE+aU9ARAmsxI8um2224L9XxKBEQC6NmzJ/xbi6jDKGp1rQrLPviEkfe+xJatO2jTaqd6KwvWtVBA8c9U2ByazEjyqWfPnqGeL1ZdAxo1IKk0bDJP1b89b948xq06KtS4Skmq5YWBxssLt5gVemwipWbevHnMmzcvtPPFKhFw99nuPrZ9e415lq8kN+nWfZsFGHnvS4lk4IYbbuA3L34USXylIPDywq3ymwgESeLSlRUpVjfccAM33HBDaOdT14AI9Zus62zZuoNh02pm+Pq/9zaGHZIEcPv81YkkpGES17A1IrmsiHxFiYAITfdvH/3SLbykZCB3BZiyuKkkTkQyUyIgsZeqf3vYtEWNvlF2sg1RhFdeCjBlcZCbFCF1wiAiSgREAvdvd7L1ocUkwQRN4hqWjYqGGUoxilUiYGZDgaEHHHBA1KFIifnVr34Fdx4WdRjSQNAkLlXZKGiYoQTxq1/9KtTzadSASABdu3ala0fNIyAihde1a1e6du0a2vli1SIgkq3Zs2fDyq0MjToQCZWa8iUKs2fPBmDo0HBqHCUCIgH8/Oc/hzVflk0iMOXZVfUuaMnj75Ob1RuWC00BRhdkQ035EoWf//zngBIBESmgoOPvIxt7X4DRBSKSmhIBkZjS+PvcFH2rikhAsbpZUES+smbyYNZMHsyscX1o06qmKmjTaidmjeuTeE3SCzI1dcNyIsUoVi0CGj4oKU3pDtUfNnjyd/Wbndd8AS1bhxpWIQUdf18MY++LmVpVpBzEKhFw99nA7KqqqjFRxyJFpPpDmNRgRcoJc+o999vz1oYcVGEFXiRITdoZNXfp5XzTqIby9Nvf/jbU86lrQCSAffbZh3322SfqMKSIBF56uYCtKuMHdkkkI2smD1YSUCbCrm9i1SIgkq2HH34YgDPOOCPiSKSelMMMf5d6NEGehxqqVUUKJez6RomAxEKuTajTpk0DlAgUnVQX9gbdOgkaaiglIuz6Rl0DEgtqQg3PlGdXJfrFKyfMYcqzqyKOSEQyUYuAiORV3Wx8IlIa1CIgIlJADVtEkiceylROJCyxSgTMbKiZTa+uTtF/KCJSAJp4SIpdrLoGNI+ANCX5psLKCXMSNxU+9thjEUcmpUwTD0lzhF3fxCoREGlKuv7tjh07RhCNRKEQawhEPfGQlJaw6xslAiIBzJw5E4BRo0ZFGoekl641JyHg0sa3b/ldXldm1HTO0lxh1zdKBEQCUCJQ/JocrRB0MqEJc/LalK+Jh6S5lAiIiERMTfkSJ7EaNSAipadugqK6rdDD7C5tMSvxc1RrCIiESYmAiBS1sGeFHN9qVr3HasqXcqdEQEREJMZ0j4BIAE8//XTUIYhITIRd38QqETCzocDQAw44IOpQpMS0bds26hAkLEGXNs7zssZBFWKeAykuYdc3seoacPfZ7j62fXstRyrNM3XqVKZOnRp1GBKG8W/ULGOcvEHj56o/jCQ8TVlc/sKub2LVIiCSrUceeQSAiy66KOJIJB+SJx8CGk8+VOQ0ZXF5C7u+iVWLgMRL4FXftg4LLaZiVDc8DwhleF4xCHskQr6tmTyYNZMHM2tcH9q0qqnG27TaiVnj+iReEwlKiYCUrcBNqNvjnQjUXRTrtlK7KMZNqimLQfMcSPbUNSBlTU2oUm40ZbHkmxIBKWuaKlbC0OSCRyJFTImAlK3Aq761mAVk7lNdsGBBgaKUctDkgkcizRB2faNEQMpW4CbUVrOAGWGFJeUi4LLGUc03IBKUEgGRAG699VYArrjiiogjkaIR9OIeJFlooNSHN0puwq5vNGpAJICnnnqKp556KuowJCZKfXij5Cbs+kaJgIiISIwpERAREYmxWCUCZjbUzKZXV1dHHYqIiEhRiFUioEWHJFsVFRVUVFREHYaIxEDY9U2sEgGRbM2dO5e5c+dGHUZsBF4nIgbrIkj8hF3faPigiBSd2+evTtwp33CdiOQJoZLLFa2U8w38rvFzmm9AIqJEQCSA66+/HoBrr7024kjio2zWiUh1cZ8wByY1uFcpi/kGpDyFXd+oa0AkgPnz5zN//vyow4iVclxqN45LPkvzhV3fqEVARIpO4HUiSmypXa1JIMVILQIiUnS01K5IeJQIiIiIxJgSAZEAdt99d3bfffeow5AcZDMkUX36EoWw6xvdIyASwKxZs6IOQXKUzZBE9elLFMKub5QIiEhslM2QRJE8UiIgEsCVV14JwM033xxxJJKLuiGHdS0CW7buoE2rneq1CKRKFkTCFHZ9o0RAJIDFixdHHYLkqOiHJGoGQqkVdn2jREBEYqHohyRqBkKJiEYNiIiIxJgSARERkRhT14BIAJ06dYo6BBGJibDrm1glAmY2FBh6wAEHRB2KlJj/+q//ijoEKXNTnl1V7/6E5AmPku9jmLJ1GONDj670THl2FbfPX514fOmAA0tmToiw65tYdQ24+2x3H9u+vW62EZHiknzRajjhUfLsh7dvHxZ6bKVo/MAuieGiayYPLpkkIAqxahEQydZll10GwG233RZxJFLONOGRQPj1jRIBkQCWL18edQgSA82d8Ci5+btywpySav6W9MKub5QIiIgUgcATHu08DyadBcB4YHybpIP8qXYDTTwkgSkREBEpAoEnPLp2CjCl6QNq4iEJSImAiIgklPLd9pIdJQIiAXTpoopQ4qFu6eXKCXMS9yxIuMKub5QISCzkelPV9OnTCxWaiEg9Ydc3SgQkFuq+5YiISH2xmlBIJFtjx45l7NixUYchIjEQdn2jFgGRAFatWhV1CCISE2HXN0oEpHxN6Q7VHzZdrv2+hY9FRKRIKRGQ8lX9IUyqjjoKkWi03zf4XAKafCjWlAiIiJSj5lzYNflQrCkREAmgZ8+eUYcgIjERdn2jREAkAK06KCJhCbu+0fBBERGRGFMiIBLA9773Pb73ve9FHYaIxEDY9Y26BkQCWLduXdQhiEhMhF3fqEVAREQkxpQIiIiIxJgSARERkRjTPQIiARx11FFRhyBSOClnIfxd4+c0A2Eowq5vlAiIBHDzzTdHHYJI4aS6uE+Y03iKbs1AGIqw6xt1DYiIiMSYEgGRAIYNG8awYcOiDkNEYiDs+kaJgEgAGzduZOPGjVGHEUtTnl1F5YQ5AFROmMOUZ8Ndq10kbGHXNyV/j4CZ7QxMBb4EFrj7AxGHJCJ5NH5gF8YP7BJ1GCJlqyhbBMxshpl9bGZvNnh+kJmtNLM/m9mE2qdPAx5z9zHASaEHKyIiRaFha9GyDz6p92+6cnFXlIkAMBMYlPyEmbUAfgmcAHQDRphZN6ATsLa22PYQYxQRkSJy+/zViZ+XffAJI+99CYCR975ULxlILidF2jXg7gvNrLLB04cDf3b39wDM7CHgZGAdNcnAcjIkNmY2FhgLsO++++Y/aClrAwYMiDoEEQmg7n6SZFu27mDYtEURRJOdsOubokwE0tibr775Q00CcARwB3CXmQ0GZqfb2d2nA9MBqqqqvIBxShm69tprow5BRAJYM3kw8FWLwJatO2jTaiceOP9IDtuvA5A6WSgmYdc3pZQIWIrn3N2/AM4NOxgRkdhJOQNhmnIRzEB46YADEz8ftl8HHjj/SIZNW1QvCWhYTkorEVgH7JP0uBPwl4hikZg54YQTAJg7d27EkYhEKOjFPaIZCBuOLqm7+CcnAanKFZuw65tSSgSWAAeaWWfgI+BM4KxoQ5K42Lx5c9QhiEhMhF3fFOWoATN7EFgMdDWzdWZ2nrtvA34APAO8Azzi7m8187hDzWx6dXV104VFRERioCgTAXcf4e5fd/dW7t7J3X9d+/zT7t7F3b/h7jdmcdzZ7j62fXstnCEiUkfj7+OtKBMBEREJj8bfx1sp3SMgEpkhQ4ZEHYJIQZXD+PtyEXZ9o0RAJIArrrgi6hBECqocxt+Xi7Drm1h1DehmQREBrWjYUKrx94DG38dErFoE3H02MLuqqmpM1LFIaTn66KMBWLBgQaRxSH5oRcP6ymX8fbkIu76JVSIgIiIhKPIZCKU+JQIiIpJfRT4DodQXq3sEREREpD4lAiIiIjEWq64BMxsKDD3ggAOiDkVKzOmnnx51CCISE2HXN7FKBDRqQLJ10UUXRR2CiMRE2PWNugZEAti0aRObNm2KOgwRiYGw65tYtQiIZOvEE08ENI+AiBRe2PWNWgRERERiTImAiIhIjCkREBERibFY3SOg4YMiIkVEUxEXhVglAho+KNkaNWpU1CGIlB9NRZxS2PVNrBIBkWwpERCRsIRd3+geAZEANmzYwIYNG6IOQ0RiIOz6Ri0CIgEMHz4c0DwCIlJ4Ydc3ahEQERGJMSUCIiJFZsqzq6icMAeAyglzmPLsqogjknKmrgERkSIzfmAXxg/sEnUYEhOxahEws6FmNr26ujrqUERERIpCrFoENI+AZGvcuHFRhyASXzGbeCjs+iZWiYBIts4444yoQxCJr5hNPBR2fROrrgGRbK1du5a1a9dGHYaIxEDY9Y1aBEQCOPvsswHNIyAihRd2faMWARERkRhTIiAiIhJjSgRERERiTImAiIhIjOlmQZEALr/88qhDEJGYCLu+iVUiYGZDgaEHHHBA1KFIiRk6dGjUIYhITIRd38Sqa8DdZ7v72Pbty2PSCQnPypUrWblyZdRhiEgMhF3fxKpFQCRbF1xwAaB5BESKWsqpiH/X+Lkin4o47PpGiYCIiJSHVBf3CXNgUoOF5spkKuJ8iVXXgIiIiNSnREBERCTGlAiIiIjEmO4REAngmmuuiToEEYmJsOsbJQIiARx77LFRhyAiMRF2faOuAZEAli9fzvLly6MOQ0RiIOz6Ri0CIgFcdtllQONxvVu3bmXdunVs2bIlgqiknLVp04ZOnTrRqlWrqEORkKWrbwpFiYBIDtatW0e7du2orKzEzKIOR8qEu7Nx40bWrVtH586dow5Hypy6BkRysGXLFnbffXclAZJXZsbuu++uliYJRaxaBLTokBSCkgApBH2uCijlVMQZyhbxdMT5EKtEwN1nA7OrqqrGRB2LSL7ccccdTJs2jUMPPZQzzjiDt99+mwkTJvDEE0/QpUsXunXrFnWIIsWlORf2GExHHKtEQCRbN910U9QhpDV16lTmzp2b6Es+6aSTAHjiiScYMmRIk4nAtm3baNlSVYFIsQi7vtH/fpEA+vTpE3UIKV144YW89957nHTSSYwePZoOHTqwdOlSzjrrLJ588kleeOEFbrjhBmbNmsU3vvGNxH6jRo1it91249VXX020JFx22WVs3ryZiooKfvOb39C1a1dmzpzJk08+yaZNm3j33Xc59dRT+dnPfgbAr3/9a/7zP/+TvfbaiwMPPJDWrVtz1113sX79ei688EI+/PBDAG677Tb69u0bye9HpBSFXd8oERAJYNGiRUDT/0GPPvroRs+dfvrpXHTRRWzatIkTTzyx0eujRo1i1KhRbNiwgeHDh9d7ranhQ3fffTd/+MMfeP755+nYsSMzZ85MxHnSSScxZMiQRsess2rVKubNm0eLFi349NNPWbhwIS1btmTevHlcddVVzJo1C6gZ0/zqq6/SunVrunbtyiWXXEKLFi24/vrreeWVV2jXrh3HHHMMhxxyCACXXnop48eP59vf/jYffvghxx9/PO+8807G9yEiXwla3+SLEgGRAK666iogvHG9Yfjud79LixYtAKiuruacc85h9erVmBlbt25NlBswYADt29f0k3br1o0PPviADRs20L9/f3bbbbfEsVatWgXAvHnzePvttxP7f/rpp3z22We0a9curLcmUtLCrm+UCIjkUab/uG3bts34eseOHUNNNHbeeefEz9deey3f+c53ePzxx1mzZk29lo3WrVsnfm7RogXbtm3D3dMed8eOHSxevJiKioqCxC0i+aV5BETKVLt27fjss88Cla2urmbvvfcGSHQvZHL44Yfzwgsv8Mknn7Bt27ZENwLAcccdx1133ZV4rKmZRYqbEgGRMnXmmWdyyy230KtXL959992MZX/84x9z5ZVX0rdvX7Zv397ksffee2+uuuoqjjjiCI499li6deuW6D644447WLp0KT169KBbt27cfffdeXk/IlIYlqmJr1xVVVX50qVLow5DCm1Se5hUnZdD1TWVN2y6f+eddzjooIPyco5S8/nnn7PLLruwbds2Tj31VEaPHs2pp54adVhlJcrPV+WEOayZPDiSc+dTzu9jSneo/rDpcnmceChdfZMLM1vm7lWpXtM9AiIB3HbbbVGHUHQmTZrEvHnz2LJlC8cddxynnHJK1CFJBlOeXcXt81cnHl864EDGD+wSYUQlIujFPY8TD4Vd3ygREAmgZ8+eUYdQdG699daoQ5BmGD+wC+MHdimbb/rlLOz6RvcIiAQwb9485s2bF3UYIhIDYdc3ahEQCeCGG24A4Nhjj404EhEpd2HXN2oREBERiTElAiIiIjGmRECkxLVo0YKePXty8MEH893vfpdNmzZlfaxRo0bx2GOPAXD++efXmyq4oQULFiTmRG+OyspKNmzYEKjsT3/604L0lTYnhnyaOXMmf/nLX0I/r0gmukdAJJ+CjjkOKsDY5IqKisTsfSNHjuTuu+/mhz/8YeL17du3J9YUaI5777034+sLFixgl112KejCKP/xH/9RsGNHYebMmRx88MHstddeUYcSmnwPW4xyGGTZDsF099hswFBg+gEHHOASAxN3zduhVqxY4StWrGj0/Ntvv12wcwY93s4775z4edq0aT5u3Dh//vnn/eijj/YRI0b4QQcd5Nu2bfMrrrjCq6qqvHv37n733Xe7u/uOHTv84osv9oMOOshPPPFEP+GEE/zRRx91d/f+/fv7kiVL3N197ty53qtXL+/Ro4cfc8wx/v777/u//uu/+l577eWHHHKIL1y40D/++GM/7bTTvKqqyquqqvzFF190d/cNGzb4wIEDvWfPnj527Fjfd999ff369fXew7Zt2/ycc87xb33rW37wwQf7L37xC3d3P+eccxLxzJkzx7t27ep9+/b1Sy65xAcPHlzzK5o40c8991zv37+/d+7c2W+//fbEcU8++WQ/9NBDvVu3bv6rX/0q8fx+++3XKIZU79PdfePGjX7yySd79+7d/YgjjvDXXnstcd5bbrklse+3vvUtf//99/3999/3b37zm37++ed7t27dfODAgb5p0yZ/9NFHfeedd/YuXbr4IYcc4ps2bfKf/OQnftBBB3n37t398ssvbxRPo89XHuz3k6fyWi7f58338fJ93rTHDKG+yQWw1NNcG2PVIuDus4HZVVVVY6KORUpL165dow6hSdu2bWPu3LkMGjQIgJdffpk333yTzp07M336dNq3b8+SJUv45z//Sd++fTnuuON49dVXWblyJW+88QZ/+9vf6NatG6NHj6533PXr1zNmzBgWLlxI586d+fvf/85uu+3GhRdeyC677MIVV1wBwFlnnZVy+eHrrruOb3/72/z0pz9lzpw5TJ8+vVHsy5cv56OPPuLNN98E4B//+Ee917ds2cIFF1yQiGHEiBH1Xl+xYgXPP/88n332GV27dmXcuHG0atWKGTNmsNtuu7F582Z69+7NsGHD2H333VP+/lK9T4CJEyfSq1cvnnjiCZ577jm+//3vN7l+wurVq3nwwQe55557OP3005k1axbf+973uOuuu7j11lupqqri73//O48//jgrVqzAzBq9Zykx7fcNNqlQgFa+sOubWCUCItmaPXs2AEOHDo04ksY2b96cmIDk3//93znvvPNYtGgRhx9+OJ07dwbgj3/8I6+//nqi/7+6uprVq1ezcOFCRowYQYsWLdhrr7045phjGh3/pZdeol+/folj1S093FC65YcXLlzIf//3fwMwePBgOnTo0Gjf/fffn/fee49LLrmEwYMHc9xxx9V7fcWKFey///6JGEaMGFEvoaFc4AMAACAASURBVBg8eDCtW7emdevW7Lnnnvztb3+jU6dO3HHHHTz++OMArF27ltWrV6dNBNK9zxdffDGxqNIxxxzDxo0bqa7OPHV1586dE3+Tww47jDVr1jQqs+uuu9KmTRvOP/98Bg8ezJAhQzIeU4pcHmcgDLu+USIgEsDPf/5zoDgTgeR7BJIlLzPs7tx5550cf/zx9co8/fTTmFnG47t7k2Ug8/LDTe3foUMHXnvtNZ555hl++ctf8sgjjzBjxox6MWSSaqnkBQsWMG/ePBYvXkzbtm05+uij2bJlS9pjpHufqc5tZrRs2ZIdO3Yknks+dsN4Nm/e3OgYLVu25OWXX2b+/Pk89NBD3HXXXTz33HMZ36fEQ9j1jUYNiMTA8ccfz7Rp09i6dSsAq1at4osvvqBfv3489NBDbN++nb/+9a88//zzjfY96qijeOGFF3j//fcBEk3mDZc5Trf8cL9+/XjggQcAmDt3Lp988kmjc2zYsIEdO3YwbNgwrr/+el555ZV6r3/zm9/kvffeS3yzfvjhh5t8z9XV1XTo0IG2bduyYsUKXnrppYzl073P5PgXLFhAx44d2XXXXamsrEzE+corryT2yyT5d/b5559TXV3NiSeeyG233ablmiUyahEQiYHzzz+fNWvWcOihh+Lu7LHHHjzxxBOceuqpPPfcc3Tv3p0uXbrQv3//RvvuscceTJ8+ndNOO40dO3aw55578uyzzzJ06FCGDx/O73//e+68807uuOMOLr74Ynr06MG2bdvo168fd999NxMnTmTEiBEceuih9O/fn3333bfROT766CPOPffcxDfsm2++ud7rFRUVTJ06lUGDBtGxY0cOP/zwJt/zoEGDuPvuu+nRowddu3blyCOPzFg+3fucNGkS5557Lj169KBt27bcd999AAwbNoz777+fnj170rt3b7p0afru8VGjRnHhhRdSUVHB3LlzOfnkk9myZQvuzpQpU5rcX6QQtAyxlK8oliGOYPhgXNQte+zuXHzxxRx44IGMHz8+6rAKqhDLEAdddCjfixNFdbxCLLKU0zED1EtahliklOmiXTD33HMP9913H19++SW9evXiggsuiDokkbKgREAkgN/+9rdRhxB748ePL/sWABEIv75RIiASwD777BN1CCISE2HXNxo1IBLAww8/HOhOdRGRXIVd36hFQCSAadOmAXDGGWdEHImIlLuw6xslAiIhKttFS0SkZKlrQCRE4wd2SQw7WjN5cF6SgLpliOu2yZMn53zMYvfEE09kXCJZRIJTi4CUnqBj9ds3nrimHKWbYricPfHEEwwZMoRu3bpFHYpI8wRZnGjNF9CydeYyeaQWASk91R/WTMjR1FYkY/qnPLuq3uNlH3xS79905XJRXV1N165dWblyJVCzSM8999wDwP3330+PHj045JBDOPvss4GalfeGDRtG79696d27N3/6058AeOGFFxItDb169eKzzz7jr3/9K/369aNnz54cfPDB/M///E+j8y9ZsoQ+ffpwyCGHcPjhh/PZZ5+xZcsWzj33XLp3706vXr0S0xnPnDmTH/zgB4l9hwwZkphIZZddduHqq6/mkEMO4cgjj+Rvf/sbixYt4sknn+RHP/oRPXv25N133+WOO+6gW7du9OjRgzPPPDNvv0eRvBv/RtN1V+W3Yds/QwtJLQIiAdSt2peN2+evTnQBLPvgE0beWzPn/ch7X+KB84/ksP06NCrXHMmrDwJceeWVnHHGGdx1112MGjWKSy+9lE8++YQxY8bw1ltvceONN/KnP/2Jjh07JubTv/TSS1MuIXzrrbfyy1/+kr59+/L555/Tpk0bpk+fzvHHH8/VV1/N9u3b2bRpU714vvzyS8444wwefvhhevfuzaeffkpFRQW33347AG+88QYrVqzguOOOY9WqzMnPF198wZFHHsmNN97Ij3/8Y+655x6uueYaTjrpJIYMGcLw4cMBmDx5Mu+//z6tW7fWcr5S8h577DH4WefQzqdEQCSAjh075rR/5YQ5jZ7bsnUHw6Ytyum4kL5rYODAgTz66KNcfPHFvPbaawA899xzDB8+PPF+6pbaTbeEcN++ffnhD3/IyJEjOe200+jUqRO9e/dm9OjRbN26lVNOOaVeEgKwcuVKvv71r9O7d2+gZrldqFnO95JLLgFqFhHab7/9mkwEvva1ryWW5z3ssMN49tlnU5br0aMHI0eO5JRTTuGUU07J/AsTKXIdO3aEtuE12KtrQCSAmTNnMnPmzKz3XzN5MGsmD2bWuD60aVXz365Nq52YNa5P4rV827FjB++88w4VFRWJb/7pltqtW0J4+fLlLF++nI8++oh27doxYcIE7r33XjZv3syRRx7JihUr6NevHwsXLmTvvffm7LPP5v777693rOYs5wtkXM63VatWiWPVLS+cypw5c7j44otZtmwZhx12WNpyIqVg5syZzFz+ZWjnUyIgEkAuicClAw5M/HzYfh144PyaVfCSuwUalsuHKVOmcNBBB/Hggw8mvsEPGDCARx55hI0bNwJfLbWbbgnhd999l+7du/OTn/yEqqoqVqxYwQcffMCee+7JmDFjOO+881IuGfyXv/yFJUuWAPDZZ58lViOsW8531apVfPjhh3Tt2pXKykqWL1/Ojh07WLt2LS+//HKT7y15Od+6/b7zne/ws5/9jH/84x98/vnnOf72RKJTkwhsDe186hoQKbCG/f51F//kJCBVuaAa3iMwaNAgRo8ezb333svLL79Mu3bt6NevHzfccAPXXXcdV199Nf3796dFixb06tWLmTNnpl1C+LbbbuP555+nRYsWdOvWjRNOOIGHHnqIW265hVatWrHLLrs0ahH42te+xsMPP8wll1zC5s2bqaioYN68eVx00UVceOGFdO/enZYtWzJz5kxat25N37596dy5M927d+fggw/m0EMPbfI9n3nmmYwZM4Y77riDhx56iPPOO4/q6mrcnfHjx/Mv//IvWf0uReJIiYBIidu+fXvK5995553Ez7/4xS8SP59zzjmcc8459cp27Ngx5ZSmd955Z6PnUu3fUO/evXnppZcaPZ+qVcXMEi0FDSV/sx8+fHji5sC+ffvWu6fhxRdfzBiPiKSnrgGREE15dlXixsHKCXPyOmRQRCQbahEQCdH4gV00pbCIFBUlAiIBPP3001GHICIx8fTTT8ON/xba+ZQIiATQtm3btK+lGy4nkot0wy2l/LVt2xZahVen6B4BkQCmTp3K1KlTGz3fpk0bNm7cqEpb8srd2bhxI23atIk6FInA1KlTmbokvHkE1CIgEsAjjzwCwEUXXVTv+U6dOrFu3TrWr18fRVhSxtq0aUOnTp2iDkMi8Mgjj8CarVzUdNG8KPlEwMz2B64G2rv78KjjkXhp1aoVnTuHNye4iEi+Rdo1YGYzzOxjM3uzwfODzGylmf3ZzCZkOoa7v+fu5xU2UhERkfIUdYvATOAuIDE1mZm1AH4JDATWAUvM7EmgBXBzg/1Hu/vH4YQqIiJSfiJNBNx9oZlVNnj6cODP7v4egJk9BJzs7jcDQ8KNUEREpLxZ1Hc71yYCT7n7wbWPhwOD3P382sdnA0e4+w/S7L87cCM1LQj31iYMqcqNBcbWPuwKrGxQpD1QnSHUTK93BDZk2LeYNfW+i/l8uRwrm32D7hOknD5vpXc+fd5Kjz5vX9nP3fdIWdrdI92ASuDNpMffpeaCXvf4bODOEOKYnu3rwNKof4+Fet/FfL5cjpXNvkH3CVJOn7fSO58+b6W36fMWbCvGeQTWAfskPe4E/CWE887O8fVSFfb7yuf5cjlWNvsG3SdIOX3eSu98+ryVHn3eAijGroGWwCpgAPARsAQ4y93fiirGppjZUnevijoOiQd93iRM+ryVv6iHDz4ILAa6mtk6MzvP3bcBPwCeAd4BHinmJKDW9KgDkFjR503CpM9bmYu8RUBERESiU4z3CIiIiEhIlAiIiIjEmBIBERGRGIt6iuFIdOzY0SsrK6MOQ0REJBTLli3b4GkmFIplIlBZWcnSpUujDkNERCQUZvZButfUNSAiIhJjSgRERERiTImAiIhIjJVMImBmM8zsYzN7M+m53czsWTNbXftvhyhjFBERKTUlkwgAM4FBDZ6bAMx39wOB+bWPRUREJKCSSQTcfSHw9wZPnwzcV/vzfcApoQYlIiJS4komEUjjX939rwC1/+4ZcTwiIiIlpdQTgcDMbKyZLTWzpevXr486HBERkaJQ6onA38zs6wC1/36crqC7T3f3Knev2mOPlJMriUgGkyZNwswS26RJk6IOSUTyoKSWITazSuApdz+49vEtwEZ3n2xmE4Dd3P3HTR2nqqrKNbOgSHbMjFKqN0QEzGyZu1eleq1kWgTM7EFgMdDVzNaZ2XnAZGCgma0GBtY+Fik4fTsWkXJRUi0C+aIWAcmXOH47juN7Fil1ZdEiICLl1RJRTu9FpJSpRUAkB1F9O47yW3m+z60WBpHCU4uAiIiIpKREQEREJMaUCIiIiMSYEgEREZEYUyIgIiISYy3zdSAz+2Gm1939F/k6l4iIiORH3hIBoF0ejyUiIiIhyFsi4O7X5etYIiIiEo683yNgZl3MbL6ZvVn7uIeZXZPv84iIiEjuCnGz4D3AlcBWAHd/HTizAOcRyTtNeysicZP3KYbNbIm79zazV929V+1zy929Z15PVP+ca4DPgO3AtnTTKNbRFMPSlKDT3mqK4eI7nog0FvYUwxvM7BuA1558OPDXApynoe+4e8+mkgARKU9RtuaoJUlKWSFaBPYHpgN9gE+A94HvufuavJ6o/jnXAFXuviFIebUISFPUIhDeuYv9eKVybpFM8toiYGY7mdmu6V539/fc/VhgD+Cb7v7tQiYBdacF/mhmy8xsbKoCZjbWzJaa2dL169cXOByR8tDwm+3ixYvr/ZuunIiUjkCJgJn9zsx2NbOdgbeBlWb2owZlfpi8ARcAY5IeF1Jfdz8UOAG42Mz6NSzg7tPdvcrdq/bYY48ChyNSHq677qtRwYsXL2bAgAEADBgwoF4ykFwuk7gnFupCkGIUqGug7mY/MxsJHAb8BFjm7j2Sykys/bEr0Bt4svbxUGChu5+f18jTxzoJ+Nzdb01XRl0D0pR8dQ1MmjSp3kVy4sSJean8w2qCNrPAZZv7+6pLLDZv3kxFRQXz58/nqKOOalQum5iLvWtAXQgStnx0DbQys1bAKcDv3X0rtTcD1nH362onFeoIHOrul7v75dQkDp2yDz8zM9vZzNrV/QwcB7xZqPNJaYrqm9ikSZMSFb67pz1vMX9TdHfcnUWLFlFRUQFARUUFixYtSrzWHHXvsU+fPmzevBmAzZs306dPn8RrEkwxf26khNT9R860Af8P+Ah4GjBgP+B/0pRdAbROetwaWBHkPNlswP7Aa7XbW8DVTe1z2GGHucRTzUc+HuUmTpzo1CTsDvjEiRMDHbvhMZItWrTIAV+0aFHacpnOmxz3okWLvKKiwgGvqKiod8ygv4dUctk3V/n+O+f7vBJfwFJPc03MetSAmbV0920pnr8aOB14vPapU4CH3f3mrE5UAOoaiK98N90We7nmls3n8VKVmzRpUr1vrYsXL6ZPnz4sWrQo0S2Qqlwh4iuEqLoG1NUgTcnUNZAxEch2RUEzOxT4d2q+EfyPu78aPNzCUyIQX8V+4S73RCCXckHl83jNvcejHN6zlKdMiUBTiw5lu6LgdmAHNYnAjiyPISISqbqWCV1opZxlvFnQa28ATLel2sfMLgUeoOamwT2B/zKzS/Ifuki4ohr6FvS8hTi3iMRAupsHkjegDXAxMBWYUbelKfs6sHPS452B14OcJ6xNNwvGFznczEUON7qFcd5czx1ULu8ll3JB5ft46Y6ZzY2UmY6Xi0K8Zykv5HqzoJk9Ss1ogLOA/wBGAu+4+6Upyr4B9Hb3LbWP2wBL3L17M/KTgtI9AvGVSx9uLmPqwzpvLucOKk73CGQ6Zi5zIpTCe5byko95BA5w92uBL9z9PmAwkO7C/hvgf81sUu3kPi8Bv25mzCI5K0RTfl0Gna8x9fk8ry4E4ctmToRcP4dxn51RCiBdU0GDyuXl2n8XAgdT0///Xobyh1Iz98ClQK8g5whzU9dAPJDnpvwwmoJzOW+u5w4ql/eSS7lUcv3dBNXUZyTo5ysfcyfk+3gSD2ToGgiaCJwPdAD6A+8BHwMXZijfAehRmxAcSs1Mg5EnAHWbEoF4IGlSm6a2ZGFd4MMo19yy+TxeGIlA0ItiPs9TJ2gSku3nMFMs+TyexEPOiUBzNuB6YC2wAHi+dnsu3+fJZVMiEA9hfGNTIhBeuXT7hnFhzOW9ZPP5inp2xnKQj5k1y0k+WgR+mmpLU3Yl8LUgx41qUyIQD2F8Y1MiEF65pvYNu0UgaLl8dynlcrw4intCVCdTIhB01MDlSQ/bAEOoGTUwOkXZWcA4d/+4yQNHRKMG4iuqu781s2Du5VIJOmVxqrLNUex/Z40aSE+/mxqZRg1k9Y2amoWEnknzWhU1CxQ9Q81SxE8CT2ZznkJtahGIL0L4xhb0vGGUa27ZfB6vEO8ljOOV++ehKaXcpK7WkvTI96JDZtaBmpEEB6Z47S3gV8AbJE0v7O4vNPtEweMZBNwOtADudffJmcqrRSC+iv0bm1oEcpfr8cq9hSioUvwmncvfrtzlPI+Amb1hZq/Xbm9Rcx/A7WmKb3D3O9z9eXd/oW7LMvYgsbUAfgmcAHQDRphZt0KdL5VU43pvvvnmRuN6jz766EDlmlM2buVA46Mld3XrB9RtDT9T2cwPIMVBf7sspGsqSN6A/ZK2vYGWGcr+ArgZOIoQhg/WnueZpMdXAldm2iffXQOkuGGpRYsWge5OT1WuEMcsl3INy2ZqxizFptt8l2tu2WyPV6rT7aYqm+qzSBmNIgmqEMcstGz/dkG7Q4LWN5nK9e/fv9nHywcydA1kXH3QzHar/fGzBi/tWtu08vcUu/Wq/ffIpOccOCbTuXKwNzXDFeusA45oWMjMxgJjAXbfffe8f6tMlWXWZaHZlCvEMculHHzVKrB27VpatmzJtm3baNmyJWvXrk28lrx8bKp9m1Iu5ZpbNpvjJf+u165dy/333w9Av379+P73v88+++zTqFyY8TWnbP/+/es9N2LECGbMmMGIESN45plneOaZZ1KWa+65G5ZbsGBBvZaxtWtrqrXzzjsv8ftLVa65500lm3MvWLCAF174qrG3f//+KeMKs1zQv91+++3X6Pc1evRoZsyYwejRNffA170+c+ZMRo0a1WS5hv8H0tVLL7zwQk71V0GkyxBqEgjep2YCofepWVp4A7Cx9uf3M+0b1gZ8l5r7Auoenw3cmWkftQgUd7nRo0fXy5JHjx7d6HcTdGso1XOplEu55pbN9njZ/E0KdVNaIX43YZRL9X+FPLdEBIknzFaQoC1Jyd+iM5UL2uLU8Ll8zTOS7y2fyMM8AncDJyY9PgH4eZB9C71RBF0DqT7MN910U6APc6pyzSkbt3Lu2VVazb3oFNMFIpdyQcrm43eTy4Uk30o5EShkMtXUJEXZXJxyvXBn87nJR6JSiIt2GO8lF/lIBJaleC7tQcPcgJbUtFp0Br4GvAZ8K9M+Gj5Y2oL2xeV7bvlSKFfIvvrm/K7zOWyrUElcc8rmK5kqhhkD85XIFfu36KC/60K851K7RyDoxfYZ4BqgkpobBq8mzTwCUWzAicAq4F3g6qbKKxGQphTbBT5ouVwuJgX4BhLJ8aJIBPIhrDHwqd5L0HMX8sId5rfohu+3GC/c+ZaPRGA3aoYLvlq73Q7sFmC/6UGOH/amRECaUmwX+KDlsqmAm3vuoJQI5C6X86a6MOWzb70QF+5C3CMgNXJOBLLdgFcKefxsNyUC0pQwLtyFaMZXi4ASgXT75uPzUMib+/LRvSLpZUoEMs4saGa3uftlZja79pdej7uflHbnmv3/4O6DMpWJgmYWlKaEMfNbIWawCzr3fqp594tthr9sj1eIWRejmokul/M2d+KcoJ+xVOcp5t+h1MhlZsHf1v57K/DzFFtGxZgEiBRaqpkmk/9Nlu9Z0Bqeu+7in2oBnrgJ+ncpp99N8re+RYsWUVFRAUBFRQWLFi1Kbr1NqJt1EUg562Jz5ft4UgDpmgrSbUAHoEdz9yumTV0D0hRyaDImi+b5Qjbj5/JechH28ZpqMg7r951PuZw313sEghw/0+9bigt5WIZ4AXASNUP1lgPrgRfc/Yc5ZyIRUNeANCXXJv+gkvctVDN+qS7+k+/jZft3yce5s6XFhCRfcl50CGjv7p8CpwG/cffDgGPzFaBIuanLtIM2x4Ka8cPQ3L9LVM3aak6XUKVrKkjeqFlS+OvAH4Hetc+9nqZsG+BiYCowo24Lcp6wNnUNSFMostEA2ZYr1DFL8XhhjdMvdvn+u0hpIEPXQNAWgf+gZlKhd919iZntD6xOU/a3wL8BxwMvAJ1ovGiRSFHKxzcx3bBXnOL+d1Erg6SVLkPIdgNe9aQWA6AV8Fy+z5PLphYByRdCnlAomxu08h1jULkcL8pZ9kTKEXm4WbALMA34V3c/2Mx6ACe5+w0pyr7s7oeb2ULgIuD/gJfdff9cEpZ80s2Cki+lMNa6FG8WzGWOhbBiFCkl+bhZ8B5qVvXbCuDurwNnpik73cw6ULM2wZPA28B/NitikSKnZtbCy/ccCyKSWtBEoK27v9zguW1pys5390/cfaG77+/ue1Jzk6FI2Zg0aVK9pjUlAvlX97ttzsgLEWm+oInABjP7BrXTDJvZcOCvacrOSvHcY1nE1iQzm2RmH5nZ8trtxEKcR6TcFPssexMnTkz8fNRRRzF//nyAet0CDcuJSHaC3iOwPzAd6AN8ArwPjHT3D5LKfBP4FvAz4EdJu+8K/Mjdv5XHuOvOOQn43N1vbc5+ukdAwlZs9wgUsg9+0qRJXHfddYnHEydOzEtCocl1RLKX6R6BQIlA0oF2pqYVYTNwhrs/kPTaycAp1MxA+GTSbp8BD7n7oixibyqeSSgRkBJQjIlAUMVyoVQiIJK9rG8WNLNdzexKM7vLzAYCm4BzgD8DpyeXdfffu/u5wBB3Pzdp+3+FSAKS/MDMXjezGbU3KaZ7L2PNbKmZLV2/fn0BwxEpDeqDFxEItvpgV2pmFhxDzU1/3wVOcfeT0+zzqpldbGZTay/OM8xsRrYBmtk8M3szxXYyNUMavwH0pOaehbQrIrr7dHevcveqPfbYI9twRMpC3PvgNepD5CsZuwbM7A137177cwtgA7Cvu6edKdDMHgVWAGdRMyPhSOAdd780n4GnOG8l8JS7H9xUWXUNSNiKrWsgrHPnU7HHJ1LMcplHYGvdD+6+HXg/UxJQ6wB3vxb4wt3vAwYD3ZsTcFBm9vWkh6cCbxbiPCIiIuWqZROvH2Jmn9b+bEBF7WOjZmrOXVPsU5c8/MPMDqZmZsHKfASbws/MrCc1wxrXABcU6DwiIiJlKWOLgLu3cPdda7d27t4y6edUSQB8NbPgtXw1s+DP8hx3XXxnu3t3d+/h7ie5e7q5DURiR/3gIhJEs4YPlgvdIyBhK4V++mLvgy/2+ESKWaZ7BJrqGmjOSX6Y6XV3/0W+ziUiIiL5EXSK4SDa1W5VwDhg79rtQqBbHs8jIjGiLg6Rwsp714CZ/REYVje6wMzaAY+6+6C8nigH6hqQsKlrQESilI9liJtjX+DLpMdfUrhRAyIiIpKDvN0jkOS3wMtm9jg1w/pOBe4rwHlEREQkR3lPBNz9RjObC/x77VPnuvur+T6PiIiI5K4QLQK4+yvAK4U4toiIiORPIe4REJFm0p3xIhIVTSgkEoJSuCO/FGIUkeyEPWpARERESoQSARERkRgriUTAzL5rZm+Z2Q4zq2rw2pVm9mczW2lmx0cVo0gqpdD3XwoxikjhlMQ9AmZ2ELAD+BVwhbsvrX2+G/AgcDiwFzAP6OLu2zMdT/cIiIhInJT8PQLu/o67r0zx0snAQ+7+T3d/H/gzNUmBiIiIBFASiUAGewNrkx6vq32uETMba2ZLzWzp+vXrQwlORESk2BVkQqFsmNk84N9SvHS1u/8+3W4pnkvZ1+Hu04HpUNM1kFWQIiIiZaZoEgF3PzaL3dYB+yQ97gT8JT8RiYiIlL9S7xp4EjjTzFqbWWfgQODliGMSEREpGSWRCJjZqWa2DjgKmGNmzwC4+1vAI8DbwB+Ai5saMSAiIiJfKZqugUzc/XHg8TSv3QjcGG5EIiIi5aEkWgRERESkMJQIiIiIxJgSARERkRhTIiAiIhJjJbHWQL6Z2XrggwZPtweqM+yW6fWOwIY8hBaFpt53MZ8vl2Nls2/QfYKU0+et9M6nz1vp0eftK/u5+x4pS7u7tppkaHq2rwNLo46/UO+7mM+Xy7Gy2TfoPkHK6fNWeufT5630Nn3egm3qGvjK7BxfL1Vhv698ni+XY2Wzb9B9gpTT5630zqfPW+nR5y2AWHYN5JuZLfU0yzuK5Js+bxImfd7Kn1oE8mN61AFIrOjzJmHS563MqUVAREQkxtQiICIiEmNKBERERGJMiYCIiEiMKREQERGJMSUCBWZm+5vZr83ssahjkfJkZjub2X1mdo+ZjYw6HilvqtPKjxKBDMxshpl9bGZvNnh+kJmtNLM/m9mETMdw9/fc/bzCRirlppmfvdOAx9x9DHBS6MFKyWvO5011WvlRIpDZTGBQ8hNm1gL4JXAC0A0YYWbdzKy7mT3VYNsz/JClTMwk4GcP6ASsrS22PcQYpXzMJPjnTcpMy6gDKGbuvtDMKhs8fTjwZ3d/D8DMHgJOdvebgSHhRijlqjmfPWAdNcnAcpTcSxaa+Xl7O9zopNBUaTTf3nz17QtqKuG90xU2s93N7G6gl5ldWejgpKyl++z9NzDMzKZRvnPGS/hSAmCtzAAAA6tJREFUft5Up5UftQg0n6V4Lu30jO6+EbiwcOFIjKT87Ln7F8C5YQcjZS/d5011WplRi0DzrQP2SXrcCfhLRLFIvOizJ2HS5y0mlAg03xLgQDPrbGZfA84Enow4JokHffYkTPq8xYQSgQzM7EFgMdDVzNaZ2Xnuvg34AfAM8A7wiLu/FWWcUn702ZMw6fMWb1p9UEREJMbUIiAiIhJjSgRERERiTImAiIhIjCkREBERiTElAiIiIjGmREBERCTGlAiIxIyZbTez5UlbxqW0w5IU115m9r+1P39oZuuTYq1ssM/RZra4wXMtzexvZvZ1M7vFzP7PzK4I872IlBKtNSASP5vdvWc+D2hmLWsnoMlFclxH1B53FFDl7j9Is89CoJOZVbr7mtrnjgXedPe/Aj8ysy9yjEukrKlFQEQAMLM1Znadmb1iZm+Y2Tdrn9/ZzGaY2RIze9XMTq59fpSZPWpms4E/mllbM3vEzF43s4drv9VXmdl5ZjYl6TxjzOwXWcT3DTP7g5ktM7P/MbNvuvsO4FHgjKSiZwIP5vTLEIkRJQIi8VPRoGsg+SK6wd0PBaYBdc3pVwPPuXtv4DvALWa2c+1rRwHnuPsxwEXAJ+7eA7geOKy2zEPASWbWqvbxucBvsoh7OnCJux9WG9vU2ucfpObij5m1Bk4EZmVxfJFYUteASPxk6hr479p/lwGn1f58HDUX8rrEoA2wb+3Pz7r732t//jZwO4C7v2lmr9f+/IWZPQcMMbN3gFbu/kZzAjazXYA+wKNmidVxW9cef4mZ7WJmXYGDgJfc/ZPmHF8kzpQIiEiyf9b+u52v6gcDhrn7yuSCZnYEkNz/nmr9+jr3AlcBK8iuNWAn4B8ZEpiHqGkVOAh1C4g0i7oGRKQpzwCXWO1XcTPrlabci8DptWW6Ad3rXnD3/6VmbfuzyOJC7e6fAu+b2Xdrj29mdkhSkQeB7wHHoKVyRZpFiYBI/DS8R2ByE+WvB1oBr5vZm7WPU5kK7FHbJfAT4HWgOun1R4A/5dBsPxI4z8xeA94CTq57wd3fBjZRcy+DRgmINIOWIRaRvDCzFtT0/28xs28A84Eu7v5l7etPAVPcfX6a/T93910KENck4HN3vzXfxxYpB2oREJF8aQu8WPuN/XFgnLt/aWb/YmarqLlJMWUSUOvTugmF8hWQmd1CTZeBWglE0lCLgIiISIypRUBERCTGlAiIiIjEmBIBERGRGFMiICIiEmNKBERERGLs/wOc7KoLS3eEDQAAAABJRU5ErkJggg==\n",
"text/plain": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"ax_spectrum, ax_residuals = datasets[0].plot_fit()\n",
"ax_spectrum.set_ylim(0.1, 40)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For more ways of assessing fit quality, please refer to the dedicated [modeling and fitting tutorial](modeling.ipynb)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Compute Flux Points\n",
"\n",
"To round up our analysis we can compute flux points by fitting the norm of the global model in energy bands. We'll use a fixed energy binning for now:"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [],
"source": [
"e_min, e_max = 0.7, 30\n",
"energy_edges = np.logspace(np.log10(e_min), np.log10(e_max), 11) * u.TeV"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we create an instance of the `~gammapy.estimators.FluxPointsEstimator`, by passing the dataset and the energy binning:"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [],
"source": [
"fpe = FluxPointsEstimator(energy_edges=energy_edges, source=\"crab\")\n",
"flux_points = fpe.run(datasets=datasets)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here is a the table of the resulting flux points:"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"Table length=10\n",
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.figure(figsize=(8, 5))\n",
"flux_points.table[\"is_ul\"] = flux_points.table[\"ts\"] < 4\n",
"ax = flux_points.plot(\n",
" energy_power=2, flux_unit=\"erg-1 cm-2 s-1\", color=\"darkorange\"\n",
")\n",
"flux_points.to_sed_type(\"e2dnde\").plot_ts_profiles(ax=ax);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The final plot with the best fit model, flux points and residuals can be quickly made like this: "
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [],
"source": [
"flux_points_dataset = FluxPointsDataset(\n",
" data=flux_points, models=model_best_joint\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAgwAAAGtCAYAAAB6LizZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdd5zdVbX38c8+Z8qZXjMhvfdCGwidBEFqQEoAwYIi4LVd21Vs2O6j6HPv9epz9QoqKApSAxJ6kyAtkACSSigJkDop03vZzx/rDHNmMpl6fmfa9/167VcyZ37nt3+hZNbZe+21nPceERERka6EBvoBREREZPBTwCAiIiLdUsAgIiIi3VLAICIiIt1SwCAiIiLdShroBxjMCgsL/eTJkwf6MURERBJizZo1e733ozr7ngKGLkyePJnVq1cP9GOIiIgkhHPu3YN9T1sSIiIi0i0FDCIiItItBQwiIiLSLQUMIiIi0i0FDCIiItItBQydcM4tdc7dWF5ePtCPIiIiMigoYOiE936F9/7qnJycgX4UERGRQUEBg4iIiHRLAYOIiIh0SwGDiIiIdEsBg4iIiHRLAUMXRkWaoKlhoB9DRERkwClg6EJ+ShOUbIDSrdBYO9CPIyIiMmDUrbJbHmpLbaRmQ2YRpGYN9EOJiIgklAKGLvn2X9ZX2EjOsMAhLXdgHktERCTBtCXRhZ8csRveeBia6tt/o7EaSrdAyUao3gfed34DERGRYcJ5/bA7qPf/rdBPyGiEtDyYex7MORfS8w+8MJQMGaMgoxBC4cQ/qIiISBw459Z474s7+55WGLpw7ZrRcNZ/wKhZsOaPcNsl8PTPYN/b7S9saYTKHbB7PZRvh+bGAXleERGRoCiHoUsOxhfbKHsP1i2HzY/A5odh7BGw4CKYeAy4aNzlm6G6BKr32EpERhEkRwb2jyAiIhIHChh6KncinPBlKP40bHoA1i+HR78N2eNg/oUw6wxITo9e7KFmn41IjgUOqZkD+vgiIiL9oYChtyLZcNhlsPBieGclrL0bnv8VrL4JZp8D88+HzNFt19eV20jJtDwHnawQEZEhSEmPnXDOLQWWHj1v6lWrnnoAmuoOfrH3Vtxp7V2w5Rl7bcrJtl0xet6B1ydFbMUhPR+cC+T5RURE+qKrpEcFDF0oLi72q1evhvpKqN5rKwUdazPEqtwF6++1LYuGaiiaAwuWwZSTINRhMUcnK0REZJBRwNBHHwQMrZob23ITmrvoMdFYA5sfhXX3QPk2W1GYdz7MOefAKpEuDOkFVggqnBzMH0RERKQHFDD00QEBQyvvbbWhZp9VfjwY3wLvvWh5Djtese2ImWdYkmTuBACuXb4WgOsvWGj1HjKLIDktiD+OiIhIl7oKGJT02BfOWfJiWq5VgazeC7X7oaWpw3UhmHScjX1vWeCw6UHYcJ8dx1ywDHxyNJfB2z1q96tnhYiIDDpaYejCQVcYOtPSAnVlFjw0Vh/8upr9sOFvNurK2Bkey3Npi7no4k9CUmr7a5PTLXCI5CpBUkREAqctiT7qVcAQq6EGavZah0vf0vk1TfXw9lPsfPbPjGneYUHB3PNg7rmW0xArnGoJkukFEFJxThERCYa2JBItJR1SJlpRp5p9turQ3KGBVVIqzDqTX64bR2jnq/xk4hvwyi3w2m0w/RSYfxEUzrBrm+uhYhtU7bKgIWOUEiRFRCShtMLQhfxJc/xp374pLvdKa6kmu6WcdN9+u2LL3mqq65uYPy6HguYSjq99muK6VaTQwNvJM3g2soRNKfPwrv3KQqXLpjycR6NL6dVz3HHNsf3+s4iIyPCkFYZBoDaUQW0og7BvJLulnNqy3eyrrP3g++u2lwOpbMw6l8fzz+Gouuc5ru4ZPll5I3tDhTyftpjVkWNocJbnkOUryGqqoMZlUB7KpS6UfpCZRURE+k8rDF3ocw5DT3gPtaV85eYneHPbbh744gkHXtPSBFv+YVUkSzZASgbMPhvmXQBZh7S/VgmSIiLST1phGIycg/R8diRN4C2fZLkJHZMkQ0kwbYmNkg12LLN1TD4xWn56vt2rsQZKt0I4JVp6WgmSIiISPwoYBoGC3Bzrhpk11uowdJYkWTQXPnQdLPqslZ/euAK2rIRRs62ew9STLcBobrAEycqdVnZaCZIiIhIH+gg6CIzPi+YfhJNsW2H0XMifagWcOsosgkXXwOV3wfFftp4VT/0Y/nopvHYr1EUrT/pmqNoNu9fbykNDTcL+PCIiMvxohWGwiuTYaK0kWbPPgoBWyWkw7yNWt+H9VZbn8NLvYM0tMPN0WHAh5E7CKkiW2kjJgsxRdl8REZFeUMAw2CWlQs44yBpjP/Sr90BT2+kKXAgmHmtj/zuW37D5Ydh4P0xYZHkO44otz6GhEvZXtrXYTstTnoOIiPSITkl0IdBTEv3RXbvt2lLYcL/1rKgthbzJFjhMP619+elQEqS35jkodhQRGelUGrqPBm3A0KqpIdpue++Bja/AEiDffsq2K/a9bVsRc861rYzY8tMuZKsNGUWQHOl22ktueAFQESgRkeFGxyqHq6QUyB5jNRlqSw9sfBVOsXbaM06Hna/ZdsWrf4F//hWmnWKrDoUz7ShnzT4bqdm24hDpJOFSRERGLAUMw0G0pgPp+XZqonoP1JbxwXaFczD2cBvl22DdcnjjIXjzMRhzqPWtmHQchMJQX2EjKS3a8Cp/WBSC0qqIiEj/KGDohHNuKbB0+vTpA/0ovZeSYSO7se10RUtj2/dzxsPxX4LiT8EbD8O6e+Dx71kNiPkXwKyzrHlWUy2UvweVO6KBQ6HyHERERjClyHfCe7/Ce391Ts4QPn4YTrbtitHz7Hhlckb776dmwcKL4dJb4dQf2ErCC/8Dty6DF35thZ/AciMqd8LudVD2HjTWdpxJRERGAH1kHO66264IJcHUxTZKNlqew7rltvIw+YRo+ekF4Pggz+GQpu2Uh3IH6k8kIiIDQAHDSHLAdkWH0xVFc+BD37NKkhvui5affgZGzbI8h6mLIZxMmq8hrbnGAoyMUZCWr3oOIiLDnAKGkah1u+KD0xV7rHlVq8wiOPpqOPzjlhi59m74+/+BVTfAvI+Q3jKdrVVhaKqD8vehYof1rUgvtJMbIiIy7ChgGMlityvqqyxwiC0GlZwGc8+DOUvh/ZetnsPLv+daknmy+VAoLYK8SW19K6pKrNZDZpGtZAwy20rVT0NEpK8UMIwwrccLDybsG8luKSe7pYIQMb0rSAc+yejc01hY9gQfCr0Kd32SN5Ln8GzaEt5Mnt3u+GW9i1AeyqXaZXZ5LDORxxy3l9UlbC4RkeFGlR67MOgrPfZBdwFDK+dbyGypIKeljGTsWGZJRR0lldZ2O5tqzgy/xNnhVeS7KnaHD+G5tMW8mnoUja5tW6KZMOWhXCpDObS4cPz/QD20YWcFlXVNLJqSH/hcqvUgIkOVSkP30XAMGPqkrsK2K+qtdfa1y9eybns5D3zxhGj56aej5afftEqRc86FeedZQmSrD8pPj7KtjqieBjB9ta20ptOVhXG5kba24nGmgEFEhiqVhpb+iWTbaKyD6j14YrYYwikw88Mw4zTY9boFDq/dGi0/vcSOZY6a3b78dEqWJUmm5Sbsh+slN7zAqi372Xr92QmZT0RkuFHAID2XHIHcCbybNAWXvcuCheYG+55zVmZ6zKF2amLdPVZ++q0n4JAFsGAZTDreyk83VNqoSI2eriiw10VEZNBSwCC95l2Y9PxxUDS386ZX2WPhuC9a+elND8P6e+Dx6+wY57wLYfaZkJIJzfVQsd0qSablR7cruu+W2VfjcoO7t4jIcKeAQfquu2OZKZmwcJn1qHj3eduuePHXsOZmmHWmvZ49LrpdES0klZptqw6R+JflDipnQURkJFDAIPGRmmmjqcECh5p9Vp8BbLthyok29rxhgcP6+6wE9eTjrYrkmEMtAGntlhlObeuWqe0KEZEBp4BB4ispBXLG2fZDzT4LHlrzHMDKTJ/yXVj0WQsaNt4PW5+FghmW5zBtiVWibK6Him22XZGeb1UkA9yuEBGRrqkBgAQjFLaKj0VzIW/ygd0yMwrh6M/A5XfCiV+zoOLpn8Btl8Art0BdmV3nmy3o2LMR9r0d3fIQEZFE0wqDBMs5q7+QlmfdMqtK2uc5JEWs9PTsc2Dby7Dublh9E7z6FzuqOf8iyJ9i17bbrhjepyta61OopoOIDBYKGKTX+vxDLCXDfvg31cfkObTY95yDCUfbKN1qxzI3PwabHoRxxbZdMeEoKwB1wOmKwnbFoEREJP4UMEjiJaVCznjIGtN5nkPeZNumOOoz1mJ7/b3wyDchdyLMvxBmnm4rE7GnK1KyIKMAIrmd9q7QJ3URkf5RwCADpzXPIWNU5222Izlw+Mdg4SXwztN2uuLZX8DLv7dtjHnnt5Wfbi0GFU6xBMn0AgjrP28RkXjR36gy8NrVc6i0PIdo3wrATk3MOA2mnwq718Lau+Gft8M/74Cpi227omi2XdvcAJU7otsV0d4VKaq/ICLSXwoYZHBJzbLRWAfVJVCznw8SJJ2DQxbaqNgJ65dbjsPbT8Lo+da3YvIJEEqy99Tut5GcEe1dkddlq20RETk4BQwyOCVHLGcha4yVnq7ZCy1Nbd/PHgPHfh6OvALeeNiSJJ/4AWSOtgqSs86ywAOsbHVZtSVKphfYlkVSSmeziojIQShgkMEtnGzBQebomATJ+rbvp2TYysK88+G952274sX/hTV/hJnR8tM54+3aliao2m1bHpFs265oDSpERKRLChhkaAiFIHOUjdpSqNrTvuFVKAyTT7Sxd7MFDhvvtxMWk46zoGLMYdEtCW+1IOrK7bRFeqFKUIuIdEMBgww9rYWg6qtsxSA2QRKgcCYs+TYsuqat/PS7z0HBdAscpp1ipykAmuraSlC3JkmqBLWIyAFUGroTzrmlzrkby8tVhnhQS82Egmkwao7lJtAhoTG9AI66Ei67E076um1JPH29lZ9e8yeoLWu71jdbnsSejbD3LVvF8D6hf5yOtpXWdH+RiEiCKGDohPd+hff+6pyc+LdYlgC0JkiOnme5Dq7D1kJSqpWevuhmOOs/oHCWtdi+bRms/Dnsf6f99Q2VVm1y93qo3AXNjQn7o8TaXlY3IPOKiHTG+QH+FDWYFRcX+9WrVw/0Y0hvtbREEyRL2leQjFX2Lqy9BzY/akmU446Ilp9eZOWn23FWRCqjsEdJkq19IPpjw84KKuuaWDQlv9/3UpVLEekp59wa731xZ99TDoMMP60JkhmF0QTJEmiqbX9N7iQ48atWfnpTa/npb0HOhLby0x/0p/DWPbOuDJLS2mo6BJAkua20pt3Kwqot+wEYlxthfJ4KUInIwNEKQxe0wjCM1FXYkcyOCZKtWpqi5afvhj2bICWzrfx0ZtGB17twNEkymMZXl9zwAqu27Gfr9WfH/d4iIgejFQaRSLaNhhrbqqgt44MKkmDVIaefCtM+ZLkLa++C1++wMeVkWLgMiua2Xd+aJNmDxlciIsOBAgYZWVLSIWUyZHXSYhui5afn26jcZeWnNz4I7/zdAoYFy2DKidHy01Gtja9CydFKkgVxqSQ5LlfHO0Vk8FDAICNTa4vtzEOigUOH0tMAWYfAMZ+DI66AzY/AurvhyR9CRpFVkJx9dvskyJZGqNpltSEi2VYQKpLd50dUzoKIDCYKGGRkCyd1KD3dycmKlHQLEOaeB++9YHkOq34bLT99hhWDai0/DbSrJBlOtTyH9AJVkhSRIS3uAYNz7v4eXLbfe39FvOcW6bOenKwIha0b5uQTYO+btuKw6UHY8DeYeKwFDmMPb5/H0FxvTa8qd1qOg9pti8gQFcQKwxzgM1183wG/DmBekf5zzvpKpOfbyYqqEstP6KhwBiz+Fhx9NWy434KGB78K+dPayk8npbZd71ti2m2nW+AQybVARURkCAgiYPiO935lVxc4534YwLwi8fXByYpqy0uo66RUeHoBFH8KDrsM3nrSTles/Bm8dCPMOde2MdI7FF9qrLHCUW5bNDgpVP8KERn04h4weO/vjMc1IoNGSgbkT4XGOstxqNlPuyOZEC0/fRbMOhO2r7Htilf+BK/dZsc1F1xkfS9i+WZLuKzeo6OZIjLoBZL06JwbD3wUOAEYC9QC64AHgYe9jz3HJjJEtPasyBpjWxU1++yHfiznYHyxjbL3Yd09dsJi88OW37BgGUw85sDy0x2OZoZ9E81OOckiMngEkfR4MzAOeAD4GVACRICZwBnAd5xz13rvn4n33CIJEU6GnHF27LJ6r606dDySCZA7AU74snXM3PiA1XR49NuQPc7KT886w/IZYkWPZk5s2kKNy7A8in4czRQRiZcgPsL8p/d+XSevrwOWO+dSgIkBzCuSWKEwZI22BMba/bbq0Fx/4HWpWXDYR61a5JZnLM/h+V/B6j9YF835F9ixzg7SfTXsf7vtaGZavh0DFREZAEHkMBwQLDjn8oAJ3vvXvfcNwFvxnldkwIRCbbUWDnYkE6w65LRTbOxeb/Uc1t5lY8pJtl0xet6B7zvgaGah5VWIiCRQYM2nnHNPA+diQclrwB5gpff+q4FMGAA1n5I+qyuPHsms6vq6qt2w7l7Y9IBdWzQH5l8EU09uX366o+R0C1DS8nU0U0TipqvmU0EGDK967w93zn0GW134vnPude/9wkAmDIACBum3ro5kxmqsgc2PWpJk+Tbb5ph3vnXMjC0/3VHAXTNFZGQZqG6VSc65McDFwHcCnEdk8Io9klm127YsOh7JBFsxmHd+tPz0i7Zd8dKN8MotVn56/gV2QqOjdl0zM6OrDnk6mikicRdkwPAj4FHgWe/9y865qcCbAc4nMnglRyBvkh3JrC45sEtmKxeCScfZ2Pe2BQ6bHoQN99lxzPkXwbgjOw8IGqpsVGy3rYqMwvbVJkVE+iGwLYnhQFsSEpjmpraiTR1rOXRUsx82RstP15ZC3hQrBDX91O4DgtRsW3WI5GjVQUS6NSA5DMOBAgYJXEuzrTZUlVgNhq401cPbT9mqw/637cTE3Nby0wVdvzdaEIqMQqsjISLSCQUMfaSAQRLG+7bAobNaDh2v3fkavH6XtdsOJcH0U2y7onBGNxM5KwSVXqiCUCJygIFKehSRnnKuZ7UcWq8de7iN8m12suKNh+2UxZjDbLti4rFWWOoA3k5s1JVbQaj0AhsqCCUi3QhkhcE5NxsrD73Ke18V8/oZ3vtH4j5hQLTCIAOqp7UcAOorLTly/b12GiN7rK04zDwDUtK7ebODtFxbdUjNjMuji8jQlNAtCefcl4DPAxuBw4B/9d7/Lfq9V7z3R8R1wgApYJBBob7KgoD6iu6vbWmCrc9a9cjd6yE5A2afDfPPtxMa3UlKi6465B9khUJEhrNEb0lcBRzpva9yzk0G7nbOTfbe/xJQmrZIb6Vm2misjdZyKKPTWg5g+QxTF9so2WB5DuvutjH5xLby0wc7MdFUCxXboHKH1XNIL+zBCoWIjARBBAzh1m0I7/1W59xiLGiYhAIGkb5LToO8yZBVb4FDzX4OGjgAFM2FU79v2xrro+Wnt6yEUbMtcOiq/LRvsSTMmn3RMtSFFkAMszLUl9zwAgB3XHPsAD+JyOAXxP/9u5xzh7V+EQ0ezgEKgQUBzCcysiSlWtXH0fMgo8iKPXUlswgWXQOX3QknfMXKVT/1Y/jrpfDardZCuyuNNVD+HuxeB2Xv20qHiIw4QawwfAJoin3Be98EfMI5d0MA84mMTOFkyBlnrbF7UgQqOc1qNsxZCu+/ZHkOL/0O1twCMz9spytyJx38/SpDLTKiBdHeelvr71vbWsfMo48mIvEWToLsMRY41OztvgiUC1mZ6YnHwP53rBDU5kdg4wqYsMgCh3HFXQcCKkMtMuIEdvjaOfdj4Argbdo2Wj1wSlBzioxooZBtP2SM6nkRqPypcPI34OirLGBYfy889G+WKzH/IphxWteBQEuT9caoLoGULMgosAqUWnUQGXaCrNZyMTDNe98Q4Bwi0lFvikC1SsuDIz4Bh14Kb//dtiv+8R/WMXPuudZJs7vy0w2VNkLJdiwzvRCSUuL35xKRARVkwLAOyAVKApxDRA7GuegP7nw7illVAo3VXb8nnAIzT4cZH4ad/7TA4dVb4Z+3w9QlsHAZFM7s+h4tjXaKo2q3ml+JDCNBBgw/BV51zq0DPlgX9d6fG+CcItKZtFwb9ZVQudtWArriHIw9zEbFdli3HN54CN56HMYcatsVk47rvrhTfYWNcEpMGWo1vxIZioIMGP4E/AxYC7QEOE/cOeeWAkunT58+0I8iEl+pWTYaqm0FoK68+/dkj4PjvgjFn4JND1nvise/Z5Uj518As86ClAyuXb4WgOsv6OT0dHMDVO6Eyl1qfiUyRAUZMOz13v8qwPsHxnu/AlhRXFx81UA/i0ggUjIs4bEn1SM/eE8mLLzYgoStz9l2xQu/htU3w6yzyGueR2m4sJuJ1fxKZKgK8v/SNc65nwL3035L4pUA5xSR3mhXPbLETld0FziEkqxK5NSToWSTlZ1efy//5u9hQ8pC2HklHLKg+5yF5norQV25M9r8qsBWP0RkUAoyYDg8+usxMa/pWKXIYJSUCrkTIOuQaOCw18pDd6doNpzyXVh0DSuX38jRdc/Bii9ZYuSCiyxRstucBW+nOWpLISnStuqg5lcig0pgAYP3fklQ9xaRgPS2emSrjFE8mnEuT6WfwY/nbLdVh7//BFbdAHM/AnOXWn2G7jTVWZJl5U67Pr1ALbdFBonAOsk4537inMuN+TrPOffvQc0nInHUWj1y9DzIGmu1FXqg0aVY3YZlf4Qzf2Z5Eqv/ALdeDM/8J5Ru7dn8vgVq98O+N23bo3ovtPQgcBGRwAS5JXGm9/7brV9470udc2cB3w1wThGJp1AYskZb9cja/ZYg2dyDWmwuZGWmJyyC/VvsZMWbj8KmFTD+KNuuGH90z2ozNNVC+fvRMtRquS0yUIIMGMLOuVTvfT2Acy4NULF5kaEoFOpQPXK3bR/0RP4UOOnrcNRnYOP9sOFv8PA3rdHV/Aut8VVSpPv7jJCW2yKDVZABw1+AJ51zN2PJjp/GajOIyFDVrnpktOx0Y03P3puWGy0//VF45+/W9OrZ/4KXf28dNOd9xFYyeqK15XbrqkNGoZ34EJHABJn0+HPn3OvAqYADfuy9fzSo+UQkwdLybNSVW+DQU+FkKz09/TTY9boFDq/dZuWnpy2x7YpRs3t2r9iW28kZFjhEcrXqIBKAuAcMzjnnvfcA3vtHgEe6ukZEhrhIDkRy2Bl+j9yW0p6/zzkrMz3mUKjYEVN++gkYPR8WLIPJJ/T8eGVjNZRVg9vW1vwqufutjm2lPVwhERnhggjD/+6c+6JzbmLsi865FOfcKc65PwGfDGBeERlAdaF0diWNsxoMkZzevTl7LBz3Bbj8Ljj285an8MT34Y7L4fU7oaGq5/fyzXYcdM9G2Psm1OyHLj6fbC/rYS7GIHDJDS9wyQ0vDPRjyAgVxJbEGVi+wl+dc1OAMiANC04eA37hvX8tgHlFZDDoS9np2PcuWAbzLoB3n7fy0y/+BtZY+WnmX2C9LeCD3hU90UKYylAWFaEcmlxby+0NOysA+v1D+I5rju3X+0WGgrgHDN77OuA3wG+cc8lAIVDrvS+L91wiMri0W95vV3Z6t33S72ngEArDlBNt7NkEa++B9ffZtsXk461bpg/1uGV2iGZyWsrIaSmjzqWxqSKFTeVhLL0KVm3ZD8C43Ajj83RkU6QzgXZ88d43AjuDnENEBo9Ol/eTUiF3ImQeAtXRfhU9KTvdatRsOOU7sOga2HAfbLgftj7L9QUzLHCYdkrfWmaHkvns8i08+V4Lb15/Xu/fLzLCqEWcyDCXqD3vni7vh3wzOS2lZLeUE6IXgQMAi0jOPJzDkl/mwpYX4emfWvnpeR+BOefa0c2eamkkt6WUWa4c9r1tNSYiOT1etRAZaRQwiEi/bCutabey0N3yfosLUxoupCyUF90mKCdEz8s+N7oUXo4cz4XnXwPbV1uew+qb4NW/wIzTrBhU/tQe368oKxXqK2yEktuaXyWldP9mkREkiGOVj2JHKR/23m+K9/1FpHcSlZB3yQ0vsGrLfrZef3bv3tjSYnUUqkqgpbF37x1/lI3SrVZ+evNjsOlBGFdsyZMTjrIy1V0oyo45etnSCFW7LOciNSta16GXJz5EhqkgVhg+iZ2U+IFzbiawCgsgnvTe9+JslIiMCKEQZBZZ3YTe9KuIlTcZTvxatPz0A7D+Xnjkm5Y7Mf9CKxTVq0qQvm3VIZzSturQl1wJkWEiiFMSu4A/An90zoWARcCZwDecc7XAY977n8d7XhEZeONye9AT4mD606+iVSQHDr8cDr0E3nnatiue/YWVn559Dsw734KT3mhusHbblbsgkm2BTSS7d/cQGQaCPiXRArwQHdc55wqB04OcU0QGTlyOJHbsV1G52zpW9kYoCaafCtM+BLvXWeDw+h1WBGrqYis/XTSnlw/mrQx2XTmEU2NWHZQKJiNDQv9L997vBW5N5JwiMoTF9quo3G3ln3vDOThkgY3KnbDuXstxePtJGD2f+fVHsyFlYe+fq7keKnfYPSM5tjKSmtX7+4gMIQqNRWTwi/aroL7SAoeGyt7fI2sMHPs5OPIK2PwwrL2Hj1XeRGkoH16/xCpJ9vqHvoe6MhtJEVtxSMvXqoMMS2rpJiJDR2oWFE6HghmQ2sc8gpR0S4S85M/cknWVBQwv/i/cugye+yWUb+vbfZvqrN327nV2aqNeOd4yvAQWBjvnvtrJy+XAGvWSEJF+Sc200VBjyZF1fag8HwqzIXUhG1IXcv1JEWuzvXGFlaCedJzlOYw5rA+FnLzlXtSWRlcdCi0fo6ddN0UGqSDXzYqjY0X067OBl4HPOufu0kkJEem3lHTInwKNdVY/oTeNrmIVzoAl34JFV8OGv9l49zkomG6Bw7RT7HhlbzXVQcU2y3eI5FquQ0pG7+8jMggEuSVRABzhvf+a9/5rWPAwCjgJuCLAecKQEFYAACAASURBVEVkpEmOWC2GojmWR0AfyzunF0Dxp+GyO+Gkf4OWJnj6erjtEljzJ1s16AvfYjUm9m6Gkk1QvRdael7dUmQwCHKFYSIQW32lEZjkva91ztUHOK+IjFT9bXQVe5/ZZ1si5PY1tl2x5mZ47S8w/TRbdehF+el2mmqh/H3Ld0jLsy2LFHXIlMEvyIDhNuBF59zfol8vBf7qnMsANgQ4r4iMdEkpkDO+LXCo3gu+D5/onYPxxTbK3rX22m88Am88BOOOsG6ZE4/ptvx0p3yLBTQ1+yA53QKHtDwrYCUyCAUSMDjnHFbt8SHgBGx98LPe+9XRSy4PYl4RkXbCSZA9FjJHQ/Ue61fRl8ABIHcSnPAVKL4SNj0A65fDo9+GnAl26mLm6b0sPx2jsQbK34tZdSjQqoMMOoEEDN5775y7z3t/JLAmiDlERHosFIasQyBjlK02VJdYfkJfRLLhsMtg4cXwzkrbrnjuv6389JylfSs/3co3WyOumr2QnBGt69B+1WFbaU3f7i3ST0FuSbzonDvKe/9ygHOIiPRcKAxZoy1wqNlHExtJoo+BQygJpn/ITlCUbIgpP30HTDnZ8hxGz+v7szZWQ3l126pDRiFAu1biIokUZMCwBDtCuRWoxrYlvPe+D3VYRUTiKBSCzFG8nzSZLF9hvSGa+5iL7RzXPtcCXEhu7hKOq1vJ0VueJ/LO33k3aTLPpS1mXcphtLj+12Go2dtMLmlc+tvn8H3Jm+ilRLVGl6EhyIDhzADvLSKDzJD84eIclS7HjmP2tUNmjLJwPg9lnM8TaWdSXL+K42qf5rLKP1IWyuP5yEm8HDmO2lDvcxNKKuooqbSAZryrouLdfZT6LCI5oyjKz+3z84r0RmABg/f+XefcCcAM7/3NzrlRQGZQ84mI9Fk/O2Ref8GCTl49Clr+Bd5fRe7auzhrx984q+FRmHmGJUnmTuj1Y167fC3rtpfz0BePa3sxJbMt16HXVSlFei7I0tDfx4o1zQJuBpKBvwDHBzWniEi/9bdDZqxQ2MpMTzoO9r1lCZKbHrRKkhOPsTyHsUf07wd9Q5WNiu3W+Cqj0OpIiMRZkFsS5wOHA68AeO93OOfU/1VEBo0ut1FaO2TWVdhxzL50yIxVMB0WXwtHXwUb7reg4cGvWQGo+RdZAmUPftAXZR3kmpamaM2JEkjJgowCK0etVQeJkyADhobo8UoPEC3YJCIytESybdRXWY5DfUX/7pdeAMWfsqOZbz9lpyue+Tm8dCPMPc9Gev5B316UHel+joZKG6Hk6FZLgVYdpN+CDBjudM7dAOQ6564CPg38LsD5RESC065D5i7bsuiPpFSYdablNOx41QKHV/4Er90G00+BBctsVaI/WhotyKnabe3A0wts1USrDtIHQSY9/odz7jSgAstjuM57/3hQ84mIJERKum0jNNbaD+K+NqRq5ZyVmR53BJS9D+vugc2PwOZHYezhtl0x6di+lZ+OVV9hI5RsgUN6gZXQFumhIFcYiAYIChJEZPhJTrMOmZmHxAQOfWitHSt3ApzwZTjqStgYLT/92HcgexzMv5AUP5EG18+thZZGWyGp2qVVB+mVuAcMzrkHvPfn9PcaEZEhITkCeZMga4wFDjX76HfgkJoFh30UFi6DLc/YdsXzv+JbLo2XIsdB5Sgrdd1frasO4ZS2VYdwcv/vK8OS876f/2F3vKFzZcAzXV0CzPPe97E3bOIUFxf71atXd3+hiEir5sa2wKEvrbUPZvd6/vnw75nf8E/CDphykm1XjJ4Xx9UBZwmerasOcXDJDS8AQ7Sw1wjknFvjvS/u7HtBbEmc14NrGgKYV0Rk4IWT49Nau6PR8/hr9qfJad7PtyZugo0r4J2nYdQcq+cw9WTrb9Ev3pI568q16iAHiHvA4L1fGe97iogMObGttatKrL12HAKH8nA+LLoGjvi4JUauuwee+jGs+q11ypx9jq0S9FdzA1TuhMpd0VWHwvjcV4asQJMeRURGvFAYssdYy+v+ttaOlZxuAcLc8+C9F2Hd3VbL4ZVbYObp0fLTE/s/T7tVh9SYVQf9+Bhp9G9cRCQR2rXW3murDi2N/b+vC8WUn37byk+/8ZBVkpwQLT897sj45Dk010PlDlt5iORYGepUFfAdKYLsJVHkvS/p8Nos7/0bQc0pIjLohUK22pAxyhIjq3bb8n88FEyDxd+ERVdHy0/fBw99HfKmWOAw/dQ4VXz0UFdmIykSbX6Vr1WHYS7Ihur/cM5d3PqFc+5rwL0BziciMnQ4Z5/Qi+ZC7iRb7o+XtDw48pNw2R1w8jdtFeKZ/wu3XQyrb4oe/YyTpjprfLV7HZRutRLaMiwFGQ4uBm50zi0DRgMbgaMDnE9EZOhpba2dlmfFn6pKetVau0vhlLby0ztfg9fvglf+bOWnp0XLTxfOiM9ceHv+2tLoqkNhlz0xZOgJsjT0TufcI8C3gBbgW957hZ4iIp1pDRzS86G2zLYqGmvid++xh9so32YnK954GN58DMYcatsVE4+zPIt4aKqDim1QuYNRTbuoCMWnpoMMrCBzGB4HdgLzgfHATc65Z7z3Xw9qThGRYSEt10ZdOVTuhsbq+N07Zzwc/69Q/GlLjly3HB77HmSNhfkXwKyzrF9GPPgWMn0lmc2VULLJtmDS8uIXmEhCBbkl8Wvv/X3R35c5547DVhtERKQnIjk26ipsqyKeUrNg4SV2/HLrs1Z++oX/gdU3w+yzYN4Fdhw0Xppqofx9y3dIy7NEyZSM+N1fAhfklsR9Hb5uAn4c1HwiIsNWJBsi2ewIv0dey/743juUBFMX2yjZaIHDuntsTDoBFl4EoxfEr/y0b7Gky5p9VksivWDQrTr84vHNfOW0mQP9GINOYKcknHOVzrmK6KhzzjU75/rZQF5EZOSqD6WxK2kcFM60TpPxVjQHPnQdfPR2W33Y+Src/yW477Pw1hPxKTgVq7HGVh12r4Oy96AhTjkb/fTLJ98c6EcYlIJcYWhXzcM59xF0SkJEpP9SMqzmQmOtlW6uK4vv/TOLouWnP2GJkWvvhqf+HV6Mlp+ec07cmlMBnaw6tOY6tP9M29rIKhESNddQasqVsCob3vv7nHPXJmo+EZFhLzkN8qdAYx1U7bLTFf1trd3x/nPPgzlL4f2XbLvi5d9Fy09/2Lpl5k2K33wQXXV4ry3XIaPQniNg20pr2F5W98HXq7bY1s+43Ajj8+KUBDrEBXlK4oKYL0NAMXH9L1lERABIjkDeZMiqj7bW3k9c/7p1IZh4jI3979iKw+ZHrGPmhKMtcBh/VBzbbGONumr22kjOgIxC7rhq0QGrDvF2yQ0vsGrLfrZef3ag8wxFQa4wLI35fROwlZ61vo4L59xU4DtAjvf+ooO9JiIyVHS7fJ2Uag2nMg+JBg77iPvntPypcPI34OirYeP9sP4+ePgbFrDMvxBmfDhO5adjNFZDWTW4bdFaFYUWJElCBZnD8Km+vtc5dxNwDlDivZ8f8/oZwC+BMPB77/31Xcz/DnClc+7url4TERl2klIgdwJkHWLHMWv2Wp5APKXlWo7DoZfC23+37Yp//Ce89DuYey7M/YhtJ8STb7Y24dV7ICWz7YRFPFc2sG0IOVDcAwbn3P+ji5DWe/+lHtzmj8D/ALfE3DcM/Bo4DdgGvOycux8LHn7a4f2f7tj4SkRkxAknQ844yBxtbbWr99oP3bjOkWLttGd8GHb+0wKHV2+Ff94OU5cwrukwtifFo812Bw1VNiq2W+Or9IK4rTooZ6FzQawwrO7vDbz3zzjnJnd4+WjgregqAc6524HzvPc/xVYjRESkM+EkyB4LGUVtn9DjHTg4B2MPs1Gx3SpIvvEQX2x8nC1JU2HLFTDp+PjXW2hpigZDJZCSBRkFEMmN+6qDBBMw3Bot0hRv44D3Y77eBiw62MXOuQLg/wCHO+e+5b3/aWevdfK+q4GrASZODCAqFhEZKOEkq96YWWSrDdUl8a+tAJA9Do77IhR/igfu+gPH166Ex6+DrDEx5acDqPLYUGkjlGQrDukF8c+nGMGCCBheAo4A257w3n8xTvftLFzsautjH/DZ7l7r5H03AjcCFBcX61SHiAw/oTBkjYaMUZbfUFUCLY3xnyclk2fTTuH5yMn85PByWHc3vPBrKz896ywLHrLHxn/eliZL+qzabQWu0vO16hAHQQQMsf9Gjo/jfbcBE2K+Hg/siOP9RURGllDIVhvSC+1ERXUJNDfEfZoWF4apJ9vYs8mOZa6/F9Yvt22KBRfBIQuD+YFeX2EjlNx2wiIpJf7zjABBBAxBfSp/GZjhnJsCbAcuBS4LaC4RkZEjFILMUXaqoWa/FYEKIHAAYNRsOOW7sOizdiRz4/2w9R9W7nrBRTB1iSVrxltLY4dVhwKrVqlVhx4LImCY7Zx7HVtpmBb9PdGvvfd+YXc3cM79FVgMFDrntgHf997/wTn3BeBR7GTETd779QE8v4jIyOScJQ2m50NtqZWdbq4PZq6MQjj6M3DEx2DzY9bs6u8/gVU32JHMuUttGyEI7VYdWnMdtOrQnSAChjn9vYH3/qMHef0h4KH+3l9ERLrgnAUNaXkWOFSVWHvqICRFrG7DnHNg28u2XbH6D/Dqn2HG6bDgQisKFYSWRltNqdpt7b7jXTdimIl7wOC9fzfe9xQRkQHQGji0rjhUlVivh0DmCsGERTZKt8Lae+DNR2HTCis7veAiGH90QFsI/oNVhwmNW6gMZUNzYzBbI0NYwppPiYjIEJaWZ6OuHCp3W7nmoORNhpO+BkdfCRtWwIb74OFvQu4kKz8988O2MhGAJJrIa9kPu9dDJNuSJCMBtBIfghQwiIhIz0VybNRV2FJ+Q1WAc+XCER+38tPv/N22K579L3j599ZBc95H7GhoILwFR3XlEE5ty3UIj9wfm4H8yZ1zC733rzvnFnjv1wYxR5Ccc0uBpdOnTx/oRxERGZwi2TbqK23FoaEyuLnCyVZ6evppsGutlZ9+7TYrPz1tiXXLLJod3PzN9VC5Ayp3WrCUUWg5DyNMUKHSp51zvwauBL4c0ByB8d6vAFYUFxdfNdDPIiIyqKVm2WiotlMV9RXBzeUcjFloo2LHB+WneesJGD3f8hwmn2CVHgPhoa7MRjjVAoe0/BGz6hD3xuLOue9H7/siEHLOXRfvOUREZJBJyYCCaVA4yz6FBy17LBz3Bbj8Ljj281Z46okfwO2Xw+t3BrtVArbqULEddq+zJM36AFdYBokgTkn80Dl3bvTeT3jv74/3HCIiMkilpEP+VGistRWHwOfLgAXLYN4F8O7zVn76xd/Amtjy0+MCfABvJ0hqSy0RM73QTpXEu8nWIBDUOsoi7/3nnHM/BhQwiIiMNMlpkD+FbUnvkdu8n2jtvuDmC4Vhyok29rxhhaA2/M22LSYdZ0HFmEODrezYVAcV2yzfIZJrSZKpmcHNl2CBBAze++9Ef/1eEPcXEZGhodGlsCfpECiaYysOtaUEGjgAjJoFS74Ni65pKz/97nNQMMPyHKYtgXCAlR19C9Tut5GUFs11yBvyqw5xz2EQERE5QFIq5E2Corm2bN9pA+I4Sy+Ao66Ey+6Ek75ulR2f/incdim8cgvUlgX/DE21UP6+5TqUvWfJoUPUyEjtFBGRwSEpBXInQOZo645Zs88+kQc6ZyrMPgdmnQ3bV0fLT99k5aenn2arDvlTg30G32J/1pp9kJxuwUxavjX+GiIUMIiISOIlpUDOeAscqkqgZm/wgYNzVmZ6/FFQ+q7lOWx+1I5mjjsSFizD+TS8C/iHeGMNlNfYKYu0fAseUtKDnTMO4h4wOOfCwGeA8cAj3vvnYr73Xe/9v8d7ThERGaLCyZAzrm3FoXov+Obg582bBCd+1bYsNj4A6++FR67lK+EinosshsbplrgZJN9igVLNXkjOiK465A3aVYcgnuoG4GRgH/Ar59x/xXzvggDmizvn3FLn3I3l5eUD/SgiIiNDOMlqK4yeB5mHgEtQgmAkBw6/HC67HU75LvUuwvnVd8JtF1ur7aqSxDxHYzWUvxfNdXjfjqUOMkEEDEd77y/z3v83sAjIdM4td86lkpAsl/7z3q/w3l+dk5OA4iMiIsPcttJedLgMhSF7jAUOWWMDrNrYcd4kmH4qv875Ov+b8xUYewS8fgf89VJ44odQsiExz+GbbcVhzybYsxlq9kNLwFs1PRTEv4kPzqp475uAq6PVHp8Chs+BVBER6ZHtZXW9f1MoDFmjrblUzV77pN/SGP+H68g53k2eCqedZ70j1t8LGx+05lej51nfiiknJiaQaayGsmpw26JtxgshOZgunT0RxJ94tXPuDO/9I60veO9/5JzbAfxvAPOJiEgALrnhhX7fY8POirjcy/kWbv/YTOuQmYjAASBrDBzzOTjiCtj8MKy9B578IWQUWQXJ2WcnpgmVb4bqPTZSMttyHYIsQtWJIEpDf+wgr/8e+H285xMRkcFnW2lNu5WFVVv2AzAuN8L4vN6fCPAuBJmjrAhSzT4LHJob4va8XUpJh/kXwtyPwHsvWLfMVb+FNX+EWWfa93LGJ+ZZGqpstJ6wyCi0Y6MJEMQpiW94738e/f0y7/1dMd/7iff+2/GeU0RE4u+Oa47t9z0uueEFVm3Zz9brz47DE2GfqjMK7VN2zf5o4FAfn3t3JxS2bpiTT4C9b1rfio0rrJrkxGOtnsPYwxPzyb+lKXqqpARSsiCjwMpRBzh3EEmPl8b8/lsdvndGAPOJiMhI45z9kCyaA7mTrPFTIhXOgMXfgsvugCM+YUmRD34Vln8G3ngYmhIUxAA0VFrHzN3rre13QHMHETC4g/y+s69FRGSYG5cb4A9z5ywh8IPAIeDaCR2lF0DxpyxwOOnf7ETDyp/Z6YrVN9sqSALcuuo9y+2o2m3By763rfS1j1/fjiCSHv1Bft/Z1yIiMsz1JWehT9LzbdSW2qmKxl4c5+yvpFRLgpx1FmxfY+WnX/kTvHYbTD/VtisKpgU2/V9feo/LF01se6G+wkYo2YKa9AKrrtkPQQQMhzrnKrDVhLTo74l+PXDnQUREZGRIy7NRVw6Vu+14YqI4B+OLbZS911Z+evPDVtthwUUw8RiIlp++dvnauE3d3b1qXAYVoRxqQxl9un8QpySGdv9OEREZHiI5NuoqbKm+oSqx8+dOhBO+AsVXwqYHYf1yePTbdqJi/oUw8/R+3b6koo6SyrZ8hXXbrTpxUVYqRdkHfj5P99WkN1fT1JxEZSibylAOza7nYYCaT4mIyPAWybZRX2krDg2ViZ//sI/CwmXwzkrbrnjul/DyH7h+9jkw/3zrpdFH1y5fy7rt5TzwxRN6+U5nz5ZeaL8Cd3724FcrYOiEc24psHT69OkD/SgiIhIvqVk26qtsxaG+ovv3xFMoCaZ/yMbu9VbPYe2dNqacbNsVo+cl8IG8bdvUlUM4FdILSOriKIQChk5471cAK4qLi68a6GcREZE4S8200VANlbsSHziABQaj59n86++DTSus/HTRHFiwDKac1Kvy00VZ/Sze1FwPlTs4JOPgexQKGEREZGRKybCTCw01ULXLPmknWtYhcMxn4chPWHLk2rvhyR9Z+el558Occ3pUfrqznIV4U8AgIiIjW0o65E+Fxlqq3RYyfIKTIwGS0y1AmHsevPeiBQ4v3WBHM2eeYUmSuRMS/1wxFDCIiIgAJKdRkjSGZF9vxzJrSxP/DC4Ek46zse9tCxw2PQgb7rPjmAuW2fHMBDeegmAqPYqIiAxZjS4V8ibDqDnW4GmgihQXTIPF37QqkkdeAXvegAe/BvdcCZseSmz5aRQwiIiIdC45AnmTLBExvYABCxzS8y1g+OjtcPI37Tme+TncdgmsvpnMlsQkbWpLQkREpCtJqVaEKfMQO45Zs48B6XSQlGrttGeeATtejZafvoVrCfPP1CNg71XWFCuo6QO7s4iIyHCSlGKJh1kxgYNvSfxzOAfjjrBRvo2X/nYDR9atguVXwZjDLM9h4jHWjjuOFDCIiIj0RjjZyjtnjrYmVzV7ByZwAMgZz/2ZF/N4+jl8f8YWWH8vPPYdyB4bLT99pp0CiQPlMIiIiPRFOBlyxkHRPAse3MC1UqoNpcOhl8Klt8KpP7BTHs//P7htGbz4G6jc2e85tMIgIiLSH+Ek+0SfUQTVe2z45oF5llASTF1so2QDvH6X5TqsvRsmnxgtPz2/T8cyFTCIiIjEQzgJssdAZjRwqCoZuMABoGgunPp9e47198LGFbBlJYyabYHD1MW9Kj+tLQkREZF4CoUtMXL0PMga26sfyoHILIJF18Dld8LxX7YeGk/9O/z1UnjtVmv/3QNaYeiEulWKiEi/hcKQNRoyRlliZNVuaGkauOdJTod5H4G558L7q6xb5ku/gzW3wMzTbdWhCwoYOqFulSIiEjehkH3KTy+0o5hVu6GlceCex4Vg4rE29r9j+Q2bH4aN93f5Nm1JiIiIJEIoBJmjLLcgZwKEUwb6iazp1snfgMvuhCM/1eWlChhEREQSKRSCjMJo4DBxcAQOaXlw5Ce7vERbEiIiIgPBOcgosF4RNfttq6I5sQ2lekMBg4iIyECKDRxqSy1waKob6Kc6gAIGERGRwcA5CxpiVxwGUeCggEFERGSwaQ0cakuhcjc01Q70EylgEBERGbTS8mzUltmKQ2PNgD2KAgYREZHBLi3XRl25rTg0Vif8ERQwiIiIDBWRHBt1Fbbi0FCVsKkVMIiIyIj3oxXr8R6OmVoAwC8e30xFbSPOOa5bOneAn64TkWwb9ZXUuTeJ+OCTIxUwiIjIiJeZmsTv/vEOtY0tH7yWlhzm6pOmDuBT9UBqFjuTJhBpqYGULGioDGwqVXoUEZER73NLppMZSW73WlYkiX9ZPG2Anqh36kLpUDgdCmZAanYgcyhgEBGRES+SHObnFy0kLTkM2OrCzy5aSCT69ZCRmgkF06BwZtwDBwUMIiIiwJJZRRw5KY+Qg+LJeSyZVTTQj9R3KRnRwGGWJUnGgQKGTjjnljrnbiwvLx/oRxERkQT66QULWDAuh5+cv2CgHyU+UtKtI+Wo2RDJ7detFDB0wnu/wnt/dU5OfKIyEREZGibkp/O3L5zAhPz0gX6U+EpOg/wp/QocdEpCRERkpGgNHBrroGqXVZDE9+itChhERERGmuQI5E2GrHqo3GU9K7oJHBQwiIiIjFRJqZA3CbIOgardtHRxqXIYRERERrqkVMidyI5K33SwSxQwiIiISLcUMIiIiEi3FDCIiIhItxQwiIiISLcUMIiIiEi3FDCIiIhItxQwiIiISLcUMIiIiEi3FDCIiIhItxQwiIiISLcUMIiIiEi31HxKREQCc8c1xw70I0icaIVBREREuqWAQURERLqlgKETzrmlzrkby8vLB/pRREREBgUFDJ3w3q/w3l+dk5Mz0I8iIiIyKChgEBERkW4pYBAREZFuKWAQERGRbilgEBERkW4pYBAREZFuKWAQERGRbilgEBERkW6pl4SIiMgQ9KMV6/EejplaAMAvHt9MRW0jzjmuWzo37vMpYBARERmCMlOT+N0/3qG2seWD19KSw1x90tRA5tOWhIiIyBD0uSXTyYwkt3stK5LEvyyeFsh8ChhERESGoEhymJ9ftJC05DBgqws/u2ghkejX8aaAQUREZIhaMquIIyflEXJQPDmPJbOKAptLAYOIiMgQ9tMLFrBgXA4/OX9BoPMo6VFERGQIm5Cfzt++cELg82iFQURERLqlgEFERES6pYBBREREuqWAQURERLqlgEFERES6pYBBREREuuW89wP9DIOWc24P8G4f354DlCfofYXA3j7MJQfq67+3wWKwPH8inyOoueJ13/7eJ1F/l+jvkfgZLP8f9sUk7/2ozr6hgCEgzrkbvfdXJ+J9zrnV3vvi3s4lB+rrv7fBYrA8fyKfI6i54nXf/t4nUX+X6O+R+Bks/x/Gm7YkgrMiwe+T+Bjq//wHy/Mn8jmCmite9+3vffR3ydAzLP/Za4VhGNAnAxHpL/09It3RCsPwcONAP4CIDHn6e0S6pBUGERER6ZZWGERERKRbcetW6ZxbC3S2XOEA771fGK+5REREJLHitiXhnJvU1fe9932tZyAiIiIDLJAchmjwMMN7/4RzLg1I8t5Xxn2igBUWFvrJkycP9GOIiIgkxJo1a/YerHBT3LYkWjnnrgKuBvKBacB44LfAh+I9V9AmT57M6tWrB/oxREREEsI5d9DdgCCSHj8PHA9UAHjv3wSKAphHRGREWbx4MYsXLx7ox5ARKoiAod5739D6hXMuic6TIUVERGSICCJgWOmc+zaQ5pw7DbiLYVomU0SGNn1iF+m5IAKGa4E9wFrgGuAh4LsBzCMiIiIJEvekR+99C/C76BCRPmj91Pv0008P6HOIiLRKROEmAFS4SUREZOiK5wrDOdFfPx/99c/RXy8HauI4j4iIiCRY3AKG1kqOzrnjvffHx3zrWufcc8CP4jWXiIiIJFYQSY8ZzrkTWr9wzh0HZAQwj0iPKRteRKR/4p70CFwJ3OScy4l+XQZ8OoB5REREJEGCOCWxBjjUOZeN9aooj/ccIiIiklhx35JwzuU45/4LeAp40jn3nzGrDSIiIjIEBZHDcBNQCVwcHRXAzQHMIyIiIgkSRA7DNO/9hTFf/9A591oA84iIiEiCBLHCUNvhlMTxQG08buycO8M594Zz7i3n3LWdfP9y59zr0fG8c+7QmO9tdc6tdc695pxTz2oREZFeCGKF4V+AP0XzFhywH7iivzd1zoWBXwOnAduAl51z93vvN8RctgU42Xtf6pw7E7gRWBTz/SXe+739fRYREZGRJohTEq/RdkoC731FnG59NPCW9/4dAOfc7cB5wAcBg/f++ZjrXwTGx2luERGRES3uAYNzLhf4/AWyPwAAFClJREFUBDAZSHLOAeC9/1I/bz0OeD/m6220Xz3o6Erg4ZivPfCYc84DN3jvb+zsTc65q4GrASZOnNivB5bBZevWrQP9CCIjjhqpDR9BbEk8hH26Xwu0xPG+rpPXOm125ZxbggUMJ8S8fLz3fodzrgh43Dm3yXv/zAE3tEDiRoDi4uKDNtOSoefdd98d6EcQERmygggYIt77rwZw323AhJivxwM7Ol7knFsI/B4403u/r/V17/2O6K8lzrl7sS2OAwIGSaxElWt+7bXXEjZfvD5JaUVERAaTIAKGPzvnrgIeAOpbX/Te7+/nfV8GZjjnpgDbgUuBy2IvcM5NBJYDH/feb455PQMIee8ro7//MGqGNSJs3bq13crCypUrAZg0aRKTJ08eoKfqGa2IiMhgEkTA0AD8X+A7tG0ZeGBqf27qvW9yzn0BeBQIAzd579c75z4b/f5vgeuAAuA30dyJJu99MTAauDf6WhJwm/f+kf48j8RHovY1Fy9ezMqVK/E++F2meKxixHNFRHvHIhIPQQQMXwWmB3F80Xv/EJYjEfvab2N+/xngM5287x3g0I6viww2Q3lFRESGtyAChvVATQD3FemXSZMmJWSeeHyiT+SKiIhITwQRMDQDrznn/k77HIb+HqsU6Rd9QhcR6bsgAob7okNE+iFRKyIiIj0RRKXHP8X7niIjkVZERGQw6XXzKedcqLXss4jIUKd6FyI906OAwTl3m3MuO1rDYAPwhnPu34J9NBGR4A21ehcKcGSg9HRLYq73vsI5dzl2rPGbwBqs3gIAzrlvAY9471+N/2OKiLQ3mOpdJLLWxVALcGT46GnAkOycSwY+AvyP974x2sQp1hbgX51zhwL/xBo/Pea9L43f48pgoGYyMtTFu97FcCxxDvp/XNrracBwA7AVCwSecc5NAtq1rfbe3w7cDuCcOxw4A1junAsDT2CrDy/F6blFZIQbTPUugv4BroJeMhj0KGDw3v8K+FXMS+9GO0Ie7PpXgVeBn0YTJE/DKjAqYJABoU9KEqThWOJcpKMuAwbnXHddJ/+ruwm89xXAPdEhIjKoqN6FSM90t8KQlZCnEBEZIENtSV8BjgyULgMG7/0PE/UgIiLSvaEW4Mjw0aMcBudcBLgSmAdEWl/33n+6k2tDWGfIsUAtsN57vzsuTysiIiIDoqenJP4MbAJOB34EXA5sjL3AOTcNq89wKvAmsAcLLmY652qwkxZ/8t63xOfRhwcdURQRkaGgp6Whp3vvvwdUR3tFnA0s6HDNvwN/AaZ570/33n/Me3+R934hcC6QA3y8Pw/rnDvDOfeGc+4t59y1nXzfOed+Ff3+6865I3r6XhHpn8WLFyesPoCIJF5PVxgao7+WOefmA7uAybEXeO8/erA3e+9LgP/uywO2itZz+DV2RHMb8LJz7n7v/YaYy84EZkTHIuB/gUU9fK/0gsrTBksrTiIy2PQ0YLjROZcHfA+4H8gErou9wDl3QVc38N4v79MTtjkaeMt7/050vtuB87DeFq3OA27xdkj5RedcrnNuDBbcdPde6QWVpxURGVl6Wrjp99HfrgSmHuSypV3dAuhvwDAOeD/m623YKkJ314zr4XsBcM5dDVwNUFBQwA9+8IN+PXR3tm7dSllZWeDztPrjH//Y73vs2rULiE+29hVXXNHve8jg0LrqlKj/luNhqD3zUHteGJrPLJ1zPakY5py7rrPXvfc/ivsTHfwZlgGne+8/E/3648DR3vsvxlzzIPBT7/2z0a+fBL6BBTldvrczxcXFfvXq1YH8eVolunJbf/aYO5anbdWf8rRaeh8+hmIC71B75qH2vDA0n3kkc86t8d4Xd/a9nm5J/P/27j3orqo+4/j34SbXSFEuAZIGGAjGC1CiFHSsIFqhFArIRVABKZHW0DoWyq0zXtLpMIJQp6NAuKYdG26BESgSrikFKSAUQgCjFFKKSQOEcBEQSHz6x96HHF7Pec95z+19z3mfz0zm7Mtaa/9Ihn1+Z+2113qtant94ACGvCVRdbEtgX8Atra9n6RpwJ62LxlBzLU8C0yq2t8WWNpkmfWaqPs7Fi9e3PVBXP24mEymp42IGH+afSTxvep9SedQjGWo5XLgMuDMcv8XwJVAuwnDA8COkrYDfgUcCRw1pMz1wMxyjMIewMu2l0l6vom6PZXFZCIiop8028Mw1IbUH8vwfttXSTodwPYqSatbvM47ynZmAvOBtYFLbT8m6cTy/AXATcD+wJPA68Bxw9VtdM2pU6d2vRutX3+tZ3raiIjxpdmZHh+lGLgIxRfu5hQTONXymqT3VcpL+kPg5TbjBMD2TRRJQfWxC6q2DXyt2brRuvSCRESML832MBxQtb0KWG57VZ2y36B4NLCDpHsokovPtx7i4Muv9YiIGOsaLW+9Wbn56pBTEyRh+8WhdWw/JOmPgKmAgMW23x5aLtbIr/WIiBjrGvUwPEjxaEHAZGBlub0p8AywXaXgMBM37VQmF+3OwxARERGjpNHy1tsBSLoAuL4cB4Ck/SgWmapWmbhpC2Av4I5yf29gAe1P3BQRERGjpNkxDB+1fWJlx/ZPJM2qLmD7OABJNwLTbC8r9ydSrOMQNWQykxgkWWMkYnA1mzC8IOnvKFajNPBFYEWdslMqyUJpObBT6yFGRL/IGiMRg6vZhOELwDeB68r9u8pjtSyQNB+YS5FcHAnc2U6QMbakV6S/9GoW0X6ctTQimtfsTI8vAn/dZNmZkg4GPlkemm37uuHqRET/yqylEeNDo9cq/9H21yXdwJqJm95h+8A6VX9KMV+DgfvbjjIiWtarX+L9OGtpeikimteoh+Ffys9zmm1Q0uHA2RRvRgj4J0mn2L6mpQgjIiJi1DV6rfLB8vPfK8ck/R4wyfbCOtXOpHir4rmy/ObAbUAShogBl1lLuys9IjGa1mqmkKQFkiaUMz8+Alwm6dx6bVaShdKKZq8TEf0tYxYiBlezb0m81/Yrkv4cuMz2NyXV62G4ueotCYAjgJ+0G2hERESMnmYThnXKCZgOp3jkUJftUyQdCnycYgxD3pKIiBjHMqHXYGg2YfgOMB+4x/YDkrYHflmvsO15km6ttC9ps1oLVUVExODLhF6Dodl5GK4Grq7afwo4tFZZSV+lSDDeAH5L0ctgYPtWgyzHTlwJTAGWAIfbXjmkzCTgn4GtyuvOtv398ty3gBOA58viZ1TWxYiIiNo6MQlXJyf0yqDP0dVUwiBpJ+B8YEvbH5L0EeBA239fo/jJwAdtv9DBOE8Dbrd9lqTTyv1Th5RZBfxNubz2JsCDkm61/Xh5/jzbTb8eGhERrcuEXoOn2UcSFwGnABcC2F4o6V+BWgnDfwOvdya8dxwEfKrcnkMxx8O7EoZy/Ypl5farkp4AtgEeJyIiRqwTv+j7cUKvqK3ZhGFD2/dLqj62qk7Z04GfSroPeLNy0PZftRYiUPRsVJKBZZK2GK6wpCnAbsB9VYdnSvoy8DOKnoiVNaoiaQYwA2Dy5MlthBwRETE4RrJa5Q6U00NL+jzlr/kaLgTuAB6lGEvQFEm3UYw/GGrYtzJqtLMxMA/4uu1XysPnA7Mo4p8FfA/4Sq36tmcDswGmT5+elDgiok2Z0GswNJswfI3iS3RnSb8CngaOrlN2le1vjDQQ2/vWOydpuaSJZe/CROC5OuXWpUgWfmT72qq2l1eVuQi4caTxRUREazJmYTA0NQOj7afKL/TNgZ0pxhN8ok7xOyXNkDRR0maVP23GeT1wTLl9DPDjoQVUPC+5BHjC9rlDzk2s2j0YWNRmPBEREeNKo9UqJ1D0LmxD8SV9W7l/MsUU0T+qUe2o8vP0qmNtvVYJnAVcJel44BngsDK+rYGLbe9PMVHUl4BHJT1c1qu8PvldSbuWcSwBvtpGLBFRQ155ixhszaxWuRK4l2Ieg78F1gP+zPbDtSrY3q6jERZtrgA+XeP4UmD/cvtuijkfatX/UqdjioiIGE8aJQzb2/4wgKSLgReAybZf7XpkERERMWY0GsPwdmXD9mrg6SQLERER40+jHoZdJFVeTRSwQbkvwLYndDW6iIiIGBOG7WGwvbbtCeWfTWyvU7XdMFko13CIiIiIPtfUa5VtOLDL7UdEREQPdDthqPnWQkRERPSXbicMu3e5/YiIiOiBriYMtpteSyIiIiLGrm73MERERMQASMIQERERDTW7WuWISPoT4IPA+pVjtr/TjWtFRERE93W8h0HSBcARwEkUb0kcBmQx9IiIiD7WjUcSe9n+MrDS9reBPYFJXbhORERE9Eg3EoY3ys/Xy+Wn3wY6voJlRERE9E43xjDcKGlT4GzgIcDAxV24TkRERPRINxKG79p+E5gn6UaKgY+/aadBSZsBVwJTgCXA4bZX1ii3BHgVWA2ssj19JPUjIiKitm48kri3smH7TdsvVx9r0WnA7bZ3BG4v9+vZ2/aulWShhfoRERExRMd6GCRtBWxDsQT2bqxZR2ICsGGbzR8EfKrcngMsAE7tYf2IiIhxrZOPJP4YOBbYFji36virwBlttr2l7WUAtpdJ2qJOOQO3SDJwoe3ZI6yPpBnADIDJkye3GXZERMRg6FjCYHsOMEfSobbnjbS+pNuArWqcOnMEzXzc9tIyIbhV0s9t3zWSOMokYzbA9OnTPZK6ERERg6rjgx5tz2tlpkfb+9Y7J2m5pIll78BE4Lk6bSwtP5+TdB3wMeAuoKn6ERERUVu/zPR4PXBMuX0M8OMa191I0iaVbeCzwKJm60dERER9/TLT41nAZyT9EvhMuY+krSXdVJbZErhb0iPA/cC/2b55uPoRERHRnG7MwzB0pscVtDnTo+0VwKdrHF8K7F9uPwXsMpL6ERER0ZzM9BgRERENdWPQ46xy852ZHsvJmyIiIqJPdXLipkOGOYftazt1rYiI6A8LFiwY7RCiQzrZw/Cn5ecWwF7AHeX+3hQzKyZhiIiI6FOdnLjpOIDyMcS0ysyK5bwHP+jUdSIiIqL3ZHd2MkNJi2x/qGp/LWBh9bF+Iel54H9arP5eoJWxG63Uez/wQgvXit/V6r/bWDFW4u9lHN26VqfabbedXt1Lch/pnLHy/2Erft/25rVOdOMtiQWS5gNzKd6QOBK4swvX6bp6f2nNkDTb9oxe1JP0syGrc0aLWv13GyvGSvy9jKNb1+pUu+2206t7Se4jnTNW/j/stG68JTFT0sHAJ8tDs21f1+nr9IEbelwvOqPf//7HSvy9jKNb1+pUu+22k3tJ/xnIv/uOPZKQJDdorJkyMXL5ZRAR7cp9JBrp5NTQd0o6SdK71oSWtJ6kfSTNYc16DtFZsxsXiYgYVu4jMaxO9jCsD3wFOJpiKuiXgA0okpJbgB/YfrgjF4uIiIie6vhbEgCS1qUYcfuG7Zc6foGIiIjoqa4kDBERETFYurG8dURERAyYJAwRERHRUBKGASNpI0lzJF0k6ejRjici+pOk7SVdIuma0Y4lxoYkDH1A0qWSnpO0aMjxz0laLOlJSaeVhw8BrrF9AnBgz4ONiDFrJPcS20/ZPn50Io2xKAlDf7gc+Fz1AUlrUyzqtR8wDfiCpGnAtsD/lsVW9zDGiBj7Lqf5e0nEuyRh6AO27wJeHHL4Y8CT5a+At4ArgIOAZymSBsi/b0RUGeG9JOJd8oXSv7ZhTU8CFInCNsC1wKGSzmdA5zOPiI6qeS+R9D5JFwC7STp9dEKLsaQbq1VGb6jGMdt+DTiu18FERN+qdy9ZAZzY62Bi7EoPQ/96FphUtb8tsHSUYomI/pV7STQlCUP/egDYUdJ2ktYDjgSuH+WYIqL/5F4STUnC0AckzQXuBaZKelbS8bZXATOB+cATwFW2HxvNOCNibMu9JNqRtSQiIiKiofQwRERERENJGCIiIqKhJAwRERHRUBKGiIiIaCgJQ0RERDSUhCEiIiIaSsIQETVJWi3p4ao/pzWu1X1VcW0t6b5y+xlJz1fFOmVInU9JunfIsXUkLZc0UdLZkv5P0sm9/G+J6CdZSyIi6nnD9q6dbFDSOuVEQe2ojmuPst1jgem2Z9apcxewraQptpeUx/YFFtleBpwi6bU244oYaOlhiIgRkbRE0rclPSTpUUk7l8c3knSppAck/Zekg8rjx0q6WtINwC2SNpR0laSFkq4sewmmSzpe0nlV1zlB0rktxLeDpJslPSjpPyTtbPu3wNXAEVVFjwTmtvWXETGOJGGIiHo2GPJIovrL9gXbfwCcD1S68c8E7rD9UWBv4GxJG5Xn9gSOsb0P8JfAStsfAWYBu5dlrgAOlLRuuX8ccFkLcc8GTrK9exnbD8vjcymSBCS9B9gfmNdC+xHjUh5JREQ9wz2SuLb8fBA4pNz+LMUXfiWBWB+YXG7favvFcvsTwPcBbC+StLDcfk3SHcABkp4A1rX96EgClrQxsBdwtfTOqs3vKdt/QNLGkqYCHwD+0/bKkbQfMZ4lYYiIVrxZfq5mzX1EwKG2F1cXlLQHUD0+QNR3MXAG8HNa611YC3hpmETnCopehg+QxxERI5JHEhHRKfOBk1T+tJe0W51ydwOHl2WmAR+unLB9HzAJOIoWvtBtvwI8Lemwsn1J2qWqyFzgi8A+ZAnniBFJwhAR9Qwdw3BWg/KzgHWBhZIWlfu1/BDYvHwUcSqwEHi56vxVwD1tPC44Gjhe0iPAY8BBlRO2HwdepxhrkbciIkYgy1tHRE9JWptifMJvJO0A3A7sZPut8vyNwHm2b69T/9e2N+5CXN8Cfm37nE63HTEI0sMQEb22IXB32QNwHfAXtt+StKmkX1AMtqyZLJReqUzc1KmAJJ1N8agivQ4RdaSHISIiIhpKD0NEREQ0lIQhIiIiGkrCEBEREQ0lYYiIiIiGkjBEREREQ/8P7JOysNQWR5EAAAAASUVORK5CYII=\n",
"text/plain": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"flux_points_dataset.plot_fit();"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Stack observations\n",
"\n",
"And alternative approach to fitting the spectrum is stacking all observations first and the fitting a model. For this we first stack the individual datasets:"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [],
"source": [
"dataset_stacked = Datasets(datasets).stack_reduce()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Again we set the model on the dataset we would like to fit (in this case it's only a single one) and pass it to the `~gammapy.modeling.Fit` object:"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [],
"source": [
"dataset_stacked.models = model\n",
"stacked_fit = Fit([dataset_stacked])\n",
"result_stacked = stacked_fit.run()\n",
"\n",
"# make a copy to compare later\n",
"model_best_stacked = model.copy()"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"OptimizeResult\n",
"\n",
"\tbackend : minuit\n",
"\tmethod : minuit\n",
"\tsuccess : True\n",
"\tmessage : Optimization terminated successfully.\n",
"\tnfev : 31\n",
"\ttotal stat : 29.88\n",
"\n"
]
}
],
"source": [
"print(result_stacked)"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"Table length=3\n",
"
\n",
"
name
value
unit
min
max
frozen
error
\n",
"
str9
float64
str14
float64
float64
bool
float64
\n",
"
index
2.5876e+00
nan
nan
False
5.867e-02
\n",
"
amplitude
2.6900e-11
cm-2 s-1 TeV-1
nan
nan
False
1.234e-12
\n",
"
reference
1.0000e+00
TeV
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 2.5876e+00 nan nan False 5.867e-02\n",
"amplitude 2.6900e-11 cm-2 s-1 TeV-1 nan nan False 1.234e-12\n",
"reference 1.0000e+00 TeV nan nan True 0.000e+00"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"model_best_joint.parameters.to_table()"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"Table length=3\n",
"
\n",
"
name
value
unit
min
max
frozen
error
\n",
"
str9
float64
str14
float64
float64
bool
float64
\n",
"
index
2.5912e+00
nan
nan
False
5.881e-02
\n",
"
amplitude
2.6933e-11
cm-2 s-1 TeV-1
nan
nan
False
1.234e-12
\n",
"
reference
1.0000e+00
TeV
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 2.5912e+00 nan nan False 5.881e-02\n",
"amplitude 2.6933e-11 cm-2 s-1 TeV-1 nan nan False 1.234e-12\n",
"reference 1.0000e+00 TeV nan nan True 0.000e+00"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"model_best_stacked.parameters.to_table()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, we compare the results of our stacked analysis to a previously published Crab Nebula Spectrum for reference. This is available in `~gammapy.modeling.models.create_crab_spectral_model`."
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 32,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEKCAYAAADq59mMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdd3hU1brH8e9O7yG9QCBAgJA6aRQpUqSpgAqoiHqxoaJ4LAgCtssREOUoctDjUY+i99gVUFCkF6WlkU5ICCS09N4zmVn3j5AxkRAyySQEWJ/n8XnIlDVrB8wv+917vUsRQiBJkiRJ+jC62hOQJEmSrj0yPCRJkiS9yfCQJEmS9CbDQ5IkSdKbDA9JkiRJbzI8JEmSJL2ZXO0JdAVnZ2fh7e19tachSZJ0TYmJiSkQQri09NwNER7e3t5ER0df7WlIkiRdUxRFybrcc7JsJUmSJOlNhockSZKkNxkekiRJkt5uiGseknS9UKvVnDt3jpqamqs9Fek6YmFhQa9evTA1NW3ze2R4SNI15Ny5c9ja2uLt7Y2iKFd7OtJ1QAhBYWEh586do2/fvm1+nyxbSdI1pKamBicnJxkcksEoioKTk5PeZ7MyPFohhODnjJ+p09Rd7alIko4MDsnQ2vNvSoZHK2JyY1j2xzJm/DyDQ+cPXe3pSFK3sGLFCvz9/QkKCkKlUnH06FEA1q5dS1VVVbvG3LBhA08//XS73rtv3z5uv/32dr23rebOncsPP/yg9/teffVVdu3a1QkzarsxY8bo1rmtXLnSYOPK8GhFuHs4/7rlX2iFlsd3Pc7z+54npzLnak9Lkq6aw4cPs3XrVmJjY0lISGDXrl14eXkBHQuP69Xy5cu55ZZb9H5ffX19J8xGhkeXGtlzJJumb2JByAJ+P/c70zZP45PET1Br1Fd7apLU5bKzs3F2dsbc3BwAZ2dnPD09WbduHRcuXGDs2LGMHTsWgCeffJLw8HD8/f157bXXdGNERUVx0003ERwczJAhQygvL2/2Gb/88gvDhw+noKCAHTt2MHz4cEJDQ5k1axYVFRUA/Pbbb/j6+jJy5Eg2btzY4lwzMzMZNWoUoaGhhIaGcuhQQ/Vg3759jBkzhpkzZ+Lr68ucOXNo3FF1+fLlREREEBAQwLx58/jrTqu7d+/mzjvv1H29c+dO7rrrLjQaDXPnziUgIIDAwEDeffddoPkZy0svvYSfnx9BQUEsXLjwkvm+/vrrzJs3j4kTJ/Lggw+Sn5/PjBkziIiIICIigoMHDwKwf/9+VCoVKpWKkJAQysvLLzn7evrpp9mwYUOz8V966SWqq6tRqVTMmTOnxe+ZPuTdVm1gZmzGvKB53NbvNt6KfIv3Yt/jp5M/sXToUoZ7Dr/a05NuUP+7JZmUC2UGHdPP047Xpvpf9vmJEyeyfPlyBg4cyC233MI999zDzTffzDPPPMM777zD3r17cXZ2BhrKW46Ojmg0GsaPH09CQgK+vr7cc889fPvtt0RERFBWVoalpaVu/E2bNvHOO+/w66+/otFoeOONN9i1axfW1tasXr2ad955h0WLFvHYY4+xZ88efHx8uOeee1qcq6urKzt37sTCwoL09HRmz56tK98cO3aM5ORkPD09GTFiBAcPHmTkyJE8/fTTvPrqqwA88MADbN26lalTp+rGHDduHE899RT5+fm4uLjw2Wef8dBDDxEXF8f58+dJSkoCoKSkpNlcioqK2LRpE6mpqSiKcsnzjWJiYvjjjz+wtLTkvvvu47nnnmPkyJGcOXOGSZMmcfz4cdasWcP777/PiBEjqKiowMLC4kp/rQC8+eabrF+/nri4uDa9/krkmYceetr05L1x7/H++PfRCA3zds5j4f6FspQl3TBsbGyIiYnho48+wsXFhXvuueeS33Abfffdd4SGhhISEkJycjIpKSmcOHECDw8PIiIiALCzs8PEpOF32L1797J69Wp++eUXHBwcOHLkCCkpKYwYMQKVSsXnn39OVlYWqamp9O3blwEDBqAoCvfff3+Ln69Wq3nssccIDAxk1qxZpKSk6J4bMmQIvXr1wsjICJVKRWZmpm4OQ4cOJTAwkD179pCcnNxsTEVReOCBB/jvf/9LSUkJhw8fZsqUKfTr149Tp06xYMECfvvtN+zs7Jq9z87ODgsLCx599FE2btyIlZVVi3OeNm2aLkx37drF008/jUqlYtq0aZSVlVFeXs6IESN4/vnnWbduHSUlJbrvX1eTZx7tMLrXaIZ6DGVD0gY+TvyYA+cO8ETwEzww+AFMjdu+yEaSOqK1M4TOZGxszJgxYxgzZgyBgYF8/vnnzJ07t9lrTp8+zZo1a4iKisLBwYG5c+dSU1ODEOKyd/Y0/gBOS0sjPDwcIQQTJkzg66+/bva6uLi4Nt0d9O677+Lm5kZ8fDxarbbZb+iNZbfG46mvr6empob58+cTHR2Nl5cXr7/+eou3rz700ENMnToVCwsLZs2ahYmJCQ4ODsTHx7N9+3bef/99vvvuOz799FPde0xMTIiMjGT37t188803rF+/nj179lwytrW1te7PWq2Ww4cPNzszg4by02233cavv/7KsGHD2LVrFyYmJmi1Wt1rumIRqTzzaCdzY3MeD36cn+74iWEew3g35l1mbJnBkewjV3tqktRpTpw4QXp6uu7ruLg4+vTpA4Ctra3u+kVZWRnW1tbY29uTm5vLtm3bAPD19eXChQtERUUBUF5errs43KdPHzZu3MiDDz5IcnIyw4YN4+DBg5w8eRKAqqoq0tLS8PX15fTp02RkZABcEi6NSktL8fDwwMjIiP/7v/9Do9G0emyNP3CdnZ2pqKi47N1Vnp6eeHp68sYbb+hCs6CgAK1Wy4wZM/j73/9ObGxss/dUVFRQWlrKrbfeytq1a9tUOpo4cSLr16/Xfd34noyMDAIDA1m8eDHh4eGkpqbSp08fUlJSqK2tpbS0lN27d7c4pqmpKWq1Ya7XyjOPDupp05N149Zx4NwB3ox8k8d2PMYk70ksDF+Iu7X71Z6eJBlURUUFCxYs0JVLfHx8+OijjwCYN28eU6ZMwcPDg7179xISEoK/vz/9+vVjxIgRAJiZmfHtt9+yYMECqqursbS0bHYr66BBg/jyyy+ZNWsWW7ZsYcOGDcyePZva2loA3njjDQYOHMhHH33EbbfdhrOzMyNHjtRda2hq/vz5zJgxg++//56xY8c2+62+JT169NCVuby9vXWltZbMmTOH/Px8/Pz8ADh//jwPPfSQ7rf/VatWNXt9eXk506dP1519NV5Qb826det46qmnCAoKor6+ntGjR/Phhx+ydu1a9u7di7GxMX5+fkyZMgVzc3PuvvtugoKCGDBgACEhIS2OOW/ePIKCgggNDeXLL7+84hxao/z1boLrUXh4uOiK/TxqNbV8mvQp/0n8D0aKEU8GP8n9g++XpSzJYI4fP87gwYOv9jRueE8//TQhISE88sgjV3sqBtPSvy1FUWKEEOEtvf66LlspijJVUZSPSktLu+TzzI3NeTL4STZP38xQj6G8E/MOM7fM5Gj20S75fEmSOl9YWBgJCQmXvVB/o7iuw0MIsUUIMc/e3r7dYxQWFlJZWanXe3rZ9uKf4/7J+nHrqdXU8uiOR1m0fxG5lbntnockSd1DTEwMBw4caHbR/UZ0XYeHIeTk5LBjxw5OnDih96rPm71uZvP0zcwPns/uM7uZtnkanyd/jlorFxhKknRtk+HRBpWVlURHR7N3717y8vL0eq+FiQVPqp5k8x2biXCPYE30Gmb9PIuonKhOmq0kSVLnk+HRRkIIcnNz2bt3L9HR0VRXV+v1fi9bL9aPX88/x/2TGk0ND29/mMUHFpNXpV8YSZIkdQcyPPSkVqtJTU1lx44dnDp1qtnCnLYY4zWGzdM380TwE+zK2sXUTVNlKUuSpGuODI92Ki8v5/Dhw+zdu5eioiK93mthYsFTqqfYPH0zYW5hrIlew91b7palLOmaYWNjc9nnLly4wMyZM684hiE7vF5Oe9u1R0dH88wzz3TCjNquaZv6zZs3N2uv0h3I8LiCEwW1XG4pjBCC7Oxsdu3axbFjx3QLmdrKy86L98e/z7qx66hSV+lKWflV+QaYuSRdHZ6enm3a+6IrwqO9wsPDWbdund7vE0LoXY1oCxke15jozCLe+S2euNNnyK42vuzr6urqSE5OZvv27WRlZen1j0dRFMb2HsvmOzYzL2geO7N2MnXzVL5I/oJ6bef09JckQxBC8OKLL+rakH/77bdAQyv0gIAAoOG357vuuovJkyczYMAAFi1aBFy5Pfjl2rl7e3vz2muvERoaSmBgIKmpqQBERkZy0003ERISwk033cSJEyeajafVahkwYAD5+fm6r318fCgoKOD7778nICCA4OBgRo8eDTQ/Y2mpBXpTmZmZDB48mPnz5xMaGsrZs2d5++23iYiIICgoSDf/yspKbrvtNoKDgwkICNB9v7y9vSkoKAAaznjGjBnTbPxDhw7x888/8+KLL6JSqXRtWa422Z6kFcFePVjnsZPQol/YeSaUn61m4ufeA1vTlk9FysrKOHjwIJ6enoSEhKDP+hJLE0sWhCxgev/prIxcydvRb7Pp5CaWDV1GuHuLCzwlCT677dLH/O+AIY9BXRV8OevS51X3QcgcqCyE7x5s/txDv7T5ozdu3EhcXBzx8fEUFBQQERGh++HbVFxcHMeOHcPc3JxBgwaxYMGCK7YHb6mde1BQENDQeyo2NpYPPviANWvW8Mknn+Dr68uBAwcwMTFh165dLF26lB9//FE3npGREffffz9ffvklzz77LLt27SI4OBhnZ2eWL1/O9u3b6dmzZ4ut0tvSAv3EiRN89tlnfPDBB+zYsYP09HQiIyMRQjBt2jQOHDhAfn4+np6e/PJLw/e4rYuXb7rpJqZNm8btt9/epnJgV5FnHq0wNTbCYsxCIh2mM9I4hXdqXsH61GYic+pRX+bkQqvVcu7cOXbs2EFiYqLeTch62/XmX+P/xdqxa6lUV/LQ9od46feXZClL6nb++OMPZs+ejbGxMW5ubtx88826hodNjR8/Hnt7eywsLPDz8yMrK+uKY7fUzr3RXXfdBTSs9G5spV5aWsqsWbMICAjgueeeu6SVOsDDDz/MF198AcCnn37KQw89BMCIESOYO3cuH3/8cYvNE9vSAr1Pnz4MGzYMgB07drBjxw5CQkIIDQ0lNTWV9PR0AgMD2bVrF4sXL+b333/X65fL7kieeVyBMDYnzWkiWfbDGZT/C3Mqd1JbYsLa0vuY4FZNkH0dRi10h66trSU+Pp4zZ86gUqnw9PRs8ybziqIwvvd4bvK8iY8TPmZD8gb2nd3HU6qnmO07GxMj+dcmXdTamYKZVevPWzvpdabxV23ti9dS+/PWXK6d+1/HazrWK6+8wtixY9m0aROZmZmXlH4AvLy8cHNzY8+ePRw9elTXGPDDDz/k6NGj/PLLL6hUqkvOhlpqge7r69vsNU2bLgohWLJkCY8//vglc4iJieHXX39lyZIlTJw4kVdffbVZO/WuaKVuKPLMo41qTWxJ8LiXX72WUtpzPNbGWjLOXyA94zhnKi//bSwuLmb//v0cPHhQt4VmW1maWPJM6DNsmr4JlYuKt6Le4u6tdxOTG9PRw5GkDhs9ejTffvstGo2G/Px8Dhw4wJAhQ9r8/su1B79cO/fWlJaW0rNnT4DLbk4F8Oijj3L//fdz9913Y2zccB0zIyODoUOHsnz5cpydnTl79myz97TUAr01kyZN4tNPP9X9/37+/Hny8vK4cOECVlZW3H///SxcuFDXtt3b25uYmIb/p5uW2ppq2u6+u5DhoadSc0/cbS14on85y2y38XfWM/H8u8Rm5lBU1/K3U6vVkpmZyfbt20lNTdW7zUkfuz7865Z/sXbMWirqKpj721yW/bGMguoCQxySJOmlvr4ec3Nz7rzzToKCgggODmbcuHG89dZbuLu3fRuCxvbgf71gHhwcrGvn/vDDD+vaubdm0aJFLFmyhBEjRrS6b8e0adOoqKjQlawAXnzxRQIDAwkICGD06NEEBwc3e8/atWt1F9QtLS2ZMmVKq3OZOHEi9913H8OHDycwMJCZM2dSXl5OYmIiQ4YMQaVSsWLFCl5++WUAXnvtNf72t78xatQoXaD91b333svbb79NSEhIt7lgLluyX0FycjLHjh1r8TlFaOlTepjAwq3YizJ+1gxnj+2dBLhZY3H5m7NwdnYmJCQENzc3vedTpa7ik8RP+Cz5MyyMLXg65GnuGXSPLGXdILpDS/b4+Hgee+wxIiMjr+o82iM6OprnnnuO33///WpPpduRLdm7kFCMyOwxgm19XyPafgqTjKNxKo3n3TR7IovM0FwmlwsKCtizZw+HDx/Wu82JlalVQylr2iaCXIJ4M/JN7tl6D8fyWg44STKkDz/8kNmzZ/PGG29c7ano7c0332TGjBmXbNQktY8887iC1s48/spKXcwptT1bc+zwrTmGt2kJtW5D6G93+e+xlZUVgYGB9OvX77KnrJcjhGDXmV28FfUWOZU5TOs/jefCnsPZ0lmvcaRrR3c485CuT/LMw8AcHBxabcXQVJWpA+5WRjzSt4InbQ/xkrKB2TkrSD51iryalr/VVVVVREZGsmfPHgoLC/Wam6IoTOgzgZ+m/8SjgY/y6+lfmbZpGl8d/0ouMJQkqVPJ8LgCT09PJk2aRHBwcJtDRFEgsdeD7HF7jB4malZp/0F41gfEniuisv7S23UbO/bu3LmT6OhovducWJla8bfQv7Fx2kYCnANYFbmKe7feS1xeywuwJEmSOkqGRxtYWloSGBjIpEmTUKlUbQsRReGCrYpdfZdxyGEmocanqCor5J00e37PN6e+hUWG9fX1pKamsm3bNk6fPq13j5y+9n3594R/s+bmNZTUlvDAtgd4+Y+XKazW74xGkiTpSuQ1j3aorq4mIyODkydPtnnthqmmivN11mzPs2J01Q7sjevIcR3DIHuFltYOKoqCh4cHoaGh9OjRQ+85Vqmr+HfCv/ki+YuG1iehC7h74N0YG+l3XUXqXuQ1D6mzyGseXcDS0pKAgADdmUjT1aWXoza2wtVS8ECfSm63zeBpox+Yn/8qmacSOFfZcinrwoULbN++nbi4OOrq6vSao5WpFc+FPceP03/Ez9mPlUdXMvuX2bKUJXVYTk4O9957L/3798fPz49bb72VtLQ0vcZoawm4LfLz8xk6dCghISHyFtwuJMOjAxpDZPLkyW0OEWi4HrLd81nqTWx5RfybKefXEJOVR6n60hBRq9UkJSXx22+/cfbs2Ta3hGjUz74fH0/4mDU3r6GwppAHtj3AqwdfpahGvz1IJAkafqm58847GTNmDBkZGaSkpLBy5Upyc3Obva61hXrt0drC2t27d+Pr68uxY8cYNWpUm8Yz9PxuRDI8DKAxRKZMmdLmEMm3GsAB74Xsc36QnsYlnKo0ZW2aPbtyzKlt4d91WVkZBw4c4MCBA3q3KVAUhUnek9hyxxYeCniILRlbuH3T7XyT+g0arfyfSGq7vXv3YmpqyhNPPKF7TKVSMWrUKPbt28fYsWO57777CAwMBOCOO+4gLCwMf39/Pvroo2ZjvfDCC4SGhjJ+/Hhdq/Sm5s6dy/PPP8/YsWNZvHgxGRkZTJ48mbCwMEaNGkVqaipxcXEsWrSIX3/9FZVKRXV1NTt27GD48OGEhoYya9YsXWnZ29ub5cuXM3LkSL7//vtWX9dS2/fGlemBgYEEBQXpWolcbpzrnbzm0Qlqamo4efIk6enpVFZWXvH1itBQpDZlR64lsyu/QGNkRqrTZAY5mrbYdNHMzAw/Pz98fX1b7PB5JadKTrHy6EqO5hxlsONglg1bRrBL8JXfKF11TevSqyNXk1rUep8lffk6+rJ4yOLLPr9u3TpOnz7Nu+++e8lz+/bt47bbbiMpKYm+ffsCUFRUhKOjI9XV1URERLB//36cnJxQFIX//ve/zJkzh+XLl5OXl8f69eubjTd37lwKCgr46aefMDY2Zvz48Xz44YcMGDCAo0ePsmTJEvbs2cOGDRuIjo5m/fr1FBQUcNddd7Ft2zasra1ZvXo1tbW1vPrqq3h7ezN//nwWLVp0xde98MILLFiwgA8++IDY2Fg++eQTFi9eTG1tLWvXrgUa+tZpNJrLjnOt0feah+xp0QksLCwICAjAx8enTSEiFGMczLTc06sc32w1wVV7KS/6nQ3F06lwG4G3bfPX19XVERcXR1ZWFqGhoXh4eOg1v349+vHxxI/ZnrWdtyPf5v5f7+dOnzt5NuxZHC0c23PIkgTAkCFDdMEBDWGzadMmAM6ePUt6ejpOTk4YGRlxzz33AHD//ffr2qz/1axZszA2NqaiooJDhw4xa9af+5O0dEv7kSNHSElJ0fXDqqurY/jw4brnGz/zSq9r2vZ948aNAOzatYtvvvlG9xoHBwe2bt3a6jjXMxkeneivIZKWlkZVVdXl36AYkeJ5L+drRjEodxPPqr/idPZu1hY8Tl9PD5zNm9+6W1xczN69e+nTp49e11ygoZQ12Xsyo3qO4sP4D/lvyn/ZfWY3z4Q8w8yBM+VdWdeA1s4QOou/v3+rW8w2/Te4b98+du3axeHDh7GysmLMmDGXbTl+ue0KGsfTarX06NHjsptHNRJCMGHCBL7++utWx7vS61pq+y6EuGSeVxrneiaveXSBxhC59dZbUalUWFlZtfr6UoueRPZ+ih3u81FMLYmtdmFduh3bL5hSrWn+j1er1XL69Gm2bdtGamqq3hcCrU2teSH8BX6Y9gO+jr68cfQN7vv1PhLyE/Q+Tun6N27cOGpra/n44491j0VFRbF///5LXltaWoqDgwNWVlakpqZy5MgR3XNarVYXQl999RUjR45s9XPt7Ozo27cv33//PdDwQzs+Pv6S1w0bNoyDBw9y8uRJoKGDQ0t3grX1dU1NnDixWWmtuLi4XeNcL2R4dCG9QkRRyLPx5w/vhcwZaExoj1rml7+La8ZXHM+vvKTpYk1NDdHR0ezYsYO8vDy959a/R38+mfgJb41+i4KqAub8OofXD71OcU2x3mNJ1y9FUdi0aRM7d+6kf//++Pv78/rrr+Pp6XnJaydPnkx9fT1BQUG88sorup32oOEMIDk5mbCwMPbs2dOmawRffvkl//nPfwgODsbf35+ffvrpkte4uLiwYcMGZs+eTVBQEMOGDWtx/422vq6pl19+meLiYl179r1797ZrnOuFvGB+FdXU1JCenk56enrr5SwaLqr75P1CSPkeNELhS6Zw3mU8/eyNL1lkaGRkRL9+/VCpVC3ut3wllepKXSnL2syaZ0KeYcaAGbKU1Q3IRYJSZ9H3grkMj26gMUTS0tKu2KLdqq4Qn7wtBNVEkS/sedn4eXw83PCwvLRcZWlpiUqlom/fvhgZ6X+SebL4JCsjVxKVE4W/kz/Lhi4j0CVQ73Ekw5HhIXUWGR4t6O7h0UifEHGoPk3PvD08UjWfIo0lo3sUMczNGFvTS/8+3dzcCAsLw9FR/zuphBBsO72NNdFrKKgu4K4Bd/G30L/hYOGg91hSx8nwkDrLdRceiqL0A5YB9kKImRcfswY+AOqAfUKIL1sb41oJj0Y1NTWkpaWRnp5+xRCp1ij8nmvC3ytf4bTwYJftDPq7OWL6lxMNY2NjBgwYQFBQEGZmZnrPqaKugg/iP+Cr419hY2bD30L/xowBMzBS5GWzriTDQ+osBgsPRVF+bsPnFQkh5l7uSUVRPgVuB/KEEAFNHp8MvAcYA58IId680gcpivJDk/B4ACgRQmxRFOVbIcQ9rb33WguPRm0NESNRj1fBfkJKt2EhavlRjOW48630c7C4ZJGhtbU1KpUKb2/vy94e2Zr04nRWHF1BTG4MAU4BLBu2jADngCu/UTKI48eP4+vr266/O0m6HCEEqampBguPdODRVj5PAd4XQvhf9gWKMhqoAL5oDA9FUYyBNGACcA6IAmbTECR/3R/yYSFE3sX3NQ2PJcA2IUScoihfCSHua2We12x4NKqpqeHEiROcPHmy1RAx11TQJ3cboVUHqBbmzFdewdfDkT7Wza+HNHbsDQsLw97eXu/5CCHYemor/4j+B0U1RcwcOJNnQp6hh4X+3X8l/Zw+fRpbW1vdKm1J6ighBIWFhZSXlzdb4AntD4+7hRDftfahbXyNN7C1SXgMB14XQky6+PWSiwfQ6sbCLZx5FAshtiqK8o0Q4t4WXj8PmAfQu3fvsKysrNaGvyY0hkh6evplF1sB2NTm0iP/MM9UPEBZvQkTbc8Q4G6Ho3nzv2sTExMGDx6Mv79/u9qclNeV80HcB3yd+jW2ZrY8G/osdw64U5ayOpFarebcuXOt/v1Lkr4sLCzo1asXpqamzR6/qtc8WgiPmcBkIcSjF79+ABgqhHj6Mu93AlbQcKbyiRBi1cVrHuuBGuCP6+2ax5XU1NSQmprKyZMnW/0hUqeF2Lx6VpQvIUN4ssX6bnq698TiL3fc2traEhoaipeXV7vmk1acxoojK4jNiyXQOZBlw5bh73TZE1JJkq4R7Q4PRVF60VBSGgl4AtVAEvALDWWjK25110J4zAIm/SU8hgghFuhxTHq53sKjUXV1te5M5HJb1ypCi0dJJKrCLThSwjbtUI46TKevsz3GTaoeiqLQq1cvwsLC2rXXwl9LWXcPupsFIQuwN9e/LCZJUvfQ3rLVZ0BPYCsQDeQBFsBAYCwQBrwkhDhwhQ/3xgBlq464XsOjUXV1te5M5HIhYqytxSt/N2HlOzESgpm8Rai7FQNsm++TYGpqip+fH35+fhgb678osLGU9VXqV9ib2fNs2LPc4XOHLGVJ0jWoveERIIRIamVQM6C3EOLkFT7cm+bhYULDBfPxwHkaLpjfJ4RIvvKhtM/1Hh6NGkMkPT39sjsPWqhLMCtM4vWyWymqM+Yu60R6u/fE2bL5xVd7e3vCwsJabDvRFieKTrDy6Epi82IJcg5i2bBl+Dn5tWssSZKuDoNd81AUxQHwEkK0qWueoihfA2MAZyAXeE0I8R9FUW4F1tJwh9WnQogVbXdHI9EAACAASURBVJ6EHhRFmQpM9fHxeSw9Pb0zPqJbqq6uJiUlhYyMjMuGSL0W0vPL+d+ypWQJN763nIWD+yCsm1wvUxSFPn36EBoaesVmji0RQrDl1Bb+Ef0PimuKZSlLkq4xHQoPRVH2AdNoaN8eB+QD+4UQzxt4np3mRjnz+KuqqipdiKjV6ktfIATO5UkE5W/CU+RySOvPbruZeLm6YtKkymRmZkZgYCCDBg1qV5uTsroy3j/2Pt+c+AZ7M3ueC3uO6T7TZSlLkrq5jobHMSFEiKIoj9Jw1vGaoigJQoigzphsZ7hRw6NRZWUlx48fv2yIKEKDZ+EfhJb8irmoZbL2PYa5m+Bvp27WdNHJyYmwsDBcXV3bNY/UolRWHFlBXH4cQS5BvDz0ZQY7ydXSktRddTQ8EoGJwOfAMiFElAyPa1NlZSXJycmcPn26xRAx1VShKT7Ne8XDyK014RHL/Vi5+eFm8+caEEVR6NevHyEhIe3q2KsVWrZkbOGdmHcoqS3h7oF383TI07KUJUndUEfDYxbwCg3rKeZf7DX1thBihuGn2jlkeDRXWVlJUlISmZmZLYaIVkB2fh6Ly/6XHOHAd6Z3YuwRir35n6chFhYWBAcH4+Pj066Vzk1LWT3Me/Bs6LOylCVJ3cw13RjREGR4tKyiokIXIo1bbTZlX5mBX+4m+mtPk6LtwxbrWTi798W8yR28zs7ORERE4OTk1K45NC1lqVxULBu2DF9H3/YekiRJBnTDhseNereVvsrLy3Uhcsk2tkLgXBJLaOFPmIhaJmvWMsJNENKjTtd00cjICB8fH4KDg3V7P+tDK7T8nPEz78a8S0ltCfcOupenQp7CzszOAEcnSVJ73bDh0UieebRNWVkZiYmJnDlz5pIQMdbWUV2Wx2eFgzlfbcQii82Uuw7Hw+7P6x6WlpaEhITQt2/fdpWySmtLWX9sPd+lfUcP8x48H/Y80/pPkw0AJekqkeEhw0MvrYWIEFBWkMUTJW9TgSXfGU+lxnMEDk0aZrm5uREREUGPHu3rsptSmMKKoytIyE8g1DWUpUOXMshxUIeOSZIk/XWkt5UvDS1KjgohKpo8PlkI8ZvBZ9pJZHi0T0lJCYmJiZw9exattnkbM+uaC/jk/ExgfSJZWlc2W87A0j0Aq4uLDI2MjBg0aBCBgYHt2nxKK7T8dPIn3o15l9K6Umb7zuYp1VPYmtka4tAkSWqD9rYneQZ4CjgOqIC/CSF+uvhcrBAitJPma3AyPDqmuLhYFyJ//ffSo/w4gXmbqNMI7tSsYrSLmiGOtbpFhlZWVrpSVnuU1pbyz2P/5LsT3+Fo4cgL4S9we7/bZSlLkrpAe8MjERguhKi42J/qB+D/hBDvNS4c7KwJG5oMD8MoKioiMTGRc+fONQsRRWipqCzju/xe5FTWs9z8S045T8C9h51ukaG7uzsRERHt2nwKILkwmZVHVpJQIEtZktRV2hseKUIIvyZf29AQICnAOCGEqjMm2xlkeBhWQUEBCQkJZGdnNwsRIUBTmMa9xR+gAD8YTaLAfTxO1g1lKyMjI3x9fQkKCmrX5lNaoWXzyc28G/Mu5XXlzPadzXzVfFnKkqRO0t7w2AM8L4SIa/KYCfApMEcIoX+/7i4mb9XtXPn5+cTHx5Obm9ssRMzVxfTJ3sqQuiMUCDt+ML0D4TkU24uXPqytrQkNDaVPnz7t+tzS2lLWxa7j+7TvZSlLkjpRe8OjF1AvhMhp4bkRQoiDhp1m55FnHp0rNzeXhIQE8vLymoWIddUZBudsIq/enEfrFzPKuZYRzjWYXbwe4uHhQUREBHZ27VvPkVyQzIqjK0gsSCTUNZRlw5Yx0GGgIQ5JkiQMcKtuYyt2GjrrAiCEiDXYDDuZDI+ukZ2dTWJiInl5eX8+KATlNXVsyXektKyUlWYbiHKchquTK0ZKwz7qvr6+BAYGtmvzKa3QsjF9I+/Fvkd5XTn3Db6P+cHzsTHTfzdESZKa62hvq78Dc4EMoPHFQggxzpCT7EwyPLqOEIILFy6QkJBAYWFhs+eMClO5vfhTrEUVvyg3c8p1Cs52DT/kbW1tCQkJoXfv3u363JKaEtYdW8cPaT/gZOnEC+EvcFvf22QpS5I6oKPhcQIIFEK0vKvQNUCGR9cTQnDu3DkSExMpKirSPW5SX4lHzg5uqt5LLaZ8bzKV8p7jcDRrWEfi6elJREQEtrbtuwieVJDEG0feILkwmTC3MJYNXcYAhwEGOSZJutF0NDx+BJ4UQuS1+sJuTIbH1SOE4MyZMyQkJFBaWqp73KI2j77ZW0ir7cHL9Y9wk1MtNztXYWGiYGJigp+fH/7+/u0qZWm0GjaebChlVdRVMGfwHJ4MflKWsiRJTx0Nj3DgJyAJqG18XAgxzZCT7EwyPK4+rVZLVlYWiYmJlJWV6R4vr9OyI88GpTSL/zX9gt32M3Bw9cZYaShlhYaG4uXl1a7PLK4p5r3Y99iYvhFnS2cWhi9kSt8pspQlSW3U0fBIBv4NJAK6HhVCiP2GnGRnkuHRfWi1Wk6fPk1iYiIVFbqON5gVpzKm8EtcKWKPCCfBeTrODo4A9OzZk/Dw8HaXshLyE1hxdAUphSlEuEewdMhSfBx8DHI8knQ962h47BdC3NwpM+tkcp1H96XRaDh58iTJyclUVVUBYKSpwyVvHyMqfsOEer4zup3cnpNxs9BiYmLC4MGDCQgIaHcp68f0H3kv9j2q1FXc73c/TwQ/gbWptaEPTZKuGx0Nj3doKFf9TPOylbxVV+owjUbDiRMnSElJoaamBgDTulJ65mzjUFVPPqyfSrhDDRNcK7E0NcbGxoawsLAOl7J+TP8RV0tXFkYsZLL3ZFnKkqQWdDQ89rbwsLxVVzIotVpNamoqqamp1NY2/I5SWa+wN88C15JYnjf9nl+tZ2LlNhgTY6XDd2XF58ez4sgKjhcdZ4j7EJYOXUr/Hv0NeUiSdM2T+3nI8Lhm1NXVkZycTFpamm5/dfPSdIbkf0cfLhAjBnHIYSY9nDwxMTHWlbLa0ytLo9XwQ9oPvHfsParV1Tzg9wBPBD+BlamVoQ9Lkq5JHT3zWAm8JYQoufi1A/CCEOJlg8+0k8jwuPbU1NSQlJTEyZMnqa+vRxEa7POPcFPpFhyVcv7L7ZzzvJ1eVhqsra0JCwtr9wLDopoi3V1ZrlauvBjxIpP6TJKlLOmG19HwuKT9utzPQ+oqVVVVJCYmkpGRgVarxVhTjXP2bn6pGMDm+mGE21UwwbUCawvzDvfKisuLY+XRlRwvOs5Qj6EsHbqUfvb9DHxEknTt6Gh4JAARQojai19bAtFCCH+Dz7STyPC49pWXlxMfH09WVhZCCGo0cCDfgqCS7TxovIOtFtPBIxwLMxMGDRpEcHBwu0tZ36d9z7pj66iuv1jKCpKlLOnG1NHwWARMAz6jobfVw8DPQoi3DD3RziLD4/pRWlpKXFycbkMqs/IsAvM2MlicJE14sdtuJtauPlhbWRISEkK/fu07cyisLmRt7Fo2n9yMm5UbL0a8yMQ+E2UpS7qhGKKr7mTgFkABdgghtht2ip1Lhsf1p6CggPj4eLKzs0EIrIviGFr8E57k8zlTyXK/nX429bi4uBAREYGjo2O7PicuL44VR1eQWpTKMI9hLBm6RJaypBtGe/fzUMQVkqUtr7ma5CLB619OTg7x8fHk5+ejaNXY5/zB1+UBHFQPZIhNIeNdq7CxtsXHx4eQkBDMzMz0/ox6bT3fnfiO9cfWU62p5kG/B3k86HFZypKue+0Nj33Aj8BPQogzTR43A0YC/wPsFUJsMPSEDU2eeVz/zp07R3x8PMXFxai1cKjQnPFFX3Gb0SG2mN1OtccobG2sCA4OZuDAge0qPxVUF7A2Zi0/ZfyEm5UbiyIWMaHPBFnKkq5b7Q0PCxqub8wB+gIlgCVgBOwA3m+6RW13JsPjxiCEIDMzk4SEBMrLyzGpyqV/zs9EaOM4L5z5zfoujN2CcHF2JDw8HDc3t3Z9zrG8Y6w4soITxScY7jGcJUOX0Ne+r4GPRpKuPkNc8zAFnIHqxvUe1xIZHjcWIQRpaWkkJSVRXV2NaUkaoQWbGMAZPhJ3kOV2K762ary9+xAWFoaVlf7lp3ptPd+e+Jb1x9ZTo6lhrv9cHgt8TJaypOuKXGEuw+OGVF9fT2pqKsnJyahra7EqiOLT0hCSat0YaXWW0S512Dm4dqjhYkF1Ae/GvMvPGT/jbu3O4ojFjO89XpaypOuCDA8ZHje0uro6EhMTSUtLo65eQ1SROXcVvs8IJZFtJrdQ6HkLjo7OhIaG0qdPn3Z9RmxuLCuOriCtOI0RniNYMnQJfezaN5YkdRcyPGR4SEB1dTXx8fFkZGRATQle2b8wUn2IUqz5xXI69R5D8fTwIDw8vF239tZr6/km9Rvej3ufWk1tQykr6DEsTSw74WgkqfO194L5duA3YJsQIrUT59fpZHhITTVdaKiUncM/bxPB4jj/0M4m23UcwT3UDPDpT0hICObm5nqPX1BdwDvR77Dl1BY8rD1YPGQx47zGyVKWdM1pb3i4A5Mv/jcQOEpDmOwWQlS0+KZuSoaH1JKCggJiY2PJy83FpCiFT4qDOVljy2SLFMKctdi5eBMQEMDgwYMxMjLSe/zonGhWRq4kvTidkT1HsmTIEnrbta95oyRdDYa428oIGApMAcYD1TSsNL8mWpTI8JBac+7cOY4dO0ZxSSnxJWbMKXgbFensNh7JGfdbcXbr2e691P9ayno44GEeCXxElrKka4LBr3koiuIMTBJCfNnRyXUFGR7SlQghyMjIID4+nuqyQpyydzG+dhf1mLDNfDLlHjfTy6vh1l4HBwe9x8+vyuedmHfYemorntaeLB6ymLFeY2UpS+rW5AVzGR5SG2k0GlJSUkhJSUFdmk2/nC3cpI3mpfrHqXQJJ8JJzUCf/qhUKiwsLPQePzonmhVHV3Cy5CSjeo5iyZAleNm1b0tdSepsN2x4yN5WUnvV1taSkJBAeno69cVZbCgYxKlqc2abH2KgoxlW7r74+/vj5+en9/UQtVbN18e/5oP4D1Br1DwU8BCPBj6KhYn+YSRJnemGDY9G8sxDaq/y8nKOHTtGVtYZjpeZ8EDeKgYrmRxSQkh1nYaDe19CQkLw9vbWe+y8qjzWRK9h2+lt9LTpyeKIxYzxGiNLWVK30dH9PJ5v4eFSIEb2tpJuFPn5+cTExJCfcwGbnN+ZUP0r5qjZZTqWXI+JuHn2Jjw8HBcXF73HjsqJYsWRFWSUZjC612heinhJlrKkbqGj4fEVEA5sufjQbUAU4At8fy3ccSXDQzKUrKwsjh07RmlhDp7Z2xirPsAj9S9h4dyP4U619PPuTWhoKDY2NnqNq9aq+er4V3wQ9wH12noeCXyEhwMelqUs6arqaHhsB2Y0ru1QFMUG+AG4k4azDz8Dz9fgZHhIhqTVaklNTSUpKYmy4ly+L+jFiXIz/mb2My497DFyUzFw0ECCgoL03j8ktzKXf0T/g22ZDaWsJUOWcLPXzZ10JJLUuo6Gx3EgWAhRd/FrcyBOCDFYUZRjQogQg8/YwGR4SJ2hrq6O+Ph40tLSyCiDu3PfwpdMEhjAMee7sHT1wd/fv12LDI9mH2Xl0ZWcKj3Fzb1uZvGQxXjZylKW1LU6Gh6v0HCW8dPFh6YCPwP/AD4SQswx4Fw7hQwPqTOVlZURExPDubNnMM2J4pbKn3FWStlvNIxT7tNwcO2JSqXS+6K6WqPmy+Nf8kH8B2i0Gh4NfJSHAh6SpSypy7Q7PJSG2z56Aa407B6oAH8IIa6pn8QyPKSukJOTQ3R0NEX5OfTI3sPY2j3cqf47vZwcGOVcg6erE6GhoXpvQpVbmcua6DX8lvkbvWx6sWToEkb3Gt1JRyFJf+romUeMECKsU2bWRWR4SF0pPT2d+Ph4CktK+SXfgYRSc942+xhh40mt+1B6eTVcVLe3t9dr3CPZR1h1dBWnSk8xxmsMiyMW08u2VycdhSR1PDzeBzYIIaI6Y3JdQYaH1NXq6+tJSEjgxIkTnCurY0L2+6g4wSl6csjhLhTnwfTr1w+VSoWlZdv7XKk1av57/L/8K/5faIVWd1eWubH+3X8l6Uo6Gh4pwCAgE6ikoXQlhBBBBp5np5HhIV0tFRUVREdHc/bMWURBMmNKf6S3kkeMEkCc+71Y27swcOBAAgMDMTU1bfO4OZU5rIlew/bM7XjZevHSkJdkKUsyuI6GR4vboQkhsgwwty4hw0O62rKzsxuuhxQVYpF9kGHV+7mjdjm+TqaMda7GwcYcPz8/ve/MOpJ9hJVHV3K69DRjvcayeMhietr07MQjkW4khmjJPhIYIIT4TFEUF8BGCHHawPPsNDI8pO5ACEFaWhoJCQkUVdSwM8+a2GITvjVbQbaVL6XuN2Nj70BQUBD9+/dv87hqjZovUr7g3wn/Riu0PBb4GHMD5spSltRhHT3zeI2GFeaDhBADFUXxpGFl+QjDT7VzyPCQupOm60OKKqtRXfiKUcSQI5w4YD+dWucQHBwdCQ4O1msPkZzKHN6OepsdWTvwsvViyZAljOo1qhOPRLredTQ84oAQILZxQaCiKAnXwjUP2VVX6s5KSkqIiooiJyeXuqJTDC/aiJ9ymlT6ss/tMWxs7XF1dSUkJESvnlmHLhxi1dFVZJZlMs5rHIuGLJKlLKldOhoekUKIIYqixAohQhVFsQYOXwvh0UieeUjdWVZWFrGxsZSVl0NOHIMqDvNg3WICemiZ5FqGtZkJPXv21Ov23jpNHV+kfMFHCR8hhOCxoMeY6z8XM2P92qVIN7aOhsdCYAAwAVgFPAx8JYT4p6En2llkeEjdnUajISkpiePHj1Neq2FfngXxRYJtZi9x3FxFnscEFDMbvL29UalUWFtbt2nc7Ips3o5+m51ZO+lt25slQ5cwsufITj4a6XphiAvmE4CJNNymu10IsdOwU+xcMjyka0VFRQUxMTGcPXuWiqpqvHN+ZoLmdyqwYq/N7VS43YSRiRkDBgwgMDCwzbsZHjp/iFWRDaWsW3rfwqKIRXjYeHTy0UjXOrkZlAwP6RqTnZ1NVFQUZWVlVJbkEFKwiSEkcRY3fnB9ARc7a0xNTRk0aBABAQGYmJhccczGUta/4/8NwOPBj/Og34OylCVdVrvCQ1GUrUKI268w8BVf0x3I8JCuRUIIUlJSSEpKora2jtr8NOxLk1hU9wgBdmqmu+RiaWmDuXnDGhFfX1+MjY2vOG52RTaro1az+8xuvO28WTJkCTf1vKkLjki61rQ3PEqAA62NC/gLIfp1fIqdS4aHdC2rqakhJiaGzMxMajWCPwosSM+v5DezF4k1DeWsx21ozR2wsrIiICCAAQMGtGkr2z/O/8Gqo6s4U36GCX0msChiEe7W7l1wRNK1or3h0ZYdaOqEEIc7MrmuIMNDuh7k5eURFRVFcXEx1bU1OF3YzW31O9GisN9yAoUe4xDGFtja2hIcHEyfPn2uGCJ1mjo2JG/g44SPURSFeUHz+B+//8HUuO2tUqTrl7zmIcNDuk4IIThx4gQJCQnU1dVRUV6MT+5WbuEIecKBj1xfx9u+4fqHg4MDwcHB9Op15c67Fyou8FbUW7pS1tKhSxnuObyzD0fq5mR4yPCQrjM1NTXExsZy+vRptFpBRdFZagszWVs3jQE2auY4pWNi23A3lYuLCyqVqk37iPx+7nfejHyTM+VnmNhnIi9GvChLWTcwGR4yPKTrVH5+PpGRkRQXF1OvhaNF5uTlZ/Oj6SvEGgWR7j4dtVXDD38PDw9UKhVOTk6tjlmrqWVD0gY+TvwYI8WIJ4Kf4IHBD8hS1g2oo4sEXYUQeX95bJAQ4oQB59ipZHhI17PGUlZ8fDxqtZoadT0WFw4yvW4LltTyh/losj0mozG1RVEUevXqhUqluuJq9fMV51kduZq9Z/fS174vS4YskaWsG0xHw+ME8IoQ4ruLX78APCKE8DP4TDuJDA/pRlBbW0t0dDSZmZkIISivqqJnzm/cqtlLEfasdVnBYHuBooCiKPTt25egoCBsbGxaHffAuQO8GfkmZ8vPMsl7EgvDF8pS1g2io+HhAXwE1ABuwHHgBSFEhaEn2llkeEg3kqZ3ZQGUlBRwoSCfr2puorelmqccolE7DAJFwcjICB8fHwIDA1vd0bBWU8tnSZ/xSeInGClGPBn8JPcPvl+Wsq5zhmhP8hSwBNACs4UQBw07xc4lw0O60QghSE1NJSEhAbVajVZATLEZtflpfGj8NieUfiS63kmNbcMyLRMTEwYMGEBAQADm5pffB+Rc+TlWR65m37l99LPvx9KhSxnqMbSrDkvqYh0989gJZAPPAL2AT4EDQoiFhp5oZ5HhId2oqquriYmJISsrCyEEdfVatDmxTK3ehJtSQqRJOJnuU6mzcAbA1NSUwYMHM3jw4Fa3xd1/dj+rIldxvuI8U7yn8EL4C7hZX/luLuna0tHwuEMIsbnJ1ybAEiHE3w07zc4jw0O60eXk5BAVFUVpaSkA5TVqemTvZXr9b+TTg7WOywlxrMfo4prCtrQ8qamv0ZWyTIxMeDL4Seb4zcHUSJayrhfyVl0ZHpKEVqslJSWF5ORk1Go1AMXlZRzPreTXGj+8LKp50W4f5S4RCKUhMCwtLfH392fgwIGX3Vv9bPlZVkeuZv+5/fS378/SoUsZ4jGky45L6jwdPfMoBxpfZAaYAhVCiLbtStMNyPCQpD9VVlbq2r4LIRACkspMMco9xiqjf3EGD44530m5vR9cbG9iY2NDQEAA/fr1u2yI7Du7jzcj32woZfWdwsLwhbhauXbloUkGZtAzD0VR7gCGCCGWGmJyXUGGhyRd6vz580RHR1NeXg6AWiOozE3ltsof6atkk2g0mBPud1Jj9ecWtnZ2dgQGBuLt7d1i36ya+ho+TfqU/yT+BxMjE+ar5nPf4PtkKesaZfCylaIoR4QQwzo8sy4iw0OSWtZ0B8P6+noAqtRajC8c4a66n8gS7nzk+BJDHGsxaXLC4eDgQFBQEL169WoxRM6WnWVV5Cp+P/87Pj18WDp0KRHuEV11WJKBdLRsdVeTL42AcOBmIUS3X2qqKMpUYKqPj89j6enpV3s6ktRtlZWVERUVRXZ2tu6xkqpqonI17K/0or9ZMS/Z7qDA9WY0xn/eyuvs7ExwcDAeHpfuSiiEYN/ZfayOWs35ivPc2vdWFoYvxMXKpUuOSeq4jobHZ02+rAcygY//2rKkO5NnHpLUNpmZmcTGxlJVVQWAEJBeYYJZTiQvK/8hnx5E9ZhGkVMEKA2nIoqi4OrqSlBQUIvNF2vqa/gk8RM+S/oMU2NT5gc3lLJMjK68+6F0dcm7rWR4SFKbqdVq4uPjSUtLQ6vVAqARUJh3hltKfyDIKINTSm8SXO+kwnag7n2KouDu7n7Z5otnys6wKnIVf5z/A58ePiwbuoxw9xZ/LkndRHs3g/onf95ldQkhxDOGmV7nk+EhSforLi4mMjKS/Px83WPV9YLa7ATurPmReK0PPzrMY4RzDWZNroc0Nl8MCgrCwcGh2ZhCCPac3cNbkW9xofICt/e7nRfCX8DZ0rmrDkvSQ3vD439aG1QI8bkB5tYlZHhIUvsIITh58iTx8fHU1NToHi+prudgnsLhcmdCTM7wvM0OLrhPQm3yZ5NFRVHo3bs3wcHB2NnZNRu3ur5aV8oyNzbnKdVT3Ot7ryxldTPtDQ8TIUR9p86si8jwkKSOqa2tJTY2llOnTtH0Z0ZmpQlm2Yd5QXxOpWLJUdsp5LuOQqv8eWuukZERffv2JTAw8JIOvlllWaw6uoqDFw4ywGEAy4YuI8wtrMuOS2pde8MjVggRevHP/xRCLOjEOXYqGR6SZBi5ublER0frOvYCaAXkFOQzvHgTo4ziycaFGOe7KOkR1Oy9xsbG9O/fn4CAAKysrHSPCyHYc2YPq6NWk12ZzdR+U3k+/HlZyuoG2hsex4QQIRf/rAuSa5EMD0kyHK1Wy/Hjx0lKStK1OQGo00Jxzklur/iefVoVv/e4izEu1Vj8pTWWqampLkQsLCx0j1epqxpKWcmfYWFswdMhT3PPoHtkKesqMsSZhwwPSZKaqaysJDo6mnPnzjUrZZXXadmXa0ZkqQ23mMTzuOU+Mj2mUmPm2Oz9ZmZmDBw4ED8/P8zMzHSPZ5ZmsipyFYcuHGKgw0BeHvYyIa4hXXZc0p/aGx5VwElAAfpf/DMXvxZCiKAW39gNyfCQpM5z9uxZYmNjdW1OGl2oNsbkwmHma79CAY5ajyfb/RbURs03nTI3N8fX1xdfX19dG3ghBLvP7GZ11GpyKnOY1n8az4U9J0tZXay94dGntUGFEFkGmFuXkOEhSZ2rvr6ehIQETpw4gUaj0T0uBGSXlONfsIXblIMUY8chxxkUO17688jS0pLBgwczcOBATEwaSlVV6io+TvyYDckbsDS25KmQp2QpqwvJRYIyPCSpS5SUlBAdHU1OTk6zx+u1kJt3gXFlP7BJM4KT9qMY51KFTQv9Eq2trfHz88PHx0e3l8jp0tOsOrqKw9mHGeQwiGXDlslSVheQ4SHDQ5K6jBCCU6dOERcXR3V1dbPnqtSwN8+co8WWPGCyi7stIkl3n06FRc9LxrG1tSUgIABvb2+MjY0RQrAzaydvRb1FblUu0/tP57mw53CyvHQ1u2QYMjxkeEhSl6utreXYsWNkZGTw158z+bVGcCGSh+u/xYZqoi1HcsZtCrWml24TZG9vj7+/P97e3hgZGVGlruKjhI/4POVzLI0tWRC6gLsH3o2xUcs7Hkrt1+7wUBQlSAiRoChKoBAisdNm2MlkeEjS1dPS2pBG58tq6JO3nRliN2rFlD097qbEeWiL4zg6OhIQEICXlxeKonCq9BSrjq7iSPYRBjsOZunQpahcVZ19+M1OewAAFr9JREFUODeUjoTHWuB94CkhxLOdNL9OJ8NDkq4ujUZDamrqJWtDoGGR4dmCIoYVb+bz+lsosfVlsmsZdmZGus69TTk7OxMQEEDPng2lrh1ZO3gr6i3yqvK4w+cOng19VpayDKS9d1u9BjgBc4AvgQIhxPJOm2UnkuEhSd1DRUUFMTExl6wNAajVwIECCw4WWPCi8ddMMEsmxe0OSqx9LhlHURRcXFwIDAzE3d2d6vpq/p3wb75I/gJLU0ueCXmGWQNnyVJWB3XkzGMaMBn4TQjxcyfNr9PJ8JCk7uXMmTP/396dR1dZnwkc/z73Zr1JWEI2IECQRRI22aW4gVIFcYNiFe2pVGtnporLdNqp1OrUVnvUo+OcVlunderUqkXb6VGnrVAdQbFFQsgCYQ0SCFlJCJDkZr3P/JEbGpJ7Q26Sm/X5nJPDve/7e988L79z75Pf8v5eMjIyqKqqarevsl6oLczijoa3GC0VZIXO5cjoG6gOa/88dBEhMTGRWbNmkZCQwJHKIzy540l2FO8gNTaVjZduZHb87N64pEGpO8njR6q6UUSeUNVHgxZhkFnyMKb/aWhoIDs7m4MHD553b0iLoqomRhRvZZ2+R5g08sdhd3A2wfd4iIgwZswYZs2aRWxsLO8ffZ9ndj5DqbuUWybfwoPzHiQ2ItbnscY/m21lycOYfquiooL09HRKS9s/nFQV8k9Vk3byj7zYsJKQ6ARuji9jWGQYHml/o6DD4SA5OZkZM2YQHhPOz7N+zq9zf40r1MWGORv40tQvWVdWACx5WPIwpl/zeDzk5eW1e25IiwYP/LU8nK1lkTzn/A/mhhwlO/5mymJmgki78k6nk+TkZGbOnEk55Ty540k+K/6M1NhUvnfp95gVP2BWV+pTljwseRgzILjdbnbt2kV+fn67AXWAqkahovAga2s3McVRyEHnVPYn3cKZyPE+zxcSEsKECROYPn0628u38+zOZyl1l7JmyhoemPsAIyNG+jzONLPkYcnDmAGlqKiI9PR0Tp8+7XN/iVtxFn7GnZ4/ECdneCtmPbUJ8301QoDmZeBTUlJImZrCa3mv8Vrua7hCXTww9wHWTFljXVl+dHWqrhO4B0imebbV9lb7vqeqPwxGsMFgycOYgaexsZHc3Fz27dvX7t4Q8I6HnGlgdOmH/LRuJcOjIlkXd4yYKBeNjggfZ/z7s0QikiN4Lus5dhbvZPqo6WxctJGZ8TODfUkDTleTxy8AF/AZ8BVgq6o+7N03oJ7vYcnDmIHr9OnT5xZb9PV91aSQXhHOB6XhvO58nAmOk2TGrqJw5KWoj5sMoXkZ+MmTJ3PMdYwXsl7gpPskq6estq6sNrqaPLJbntkhIiHAi0AccDvwt5anDA4EljyMGdhUlc8//5zMzExqamp8lnE3CQXFJ7ipehNzHYc47khmb8LNnIxO9XveiIgIxl40lq31W3nz4JtEh0WzYc4G68ry6mry2K+q09ps+z5wLZCgqlN6PNIgseRhzOBQV1dHZmYmeXl5eDwen2XKa4Waoj3c0fAW4x1lvBr9DSRxFg4/4yEALpcL13gXb5a9ya7SXcwYNYONl25kRtyMIF3JwNDV5PEa8Jqq/rnN9nuAl1TVx0r8/ZMlD2MGl9LSUtLT06moqPBb5niVElr8GS/XXk18pHB37F5iYkZSFxLj9xiXy0V5Qjmvn3id8tpy1kxdwwNzHmBExIhgXEa/Z7OtLHkYM+h0tNhiC49CVmUYH5SE847jnxnlqCJzxLXkx16Fx+H/79/Q6FAyXZn8qfhPRIdF8+DcB1k9ZTUOP2Mog1VXWx7fVtWnva/XqupbrfY9qaqPBCXaILDkYczgdfbsWXbt2sWJEyd8DqgD1Hsgr6SCq8/8nmXO3ZTJKHLibqRw2DyfNxm2cEe5+Qt/IfdMLjPjZrJx0Uamx00P1qX0O11NHudmVLWdXWWzrYwx/Ymqkp+fT2Zmps/FFlucaRBKCj/n1rpNTHfk85Poh4hJnIyzg/EQgKOuo7xf+z6nG06zdupaNszdwPDw9g+uGmy6mjx2t8yoav3a1/v+zpKHMUNDXV0d2dnZHD582Odiiy2KaoTK4gP8pmYh8eEeHh6xnYgRo6kOjfN/bupID0vn09pPGRY2jIfmPcTNk28e1F1ZHSWPjq5a/bz29T5oROQiEfmliLzd0TZjjAkPD2fBggUsX76cuDj/iWC0S5k2cSrrxlcTovVcU/k61+c/wcXFfyC0yfdU4HDCWVK/hLscdzHSM5LHPn2MO/94J7nlucG6nH6to5ZHE1ANCBAJtPyPChDRmdlWIvIKsAooVdUZrbZfB7wAOIFfqOqPO3Gut1X1Sxfa5ou1PIwZepqamjh06BA5OTnU1dX5Ldfogf0na1h06h1ucnxClUSRFXs9x0YuQcX3vR6qSi65fMRH1HhqWHvxWu6fc/+g68rqs9lWInIFUAX8d0vy8C57chBYDhQAO2m+8dAJPNXmFF9T1VLvcZY8jDEBq6qqIiMjg+PHj/sdUAeoaRQ+Ly7hxupNLHbu44no75OckEhoB/0ztVrLdrazy7OLYaHDeHj+w9w8ZfB0ZfXpVF0RSQHea5U8FgOPq+q13vffBVDVtomj7XkCSh4ici9wL8D48ePn5efnd/NKjDEDlapSUFBARkYGZ8+e7bBsWa2QV1zKe1VTGRHaxBPD3kVGTqIyItnvMaVayhbdQoEWkDYijccue4y0UWk9fRm9rqtjHsEyFjje6n2Bd5tPIjJKRH4GzGlJNL62taWqL6vqfFWdHx8f34PhG2MGGhFh3LhxrFixgrS0NEJC2j9IqkV8hHJpSjzrU84yylnN4qo/s7Lgx8w88RsiGyt9HpMgCayTdayUleRX5nPbe7fx/a3f53Sd71WBB4O+aHmsBa5V1Xu8778CLFTV+4MVg3VbGWNaO3nyJBkZGZSVlXXYleVR2FfewMUV77NONuMRJ1nDv8iRUctocoT7PKZWa/lEPyFDM4h2RnPfzPu4bdZtA7Irq7+1PAqAca3eJwOFfRCHMWaIiouL4+qrr2b+/PlERPhevh3AITA9LpSmSavYGP1D/q/pEmZWbmFHSRN1fmYCR0gE1ziu4auOrzKsaRhPZT7F6k2rST82uP6A7YuWRwjNA+ZXAydoHjBfp6p7gxWDtTyMMf5UV1efG1D3t9hii8p6Ib24jo/OJBId0sRLUa9QM2oWZa6pPsurKnt0Dx/pR7hxs2zkMr695NuMGTUmGJfS4/pyttUbwFU0L+VeAjymqr8UkZXAv9M8w+oVVf1R0ILAkocxpmOqyokTJ9i9e7ffpxe2drzGyY7Cen6iPyBZTrI/fDYHE2/iTFiiz/K1WsvH+jG7dTcuXKxLXsf6ReuJifa/SGN/MGQXRhSRG4AbJk+e/PVDhw71dTjGmH6uvr6evXv3cvDgQb+LLbZQhf2nlYTSbXxN3iVS6smJvpyD8atocEb6PKZES9js2UwhhYxzjOMbk7/B8tnLcblcwbicbhuyyaOFtTyMMYEoLy8nIyOD0tLSDgfUARo8sKeslgWV/8sVjiweifohX0hQXCG+j1NVcjSHj/QjaqllUdgi7km9h1kXzyIy0nfS6SuWPCx5GGMC1NTURF5eHjk5Objd7guWr2oUPil2sL0yhhhnA29EPktp7CJOxMzxuXKvW93nurKiiOI613XcOv1WJk2aRHi475lcvc2ShyUPY0wXVVdXs3v3bo4dO3bBAXWAkloHGUXVPNrwAtMcxzkaOonchNVURKb4LF+sxWz2bKaIIpJJ5paYW7hy+pVMnDiRsLCwHr6awFjysORhjOmGlgH1zMxMKit93yjY1qEzDsJKdnIvvyNeTpMbuYCc0V+mwdG+a0pVydZstupWaqllnszj+hHXc0naJUyYMIHQ0L55cKslD0sexpgeEMiAOkCTQk65Mrl8C/NkH09FfodrkuoYFtIEPm4abNuVtdSxlMtiLyM1NZVx48b1ehIZssnDZlsZY4IhkAF1AHeTsK00jE8rIhlJFe9EPs7nsVdydPgXfK7cW6RFbPFsoYgixjGO5c7lpMalkpqaytixYztcXqUnDdnk0cJaHsaYnhbogDpAeZ2D7OKz3Fv3Cosc+ylyjmFvws0UR7V/tK2qkqVZbNNt57qyLndcTnJCMqmpqYwZMwan0/eS8T3FkoclD2NMkFRVVZGVlcWxY8c6fHpha0ernNQU7+VezyZSHCXkhaexY+y9eBztu6Xc6mabbiNTM4kiimWyjOnO6SQmJpKWlkZSUhIOR3BWmrLkYcnDGBNELQPqWVlZnDp1qlPHeBT2nHIQV/YJKVrAq671XJvkJj6kliZH+1lW7bqyHMtJCkkiMTGR1NRUEhMTezyJWPKw5GGM6QV1dXXk5uZ2ekAdoN4D209GsK0sgglSyB/CHmffyKs5PHJpuyTiUc+5WVl11DFf5rNElhAVGkVSUhJpaWnEx8cjPu4r6QpLHpY8jDG9KNABdYAzDUJWUTWr3b/lWmc6FRLLnvgbOBYzv93MLLe62apbydIsoolmqSwlVVIJCwtj/PjxLF68uEeuw5KHJQ9jTC9rbGwkLy+PPXv2dHpAHaDQ7aSwKJ+7G99gpuMo+aGT+Hjcg+CjS6pQC9ni2UIxxee6ssa7xrNmzZoeuYYhmzxsqq4xpq+dPXuWrKwsjh8/3ukBdVU4cMaJlOwmuuk0WyJXsCKphpSQCtwhI84r61HPuVlZ9dSzKHQRz699nqjQqG7HPmSTRwtreRhj+lLLM9SzsrI6fYc6QKMHPqsI58PSCOaxl1fDnmbvsCs5MGoF9c7zV+Kt0Rq26layNZuEyAS+teBbXJdyXbfGPyx5WPIwxvQDtbW15ObmcvjwYerr6zt9XE2jsLuklqvOvsNa51ZqHC72xq4gb8TleOT8GwbLw8r51PUpueW5LExayCOLHmHSiEldire/PYbWGGOGpIiICObOncvSpUtJTEzsdKvAFaIsGRtO4cQv883QH7C7cSKLyt9mSf4LqOf8BkCyM5nXV77Oo5c+yv6K/Ty98+lgXIq1PIwxpi80NjZy+PBhcnNzqampCejYvLNOKksO0VhfR2b4Qq5PqmJGSAGV4WOJjIw8N2BeUVuBu9HN2OixXYqxo5ZH7yyQYowx5jwhISFMmzaNMWPGkJ2dHdCA+qSYJjzRF5FxKoxTpQ7Kj2WyKvRl9kct4si4tefKxUbEBit8a3kYY0xf83g8HD9+nJycnIAG1AHqmmBnmYc5lZtZ7/wTKg5kyQbCrngIwqO7FdeQHTC3qbrGmIHE7XafG1Dv7B3qLSrrhcziam5y/44bnH+jZNRC4r65GafDZlt1mbU8jDEDSWlpKZmZmZSVlXX6DvVzxza5KK08xeHiSqoSF/DoqjSWTI7rUhw228oYYwaQhIQEli5dypw5c4iMbP/kwY5MiFae2vA17rr9dqrqGvnNjvygxGgD5sYY0w+FhoaSlpZ2bkC9oKCgU89QBxARVs0awzWpidTUd24QPlDW8jDGmH5sxIgRXHbZZSxevJhhw4YFdGxEqJPYqPbLu/cEa3kYY0w/53A4mDhxIomJiezdu5cjR44EPKDe4zH16W83xhjTaS6Xi/nz53PFFVeQkJDQY8/t6ApreRhjzAAiIowePZrY2FgOHTrEgQMHAlryvadYy8MYYwag8PBwZsyYwdKlS0lOTg7ac8z9seRhjDEDWGxsLJdffjkLFy4kJiam137voO62anWHeV+HYowxQeN0Opk8eTJJSUnk5eX1yu8c1C0PVX1XVe8dPnx4X4dijDFBFx0dzezZs3vldw3q5GGMMSY4LHkYY4wJmCUPY4wxAbPkYYwxJmCWPIwxxgTMkocxxpiAWfIwxhgTMEsexhhjAmbJwxhjTMCGxPIkwBkROdRq13DgtI9D/G2PA072fITd4i/WvjxnIMd3tuyFynV1v6/tQ6Weu3veQI/tjboOdJ/VdedM8LtHVYfcD/BygNvT+zrmzsbal+cM5PjOlr1Qua7u97V9qNRzd88b6LG9UdeB7rO67v7PUO22ejfA7f1RMGLt7jkDOb6zZS9Urqv7B0pdByvO7pw30GN7o667uq8/6Y917Zd4M5PpgIikq+r8vo7DBJfV89Bhdd19Q7XlEaiX+zoA0yusnocOq+tuspaHMcaYgFnLwxhjTMAseRhjjAmYJQ9jjDEBs+TRDSJykYj8UkTe7utYTM8TkSgReVVE/lNE7ujreEzw2Gc5cEM2eYjIKyJSKiJ72my/TkQOiMhhEfnXjs6hqkdU9e7gRmp6UoD1vhp4W1W/DtzY68Gabgmkru2zHLghmzyAXwHXtd4gIk7gp8AKIA24XUTSRGSmiLzX5ieh90M2PeBXdLLegWTguLdYUy/GaHrGr+h8XZsADeq1rTqiqttEJKXN5oXAYVU9AiAibwI3qepTwKrejdAEQyD1DhTQnEAyGdp/aA1IAdZ1bu9GN/DZB+J8Y/n7X5rQ/OUx1l9hERklIj8D5ojId4MdnAkaf/X+e2CNiLzEwFniwnTMZ13bZzlwQ7bl4Yf42Ob3LkpVLQf+IXjhmF7is95VtRpY39vBmKDyV9f2WQ6QtTzOVwCMa/U+GSjso1hM77F6HzqsrnuIJY/z7QSmiMhEEQkDbgPe6eOYTPBZvQ8dVtc9ZMgmDxF5A/grcLGIFIjI3araCNwHvA/sAzap6t6+jNP0LKv3ocPqOrhsYURjjDEBG7ItD2OMMV1nycMYY0zALHkYY4wJmCUPY4wxAbPkYYwxJmCWPIwxxgTMkocxgIg0iUhmq58Ol+PvLa3iGiMiO7yvj4lIWatYU9occ5WI/LXNthARKRGR0SLyjIgUi8i3evNazOBia1sZ08ytqpf05AlFJMR7U1p3tI5rkfe8dwHzVfU+P8dsA5JFJEVVj3q3XQPsUdUi4F9EpLqbcZkhzloexnRARI6KyL+JSIaI5IjINO/2KO/DhnaKyG4Rucm7/S4ReUtE3gU2i4hLRDaJSLaI/NbbepgvIneLyPOtfs/XReS5LsQ3SUT+LCK7RORjEZmmqh7gLeDLrYreBrzRrf8MY1qx5GFMs8g23Vatv3hPqupc4CWgpatnI/Chqi4AlgLPiEiUd99i4Kuqugz4J+CUqs4CngDmecu8CdwoIqHe9+uB/+pC3C8D96vqPG9sL3q3v0FzwkBEwoGVwO+6cH5jfLJuK2OaddRt9Xvvv7tofjQtwBdp/vJvSSYRwHjv6y2qWuF9fRnwAoCq7hGRbO/rahH5EFglIvuAUFXNCSRgEYkGvgC8JXJupfFw7/l3iki0iFwMpAJ/U9VTgZzfmI5Y8jDmwuq8/zbx98+MAGtU9UDrgiKyCGg9nuDr+REtfgE8Auyna60OB1DZQdJ7k+bWRyrWZWV6mHVbGdM17wP3i/dPfhGZ46fcJ8Ct3jJpwMyWHaq6g+ZnS6yjC1/uqnoG+FxE1nrPLyIyu1WRN4A7gWXYsuOmh1nyMKZZ2zGPH1+g/BNAKJAtInu87315EYj3dld9B8gGTrfavwnY3o0upTuAu0UkC9hL8/O4AVDVXKCG5rEZm11lepQtyW5MEImIk+bxjFoRmQR8AExV1Xrv/veA51X1Az/HV6lqdBDiehyoUtVne/rcZmiwlocxweUCPvG2DP4H+EdVrReRESJykOaBep+Jw+tMy02CPRWQiDxDc3eWtUZMl1nLwxhjTMCs5WGMMSZgljyMMcYEzJKHMcaYgFnyMMYYEzBLHsYYYwJmycMYY0zA/h+/pGccHzGr8QAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plot_kwargs = {\n",
" \"energy_range\": [0.1, 30] * u.TeV,\n",
" \"energy_power\": 2,\n",
" \"flux_unit\": \"erg-1 cm-2 s-1\",\n",
"}\n",
"\n",
"# plot stacked model\n",
"model_best_stacked.spectral_model.plot(\n",
" **plot_kwargs, label=\"Stacked analysis result\"\n",
")\n",
"model_best_stacked.spectral_model.plot_error(**plot_kwargs)\n",
"\n",
"# plot joint model\n",
"model_best_joint.spectral_model.plot(\n",
" **plot_kwargs, label=\"Joint analysis result\", ls=\"--\"\n",
")\n",
"model_best_joint.spectral_model.plot_error(**plot_kwargs)\n",
"\n",
"create_crab_spectral_model(\"hess_pl\").plot(\n",
" **plot_kwargs, label=\"Crab reference\"\n",
")\n",
"plt.legend()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exercises\n",
"\n",
"Now you have learned the basics of a spectral analysis with Gammapy. To practice you can continue with the following exercises:\n",
"\n",
"- Fit a different spectral model to the data.\n",
" You could try `~gammapy.modeling.models.ExpCutoffPowerLawSpectralModel` or `~gammapy.modeling.models.LogParabolaSpectralModel`.\n",
"- Compute flux points for the stacked dataset.\n",
"- Create a `~gammapy.estimators.FluxPointsDataset` with the flux points you have computed for the stacked dataset and fit the flux points again with obe of the spectral models. How does the result compare to the best fit model, that was directly fitted to the counts data?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## What next?\n",
"\n",
"The methods shown in this tutorial is valid for point-like or midly extended sources where we can assume that the IRF taken at the region center is valid over the whole region. If one wants to extract the 1D spectrum of a large source and properly average the response over the extraction region, one has to use a different approach explained in [the extended source spectral analysis tutorial](extended_source_spectral_analysis.ipynb)."
]
},
{
"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
}