\n",
"\n",
"**This is a fixed-text formatted version of a Jupyter notebook**\n",
"\n",
"- Try online [![Binder](https://static.mybinder.org/badge.svg)](https://mybinder.org/v2/gh/gammapy/gammapy-webpage/v0.18.2?urlpath=lab/tree/pulsar_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",
"[pulsar_analysis.ipynb](../_static/notebooks/pulsar_analysis.ipynb) |\n",
"[pulsar_analysis.py](../_static/notebooks/pulsar_analysis.py)\n",
"
\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Pulsar analysis with Gammapy"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Introduction"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This notebook shows how to do a pulsar analysis with Gammapy. It's based on a Vela simulation file from the CTA DC1, which already contains a column of phases. We will produce a phasogram, a phase-resolved map and a phase-resolved spectrum of the Vela pulsar using the class PhaseBackgroundEstimator. \n",
"\n",
"The phasing in itself is not done here, and it requires specific packages like Tempo2 or [PINT](https://nanograv-pint.readthedocs.io)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Opening the data"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's first do the imports and load the only observation containing Vela in the CTA 1DC dataset shipped with Gammapy."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"from gammapy.utils.regions import SphericalCircleSkyRegion\n",
"from astropy.coordinates import SkyCoord\n",
"import astropy.units as u\n",
"\n",
"from gammapy.makers import (\n",
" SafeMaskMaker,\n",
" PhaseBackgroundMaker,\n",
" SpectrumDatasetMaker,\n",
")\n",
"from gammapy.maps import Map, WcsGeom, MapAxis\n",
"from gammapy.data import DataStore\n",
"from gammapy.datasets import Datasets, SpectrumDataset, FluxPointsDataset\n",
"from gammapy.modeling.models import PowerLawSpectralModel, SkyModel\n",
"from gammapy.modeling import Fit\n",
"from gammapy.estimators import FluxPointsEstimator"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Load the data store (which is a subset of CTA-DC1 data):"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"data_store = DataStore.from_dir(\"$GAMMAPY_DATA/cta-1dc/index/gps\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Define obsevation ID and print events:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"EventList\n",
"---------\n",
"\n",
" Instrument : None\n",
" Telescope : CTA\n",
" Obs. ID : 111630\n",
"\n",
" Number of events : 101430\n",
" Event rate : 56.350 1 / s\n",
"\n",
" Time start : 59300.833333333336\n",
" Time stop : 59300.854166666664\n",
"\n",
" Min. energy : 3.00e-02 TeV\n",
" Max. energy : 1.52e+02 TeV\n",
" Median energy : 1.00e-01 TeV\n",
"\n",
" Max. offset : 5.0 deg\n",
"\n"
]
}
],
"source": [
"id_obs_vela = [111630]\n",
"obs_list_vela = data_store.get_observations(id_obs_vela)\n",
"print(obs_list_vela[0].events)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now that we have our observation, let's select the events in 0.2° radius around the pulsar position."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"EventList\n",
"---------\n",
"\n",
" Instrument : None\n",
" Telescope : CTA\n",
" Obs. ID : 111630\n",
"\n",
" Number of events : 843\n",
" Event rate : 0.468 1 / s\n",
"\n",
" Time start : 59300.833333333336\n",
" Time stop : 59300.854166666664\n",
"\n",
" Min. energy : 3.00e-02 TeV\n",
" Max. energy : 4.33e+01 TeV\n",
" Median energy : 1.07e-01 TeV\n",
"\n",
" Max. offset : 1.7 deg\n",
"\n"
]
}
],
"source": [
"pos_target = SkyCoord(ra=128.836 * u.deg, dec=-45.176 * u.deg, frame=\"icrs\")\n",
"on_radius = 0.2 * u.deg\n",
"on_region = SphericalCircleSkyRegion(pos_target, on_radius)\n",
"\n",
"# Apply angular selection\n",
"events_vela = obs_list_vela[0].events.select_region(on_region)\n",
"print(events_vela)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's load the phases of the selected events in a dedicated array."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<Column name='PHASE' dtype='float32' length=10>\n",
"
\n",
"
0.81847286
\n",
"
0.45646095
\n",
"
0.111507416
\n",
"
0.43416595
\n",
"
0.76837444
\n",
"
0.3639946
\n",
"
0.58693695
\n",
"
0.51095676
\n",
"
0.5606985
\n",
"
0.2505703
\n",
"
"
],
"text/plain": [
"\n",
" 0.81847286\n",
" 0.45646095\n",
"0.111507416\n",
" 0.43416595\n",
" 0.76837444\n",
" 0.3639946\n",
" 0.58693695\n",
" 0.51095676\n",
" 0.5606985\n",
" 0.2505703"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"phases = events_vela.table[\"PHASE\"]\n",
"\n",
"# Let's take a look at the first 10 phases\n",
"phases[:10]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Phasogram"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Once we have the phases, we can make a phasogram. A phasogram is a histogram of phases and it works exactly like any other histogram (you can set the binning, evaluate the errors based on the counts in each bin, etc)."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"nbins = 30\n",
"phase_min, phase_max = (0, 1)\n",
"values, bin_edges = np.histogram(\n",
" phases, range=(phase_min, phase_max), bins=nbins\n",
")\n",
"bin_width = (phase_max - phase_min) / nbins\n",
"\n",
"bin_center = (bin_edges[:-1] + bin_edges[1:]) / 2\n",
"\n",
"\n",
"# Poissonian uncertainty on each bin\n",
"values_err = np.sqrt(values)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEWCAYAAABi5jCmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAcD0lEQVR4nO3deZwdZZ3v8c+XJGyyBAiJAYmNEBDkyjJBxLWRRcCFKIIKQmSA6DhwnTsjyjCO4mtmHNQZde5cGG8GkDDsqwSGdQIJooAkyCoiiyGAnYVA2GSEwG/+qKdN1Umf7jqdrrP19/169atPVZ2q8zvPOX2+p56qeloRgZmZWb91Wl2AmZm1FweDmZkVOBjMzKzAwWBmZgUOBjMzK3AwmJlZgYOhjUiaJ+m4VtfR7iQdKenGQZb3SnqqmTUNh6TPS7qt1XWsLUkbSLpa0vOSLm1xLV3Rpq3mYGgySYskvSLpJUlLJf1Y0katrquTRMT5EXFA/7SkkLR9K2vqZJJOlXTeWmziU8AkYIuIOKzOY/wfSUtSeJwtab0699tB0lWSlkt6VtINknZci9psGBwMrfGxiNgI2APYE/h6i+spRdLYVtdgmTZ7Ld4K/CYiVg20UNKHgZOBfYEe4G3At+psazwwB9iRLGx+AVw1wvXaEBwMLRQRTwPXAbvkZr9V0s8kvSjpRkkT+hdIujT3retWSe/ILTtY0q/Sek9L+kpu2fGSHk3fwOZI2iq37ABJD6dtniFpfn93Vtot/5mkH0h6FjhV0naSbpa0QtIzks6XND63vUWSTpJ0n6SXJZ0laZKk61Jt/yVps4HaIz32oen2+9KewMFpej9J9+Tqui3dvjWtfm/aC/t0bnt/JWmZpD5Jx9R7HSQdI+mhVN/jkr6QW9Yr6al625K0RepGeUHSXZL+PldbT3oOY3P3r9tdKOlfJD2ZtrVQ0vtzy06VdJmk8yS9AHx+gPU3kPTPkp5Ir+dtad4aXWvpddpP0oHAKcCnU/vdW6e2nVLtKyU9KOnjaf63gG/k1j92gNVnAGdFxIMR8RzwdwPVDxARv4iIsyLi2Yh4DfgBsKOkLerUtUV6T78g6RfAdjXL3y7ppvTef1jS4TXrDvjajXYOhhaStA1wMPDL3OwjgGOAicC6wFdyy64DpqZldwPn55adBXwhIjYmC5qb02N8CPhH4HBgMvAEcFFaNgG4DPhrYAvgYeA9NWXuBTyeHvMfAKXtbQXsBGwDnFqzzqHA/sAOwMdS3acAE8jec/+7TpPMB3rT7Q+kx/1gbnp+7QoR8YF0c9eI2CgiLk7TbwY2BbYGjgVOrxdIwDLgo8AmZG3/A0l75JYPtq3TgZfTfWakn+G6C9gN2By4ALhU0vq55YeQvV7jKb72/f4J+BOy13Bz4KvAG4M9YERcD3wbuDi1366195E0DrgauJHsfXAicL6kHSPimzXrnzXAw7wDyAfOvcCkeh/2NT4ALImIFXWWnw78N9l7+0/TT3/dbwJuImvLicBngTO0+gvVSL523SUi/NPEH2AR8BKwkuxD+gxgg7RsHvD13H2/BFxfZzvjgQA2TdOLgS8Am9Tc7yzgu7npjYDXyHbpjwZuzy0T8CRwXJr+PLB4iOczHfhlzfM7Mjd9OfBvuekTgZ/U2da+wH3p9vXAccAdaXo+8MlcXbfl1gtg+9x0L/AKMDY3bxnw7pKv0U+ALw+1LWBMassdc8v+vr+21MZRs+68mva9bZA6niMLPMjC99ZB7rtOqnPXAZb1Ak8N8D7cL7ft8wbZ9vuBJcA6uXkXAqeWXP8x4MDc9LjULj1DvA5vAZ4GPltneX/7vz0379u59v808NOadf4/8M2hXrvR/uM9htaYHhHjI+KtEfGliHglt2xJ7vbvyT7IkTRG0mmSHktdCYvSffq7mg4l2/t4InXJ7J3mb0UWQABExEvACrJvv1uRBUH/sgBqz+Z5Mj8haaKki5R1V70AnJerod/S3O1XBpiud7D9dmAHSZPIvjmfC2yT9mzeBdxaZ72BrIhin/cf27KWpIMk3ZG6G1aStWP+OdXb1pbAWIptVGivRqTuqodSN9BKsr2UfB2DbXsCsD7Zh/BI2wp4MiLyex9PkL2HyniJbG+sX//tF+utIGlLsj2UMyLiwjp3G6j9n8jdfiuwV+r+Wpna9EiyPYQRfe26jYOhcxxB1pWwH9kHRk+aL4CIuCsiDiHbZf4JcEla/juyP5Dsztnu9RZk38T6yL6V9S9TfjqpHX73H9O8d0bEJsDn+mtYWxHxe2Ah8GXggYh4Ffg58JfAYxHxzEg8Tp6ys2MuJ+uGmRQR44FrKfeclgOrKLbZNrnbL6ffG+bmvblOHe8HvkbW5bdZquP5mjoGGwr5GbIule0GWPZyvgZJY8g+GMtsF7L30DaS8p8XU8jeQ2U8COS7qHYFlkad7qHUTXcjMCci/mGQ7fa3f77Np+RuPwnMT1/C+n82iog/Y+jXblRzMHSOjYE/kH3b35BslxkASesqO7d/08gO2L0AvJ4WXwAcI2m39CH4beDOiFgE/CfwvyRNTwdI/5w6H1w1dbwErJS0NXDSiD3DzHzgBFYfT5hXMz2QpWRnugzHusB6pA8KSQcBBwy+SiYiXgeuIDsov6Gkt5N1z/UvX0724fm5tMf3pwz8wQ1Zu65KdYyV9A2K37KHquUN4Gzg+5K2So+3d3rNfwOsL+kj6XjB19Nz7rcU6Kn54M+7kyxcvippnKResmNHF5Us71zgWEk7pw/9rwPnDHRHSZsANwA/i4iTB9voAO2/M8XjBNeQ7YEeleoeJ2lPSTsN9dqNdg6GznEu2W7y08CvgDtqlh8FLErdO18k+yZPRMwF/pbsW3Ef2QfTZ9KyZ4DDgO+SBc7OwAKyAKrnW2Sn2T5PFixXrP1TK5hP9iF5a53pgZwKzE7dBYcPcr81RMSLZAfDLyHr0z+C7HTJsk4g24NbAvwHWd97vv2OJwvPFWQHYX9eZzs3kB2k/w3Z6/zfNN618RXgfrKD2M8C3yE7LvA82fGqM8nePy9T7DLsvyhthaS7azea9tw+DhxEtmdyBnB0RPy6TFGRHeD+LnBLem5PkPXzA6DsjLVT0uQnyE7hPiad5dT/M6V2u8kJZN16S8jC5se5x32RLOQ/Q7bXsyS1yXq5dQd77UYtpYMuZqRvjE+RHTy+pdX1dCJJ3wHeHBE+w6XD+LVbzXsMo5ykD0san7ocTiHr067dG7E60nny71TmXWSns17Z6rpsaH7t6munqyetNfYmOw6xLlkX1fSas6RscBuTdUFsRXYa6z/jK3U7hV+7OtyVZGZmBe5KMjOzgo7oSpowYUL09PS0ugwzs46ycOHCZyJiy6HvWdQRwdDT08OCBQtaXYaZWUeR9MTQ91qTu5LMzKzAwWBmZgUOBjMzK3AwmJlZgYPBzMwKHAxmZlbgYDAzswIHg5mZFTgYzMyswMFg1iF6e3vp7e1tdRk2CjgYzMyswMFgZmYFDgYzMytwMJiZWYGDwczMChwMZmZW4GAwM7MCB4OZmRU4GMzMrMDBYGZmBQ4GMzMrcDCYmVmBg8HMzAocDGZmVuBgMDOzAgeDmZkVOBjMzKzAwWBmZgVjq9y4pEXAi8DrwKqImCZpc+BioAdYBBweEc9VWYeZmZXXjD2GfSJit4iYlqZPBuZGxFRgbpo2M7M20YqupEOA2en2bGB6C2owM7M6qg6GAG6UtFDSzDRvUkT0AaTfEwdaUdJMSQskLVi+fHnFZZqZWb9KjzEA742I30maCNwk6ddlV4yIWcAsgGnTpkVVBZqZWVGlewwR8bv0exlwJfAuYKmkyQDp97IqazAzs8ZUFgyS3iRp4/7bwAHAA8AcYEa62wzgqqpqMDOzxlXZlTQJuFJS/+NcEBHXS7oLuETSscBi4LAKazDrWAfv80H6nlz8x+lHnu4DYPftt13jvpO3mcK1t8xvWm3W3SoLhoh4HNh1gPkrgH2relyzbtH35GIu32d1CBxxzUoALthnzWA49JbfNq0u636+8tnMzAocDGZmVuBgMDOzAgeDmZkVOBjMzKzAwWBmZgUOBjMzK3AwmJlZgYPBzMwKHAxmZlbgYDAzswIHg5mZFTgYzMyswMFgZmYFDgYzMytwMJiNsN7eXnp7e1tdhtmwORjMzKzAwWBmZgUOBjMzK3AwmJlZgYPBzMwKHAxmZlbgYDAzs4KxrS7AzMq54KO7t7qEttF/nci8efNaWke38h6DmZkVOBjMzKzAwWBmZgUOBrMSqhr/yOMqWTuqPBgkjZH0S0nXpOnNJd0k6ZH0e7OqazAzs/KascfwZeCh3PTJwNyImArMTdNmZtYmKg0GSW8BPgKcmZt9CDA73Z4NTK+yBjMza0zVeww/BL4KvJGbNyki+gDS74kDrShppqQFkhYsX7684jLNzKxfZcEg6aPAsohYOJz1I2JWREyLiGlbbrnlCFdnZmb1VHnl83uBj0s6GFgf2ETSecBSSZMjok/SZGBZhTWYmVmDKttjiIi/joi3REQP8Bng5oj4HDAHmJHuNgO4qqoazMysca24juE0YH9JjwD7p2kzM2sTTRlELyLmAfPS7RXAvs143H4ecMvMrDxf+WxmZgUedttsLRy8zwfpe3JxYd4jT/cBsPv2265x/8nbTOHaW+Y3pTaz4XIwmK2FvicXc/k+xQA44pqVAFywz5rBcOgtv21KXWZrw11JbajbB1br9udn1ukcDGZmVuBgMDOzAgeDmZkVOBisrVV5PMLHOswG5mAwM7MCn65qZm2v9noRXytSLQeDmbW92utFfK1ItdyVZGYN8bGZ7uc9BrMajQxzsXTpEmDNb6319C1ZUtjGYF0ijW7bbKQ4GMxqNDLMxd7nP93QtuON10t3iTS6bbOR4q4kMzMrcDCYmVmBg8HMzAocDGZmVtBwMEjaTNI7qyjGzMxar1QwSJonaRNJmwP3Aj+W9P1qSxucz6XOuB3MBtcufyPtUkcZZU9X3TQiXpB0HPDjiPimpPuqLGy4GjkH3ZfOm5mtqWwwjJU0GTgc+JsK61lrjZyD7kvnzczWVPYYw7eAG4BHI+IuSW8DHqmuLDMza5Wyewx9EfHHA84R8XirjzF0mv6+xXnz5q2xrBtGjuz252c2mpQNhn8F9igxz4ah20eO7PbnV+uCj+7e6hLM1sqgwSBpb+A9wJaS/jK3aBNgTJWFmZlZawy1x7AusFG638a5+S8An6qqKDMza51BgyEi5gPzJZ0TEU80qSYbxRoZlrqR4xFVDqXdjgY75mMjpx3auYoayh5jWE/SLKAnv05EfKjeCpLWB24F1kvrXJauf9gcuDhtaxFweEQ8N5zirfs0Mix1I8cjqhxK26zblA2GS4EfAWcCr5dc5w/AhyLiJUnjgNskXQd8EpgbEadJOhk4Gfhag3WbmVlFygbDqoj4t0Y2HBEBvJQmx6WfAA4BetP82cA8HAxmZm2jbDBcLelLwJVkewIARMSzg60kaQywENgeOD0i7pQ0KSL60vp9kibWWXcmMBNgypQpJcu0obRDn6h1jmYOMdPIe7PbTwlu9d9p2WCYkX6flJsXwNsGWykiXgd2kzQeuFLSLmULi4hZwCyAadOmRdn1zGzkeIiZ0alUMETEWp2iERErJc0DDgSWSpqc9hYmA8vWZttmZjaySgWDpKMHmh8R5w6yzpbAaykUNgD2A74DzCHbAzkt/b6q0aLNmq2qrotu7xKxzlS2K2nP3O31gX2Bu4G6wQBMBman4wzrAJdExDWSbgcukXQssBg4rPGyzWwktbpP29pL2a6kE/PTkjYF/mOIde4D1vg6FBEryILFzMza0HD/5/PvgakjWYiZmbWHsscYriY7CwmywfN2Ai6pqqhajzz88BqnxnXzcAaDqWrICLNu0Q5/I53+nyTLHmP4p9ztVcATEfFUBfUM6LXXXvVwBklVQ0aYdYt2+Bvp9NN8S3UlpcH0fk02wupmwKtVFmVmZq1TtivpcOB7ZMNXCPhXSSdFxGUV1jZifEqgmVl5ZbuS/gbYMyKWwR+vUfgvoCOCoRUa+XeWnXZcpJH+U+i852c22pUNhnX6QyFZwfDPaBoVGvl3lp12XKSR/lPovOdno5uv6SgfDNdLugG4ME1/Gri2mpKs27u+uv35mXW6of7n8/bApIg4SdIngfeRHWO4HTi/CfWZmVmTDbXH8EPgFICIuAK4AkDStLTsY5VWZ2altMO5+92uXf49bCPHL4drqGDoSUNbFETEAkk9I1ZFh3JfpLWLdjh3v9u1y7+HbeT45XaPLRrWYwwVDOsPsmyDYT2iWQOqPB7hYx1mAxsqGO6SdHxE/Ht+ZhoZdWF1ZdlIqO1egM66LN+6l095Xq2Rv1NoTlsMFQx/Qfaf145kdRBMA9YFPlFlYbb2arsXoLMuy7fu5VOeV2vk7xSa0xaDBkNELAXeI2kfoP/fcv5nRNxceWXWcp12DKXT6u1U7oLrfmX/H8MtwC0V12JmZm2g7AVuZtYlGunT7ua+favPwWA2yjTSp93NfftWn8c7MjOzAgeDmZkVOBjMzKygY48x+JS5jNvBbHDt8jfSLnWU4T0GMzMr6Ng9BjNrf+1+0aFP3R2Yg6FJOmk3cji6/flZd/KpuwNzV5KZmRU4GMzMrMBdSTXq9Yl6mGAzGy0qCwZJ2wDnAm8G3gBmRcS/SNocuBjoARYBh0fEc1XVMVI8TLB1Eh/zaY52aOcqaqiyK2kV8FcRsRPwbuDPJe0MnAzMjYipwNw0bWZmbaKyYIiIvoi4O91+EXgI2Bo4BJid7jYbmF5VDWZm1rimHGOQ1APsDtwJTIqIPsjCQ9LEOuvMBGYCjBs7phlljgr1djs78Xzu2pp9zKeztUO3TLtodVtUHgySNgIuB/4iIl6QVGq9iJgFzALYcP31oroKDTrzfO7amn3Mx2xkVHq6qqRxZKFwfkRckWYvlTQ5LZ8MLKuyBjMza0yVZyUJOAt4KCK+n1s0B5gBnJZ+X1VVDWZWTqu7Lqy9VNmV9F7gKOB+SfekeaeQBcIlko4FFgOHVViDmZk1qLJgiIjbgHoHFPat6nHNzGzteEgMMzMr8JAYZjYiOvGUZxuYg8HMRkQnnvJsA3NXkpmZFXiPwcwsx6fujvJgcJ+omdmaRnUwuE/UzGxNozoY1pZ3Oc2sG/ngs5mZFTgYzMyswF1JVlendZV1Wr1m7cp7DGZmVuBgMDOzAgeDmZkVOBjMzKzAwWBmZgUOBjMzK/DpqmZWGZ9C3Jm8x2BmZgUOBjMzK3BXUg3v+prZaOc9BjMzK3AwmJlZgYPBzMwKHAxmZlbgYDAzswIHg5mZFTgYzMysoLJgkHS2pGWSHsjN21zSTZIeSb83q+rxzcxseKrcYzgHOLBm3snA3IiYCsxN02Zm1kYqC4aIuBV4tmb2IcDsdHs2ML2qxzczs+Fp9pAYkyKiDyAi+iRNrHdHSTOBmQDjxo5pUnlmZta2B58jYlZETIuIaWPHOBjMzJql2cGwVNJkgPR7WZMf38zMhtDsYJgDzEi3ZwBXNfnxzcxsCFWernohcDuwo6SnJB0LnAbsL+kRYP80bWZmbaSyg88R8dk6i/at6jHNzGztte3BZzMzaw0Hg5mZFTgYzMyswMFgZmYFDgYzMytwMJiZWYGDwczMChwMZmZW4GAwM7MCB4OZmRU4GMzMrMDBYGZmBQ4GMzMrcDCYmVmBg8HMzAocDGZmVuBgMDOzAgeDmZkVOBjMzKzAwWBmZgUOBjMzK3AwmJlZgYPBzMwKHAxmZlbgYDAzswIHg5mZFTgYzMyswMFgZmYFLQkGSQdKeljSo5JObkUNZmY2sKYHg6QxwOnAQcDOwGcl7dzsOszMbGCt2GN4F/BoRDweEa8CFwGHtKAOMzMbgCKiuQ8ofQo4MCKOS9NHAXtFxAk195sJzEyTuwAPNLXQ9jUBeKbVRbQJt8VqbovV3Bar7RgRGze60tgqKhmCBpi3RjpFxCxgFoCkBRExrerCOoHbYjW3xWpui9XcFqtJWjCc9VrRlfQUsE1u+i3A71pQh5mZDaAVwXAXMFXStpLWBT4DzGlBHWZmNoCmdyVFxCpJJwA3AGOAsyPiwSFWm1V9ZR3DbbGa22I1t8VqbovVhtUWTT/4bGZm7c1XPpuZWYGDwczMCtoqGIYaKkOZ/5uW3ydpj1bU2Qwl2uLI1Ab3Sfq5pF1bUWfVyg6fImlPSa+n62S6Upm2kNQr6R5JD0qa3+wam6XE38emkq6WdG9qi2NaUWczSDpb0jJJA17rNazPzYhoix+yA9GPAW8D1gXuBXauuc/BwHVk10K8G7iz1XW3sC3eA2yWbh/UjW1Rph1y97sZuBb4VKvrbuF7YjzwK2BKmp7Y6rpb2BanAN9Jt7cEngXWbXXtFbXHB4A9gAfqLG/4c7Od9hjKDJVxCHBuZO4Axkua3OxCm2DItoiIn0fEc2nyDrLrQbpN2eFTTgQuB5Y1s7gmK9MWRwBXRMRigIjo1vYo0xYBbCxJwEZkwbCquWU2R0TcSvb86mn4c7OdgmFr4Mnc9FNpXqP36QaNPs9jyb4RdJsh20HS1sAngB81sa5WKPOe2AHYTNI8SQslHd206pqrTFv8P2Ansotn7we+HBFvNKe8ttPw52YrhsSop8xQGaWG0+gCpZ+npH3IguF9lVbUGmXa4YfA1yLi9ezLYdcq0xZjgT8B9gU2AG6XdEdE/Kbq4pqsTFt8GLgH+BCwHXCTpJ9GxAtVF9eGGv7cbKdgKDNUxmgZTqPU85T0TuBM4KCIWNGk2pqpTDtMAy5KoTABOFjSqoj4SXNKbJqyfx/PRMTLwMuSbgV2BbotGMq0xTHAaZF1sj8q6bfA24FfNKfEttLw52Y7dSWVGSpjDnB0Osr+buD5iOhrdqFNMGRbSJoCXAEc1YXfCPsN2Q4RsW1E9ERED3AZ8KUuDAUo9/dxFfB+SWMlbQjsBTzU5DqboUxbLCbbc0LSJGBH4PGmVtk+Gv7cbJs9hqgzVIakL6blPyI76+Rg4FHg92TfCrpOybb4BrAFcEb6trwqumxEyZLtMCqUaYuIeEjS9cB9wBvAmRHRdcPVl3xf/B1wjqT7ybpSvhYRXTkUt6QLgV5ggqSngG8C42D4n5seEsPMzAraqSvJzMzagIPBzMwKHAxmZlbgYDAzswIHg5mZFTgYzIA0Mus9kh6QdKmkDSX11Bux0qybORjMMq9ExG4RsQvwKvDFVhdk1ioOBrM1/RTYPt0eI+nf05j+N0raAEDS8ZLuSuP9X56uNEbSYWmv4940JAWSxkj6Xrr/fZK+0JqnZVaOg8EsR9JYsv9vcX+aNRU4PSLeAawEDk3zr4iIPSNiV7JhJ45N878BfDjN/3iadyzZMAR7AnsCx0vatvpnYzY8DgazzAaS7gEWkI2zc1aa/9uIuCfdXgj0pNu7SPppGnLhSOAdaf7PyIZiOJ5suAaAA8jGqrkHuJNsKJOpVT4Zs7XRNmMlmbXYKxGxW35GGoPqD7lZr5MNZw1wDjA9Iu6V9HmysWqIiC9K2gv4CHCPpN3Ixuo5MSJuqPIJmI0U7zGYDc/GQJ+kcWR7DABI2i4i7oyIbwDPkA13fAPwZ+m+SNpB0ptaUbRZGd5jMBuevyXrFnqC7HjExmn+9yRNJdtLmEv2/4jvI+uCujv9q8nlwPRmF2xWlkdXNTOzAnclmZlZgYPBzMwKHAxmZlbgYDAzswIHg5mZFTgYzMyswMFgZmYF/wOa4Lb7pBHpSwAAAABJRU5ErkJggg==\n",
"text/plain": [
"