\n",
"\n",
"**This is a fixed-text formatted version of a Jupyter notebook**\n",
"\n",
"- Try online[![Binder](https://static.mybinder.org/badge.svg)](https://mybinder.org/v2/gh/gammapy/gammapy-webpage/v0.19?urlpath=lab/tree/tutorials/analysis/time/pulsar_analysis.ipynb)\n",
"- You may download all the notebooks in the documentation as a\n",
"[tar file](../../../_downloads/notebooks-0.19.tar).\n",
"- **Source files:**\n",
"[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"
]
},
{
"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": {
"execution": {
"iopub.execute_input": "2021-11-22T21:06:28.132270Z",
"iopub.status.busy": "2021-11-22T21:06:28.131270Z",
"iopub.status.idle": "2021-11-22T21:06:28.295760Z",
"shell.execute_reply": "2021-11-22T21:06:28.295986Z"
}
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"execution": {
"iopub.execute_input": "2021-11-22T21:06:28.298595Z",
"iopub.status.busy": "2021-11-22T21:06:28.298264Z",
"iopub.status.idle": "2021-11-22T21:06:28.816236Z",
"shell.execute_reply": "2021-11-22T21:06:28.816402Z"
}
},
"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, RegionGeom\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": {
"execution": {
"iopub.execute_input": "2021-11-22T21:06:28.818605Z",
"iopub.status.busy": "2021-11-22T21:06:28.818318Z",
"iopub.status.idle": "2021-11-22T21:06:28.830699Z",
"shell.execute_reply": "2021-11-22T21:06:28.830864Z"
}
},
"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": {
"execution": {
"iopub.execute_input": "2021-11-22T21:06:28.833368Z",
"iopub.status.busy": "2021-11-22T21:06:28.833036Z",
"iopub.status.idle": "2021-11-22T21:06:28.853294Z",
"shell.execute_reply": "2021-11-22T21:06:28.853516Z"
}
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"No HDU found matching: OBS_ID = 111630, HDU_TYPE = rad_max, HDU_CLASS = None\n"
]
},
{
"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": {
"execution": {
"iopub.execute_input": "2021-11-22T21:06:28.856026Z",
"iopub.status.busy": "2021-11-22T21:06:28.855744Z",
"iopub.status.idle": "2021-11-22T21:06:28.874456Z",
"shell.execute_reply": "2021-11-22T21:06:28.874680Z"
}
},
"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": {
"execution": {
"iopub.execute_input": "2021-11-22T21:06:28.877542Z",
"iopub.status.busy": "2021-11-22T21:06:28.877216Z",
"iopub.status.idle": "2021-11-22T21:06:28.879051Z",
"shell.execute_reply": "2021-11-22T21:06:28.879226Z"
}
},
"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": {
"execution": {
"iopub.execute_input": "2021-11-22T21:06:28.881702Z",
"iopub.status.busy": "2021-11-22T21:06:28.881424Z",
"iopub.status.idle": "2021-11-22T21:06:28.882514Z",
"shell.execute_reply": "2021-11-22T21:06:28.882836Z"
}
},
"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": {
"execution": {
"iopub.execute_input": "2021-11-22T21:06:28.902193Z",
"iopub.status.busy": "2021-11-22T21:06:28.901886Z",
"iopub.status.idle": "2021-11-22T21:06:28.951129Z",
"shell.execute_reply": "2021-11-22T21:06:28.951409Z"
},
"nbsphinx-thumbnail": {
"tooltip": "Produce a phasogram, phased-resolved maps and spectra in pulsar analysis."
}
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEWCAYAAABi5jCmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAcD0lEQVR4nO3deZwdZZ3v8c+XJGyyBAiJAYmNEBDkyjJBxLWRRcCFKIIKQmSA6DhwnTsjyjCO4mtmHNQZde5cGG8GkDDsqwSGdQIJooAkyCoiiyGAnYVA2GSEwG/+qKdN1Umf7jqdrrP19/169atPVZ2q8zvPOX2+p56qeloRgZmZWb91Wl2AmZm1FweDmZkVOBjMzKzAwWBmZgUOBjMzK3AwmJlZgYOhjUiaJ+m4VtfR7iQdKenGQZb3SnqqmTUNh6TPS7qt1XWsLUkbSLpa0vOSLm1xLV3Rpq3mYGgySYskvSLpJUlLJf1Y0katrquTRMT5EXFA/7SkkLR9K2vqZJJOlXTeWmziU8AkYIuIOKzOY/wfSUtSeJwtab0699tB0lWSlkt6VtINknZci9psGBwMrfGxiNgI2APYE/h6i+spRdLYVtdgmTZ7Ld4K/CYiVg20UNKHgZOBfYEe4G3At+psazwwB9iRLGx+AVw1suXaUBwMLRQRTwPXAbvkZr9V0s8kvSjpRkkT+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+lHnu4DYPftt13jvpO3mcK1t8xvWm3W3SoLhoh4HNh1gPkrgH2relyzbtH35GIu32d1CBxxzUoALthnzWA49JbfNqssGwV85bOZmRU4GMzMrMDBYGZmBQ4GMzMrcDCYmVmBg8HMzAocDGZmVuBgMDOzAgeDmZkVOBjMzKzAwWBmZgUOBjMzK3AwmJlZgYPBzMwKHAxmZlbgYDAbYb29vfT29ra6DLNhczCYmVmBg8HMzAocDGZmVuBgMDOzAgeDmZkVOBjMzKzAwWBmZgVjW12AmZVzwUd3b3UJbaP/OpF58+a1tI5u5T0GMzMrcDCYmVmBg8HMzAocDGYlVDX+kcdVsnZUeTBIGiPpl5KuSdObS7pJ0iPp92ZV12BmZuU1Y4/hy8BDuemTgbkRMRWYm6bNzKxNVBoMkt4CfAQ4Mzf7EGB2uj0bmF5lDWZm1piq9xh+CHwVeCM3b1JE9AGk3xMHWlHSTEkLJC1Yvnx5xWWamVm/yoJB0keBZRGxcDjrR8SsiJgWEdO23HLLEa7OzMzqqfLK5/cCH5d0MLA+sImk84ClkiZHRJ+kycCyCmswM7MGVbbHEBF/HRFviYge4DPAzRHxOWAOMCPdbQZwVVU1mJlZ41pxHcNpwP6SHgH2T9NmZtYmmjKIXkTMA+al2yuAfZvxuP084JaZWXm+8tnMzAo87LbZWjh4nw/S9+TiwrxHnu4DYPftt13j/pO3mcK1t8xvSm1mw+VgMFsLfU8u5vJ9igFwxDUrAbhgnzWD4dBbftuMsszWiruS2lC3D6zW7c/PrNM5GMzMrMDBYGZmBQ4GMzMrcDBYW6vyeISPdZgNzMFgZmYFPl3VzNpe7fUivlakWg4GM2t7tdeL+FqRarkrycwa4mMz3c97DGY1GhnmYunSJcCa31rr6VuypLCNwbpEGt222UhxMJjVaGSYi73Pf7qhbccbr5fuEml022YjxV1JZmZW4GAwM7MCB4OZmRU4GMzMrKDhYJC0maR3VlGMmZm1XqlgkDRP0iaSNgfuBX4s6fvVljY4n0udcTuYDa5d/kbapY4yyp6uumlEvCDpOODHEfFNSfdVWdhwNXIOui+dNzNbU9lgGCtpMnA48DcV1rPWGjkH3ZfOm5mtqewxhm8BNwCPRsRdkt4GPFJdWWZm1ipl9xj6IuKPB5wj4vFWH2PoNP19i/PmzVtjWTeMHNntz89sNCkbDP8K7FFing1Dt48c2e3Pr9YFH9291SWYrZVBg0HS3sB7gC0l/WVu0SbAmCoLMzOz1hhqj2FdYKN0v41z818APlVVUWZm1jqDBkNEzAfmSzonIp5oUk02ijUyLHUjxyOqHEq7HQ12zMdGTju0cxU1lD3GsJ6kWUBPfp2I+FC9FSStD9wKrJfWuSxd/7A5cHHa1iLg8Ih4bjjFW/dpZFjqRo5HVDmUtlm3KRsMlwI/As4EXi+5zh+AD0XES5LGAbdJug74JDA3Ik6TdDJwMvC1Bus2M7OKlA2GVRHxb41sOCICeClNjks/ARwC9Kb5s4F5OBjMzNpG2WC4WtKXgCvJ9gQAiIhnB1tJ0hhgIbA9cHpE3ClpUkT0pfX7JE2ss+5MYCbAlClTSpZpQ2mHPlHrHM0cYqaR92a3nxLc6r/TssEwI/0+KTcvgLcNtlJEvA7sJmk8cKWkXcoWFhGzgFkA06ZNi7LrmdnI8RAzo1OpYIiItTpFIyJWSpoHHAgslTQ57S1MBpatzbbNzGxklQoGSUcPND8izh1knS2B11IobADsB3wHmEO2B3Ja+n1Vo0WbNVtVXRfd3iVinalsV9KeudvrA/sCdwN1gwGYDMxOxxnWAS6JiGsk3Q5cIulYYDFwWONlm9lIanWftrWXsl1JJ+anJW0K/McQ69wHrPF1KCJWkAWLmZm1oeH+z+ffA1NHshAzM2sPZY8xXE12FhJkg+ftBFxSVVG1Hnn44TVOjevm4QwGU9WQEWbdoh3+Rjr9P0mWPcbwT7nbq4AnIuKpCuoZ0GuvverhDJKqhoww6xbt8DfS6af5lupKSoPp/ZpshNXNgFerLMrMzFqnbFfS4cD3yIavEPCvkk6KiMsqrG3E+JRAM7PyynYl/Q2wZ0Qsgz9eo/BfQEcEQys08u8sO+24SCP9p9B5z89stCsbDOv0h0KyguGf0TQqNPLvLDvtuEgj/afQec/PRjdf01E+GK6XdANwYZr+NHBtNSVZt3d9dfvzM+t0Q/3P5+2BSRFxkqRPAu8jO8ZwO3B+E+ozM7MmG2qP4YfAKQARcQVwBYCkaWnZxyqszcxKaodz97tdu/x72EaOXw7XUMHQk4a2KIiIBZJ6RqyKDuW+SGsX7XDufrdrl38P28jxy+0eWzSsxxgqGNYfZNkGw3pEswZUeTzCxzrMBjZUMNwl6fiI+Pf8zDQy6sLqyrKRUNu9AJ11Wb51L5/yvFojf6fQnLYYKhj+guw/rx3J6iCYBqwLfKLCumwE1HYvQGddlm/dy6c8r9bI3yk0py0GDYaIWAq8R9I+QP+/5fzPiLi58sqs5TrtGEqn1dup3AXX/cr+P4ZbgFsqrsXMzNpA2QvczKxLNNKn3c19+1afg8FslGmkT7ub+/atPo93ZGZmBQ4GMzMrcDCYmVlBxx5j8ClzGbeD2eDa5W+kXeoow3sMZmZW0LF7DGbW/tr9okOfujswB0OTdNJu5HB0+/Oz7uRTdwfmriQzMytwMJiZWYG7kmrU6xP1MMFmNlpUFgyStgHOBd4MvAHMioh/kbQ5cDHQAywCDo+I56qqY6R4mGDrJD7m0xzt0M5V1FBlV9Iq4K8iYifg3cCfS9oZOBmYGxFTgblp2szM2kRlwRARfRFxd7r9IvAQsDVwCDA73W02ML2qGszMrHFNOcYgqQfYHbgTmBQRfZCFh6SJddaZCcwEGDd2TDPKHBXq7XZ24vnctTX7mE9na4dumXbR6raoPBgkbQRcDvxFRLwgqdR6ETELmAWw4frrRXUVGnTm+dy1NfuYj9nIqPR0VUnjyELh/Ii4Is1eKmlyWj4ZWFZlDWZm1pgqz0oScBbwUER8P7doDjADOC39vqqqGsysnFZ3XVh7qbIr6b3AUcD9ku5J804hC4RLJB0LLAYOq7AGMzNrUGXBEBG3AfUOKOxb1eOamdna8ZAYZmZW4CExzGxEdOIpzzYwB4OZjYhOPOXZBuauJDMzK/Aeg5lZjk/dHeXB4D5RM7M1jepgcJ+omdmaRnUwrC3vcppZN/LBZzMzK3AwmJlZgbuSrK5O6yrrtHrN2pX3GMzMrMDBYGZmBQ4GMzMrcDCYmVmBg8HMzAocDGZmVuDTVc2sMj6FuDN5j8HMzAocDGZmVuCupBre9TWz0c57DGZmVuBgMDOzAgeDmZkVOBjMzKzAwWBmZgUOBjMzK3AwmJlZQWXBIOlsScskPZCbt7mkmyQ9kn5vVtXjm5nZ8FS5x3AOcGDNvJOBuRExFZibps3MrI1UFgwRcSvwbM3sQ4DZ6fZsYHpVj29mZsPT7CExJkVEH0BE9EmaWO+OkmYCMwHGjR3TpPLMzKxtDz5HxKyImBYR08aOcTCYmTVLs4NhqaTJAOn3siY/vpmZDaHZwTAHmJFuzwCuavLjm5nZEKo8XfVC4HZgR0lPSToWOA3YX9IjwP5p2szM2khlB58j4rN1Fu1b1WOamdnaa9uDz2Zm1hoOBjMzK3AwmJlZgYPBzMwKHAxmZlbgYDAzswIHg5mZFTgYzMyswMFgZmYFDgYzMytwMJiZWYGDwczMChwMZmZW4GAwM7MCB4OZmRU4GMzMrMDBYGZmBQ4GMzMrcDCYmVmBg8HMzAocDGZmVuBgMDOzAgeDmZkVOBjMzKzAwWBmZgUOBjMzK3AwmJlZgYPBzMwKWhIMkg6U9LCkRyWd3IoazMxsYE0PBkljgNOBg4Cdgc9K2rnZdZiZ2cBascfwLuDRiHg8Il4FLgIOaUEdZmY2AEVEcx9Q+hRwYEQcl6aPAvaKiBNq7jcTmJkmdwEeaGqh7WsC8Eyri2gTbovV3BaruS1W2zEiNm50pbFVVDIEDTBvjXSKiFnALABJCyJiWtWFdQK3xWpui9XcFqu5LVaTtGA467WiK+kpYJvc9FuA37WgDjMzG0ArguEuYKqkbSWtC3wGmNOCOszMbABN70qKiFWSTgBuAMYAZ0fEg0OsNqv6yjqG22I1t8VqbovV3BarDastmn7w2czM2puvfDYzswIHg5mZFbRVMAw1VIYy/zctv0/SHq2osxlKtMWRqQ3uk/RzSbu2os6qlR0+RdKekl5P18l0pTJtIalX0j2SHpQ0v9k1NkuJv49NJV0t6d7UFse0os5mkHS2pGWSBrzWa1ifmxHRFj9kB6IfA94GrAvcC+xcc5+DgevIroV4N3Bnq+tuYVu8B9gs3T6oG9uiTDvk7nczcC3wqVbX3cL3xHjgV8CUND2x1XW3sC1OAb6Tbm8JPAus2+raK2qPDwB7AA/UWd7w52Y77TGUGSrjEODcyNwBjJc0udmFNsGQbRERP4+I59LkHWTXg3SbssOnnAhcDixrZnFNVqYtjgCuiIjFABHRre1Rpi0C2FiSgI3IgmFVc8tsjoi4lez51dPw52Y7BcPWwJO56afSvEbv0w0afZ7Hkn0j6DZDtoOkrYFPAD9qYl2tUOY9sQOwmaR5khZKOrpp1TVXmbb4f8BOZBfP3g98OSLeaE55bafhz81WDIlRT5mhMkoNp9EFSj9PSfuQBcP7Kq2oNcq0ww+Br0XE69mXw65Vpi3GAn8C7AtsANwu6Y6I+E3VxTVZmbb4MHAP8CFgO+AmST+NiBcqrq0dNfy52U7BUGaojNEynEap5ynpncCZwEERsaJJtTVTmXaYBlyUQmECcLCkVRHxk6ZU2Dxl/z6eiYiXgZcl3QrsCnRbMJRpi2OA0yLrZH9U0m+BtwO/aE6JbaXhz8126koqM1TGHODodJT93cDzEdHX7EKbYMi2kDQFuAI4qgu/EfYbsh0iYtuI6ImIHuAy4EtdGApQ7u/jKuD9ksZK2hDYC3ioyXU2Q5m2WEy254SkScCOwONNrbJ9NPy52TZ7DFFnqAxJX0zLf0R21snBwKPA78m+FXSdkm3xDWAL4Iz0bXlVdNmIkiXbYVQo0xYR8ZCk64H7gDeAMyOi64arL/m++DvgHEn3k3WlfC0iunIobkkXAr3ABElPAd8ExsHwPzc9JIaZmRW0U1eSmZm1AQeDmZkVOBjMzKzAwWBmZgUOBjMzK3AwmAFpZNZ7JD0g6VJJG0rqqTdipVk3czCYZV6JiN0iYhfgVeCLrS7IrFUcDGZr+imwfbo9RtK/pzH9b5S0AYCk4yXdlcb7vzxdaYykw9Jex71pSAokjZH0vXT/+yR9oTVPy6wcB4NZjqSxZP/f4v40aypwekS8A1gJHJrmXxERe0bErmTDThyb5n8D+HCa//E071iyYQj2BPYEjpe0beVPxmyYHAxmmQ0k3QMsIBtn56w0/7cRcU+6vRDoSbd3kfTTNOTCkcA70vyfkQ3FcDzZcA0AB5CNVXMPcCfZUCZTK3smZmupbcZKMmuxVyJit/yMNAbVH3KzXicbzhrgHGB6RNwr6fNkY9UQEV+UtBfwEeAeSbuRjdVzYkTcUF35ZiPHewxmw7Mx0CdpHNkeAwCStouIOyPiG8AzZMMd3wD8WbovknaQ9KZWFG1WhvcYzIbnb8m6hZ4gOx6xcZr/PUlTyfYS5pL9P+L7yLqg7k7/anI5ML3J9ZqV5tFVzcyswF1JZmZW4GAwM7MCB4OZmRU4GMzMrMDBYGZmBQ4GMzMrcDCYmVnB/wCa4Lb7uGjg3gAAAABJRU5ErkJggg==\n",
"text/plain": [
"