\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.16?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/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.cube import SafeMaskMaker\n",
"from gammapy.maps import Map, WcsGeom\n",
"from gammapy.data import DataStore\n",
"from gammapy.modeling.models import PowerLawSpectralModel, SkyModel\n",
"from gammapy.modeling import Fit, Datasets\n",
"from gammapy.spectrum import (\n",
" PhaseBackgroundMaker,\n",
" SpectrumDataset,\n",
" SpectrumDatasetMaker,\n",
" FluxPointsEstimator,\n",
" FluxPointsDataset,\n",
")"
]
},
{
"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 info:\n",
"- Number of events: 101430\n",
"- Median energy: 0.1 TeV\n",
"- OBS_ID = 111630\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 info:\n",
"- Number of events: 843\n",
"- Median energy: 0.107 TeV\n",
"- OBS_ID = 111630\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+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAcD0lEQVR4nO3deZwdZZ3v8c+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": [
"