{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Results\n",
"\n",
"This notebook shows how to work with the results from the joint Crab paper.\n",
"\n",
"We will give a quick overview of the files in the `results` folder, and then show how to load the spectral fit results (best-fit parameters and covariance matrix) and make a plot of the parameter constraints.\n",
"\n",
"As explained in the [README.md](README.md), the main analysis and the figures for this paper were produced by Python scripts (`make.py` and `joint_crab` package). For further examples, please look there."
]
},
{
"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": [],
"source": [
"import yaml\n",
"import numpy as np\n",
"import pandas as pd\n",
"from multinorm import MultiNorm\n",
"from astropy.table import Table"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Overview\n",
"\n",
"All results are stored in standard machine-readable formats that are easy to read from Python. We used the YAML format for structured data (using the Python `yaml` package), ECSV and FITS for table (using the Python `astropy` package)."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"README.md\n",
"\u001b[1m\u001b[34merrorbands\u001b[m\u001b[m\n",
"\u001b[1m\u001b[34mfigures\u001b[m\u001b[m\n",
"\u001b[1m\u001b[34mfit\u001b[m\u001b[m\n",
"\u001b[1m\u001b[34mmaps\u001b[m\u001b[m\n",
"provenance.yaml\n",
"\u001b[1m\u001b[34mspectra\u001b[m\u001b[m\n",
"\u001b[1m\u001b[34msummary\u001b[m\u001b[m\n"
]
}
],
"source": [
"!ls -1 results"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"contours_fact.yaml\n",
"contours_fermi.yaml\n",
"contours_hess.yaml\n",
"contours_joint.yaml\n",
"contours_joint_systematics.yaml\n",
"contours_magic.yaml\n",
"contours_veritas.yaml\n",
"fit_fact.yaml\n",
"fit_fermi.yaml\n",
"fit_hess.yaml\n",
"fit_joint.yaml\n",
"fit_joint_systematics.yaml\n",
"fit_magic.yaml\n",
"fit_veritas.yaml\n"
]
}
],
"source": [
"!ls -1 results/fit"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"samples_fact.fits.gz\n",
"samples_fermi.fits.gz\n",
"samples_hess.fits.gz\n",
"samples_joint.fits.gz\n",
"samples_magic.fits.gz\n",
"samples_veritas.fits.gz\n",
"sed_fact.ecsv\n",
"sed_fermi.ecsv\n",
"sed_hess.ecsv\n",
"sed_joint.ecsv\n",
"sed_magic.ecsv\n",
"sed_veritas.ecsv\n"
]
}
],
"source": [
"!ls -1 results/errorbands"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"data.md\n",
"data.tex\n",
"results.md\n",
"results.tex\n"
]
}
],
"source": [
"!ls -1 results/summary/"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Read YAML fit results\n",
"\n",
"The fit results are given in the files `results/fit/fit_{name}.yaml`, where `name = {\"fermi\", \"magic\", \"veritas\", \"fact\", \"hess\", \"joint\"}`. As an example, you can look at [results/fit/fit_fermi.yaml](results/fit/fit_fermi.yaml) to see the structure. For each fit parameter, there's the name, value and error, and there's also the covariance matrix."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"{'parameters': [{'name': 'amplitude',\n",
" 'value': 3.8491280706505666e-11,\n",
" 'factor': 3.8491280706505666,\n",
" 'scale': '1e-11',\n",
" 'unit': 'cm-2 s-1 TeV-1',\n",
" 'min': nan,\n",
" 'max': nan,\n",
" 'frozen': False,\n",
" 'error': 1.1435170001020147e-12},\n",
" {'name': 'reference',\n",
" 'value': 1.0,\n",
" 'factor': 1.0,\n",
" 'scale': 1.0,\n",
" 'unit': 'TeV',\n",
" 'min': nan,\n",
" 'max': nan,\n",
" 'frozen': True,\n",
" 'error': 0.0},\n",
" {'name': 'alpha',\n",
" 'value': 2.5053735278333207,\n",
" 'factor': 2.5053735278333207,\n",
" 'scale': 1.0,\n",
" 'unit': '',\n",
" 'min': nan,\n",
" 'max': nan,\n",
" 'frozen': False,\n",
" 'error': 0.02549508100457278},\n",
" {'name': 'beta',\n",
" 'value': 0.2362106123242955,\n",
" 'factor': 2.3621061232429548,\n",
" 'scale': 0.1,\n",
" 'unit': '',\n",
" 'min': nan,\n",
" 'max': nan,\n",
" 'frozen': False,\n",
" 'error': 0.02450886256203386}],\n",
" 'covariance': [[1.3076311295223112e-24,\n",
" 0.0,\n",
" 1.2371922299026042e-16,\n",
" 1.1603842009321198e-14],\n",
" [0.0, 0.0, 0.0, 0.0],\n",
" [1.2371922299026042e-16, 0.0, 0.0006499991554297277, 0.000492142818491591],\n",
" [1.1603842009321198e-14, 0.0, 0.000492142818491591, 0.0006006843440846649]],\n",
" 'statname': 'wstat',\n",
" 'statval': 41.12925420683606,\n",
" 'fit_range': {'min': 891250938.1337459,\n",
" 'max': 28183829312.64455,\n",
" 'unit': 'keV'}}"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"info = yaml.load(open(\"results/fit/fit_joint.yaml\"))\n",
"info"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Usually you'll have to write a few lines of code to get the data in the form you want.\n",
"\n",
"For example, here's how to put the parameter information into a `pandas.DataFrame`:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
name
\n",
"
value
\n",
"
error
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
amplitude
\n",
"
3.849128e-11
\n",
"
1.143517e-12
\n",
"
\n",
"
\n",
"
1
\n",
"
reference
\n",
"
1.000000e+00
\n",
"
0.000000e+00
\n",
"
\n",
"
\n",
"
2
\n",
"
alpha
\n",
"
2.505374e+00
\n",
"
2.549508e-02
\n",
"
\n",
"
\n",
"
3
\n",
"
beta
\n",
"
2.362106e-01
\n",
"
2.450886e-02
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" name value error\n",
"0 amplitude 3.849128e-11 1.143517e-12\n",
"1 reference 1.000000e+00 0.000000e+00\n",
"2 alpha 2.505374e+00 2.549508e-02\n",
"3 beta 2.362106e-01 2.450886e-02"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pd.DataFrame(\n",
" [\n",
" {\"name\": _[\"name\"], \"value\": _[\"value\"], \"error\": _[\"error\"]}\n",
" for _ in info[\"parameters\"]\n",
" ],\n",
" columns=[\"name\", \"value\", \"error\"],\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The `yaml.load` call returns Python objects (str, float, list, dict). Often for computations you want a Numpy array. It's easy to make:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[1.30763113e-24, 0.00000000e+00, 1.23719223e-16, 1.16038420e-14],\n",
" [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00],\n",
" [1.23719223e-16, 0.00000000e+00, 6.49999155e-04, 4.92142818e-04],\n",
" [1.16038420e-14, 0.00000000e+00, 4.92142818e-04, 6.00684344e-04]])"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"cov = np.array(info[\"covariance\"])\n",
"cov"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Check parameter correlations\n",
"\n",
"If you look at Table 2 or Figure 4 in the paper, you will find that the parameter errors in the joint fit are much smaller than the parameter errors from each individual instrument dataset.\n",
"\n",
"Why is that?\n",
"\n",
"The reason is that there are strong correlations between the three fit parameters, especially for the Fermi-LAT dataset. The likelihood for each dataset is approximately a three-dimensional multivariate normal, and the joint likelihood is approximately the product of those likelihood, and the parameter correlations lead to a well-constrained joint likelihood.\n",
"\n",
"Let's load all fit results and compute the product multivariate normal, and then create a Figure similar to Fig. 4 in the paper. We will use the `MultiNorm` class from the [multinorm](https://multinorm.readthedocs.io/) module to help a bit with the covariance matrix related computations.\n",
"\n",
"### Load all results\n",
"\n",
"First let's load all fit results, and store each one in a `MultiNorm` object:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"def read_fit_result(name):\n",
" data = yaml.load(open(f\"results/fit/fit_{name}.yaml\"))\n",
" mean = [_[\"value\"] for _ in data[\"parameters\"]]\n",
" cov = data[\"covariance\"]\n",
" # Use parameter names that are consistent with the paper\n",
" names = [\"phi\", \"reference\", \"gamma\", \"beta\"]\n",
" mn = MultiNorm(mean, cov, names)\n",
" mn = mn.drop([\"reference\"])\n",
" return mn\n",
"\n",
"\n",
"results = {\n",
" name: read_fit_result(name)\n",
" for name in (\"fermi\", \"magic\", \"veritas\", \"fact\", \"hess\", \"joint\")\n",
"}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Check results\n",
"\n",
"We can now do any covariance-matrix related computations, e.g. to compute parameter errors or correlations:"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
mean
\n",
"
err
\n",
"
\n",
"
\n",
"
name
\n",
"
\n",
"
\n",
"
\n",
" \n",
" \n",
"
\n",
"
phi
\n",
"
3.849128e-11
\n",
"
1.143517e-12
\n",
"
\n",
"
\n",
"
gamma
\n",
"
2.505374e+00
\n",
"
2.549508e-02
\n",
"
\n",
"
\n",
"
beta
\n",
"
2.362106e-01
\n",
"
2.450886e-02
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" mean err\n",
"name \n",
"phi 3.849128e-11 1.143517e-12\n",
"gamma 2.505374e+00 2.549508e-02\n",
"beta 2.362106e-01 2.450886e-02"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# These are the numbers in the last row in Table 2 in the paper\n",
"results[\"joint\"].parameters"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"def plot_ellipses_panel(results, names, ax, legend):\n",
" colors = [\n",
" \"#21ABCD\",\n",
" \"#FF9933\",\n",
" \"#893F45\",\n",
" \"#3EB489\",\n",
" \"#002E63\",\n",
" \"crimson\",\n",
" \"black\",\n",
" ]\n",
" for name, color in zip(results, colors):\n",
" m = results[name].marginal(names)\n",
" ellipse = m.to_matplotlib_ellipse(ec=color, fc=\"none\", label=name)\n",
" ax.add_patch(ellipse)\n",
" ax.plot(m.mean[0], m.mean[1], \".\", color=color)\n",
" ax.set_xlabel(names[0])\n",
" ax.set_ylabel(names[1])\n",
" if legend:\n",
" ax.legend(loc=1)\n",
"\n",
"\n",
"def plot_ellipses(results):\n",
" fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(10, 3))\n",
" plot_ellipses_panel(results, [\"phi\", \"gamma\"], axes[0], legend=False)\n",
" plot_ellipses_panel(results, [\"phi\", \"beta\"], axes[1], legend=False)\n",
" plot_ellipses_panel(results, [\"gamma\", \"beta\"], axes[2], legend=True)\n",
" fig.tight_layout()\n",
"\n",
"\n",
"plot_ellipses(results)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As mentioned at the start, we can see that the joint fit result is approximately given by the product of the multivariate normal approximation for the fit result of each individual dataset in this analysis."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## SED\n",
"\n",
"In `results/errorbands`, the SEDs corresponding to the best-fit model are pre-computed, including 1 sigma error bands. Here we show how to upen up such a table and plot the SED."
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"sed = Table.read(\"results/errorbands/sed_veritas.ecsv\")"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"
\n",
" name dtype unit format\n",
"------------- ------- --------------- ------\n",
" energy float64 TeV 3g\n",
" dnde_mean float64 1 / (cm2 s TeV) 3g\n",
" dnde_median float64 1 / (cm2 s TeV) 3g\n",
" dnde_lo float64 1 / (cm2 s TeV) 3g\n",
" dnde_hi float64 1 / (cm2 s TeV) 3g\n",
" dnde_fit float64 1 / (cm2 s TeV) 3g\n",
" e2dnde_fit float64 erg / (cm2 s) 3g\n",
" e2dnde_mean float64 erg / (cm2 s) 3g\n",
"e2dnde_median float64 erg / (cm2 s) 3g\n",
" e2dnde_lo float64 erg / (cm2 s) 3g\n",
" e2dnde_hi float64 erg / (cm2 s) 3g\n"
]
}
],
"source": [
"sed.info()"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Text(0,0.5,'E^2 (dN/dE) (erg cm-2 s-1)')"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEOCAYAAACjJpHCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3Xd8HNXV8PHfmdmmLsumBUOMgQAOphoTSugkJsH0XkKAYMgTEp5ACgSeh4Q08qQBCYQSUxI6ptnGtGAbAqGZZloIHQwvYFvNKtvP+8fsmrWQVitpZ4t0vp/P4t3ZnZnjQd6jO/fec0VVMcYYY4bLKXcAxhhjqpslEmOMMSNiicQYY8yIWCIxxhgzIpZIjDHGjIglEmOMMSNiicQYY8yIWCIxxhgzIpZIjDHGjIglEmOMMSMSKHcApTBhwgSdNGlSucMwxpiq8swzz6xQ1bUG+9yoTiQiMhOYuckmm7BkyZJyh2OMMVVFRN4t5HOj+taWqs5T1VlNTU3lDsUYY0atUZ1IjDHG+M8SiTHGmBGxRGKMMWZELJEYY4wZEUskxhhjRsQSiTHGmBEZ1fNIiiGeTBNLptbYJiI4Ao4I4P2Zfe04Uo4wjTGmbCyRDCKZTtMTTw3+wRzCmslGZM3X2W2uI6uTkIglIGNMdbJE4gMFVJW0Zl8NLpt81tgmn33fkTX/FD5NVo5IJjlZYjLGlI4lkgqRTT6f2Zh/w4BEIOA4BFwh4MgaySabfOw2nDGmGCyRjFKqkEilSeS5K5dt5biO9wg4QsAVgo5jScYYUzBLJGPY6ltwKf1Mwsm2aLIJJjfZ2G0zY0wuSySmX/laNLlJZnVLJvPckowxY0/FJxIRmQycCzSp6mGZbXXAZUAcWKyqN5QxxDFnoCQjgNOnBeM6gitCwLUpS8aMVr7+6xaRq0XkExF5qc/2GSLymoi8ISJn5zuGqr6lqif32XwIMEdVTwEOKHLYZpgUSKWVWNIbMr0qmqS9J8HK7jgfd0ZZ0RWjvSdOVyxJLJn67OACY0xV8rtFci3wZ+Bv2Q0i4gKXAvsCy4CnRWQu4AK/7rP/Sar6ST/HnQi8mHk+tEkepmxSaV2daLJcRwi6DiHXIehay8WYauRrIlHVR0RkUp/N04E3VPUtABG5GThQVX8N7F/goZfhJZPnsTIvVc1LLimimftkIhDMDFsOuk5mJJn9LzamkpXjX+j6wPs5r5dltvVLRMaLyOXAtiJyTmbzHcChIvIXYN4A+80SkSUismT58uVFCt34TRXiKe/WWEevd1vsk84ord1xOqMJeuMp4sm03RYzpoKUo7O9v2E9A34rqOpK4LQ+27qBE/OdRFWvBK4EmDZtmn3r+CWZRNpacVpbcVauwGlrQzo7kVWdyKounK5VEIsi0RgS7YV4HEmlIJWCdNrLHIEAuC7quhAKoZEaNBLxHvX1aGMT2thAqrGJZEsL6ZYJpFtacMY1EwwEvEmXNv/FmLIpRyJZBmyQ83oi8GEZ4jCDSadxPliG+847uO+/h7vsfZz338f98AOcjz/C+egjnBXLkTK1DjQQIL3OOqTXWZfUuusRXe9z6IYbIhtNwpk0CXeTjQmsszbi2K0xY/xUjkTyNLCpiGwEfAAcBRzjx4lEZCYwc5NNNvHj8KOGdHbi/uffBF57jcC/XyXw+n9w33oT9523kVgs774qQrplPOmWFtLjx6PjWkg3NaH1DWhDg9eiqKlBw14Lg1AIXAd1XHBd7yCpFJJKQirtnS8WRXp7kWgUWbUKWbUKp7MD6ezEaV2JtLbitK7E6ezE/eAD3A8+IDhAfOmGBlIbb0J6441h8y1gyhScL04hsPlmEA4X90IaM0aJn/eaReQmYA9gAvAxcL6qzhaRrwEX4Y3UulpVf+lbEHi3tpYsWTKsfXviSVZFk0WOqExUcT78gOBzzxJ4cSmBl5YSXLoU9713B9wltc66pCZNIrXh50lP3IDUxA1Irb8+6XXWJb3uuqTXWhuCA32N+0Oy/4lGcT/+GPfjj3A/+n/IBx/gvPcezvvv4r77Hu47b+F0dPR7DHVdUl/YjPRWW6PbbINsuzXutB1wx7eU8q9iTEUTkWdUddqgnxsLnZZjNpF0dxN8ZgmhJ58g+MzTBJ59Bvfjjz7zMQ2FSH5hM1KbbU5y8y1IbvoFUptsSmqjyWh9fcnCzU5odEVw3cyfTqbkfqb8/lD6QDSdJrV8BenXX0f/8x/4979xXn3Ve7z9FpJOf2af5MabkNp+Gjp9OrLLzrjbbEMgYi0XMzZZImGNW1unvP7668M6RjUlEmltJfT4YwQffYTQY48ReGmp17GdI93UTGLb7UhutTXJrbYiseVWpDbZtOStCkeEkOsN8/20YGTp+jK0u5vUiy+Reu45eP4FnOeeI7D0+c/cytPaWhLbTyO5y67oHnvg7LQzgboagjYk2YwBlkhyjNoWSU8PoX89RmjRQ4QeXuwljpz/n+q6JKduTWLHHUnsMJ3EdtNIbTQZytD57IgQCjiEAw5B16vTVWlSsTjJF14g/eTTyJNPEHjicQJvvrHGZzQcJrHDdOJ77EV6n31xd5hGMBQkFLDEYkYfSyQ5Rk0iUcV943XC999L6MEHCD3+2Bq/QWsoRGL6jsR33Y34zruS3H5aSW9N5XIdbzhuMJBteVTfF20ylSb+/z5CH30M+ecjBP/5CMEXl67xmXTzOOK770H8KzPQ/fYjuP56Vfv3NaYvSyQ5qjqRJJMEH/8X4XvmEr7vXgJvv7X6LRUhufW2xPfci9iee5OYviPU1JQlTEeEcNArdRJyR998DlUlnkoT/3g5uvhhAosWEnroQQLvvL3G5xLbbU9sxtdIzDyAwNZbEQ4GrLViqpYlEqq4jyQWI7TwH0Tm3U343gU4rStXv5Ue10Jsn32Jf3UGsb32QcdPKG1sORwRIkGHcMAdc1+WyVSaWDJN4j+v4zxwP+H77yP0yOI1WojJjSYT2/8A4gcchLvzlwgFA4QDjpXaN1XDEkmOqmiRxOOEFi8kcufthO+Zt8aw1eTGmxDb/wBiX/s6iR12/HT+RRmM5eQxkHSmEGWsoxMWLiS04B4iC+bjrPi0NE9q4kSiBx1K7JBDcafvSDjkWlIxFc8SSY6KTSSqBJ96gsgtNxO5Yw5OW+vqtxJTtyJ24CFEZx5AarPNvWqGZRQOOESCLpFg+ZJYNVDNJJVogvRjjxKadzeRuXfhLlu2+jPJSRsRPfJoYkcchbv5ZkSCllRMZbJEkqPSEonz7rvU3Ph3IjfdQODdd1ZvT24xheghhxE9+FBSm36hqOccDgFqQi51ocCo6/MohWxSicYSpB9/nMgdcwjfdQfuR5/O5UlsP43eo48jdvgRRNaeQE3QtY56UzEskVBhfSS9vUTm3UXk738j/Mji1ZtTn/sc0cOOJHrEUSS3nFr2lgdYAvFDOq1Ekymi0QSyaBGRW28iPPcunK4uwBtWHJ15INHjTiC9557URIJEAq5df1NWlkhylLNF4r72b2qunU3NjTfgtLcBoJGI96Vx7DeI77Z7Wfs8crmOUBN0qQnaF5ifEqk0vYkU0Y5VhOfPJXL93wk9vGj1HKDkpI3o/eZJRI/7BqH11qMmZP1RpjwskeQoeSJJJAjPn0vtX68g9Og/P928zXb0nnAi0UMOQ5ubhxWPH8IBh5qQSzhQGQltrFBVogkvqaTefpuam26g5u/X4r7vLdejwSCxmQfSc8qp6C67UhMOUBN0rS/FlIwlkhylSiTOJx9Tc9011Fx9Fe6HXmX8dF0d0cOPovebJ5HcdrthxeAHEagJutSGAhU5y3ysSWZaKb3ROMF/PEjNNbMJ37dgdT2wxNSt6Jn1bWKHH0mkqZ6aoGtlWozvLJHk8DuRuK+8TN2lfyJy601IPA5A8gub0TPrNKJHHoM2Ng7r3H4QgbpQgNqQ/WZbibKtlO54En3vPWquvZraa2avHkqcbhlPz8mn0HvKqbif+xw1QZdI0EZ8GX9YIsnhSyJRJbRoIbV/uojwwn94m0SI7fd1emd9m/gee1ZEx3mWdaBXn1gyRU8sRby7h8idt1N7xWUEn3sW8MrhRI84iu7TzyA9Zcrqvi0b8WWKyRIJPo3aSqUI330ndX/8PcGlzwNehdje475Bz2mnk9p44yJEXlyRgEt9xG5hVat4Mk13LEk8mSL4xL+o/fMlhO+Zt7pzPrbf1+n+/lkkdtyJoOtQa5MdTZEUPZGISB0QVdXUoB+uMEVpkcTj1Nx0A7UX/Z7AW28CkFp7HXpO/S96T/oW2lJ5CyIFHKEhYpVpR4tEyksosWQa9623qL3sEmr+fh0SjQIQ33kXus/6EfG990VyRuBZK8UM14gTiYg4eMvgHgvsAMSAMLAcWABcqarD+zW/xEaUSDq7SF41m7qLf796NE1y0kb0nHEmvcccB5FIMUMtChFoCAepCdkorNFodQsllUaWf0Lt5ZdRe9UVOB3tgFc4suvHPyH+1f1AxFopZtiKkUgeBv4B3A28pKrpzPYWYE+8ddbvVNXrixa1T4adSG64gfSPfoSTGYGV3HwLun/4Y6IHHQqBcix3P7hIwKUhYv0gY0EsmaI7liKRSiOdndRcfRV1f7p4dcd8Yqtt6D7nXGL7fR1EbKSeGbJiJJKgqiYGOcmgn6kEw04kV18NJ59MYsupdP/wbGIHHFSWRaEK4YjQWBOwuSBjUDSRoiuWJJVW6Omh9pq/UnvxH1cvq5zYbnu6zj2f+N77rB4AEgm4NtHRDMrXznYRqVfVrmFFVgbDTiSJBNH5C+jYc9+KTSDgjcZqCAfstsUY1xv3EkpaFXp7qbl2NnV/+B3uJx8DEN9pZ7r+92ckdt519T6uI0QyfSnWSjF9+Z1I3lPVDYcVWQlVVK0tH7iO0Gid6SaHqtIdT9ETS6IA3d3U/vUK6i76w+p1bWJf3Y+u8y8g+cUt19jXKhyYvopxa+vMgfYBzlXVyhumNIBKq/5bDNYKMfmk0kpXNEk06Q2ylFWrqL30Emr/dBFOVxcqQvTIo+n6n5+SnrjBGvsGHKE2FLCJjqbgRJLvV9lfAeOAhj6P+kH2Mz5yRGiuDdIYCdo/cjMg1xGaaoM01QRxRNCGBrrPPpcVz79Mz2n/BYEANTffyITtt6L+Z/+LdHau3jeZVjqjCZZ3xVgVTZBMpcv4NzHVIF+L5F/Ad1X1mX7ee19VN+hnt4o0Wlok4YBDYyRoI7LMkKgqq2JJeuOfTgFz3nmHhgv+l8jttwGQnrAWXeecR+83T+p3RGLIdagN222vsaYYLZITgXcHeG/QA5viEaAxEqS5NmRJxAyZiNeXNq42tLpDPT1pEh1X/42VDz1MfMcv4axYTuNZZ9Cy+84EH3n4M8eIp9K09yRY0RWjN55iNFfEMEM3YCJR1ddUdUXuNhFZN/Pex34HZjwBR2ipC9nkQjNioYDD+LoQdeEA2V9HktOm03b/QtqvvZ7UhhsSfOlFWmbOoOn4o3He/ezvkanMba8VXXFvhFjaEooZel/HAl+iMP2qCbm01IWsxIUpGhGhPhzwfq6yrVsRYgcfyoqnnqfrvPPR2loic+9iwvRtqPu/X0OmBEuutCrdsSQrumJ09Fo/ylg31G8ou69SAiLQVGMd6sY/AddhfH2Y2tyWbk0N3T88mxVLXqD3sCOQaJT6X17A+C9tT+iB+/o9juJNiFzZHaetO04sWXWl+EwRDDWRXOVLFD4RkZkicmVHR0e5QylYyHUYXxcmErRbWcZ/DZm+EyfnF5b0+hPpnH0drfPvJ7n5FgTefotxhx9M07FH4nywbMBjZftRVnbFiCYsoYwlgyYSEdlYRMKZl6+IyPdEpHLWic1DVeep6qympqZyh1KQ+nCAcXUhm2FsSirbdxLuM7E18eXdWPnok6z65YWk6+uJzJ/L+OnbUnvZnyA58EjGZFrp6PU65i2hjA2FtEhuB1IisgkwG9gIuNHXqMYYR4RxtV4nqDHl4DhCc22I+r4/g8EgPaefwcqnniO6/wE4XV00nPMjWvbajcDzz+U9ZiqTUKyFMvoVkkjSqpoEDgYuUtXvA+v5G9bY4d3KClmZE1MR6sIBmmuDn1ncM73+RDpuuIW2m+eQ2mADgi88R8teX6b+f8+F3t68x0xaQhn1Cvn2SojI0cAJwPzMtqB/IY0ddZlbWTY3xFSScMBlfF3401FdOeL7fZ2VTzxL93+dDuk0dRf/gfE770Dwn48MetzchGKd8qNLIYnkRGAn4Jeq+raIbARU/BoklUzwRmV95jaCMRXCzcxfivQzk13r6+n69W9p/cdikltMIfDWm7Ts/1UazjwDWbVq0GMn00p7T4K27jgJGzY8KozqNduzKqlESrZWVtDmhpgq0R1L0hUb4N9APE7dH39H3W8vRBIJUhtuSMefryCx+x4FHz8ScKkL25LAlagYJVJMkQUz/SGWREw1GajfBIBQiO4f/4TWxY+R2Hpb3Pfeo+WA/Wj4/veQrsKWLIomvXko7T3WQqlW9o1WIuGAw7haK7hoqlM44NJSO/DQ9OSWU2l96GFvZnwwSO3VV9Gy644En3y84HPEkmlaLaFUJUskJVATcmmuDdksdVPVAm7/801WCwbp/uHZXutky6neRMYZ+1D/0/+BWKzg81hCqT4DJhIRcUXkVBH5uYjs0ue98/wPbXSoDwdojNggNzM6iHjzTfLNeUpuOZXWhf+k+/s/AKDuj7+jZe/dcP/96pDOZQmleuRrkVwB7A6sBC4RkT/kvHeIr1EVSblLpDRGgjbJ0IxKg/6CFA7T9dOf03bvgyQnbUTwxaWM331naq66HIY4wMcSSuXLl0imq+oxqnoRsCNQLyJ3ZMqlVMU9mnKVSMkO77XS72Y0827ZBvN+GSS+tDOtjz5J73HfQKJRGn/wfZqPOATnk6GvRJFNKB09Vm240uRLJKHsE1VNquos4HlgId5yu6YfAjTVBq3oohkTwgHXm1Sbp/9PGxrovPQK2q+7gXTzOMIP3EfLLjsSWvjQsM6ZHeVlCaVy5EskS0RkRu4GVb0AuAaY5GdQ1UoExtWFbDlSM6YEXYeWAoqNxg46hJX/eor4l3fD/eRjxh28v1diJZEY1nlXJ5TeBClbYKus8q2QeJyqfmYRAlX9q6pa73EfjggttTZHxIxNruP9/PdXViVXev2JtN29wBsm7DjUXfwHWr66F8477wz73NFEipVdMVZFE7YEcJkM6VtPRK70K5Bqli0nYTNzzVjmZP4dhAb7d+C6dP/wbNoWPEBq4kSCzyxh/G47Eb5n3rDPrUBPPMXyrpgtAVwGQ/3mG3Sq/FjjOl4JeFtDxJjs8OBgvzW6+krstAsr//kk0a/tj9PRTvMxR1B/7o+HfasLvAFh2SWAV0XtllepDDWRfOJLFFUq6Dp5Z/saMxaJCE21wTWX8R2AtrTQceOtrPrlhWggQN2fL2HcfvvgLHt/RDFkWygru2J0WkLx3ZASiarOGPxTY0Mk4FrJE2PyaIgEC5uMK0LP6WfQdu+DpCZOJPT0U4zfbWdCixaOOAYFei2h+K6QpXanicidIvKsiCwVkRdFZGkpgqtUdeEATbVBK3lizCAKmWuSlZj+JVb+8wlie++Ls3IFzQfvT91vfwPpkQ/xzU0oq6IJ60MpskJaJDfgDfk9FJgJ7J/5c8yxdUSMGbpwIFtrbvDPast42m+7k66zzwWg/hc/pfmow5D29qLEkr3ltSLTKW+jvIqjkESyXFXnqurbqvpu9uF7ZBXGJhoaM3yhgNefmG/i4mquS/c559F+252kx7UQvv9eWvbcFffVV4oWj5LtlI/b8r9FUEgiOV9E/ioiR4vIIdmH75FVEJtoaMzIBQqcuJgV3/errFz8GImpW3urMO69G+G77ihqTGn1lv9t7Y4TT9os+eEqdKndbYAZeLe0sre3xgSbaGhM8WQnLhaaTNKTJtH6wEJ6Dz8Sp7ub5hOO9crSp4rbikik0rT1eLPkrf9k6Aq52b+1qk71PZIKZbeyjCkuJzP3qrU7TrqQPoraWjqvuobktttR/z8/oe6PvyPw6it0XHUN2thY1NiiiRSxZIrGiN3GHopCfs1+QkSm+B7JAERksojMFpE5+bYZY6pHthpEQX0m4A0R/s73aL9jnlf48b4FtOyzO+6bbxY9NlXo6E3Q1h234cIFKiSR7Ao8LyKvDXX4r4hcLSKfiMhLfbbPyBzvDRE5O98xVPUtVT15sG3GmOriVYUIFp5MgPgee9K66FGSm29B4LV/07LXlwktXuRLfPFUmpU2uqsghSSSGcCmwFcY+vDfazP7ryYiLnApsB8wBThaRKaIyFQRmd/nsXaB5zHGVKGA6zCuNljQ0OCs1OTJtD64mNh+X8dpb6P5kJnUXP1XX+LLju5a2R0nlrTRXQMpJJGsB7TmDPttBdYt5OCq+kjm87mmA29kWhVx4GbgQFV9UVX37/OwkizGjHJeMilsnkmWNjbSfuOtdP/3WUgqReP3v0v92T+AZNKXGFNppb0nQXtP3NZA6UchieQvQFfO6+7MtuFaH8gtpLMss61fIjJeRC4HthWRcwba1s9+s0RkiYgsWb58+QjCNcb4Leg6NNeEhrb0quPQ9bNf0HHZlWgwSN1fLvUmL3Z2+hXm6lUarWT9mgpJJKI5V0xV0xQ22mvA4/WzbcD/I6q6UlVPU9WNVfXXA23rZ78rVXWaqk5ba621RhCuMaYUQgHHKz00xP2ixx5P29wFpFvGE37wfsZ9dS+c99/zJUbInR0fpzdut7ugsETyloh8T0SCmccZwFsjOOcyYIOc1xOBD0dwPGPMKBEOuDTWDH3dvMTOu9K68J8kv7AZwVdepmWf3Qk896wPEX4qrUpn1CYzQmGJ5DRgZ+ADvCSwIzBrBOd8GthURDYSkRBwFDB3BMcbkIjMFJErOzo6/Di8McYHkaBbWNXgPlIbbUTrA4u8pXw/+oiWr+1L6N57fIhwTTaZsYBEoqqfqOpRqrq2qq6jqscU2gkuIjcBjwObicgyETlZVZPA6cD9wKvArar68kj+Enlin6eqs5qamvw4vDHGJzUhl7phFEfVceNou2MevUcfi/T00HzMEdRcdYUPEX5WNOEVg+yO+dPhX8l8LWOrqkcPsH0BsMDPcxtjqlt9OEAqrUMvqhgK0fmXq0htNJn6X/2cxh/8N+6y9+k6/wJw/C11pEBXLElvIkVDJDBm6vNZASljTMVqqgkSDgzja0qE7h//xBvRFQhQd9HvaZx1EsRixQ+yH2NtuPCoTiTWR2JM9WuqCQ67aGr02ONpv+UO0vX11Nx2C+MOPRAp4ffBWBkuLIP95UTkzH42dwDPqOrzvkRVZNOmTdMlS5aUOwxjzDCl00prz/BrXwVeeJ7mww/G/fgjElO3pv32u0ivU9C86qIR8W7X1YaqZ2E8EXlGVacN9rlC0vw0vJFb62ces4A9gKtE5EcjCdIYYwqRrRg83NWtk1tvQ+sDi0huvAnBF19g3L574r75RnGDHIQqrIomWdEVG3XlVgpJJOOB7VT1LFU9Cy+xrAXsBnzTx9hGzG5tGTN6uNlkMsz9s2ubJLbdjsC779Dylb0IPPtMUWMsRLb/pG0UzT8pJJFsCMRzXieAz6tqL1CanqthsuG/xowuQdcZ1oTFLJ2wFm3z7ye21z44K5YzbuYMgo88XMQICxfPzD8ZDR3yhSSSG/HWJDlfRM4HHgNuEpE6oHiLKBtjTAEiweHNMcnS+nrab7md3sOOwOnqYtxhBxK+Z14RIxyaWDLNym5vQmO1rn9SyITEnwOnAO14neynqeoFqtqtqsf6HaAxxvRVHw4QGckcjVCIzquuoedbpyKxGE3HH03kxuuLF+AwRBMpVnbFWBWtvhnyedO6iDjAUlXdEij9zURjjBlAY02AVI+SGO5tIcdh1e/+SHrcOOp/eyFN3z4F6Win99unFzfQIcgWhOxNpKgNBagLuchwRxiUUN4WSabS7wsismGJ4ikq62w3ZvQSEZprhrbCYj8Hofu881n1q98A0Hj2D6n77YXeEKsyUvUW1FreFaMnXvkrNBYyj2QhsAPwFN5aJACo6gH+hlY8No/EmNErmUrT2hMf8Xd/5O/X0fjdbyOqdJ9xJl0/+wXDHm9cZK4j3u28YGlLrhQ6j6SQHqufFSEeY4zxRcB1aKoJ0tGTGHhhowJEjz8Bra2ladZJ1F38B6Sri1W/+6Pv9bkKkUorHb0JeuIp6sMBQsMpG+OjQjrbHwbeAYKZ508D/hb6N8aYIQgHXBqGUXq+r9ihh9N+wy1oOEzt7Ctp/M6pkKqcyYPZkvWVVnJl0EQiIqcAc4BsLeb1gbv8DMoYY4ZquKXn+4rP+Bptt92F1tZSc+P1XrFHn9aCH66eeIqV3fHhDzQoskLaR98BdgE6AVT1dWBtP4MqFutsN2ZsKVY/QmL3PWi7Yy7phgZq5txK04nHQzw++I4llEprxRSELCSRxFR19RUUkQB51livJDaz3ZixpzESIDTMasG5EjvtQttd80k3NRGZexfN3zimZGXohyK7fnw563cVcrUfFpGfADUisi9wG1C+aaDGGJOHiNBcGyTgjHzEVXLadNrm3kt6XAvhe++h+bijIBotQpTFlVavfldHT3kmMxaSSM4GlgMvAqfirWx4np9BGWPMSHjJJDSyOSYZyW22pW3+faTHTyD8wH00H3tkRSYTgGgyxYru2NBXlRyhQkZtpVX1KlU9XFUPyzyviltbxpixy3W8lkkxZoIkt5xK27x7vWTyjwdoPvpw6O0twpGLTxU6ekvbOhkwkYjIvExn9WfG1InIZBG5QERO8jc8Y4wZvpFWC86V/OKWtM6/j/SEtQgv/EdFJxPwWietPaUZIJCvRXIK8GXg3yLytIgsEJGFIvI23lDgZ1T16pJEOUw2assYEwm61IaKMyM8NeWLtM6/j9RaaxNe9FBF3+YCr++kFAYtkQIgIpOA9YBe4D+q2uNvWMVlJVKMMe09cWJFWkjKffUVWvafgbNiObGvzKD9+pshHC7KsYtJBNZuiIxg/+IttYuqvqOqj6vq89WWRIwxBqCpJohbhJFcAKktptA2dwHplvFeB/w3jqnHgS4eAAAUoUlEQVS4eSalVFkFW4wxxifZasHFKsOY/OKWXjIZ10L4vgU0nXgcJBJFOnp1sURijBkzAkXsfAdITt0qM2mxmcj8eTTNOqmianOViiUSY8yYEgm61BSp8x0y80zunEe6oYHIHXO8Qo/pyqiBVSp5E4mI7CQil4rIUhFZLiLvZUZvfUdErO6IMaYqNYQDBItQRiUruf002m+70yv0eNMNNJx5RtkXxyqlfPNI7gW+BdwPzMAbtTUFb1Z7BLhbRKpmcStjjMkSEZpqgkVdtyqx0y603Xw7GolQe81fqf/Jj8ZMMhlw+K+ITFDVFXl3LuAz5SQiM4GZm2yyySmvv/56ucMxxlSYaCJFR29xO8hDD95P89GHI4kEXWefS/c55asoVQnDfyfkHGyNAdIi8iWASk4iYNV/jTH5Fbu/BCC+71fpmH0d6jjUX/hLav90UVGPX4nyJZIbc54/3ue9y3yIxRhjSq4hHChKpeBcsQMPpvNSby3AhvPOoeaa2UU9fqXJl0hkgOf9vTbGmKq0ur+kyMeNHnMcnb/9AwAN3/8u4Tm3FvkMlSNfItEBnvf32hhjqlbAdYqy5ntfvbO+zar/vQBRpenUkwk9eH/Rz1EJ8i1wPFFELsFrfWSfk3m9vu+RGWNMCdWEXGLJVNHqcWX1nPkDnLZW6v50Ec3HH03bXfNJfGnnop6j3PIlkh/mPO9b8dAqIBpjRp3GSJAV3bHijtoVoevnv8Jpa6Pm+utoPuJQ2hY8QHLLqUU8SXkNmEhU9bpSBmKMMeXmOEJjJFj0IcGI0Hnxn5H2NiLz59J88Eza7l9IavLk4p6nTAZMJCIyjzx9IapqkxGNMaNOJOgSS6SJJotcMysQoGP2dcgRhxB+eBHNh8yk7YGFpNdep7jnKYN8ne2/A34PvI23DslVmUcX8JL/oRljTHk0RAJFWe/9MyIROm64hcQ22xF4+y2aDz0I6ews/nlKbMBEoqoPq+rDwLaqemRmct88VT0G2LV0IRpjTGk5jtBYk68Lefi0oYG2OXeSnLwxwaXP03TskRCL+XKuUimkatlaIrL6Rp6IbASs5V9IxWNL7RpjhisccIkEizvrPUvXWpu2O+eTWmddwo8s9srPV3HF4EISyfeBxSKyWEQWA4uAM3yNqkisRIoxZiQa/brFBaQnTaJ9zl2kGxuJ3HUHDef8sGqLPOar/rsegKreB2yKlzzOADZT1QdKE54xxpSPiH+3uACSW21N+w23oqEQtZdfVrV1ufK1SK4WkSdE5EJgJ+BlVX1BVav7Zp4xxgxBOFD8wo65ErvtTsdfrgKg4X9+QmTOLb6dyy/5Otv3A/YAFgMHA0+IyB0iMktENixNeMYYU34N4QBukQs75ooddgSrfnkhAI2nnUJo8SLfzuWHvH0kqhpV1ftU9YxMTfqz8Oae/FlEnipJhMYYU2Z+FXbM1XP6GXR/57tIIkHT8UcReLl6ZlkMaa1JVX1bVS/LTEa0IcDGmDEj6DrUhf3rLwHo+sWFRA86BKezk+bDD8L58ANfz1cs+Wa2v82aM9sl57Wq6sZ+BmaMMZWmLhwgnkwTT/k0VNdx6LhiNs5H/4/QE4/TfMQhtN37D7ShwZ/zFUm+Fsk0YIecx3S8me4CPO9/aMYYU3kai7zW+2dEIrTfdBvJjTch+OJSmk44FhJFrv1VZPk621eq6kqgDdgfb/7ITsDXVfXQEsVnjDEVxc0UdvSTtoynfc7dpCesRfihB2k484yKnmOSbx5JUEROBV4BvgwcqKrHqeorJYvOGGMqUCTo36z3rNTkybTdMgetqaH2b9dQe/EffD3fSOS7tfU2cA5wObAA2FpEDsk+ShKdMcZUKD9nvWclp02n4wpvvfeG888jfNcdvp5vuPIlkn/g3c7aGpjZ57G//6EZY0zlyg4J9lvswINZ9bNfANB06skEllTezIt8C1t9s4RxGGNM1QkFvCHB3bGkr+fpOeNM3DffpPZv19B81OG0PvQI6c9/3tdzDkW+4b9n5ttRVSv3hp0xxpRIXcgllkiRTPvYGS7Cqj9cjPveu4QXL2TckYfQ+sAitLHRv3MOQb5bWw2ZxzTg28D6mcdpwBT/QzPGmMpXilnvAASDdFx3A8kvbEbg1VdoOukbkPS3JVSofMN/f6aqPwMmANup6lmqehawPTCxVAECiMhkEZktInNyth0kIleJyN0i8pVSxmOMMbkCJZj1DqDNzbTfcgfplvGEH7yf+vPO9v2chSikRMqGQDzndRyYVOgJRORqEflERF7qs32GiLwmIm+ISN6roapvqerJfbbdpaqnAN8Ejiw0HmOM8UNdOEDQHVLVqWFJTZ5M+w03o8EgdX+5lJqr/+r7OQdTyN/678BTIvJTETkfeBK4bgjnuBaYkbtBRFzgUmA/vNtkR4vIFBGZKiLz+zzWHuT452WOZYwxZdUYCfh/iwtI7LwrnRd7X3sNP/hvyl0teNBEoqq/BE7Em+HeDpyoqr8u9ASq+gjQ2mfzdOCNTEsjDtyMN+HxRVXdv8/jk/6OK57fAPeq6rP9vD9LRJaIyJLly5cXGq4xxgxbqW5xAUSPPZ7u/z4LSaVoOuFY3DffLMl5+5NvZnt99rmqPquqF2cez/X3mSFaH3g/5/WyzLaBYhkvIpcD24rIOZnN3wX2AQ4TkdP67qOqV6rqNFWdttZaVbHEvDFmFKgLBwj4uHZJrq7zLyC239dx2ttoPvowpKOjJOftK1+L5G4R+b2I7CYiddmNmY7vk0XkfvrcshqC/q7ygGPnMnW/TlPVjbOtIVW9RFW3z2y/fJhxGGNM0ZVkFBd41YKvuobkFlMIvPZvmk4+AVKpUpx5zTAGekNV9wYeAk4FXhaRThFZCVwPrAucoKpzBtp/EMuADXJeTwQ+HOaxBiQiM0Xkyo4yZWljzNhUyltc2tBA+01zPh3Jdf55JTlvLtESVJQUkUnAfFXdMvM6APwH2Bv4AHgaOEZVX/bj/NOmTdMlS5b4cWhjjBlQW3fcv7VL+gg++k/GHfg1JJmk44rZRI86BhFYuyEy7GOKyDOZ1XHz8n2smojcBDwObCYiy0TkZFVNAqcD9wOvArf6lUSMMaZcGkt1iwtI7PplVv32j955v/dfBJ4p3S/Pvre9VPXoAbYvwKsqbIwxo5LrCA2RIJ3R0ixM1XvStwgsfYHaa/5K83FH0vbwY9Awyffz+j97poysj8QYU241IZdwoHRftav+7/fEd9oZ98MPaTzuaIjFfD/nqE4kqjpPVWc1NTWVOxRjzBjWECndLS5CITr+diOpiRMJPfkEfOc7vq+umG8eyVQReUJE3heRK0VkXM57lVcQ3xhjKpTrCPWR0oziAkivvQ7t19+CRiJw++3w/vuD7zQC+VokfwF+CkzFG2H1qIhsnHnP/9VcisBubRljKkVtqDS1uLKS225Hx7V/h6eegg039PVc+f5W9ap6n6q2q+rv8EZZ3SciXyLP5MFKYre2jDGVpFS1uLLiX9sfNt3U9/PkSyQiIqu/gVV1EXAoXhHHylmayxhjqkTAdagt0UTFUsqXSH4DbJG7QVWX4k0irMwV6I0xpsLVhVzcEtXiKpV8a7bfOMD294BTfIvIGGNGMRGhMRKkrSc++IerxKA9PyKyVSkC8YN1thtjKlEoULpaXKWQN5GIyD7AZSWKpeiss90YU6nqS1hu3m/55pEcC/wfcHDpwjHGmLGjZOXmfZavbTUbmKKqtrygMcb4IOA61EcCrIomyx3KiOS7tXUBMFtEakoVjDHGjDW1oQChEk5U9EO+ha1+hdcquat04RSXdbYbY6pBQ4knKhZb3jSoqtfj9ZNUJetsN8ZUg2qfqDhoe0pVHypFIMYYM5ZV80TFfKO2fpTz/PA+7/3Kz6CMMWasyU5UrEb5WiRH5Tw/p897M3yIxRhjxrRQwCESdMsdxpDlLdo4wPP+XhtjjCmChnAAqbJv2HyJRAd43t/rimSjtowx1cZxhPoq63jPl0i2FpFOEVkFbJV5nn09tUTxjYiN2jLGVKPaUHWVT8lX/bf6btQZY8wo0VBFFYKrezqlMcaMUqGAQyRQHb/PWyIxxpgKVV8lM94tkRhjTIVyHamKdUsskRhjTAWrDbk4FT4e2BKJMcZUMBGhIVLZrZJRnUhsHokxZjSIBF2CFVxqvnIjKwKbR2KMGS0qeZLiqE4kxhgzWoQCDuFAZX5lV2ZUxhhjPqM+XJnDgS2RGGNMlQi4DpFQ5U1StERijDFVpD5UedWBLZEYY0wVqcTqwJZIjDGmylRadWBLJMYYU4UaKmhZXkskxhhThSqpOrAlEmOMqVKVUh14VCcSK5FijBnNXEeorYCO91GdSKxEijFmtKurgOrAozqRGGPMaCdS/uHAlkiMMabK1YTcsg4HtkRijDGjQH0Z1yyxRGKMMaNAOOASKtOaJZZIjDFmlChXq8QSiTHGjBJB1yESLP0kRUskxhgzipRjzRJLJMYYM4q4jlBT4jVLLJEYY8woU1fiNUsskRhjzCjjOEJdqHQd75ZIjDFmFKotYekUSyTGGDMKlbJ0iiUSY4wZpUo1FNgSiTHGmBGp+EQiIpNFZLaIzMnZtoWIXC4ic0Tk2+WMzxhjxjpfE4mIXC0in4jIS322zxCR10TkDRE5O98xVPUtVT25z7ZXVfU04AhgWvEjN8YYUyi/WyTXAjNyN4iIC1wK7AdMAY4WkSkiMlVE5vd5rD3QgUXkAOBR4CH/wjfGGDMYX7v0VfUREZnUZ/N04A1VfQtARG4GDlTVXwP7D+HYc4G5InIPcGNxIjbGGDNU5egjWR94P+f1ssy2fonIeBG5HNhWRM7JbNtDRC4RkSuABQPsN0tElojIkuXLlxcxfGOMMbnKUXO4vxkyOtCHVXUlcFqfbYuBxflOoqpXAlcCTJs2bcDjG2OMGZlytEiWARvkvJ4IfFiGOIwxxhRBORLJ08CmIrKRiISAo4C5fpxIRGaKyJUdHR1+HN4YYwwgqv7d9RGRm4A9gAnAx8D5qjpbRL4GXAS4wNWq+kvfgvDiWA68m7OpCegYwusJwAofQut7nmLuN9hnBnq/v+2Vcr36O1ex9rHrNfR98n3OrtfQPjeS69V3WzGv1+dVda1BP6WqY+4BXDnE10tKEUcx9xvsMwO939/2Srlew71mdr382Sff5+x6le569d1WyuuVfVT8zHafzBvi61LFUcz9BvvMQO/3t71Srtdwz2XXy5998n3OrtfQPjeS69V3WymvF+Dzra3RQkSWqKrNoC+QXa+hses1NHa9hqYU12ustkiG6spyB1Bl7HoNjV2vobHrNTS+Xy9rkRhjjBkRa5EYY4wZEUskxhhjRsQSiTHGmBGxRDICInKQiFwlIneLyFfKHU+l62+RMrMmEakTkesyP1fHljueamA/V0Pjx/fWmE0kRVp06y5VPQX4JnCkj+GWnV+LlI0FQ7x2hwBzMj9XB5Q82AoxlGs2Vn+ucg3xehX9e2vMJhKKu+jWeZn9RrNr8WmRsjHgWgq8dnhFTLPLLKRKGGOluZbCr5kZ3vUq2vdWOcrIVwQtwqJbIiLAhcC9qvqsvxGXVzGu11g1lGuHVx17IvA8Y/gXvSFes1dKG13lGcr1EpFXKfL31pj9QR3AkBbdAr4L7AMcJiKn5fncaDXiRcrGsIGu3R3AoSLyF8pQ6qLC9XvN7OdqQAP9jBX9e2vMtkgGMNRFty4BLvEvnIo34kXKxrB+r52qdgMnljqYKjHQNbOfq/4NdL2K/r1lLZI12aJbQ2PXa/js2g2dXbOhKdn1skSyppItujVK2PUaPrt2Q2fXbGhKdr3GbCLJLLr1OLCZiCwTkZNVNQmcDtwPvArcqqovlzPOSmHXa/js2g2dXbOhKff1sqKNxhhjRmTMtkiMMcYUhyUSY4wxI2KJxBhjzIhYIjHGGDMilkiMMcaMiCUSY4wxI2KJxIx5IpISkedzHnnL4ZeSiMzJrLfxZCa290RkeU6skwbY7xci8vM+26aJyNLM84dEpMn/v4EZC2weiRnzRKRLVeuLfMxAZkLYSI7xReAXqnpwzrZvAtNU9fQC9r1TVb+Qs+13wEpV/bWInAxMUNXfjCRGY8BaJMYMSETeEZGficizIvKiiGye2V6XWUjoaRF5TkQOzGz/pojcJiLzgAdExBGRy0Tk5cyaLAtE5DAR2VtE7sw5z74ickc/IRwL3F1AnPuJyOOZOG8RkbrMDOaoiGyf+YwAhwM3Z3a7GzhmJNfHmCxLJMZATZ9bW7mrxq1Q1e2AvwA/yGw7F1ioqjsAewK/FZG6zHs7ASeo6l54qx1OAqYC38q8B7AQ2EJE1sq8PhG4pp+4dgGeyRd4ZsGws4G9M3EuBc7IvH0TXn2l7LE+VNW3AVR1BdAgIs35jm9MIayMvDHQq6rbDPBetqXwDF5iAPgKcICIZBNLBNgw8/xBVW3NPN8VuE1V08BHIrIIvDreIvJ34DgRuQYvwXyjn3OvBywfJPad8Va/+5fX6CAEPJp57ybgYRH5EV5CuanPvssz52gf5BzG5GWJxJj8Ypk/U3z670WAQ1X1tdwPisiOQHfupjzHvQZv4aooXrLprz+lFy9J5SPAfap6fN83VPUdEfkQ+DJwMLB9n49EMucwZkTs1pYxQ3c/8N1MvwMisu0An3sUb7VDR0TWAfbIvqGqH+KtDXEe3nrb/XkV2GSQWP4F7C4ikzOx1InIpjnv34S3iNGrqvpRdqOIOMAE1lxBz5hhsURizGf7SC4c5PM/B4LAUhF5KfO6P7fjLS70EnAF8CTQkfP+DcD7qjrQmuP3kJN8+qOqHwMnA7eIyAt4ieULOR+5FdiSTzvZs6YDj6pqKt/xjSmEDf81xkciUq+qXSIyHngK2CXbMhCRPwPPqersAfatARZl9inqF76IXIq3PsXDxTyuGZusj8QYf83PjIwKAT/PSSLP4PWnnDXQjqraKyLnA+sD7xU5rucsiZhisRaJMcaYEbE+EmOMMSNiicQYY8yIWCIxxhgzIpZIjDHGjIglEmOMMSNiicQYY8yI/H+XzkJXeJRoXQAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots()\n",
"ax.plot(sed[\"energy\"], sed[\"e2dnde_fit\"], lw=2, c=\"r\")\n",
"ax.fill_between(sed[\"energy\"], sed[\"e2dnde_lo\"], sed[\"e2dnde_hi\"], alpha=0.1)\n",
"ax.loglog()\n",
"ax.set_xlabel(\"Energy (TeV)\")\n",
"ax.set_ylabel(\"E^2 (dN/dE) (erg cm-2 s-1)\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Some key points to take away:\n",
"\n",
"- YAML is great to store structured data, and easy to read / write with the Python `yaml` module.\n",
"- ECSV is nice to store tabular data if you use `astropy.table.Table`.\n",
"- FITS is also a good format for tabular data. It supports columns which contain multidimensional arrays, which ECSV does not.\n",
"- Separate computation from plotting, for easier reproducibility, debugging and re-use."
]
}
],
"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.6.0"
}
},
"nbformat": 4,
"nbformat_minor": 2
}