\n",
"\n",
"**This is a fixed-text formatted version of a Jupyter notebook**\n",
"\n",
"- Try online [![Binder](https://mybinder.org/badge.svg)](https://mybinder.org/v2/gh/gammapy/gammapy-webpage/v0.14?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/en/latest/readme.html)."
]
},
{
"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.maps import Map, WcsGeom\n",
"from gammapy.cube import fill_map_counts\n",
"from gammapy.data import DataStore\n",
"from gammapy.modeling.models import PowerLawSpectralModel\n",
"from gammapy.modeling import Fit, Datasets\n",
"from gammapy.spectrum import (\n",
" PhaseBackgroundEstimator,\n",
" SpectrumExtraction,\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+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAcFUlEQVR4nO3deZwdZZ3v8c+XJGwSIEiCEYkNhJ3LdgO4DGPCIohLcOM6okQEIjJ6Z+5lEERFnGFm8OpVdMYZblQkDJsoIMsIyAQSUAFJNGwChiUkQJMEJIRNIfC7f9TTpKpzTnedTtfZ+vt+vfrVtZ9fPeec+lU9T9VzFBGYmZn1Wa/VAZiZWXtxYjAzswInBjMzK3BiMDOzAicGMzMrcGIwM7MCJ4YWkjRX0nGtjqPTSXpe0nYDzF8s6eBmxjQUkkLS5FbHsa4kfVbSsvS+vLHFsXRFmTabE0PF0kHppfQlWSbpR5I2aXVc3SQiNomIhwEknSfpzFbH1Kkk9aSD6eghrj8G+Bbw7vS+PF1jmb0kLZD0Yvq/1wDb+6akRZKek3S/pKOHEpc1xomhOd4fEZsA+wD7Al9ucTzrbKgHDht+bfZebAVsCNxba6ak9YErgQuAccBs4Mo0vZYXgPcDmwEzgO9IesdwB21FTgxNFBGPA9cCu+cmv1XSr9IZ0S8kbdk3Q9JPJD0p6VlJN0vaLTfvcEm/T+s9LunvcvPeJ2mhpJWSfi1pj9y8XVIV1kpJ90r6QG7eGyVdLWmVpDsknSnpl7n5IemvJS0CFqVp35G0NK2zQNIBueXPSPtwQYrzbkk7SvqipOVpvXfXKitJx0i6Ojf+oKRLc+NL+840+6oLJM0EjgK+kK7Qrs5tci9Jd6Wy/LGkDeu87vaSbpT0tKSnJF0oafPc/MWS/q7etiR9QVKvpCckHZevyuhfdSjpU/ny7RfHeyX9LpXrUkln5Ob1ndUfK2kJcGOdbUxPn4NVkh6SdFhuHw7OLXeGpAvS6M3p/8pUhm+vsd0NJJ2d9vGJNLyBpB2BB3Lr14prKjAaODsi/hwR3wUEHFhrHyLiqxFxf0S8FhG3A7cAa8WUi+3kXPl/ukbc35S0RNnV+zmSNsrNr/vejTgR4b8K/4DFwMFpeBuyM6l/SONzgYeAHYGN0vhZuXU/DYwFNgDOBhbm5vUCB6ThccA+aXgfYDmwPzCK7CxrcdrGGOBB4DRgfbIv43PATmndS9LfxsCuwFLgl7nXDOAGYAtgozTtE8Abyb7sJwFPAhumeWcAfwIOTfPPBx4BvpRiOR54pE65bQesJDt5mQg8Cjyem/cMsF4urslp+DzgzBrvwW+AN6fY7wNOqPO6k4FDUnmNJztQnl1mW8Bhaf93S2X4H/1imwscl9vWp2qUb9+yU4H/lvZ/D2AZcESa15OWPR94Q9970W8/9gOeTfuyHrA1sHP/z2Tufbqg37ZHD/CZ/nvgNmBCKqNfs+YzPeD6wP8Cru037RrgpBLfpY3IPveH1Zl/WCqn3VO5XNSvTM8Grkrv21jgauCfy7x3I+2v5QF0+1/6Ej5PdpB7FPg31hxU5wJfzi17InBdne1snj6om6XxJcBngE37LffvfV/S3LQHgHcBB6QP/3q5eRenA8Mo4BVSkkjzzqxx4DpwkP19BtgzDZ8B3JCb9/5UFqPS+Ni0zc3rbGspWaL7GDCL7IC8M3AMcFW/uAZLDJ/Ijf8f4JyS798RwO/KbAs4t+9Ak8YnM8TEUCOOs4Fvp+GetOx2A8T9//qWr/OZXJfE8BBweG78UGBxmfWBrwCX9Jt2IXBGifdiNnAdoDrzz6V4YrVjX5mSXZW8AGyfm/920onJYO/dSPtzVVJzHBERm0fEWyPixIh4KTfvydzwi8AmAJJGSTorVQGsIvsyA/RVNX0YOBx4VNK83CX/W4GTUlXRSkkrya5U3pz+lkbEa7nXfJTsbHI82Vn90ty8/HDNaZJOknRfqlZZSVYXvGVukWW54ZeApyLi1dw4fftcwzyyM+e/TMNzyRLcu9J4I2qWc3+SJki6RFn13CqyuvAt+y1Wb1tvZvDyK0XS/pJukrRC0rPACTXiGGj725AdwKvwZrLPTZ9H07Qyngc27TdtU7Ir17okfYPsSuDISEfuOnHlyyQf43iyK4EFue/FdWl6rXWH/N51AyeG9vVxYDpwMNnBtidNF0BE3BER08ku538G9NW/LwX+MSWivr+NI+Ji4AlgG0n5930S8DiwAlgNvCU3b5sacb3+pUztCacARwLjImJzsuoLDXmvi/oSwwFpeB6DJ4Z17S74n9M29oiITcmqysruTy8Dl98LZAenPm8aYFsXkVV7bBMRmwHn1IhjoH1dCmxfZ95AcZQpvyfITkD6TErTyrgX2ENSfl/2oE5jNYCkrwHvIbvTadUA2+6lWOaTcsNPkZ2I7Jb7XmwW2U0hfesO9tkfMZwY2tdY4M/A02Rf4n/qmyFpfUlHSdosIl4BVgF9Z+HfB05IZ5yS9IbUkDkWuJ3soPAFSWMkTSWr3rkkncVfDpwhaWNJOwOD3Ro4liyZrABGSzqdtc8G18U8YBpZ1dtjZA2Ph5G1afyuzjrLyNoghmosqepP0tbAyQ2seylwjLIG/o2B0/vNXwh8KJXvZODYQeL4Y0T8SdJ+ZCcKjfhhiuUgSetJ2jq9p31xfCx9BqYAH8mttwJ4jYHL8GLgy5LGK7tZ4nSyK6sy5pJ9Vv9nagz+XJperwH9i2T7fkjUuPW1n0uBT0naNZX/V/tmpKvk7wPfljQhbXtrSYfm1h3ovRtRnBja1/mkBlfg92SNfXmfBBan6o4TyM5siYj5ZI26/0pW3/8gWV02EfEy8AGys6+nyNo7jo6I+9M2P0d2dfIkWePbxWTJqZ7rye6y+kOK9U8M4yV4RPyB7CB9SxpfBTwM/CpXHdXfD4FdU3XBz4bwsl8ja9d4FvhPsmRZNt5rge8CN5GV+61pVl8Zfht4mSx5zSarW6/nRODvJT1HdpC6dIBla8XyG7K2mG+T7cs81pzlf4XsauIZsv29KLfei8A/Ar9KZfi2Gps/E5gP3AXcDfw2TSsT18tk7TZHk7W7fZqsqvVlgHTCk796+CeyM/9F6S6p5yWdVmfb15K1xdxIVv79k80pafpt6XvzX8BOuXUHeu9GFNWvrrORTtLXgTdFxIxWx9KJJO0C3ANsEBGrWx2PlTfS3ztfMdjrJO0saY9UBbUfWVXHFa2Oq5NI+mCq6hsHfB24eiQeWDqR37s1nBgsbyxZ1ckLZFUX/5fsKVUr7zNk9fQPkdWlf7a14VgD/N4lrkoyM7MCXzGYmVlBO3W+VdeWW24ZPT09rQ7DzKyjLFiw4KmIGD/4kkUdkRh6enqYP39+q8MwM+sokh4dfKm1uSrJzMwKnBjMzKzAicHMzAqcGMzMrMCJwczMCpwYzMyswInBzMwKnBjMzKzAicHMzAqcGMw6xNSpU5k6dWqrw7ARwInBzMwKnBjMzKzAicHMzAqcGMzMrMCJwczMCpwYzMyswInBzMwKnBjMzKzAicHMzAqcGMzMrMCJwczMCpwYzMyswInBzMwKnBjMzKzAicHMzAqcGMzMrMCJwczMCpwYzMysYHSVG5e0GHgOeBVYHRFTJG0B/BjoARYDR0bEM1XGYWZm5TXjimFaROwVEVPS+KnAnIjYAZiTxs3MrE20oippOjA7Dc8GjmhBDGZmVkfViSGAX0haIGlmmrZVRPQCpP8Taq0oaaak+ZLmr1ixouIwzcysT6VtDMA7I+IJSROAGyTdX3bFiJgFzAKYMmVKVBWgmZkVVXrFEBFPpP/LgSuA/YBlkiYCpP/Lq4zBzMwaU1likPQGSWP7hoF3A/cAVwEz0mIzgCurisHMzBpXZVXSVsAVkvpe56KIuE7SHcClko4FlgAfrTAGs451+LR30bt0yevjix7vBWDvyduutezEbSbx85vmNS02626VJYaIeBjYs8b0p4GDqnpds27Ru3QJl01bkwQ+fs1KAC6atnZi+PBNjzQtLut+fvLZzMwKnBjMzKzAicHMzAqcGMzMrMCJwczMCpwYzMyswInBzMwKnBjMzKzAicHMzAqcGMzMrMCJwczMCpwYzMyswInBzMwKnBjMzKzAicHMzAqcGMyG2dSpU5k6dWqrwzAbMicGMzMrcGIwM7MCJwYzMytwYjAzswInBjMzK3BiMDOzAicGMzMrGN3qAMysnIvet3erQ2gbfc+JzJ07t6VxdCtfMZiZWYETg5mZFTgxmJlZgRODWQlV9X/kfpWsHVWeGCSNkvQ7Sdek8S0k3SBpUfo/ruoYzMysvGZcMfwNcF9u/FRgTkTsAMxJ42Zm1iYqTQyS3gK8F/hBbvJ0YHYang0cUWUMZmbWmKqvGM4GvgC8lpu2VUT0AqT/E2qtKGmmpPmS5q9YsaLiMM3MrE9liUHS+4DlEbFgKOtHxKyImBIRU8aPHz/M0ZmZWT1VPvn8TuADkg4HNgQ2lXQBsEzSxIjolTQRWF5hDGZm1qDKrhgi4osR8ZaI6AE+BtwYEZ8ArgJmpMVmAFdWFYOZmTWuFc8xnAUcImkRcEgaNzOzNtGUTvQiYi4wNw0/DRzUjNft4w63zMzK85PPZmZW4G63zdbB4dPeRe/SJYVpix7vBWDvyduutfzEbSbx85vmNSU2s6FyYjBbB71Ll3DZtGIC+Pg1KwG4aNraieHDNz3SlLjM1oWrktpQt3es1u37Z9bpnBjMzKzAicHMzAqcGMzMrMCJwdpale0Rbuswq82JwczMCny7qpm1vf7Pi/hZkWo5MZhZ2+v/vIifFamWq5LMrCFum+l+vmIw66eRbi6WLXsSWPustZ7eJ58sbGOgKpFGt202XJwYzPpppJuLt1/4eEPbjtdeLV0l0ui2zYaLq5LMzKzAicHMzAqcGMzMrMCJwczMChpODJLGSdqjimDMzKz1SiUGSXMlbSppC+BO4EeSvlVtaAPzvdQZl4PZwNrlO9IucZRR9nbVzSJilaTjgB9FxFcl3VVlYEPVyD3ofnTezGxtZRPDaEkTgSOBL1UYzzpr5B50PzpvZra2sm0MXwOuBx6MiDskbQcsqi4sMzNrlbJXDL0R8XqDc0Q83Oo2hk7TV7c4d+7cteZ1Q8+R3b5/ZiNJ2cTwL8A+JabZEHR7z5Hdvn/9XfS+vVsdgtk6GTAxSHo78A5gvKT/nZu1KTCqysDMzKw1BrtiWB/YJC03Njd9FfCRqoIyM7PWGTAxRMQ8YJ6k8yLi0SbFZCNYI91SN9IeUWVX2u1ooDYfGz7tUM5VxFC2jWEDSbOAnvw6EXFgvRUkbQjcDGyQ1vlpev5hC+DHaVuLgSMj4pmhBG/dp5FuqRtpj6iyK22zblM2MfwEOAf4AfBqyXX+DBwYEc9LGgP8UtK1wIeAORFxlqRTgVOBUxqM28zMKlI2MayOiH9vZMMREcDzaXRM+gtgOjA1TZ8NzMWJwcysbZRNDFdLOhG4guxKAICI+ONAK0kaBSwAJgPfi4jbJW0VEb1p/V5JE+qsOxOYCTBp0qSSYdpg2qFO1DpHM7uYaeSz2e23BLf6e1o2McxI/0/OTQtgu4FWiohXgb0kbQ5cIWn3soFFxCxgFsCUKVOi7HpmNnzcxczIVCoxRMQ63aIRESslzQUOA5ZJmpiuFiYCy9dl22ZmNrxKJQZJR9eaHhHnD7DOeOCVlBQ2Ag4Gvg5cRXYFclb6f2WjQZs1W1VVF91eJWKdqWxV0r654Q2Bg4DfAnUTAzARmJ3aGdYDLo2IayTdClwq6VhgCfDRxsM2s+HU6jptay9lq5I+nx+XtBnwH4Oscxew1ulQRDxNlljMzKwNDfU3n18EdhjOQMzMrD2UbWO4muwuJMg6z9sFuLSqoPpb9MADa90a183dGQykqi4jzLpFO3xHOv2XJMu2MXwzN7waeDQiHqsgnppeeeVld2eQVNVlhFm3aIfvSKff5luqKil1pnc/WQ+r44CXqwzKzMxap2xV0pHAN8i6rxDwL5JOjoifVhjbsPEtgWZm5ZWtSvoSsG9ELIfXn1H4L6AjEkMrNPJzlp3WLtJI/Sl03v6ZjXRlE8N6fUkheZqh39E0IjTyc5ad1i7SSP0pdN7+2cjmZzrKJ4brJF0PXJzG/wfw82pCsm6v+ur2/TPrdIP95vNkYKuIOFnSh4C/IGtjuBW4sAnxmZlZkw12xXA2cBpARFwOXA4gaUqa9/5KozOzUtrh3v1u1y4/D9tI++VQDZYYelLXFgURMV9Sz7BF0aFcF2ntoh3u3e927fLzsI20X27/0OIhvcZgiWHDAeZtNKRXNGtAle0Rbuswq22wxHCHpOMj4vv5ialn1AXVhWXDoX/1AnTWY/nWvXzL8xqNfE+hOWUxWGL4W7JfXjuKNYlgCrA+8MEqA7N11796ATrrsXzrXr7leY1GvqfQnLIYMDFExDLgHZKmAX0/y/mfEXFj5ZFZy3VaG0qnxdupXAXX/cr+HsNNwE0Vx2JmZm2g7ANuZtYlGqnT7ua6favPicFshGmkTrub6/atPvd3ZGZmBU4MZmZW4MRgZmYFHdvG4FvmMi4Hs4G1y3ekXeIow1cMZmZW0LFXDGbW/tr9oUPfulubE0OTdNJl5FB0+/5Zd/Ktu7W5KsnMzAqcGMzMrMBVSf3UqxN1N8FmNlJUlhgkbQOcD7wJeA2YFRHfkbQF8GOgB1gMHBkRz1QVx3BxN8HWSdzm0xztUM5VxFBlVdJq4KSI2AV4G/DXknYFTgXmRMQOwJw0bmZmbaKyxBARvRHx2zT8HHAfsDUwHZidFpsNHFFVDGZm1rimtDFI6gH2Bm4HtoqIXsiSh6QJddaZCcwEGDN6VDPCHBHqXXZ24v3c/WN2m09na4dqmXbR6rKoPDFI2gS4DPjbiFglqdR6ETELmAWw8YYbRHURGnTm/dz9Y3abj9nwqPR2VUljyJLChRFxeZq8TNLENH8isLzKGMzMrDFV3pUk4IfAfRHxrdysq4AZwFnp/5VVxWBm5bS66sLaS5VVSe8EPgncLWlhmnYaWUK4VNKxwBLgoxXGYGZmDaosMUTEL4F6DQoHVfW6Zma2btwlhpmZFbhLDDMbFp14y7PV5sRgZsOiE295ttpclWRmZgW+YjAzy/GtuyM8MbhO1MxsbSM6MbhO1MxsbSM6MawrX3KaWTdy47OZmRU4MZiZWYGrkqyuTqsq67R4zdqVrxjMzKzAicHMzAqcGMzMrMCJwczMCpwYzMyswInBzMwKfLuqmVXGtxB3Jl8xmJlZgRODmZkVuCqpH1/6mtlI5ysGMzMrcGIwM7MCJwYzMytwYjAzswInBjMzK3BiMDOzAicGMzMrqCwxSDpX0nJJ9+SmbSHpBkmL0v9xVb2+mZkNTZVXDOcBh/WbdiowJyJ2AOakcTMzayOVJYaIuBn4Y7/J04HZaXg2cERVr29mZkPT7C4xtoqIXoCI6JU0od6CkmYCMwHGjB7VpPDMzKxtG58jYlZETImIKaNHOTGYmTVLsxPDMkkTAdL/5U1+fTMzG0SzE8NVwIw0PAO4ssmvb2Zmg6jydtWLgVuBnSQ9JulY4CzgEEmLgEPSuJmZtZHKGp8j4q/qzDqoqtc0M7N117aNz2Zm1hpODGZmVuDEYGZmBU4MZmZW4MRgZmYFTgxmZlbgxGBmZgVODGZmVuDEYGZmBU4MZmZW4MRgZmYFTgxmZlbgxGBmZgVODGZmVuDEYGZmBU4MZmZW4MRgZmYFTgxmZlbgxGBmZgVODGZmVuDEYGZmBU4MZmZW4MRgZmYFTgxmZlbgxGBmZgVODGZmVuDEYGZmBU4MZmZW0JLEIOkwSQ9IelDSqa2IwczMamt6YpA0Cvge8B5gV+CvJO3a7DjMzKy2Vlwx7Ac8GBEPR8TLwCXA9BbEYWZmNSgimvuC0keAwyLiuDT+SWD/iPhcv+VmAjPT6O7APU0NtH1tCTzV6iDahMtiDZfFGi6LNXaKiLGNrjS6ikgGoRrT1spOETELmAUgaX5ETKk6sE7gsljDZbGGy2INl8UakuYPZb1WVCU9BmyTG38L8EQL4jAzsxpakRjuAHaQtK2k9YGPAVe1IA4zM6uh6VVJEbFa0ueA64FRwLkRce8gq82qPrKO4bJYw2WxhstiDZfFGkMqi6Y3PpuZWXvzk89mZlbgxGBmZgVtlRgG6ypDme+m+XdJ2qcVcTZDibI4KpXBXZJ+LWnPVsRZtbLdp0jaV9Kr6TmZrlSmLCRNlbRQ0r2S5jU7xmYp8f3YTNLVku5MZXFMK+JsBknnSlouqeazXkM6bkZEW/yRNUQ/BGwHrA/cCezab5nDgWvJnoV4G3B7q+NuYVm8AxiXht/TjWVRphxyy90I/Bz4SKvjbuFnYnPg98CkND6h1XG3sCxOA76ehscDfwTWb3XsFZXHXwL7APfUmd/wcbOdrhjKdJUxHTg/MrcBm0ua2OxAm2DQsoiIX0fEM2n0NrLnQbpN2e5TPg9cBixvZnBNVqYsPg5cHhFLACKiW8ujTFkEMFaSgE3IEsPq5obZHBFxM9n+1dPwcbOdEsPWwNLc+GNpWqPLdING9/NYsjOCbjNoOUjaGvggcE4T42qFMp+JHYFxkuZKWiDp6KZF11xlyuJfgV3IHp69G/ibiHitOeG1nYaPm63oEqOeMl1llOpOowuU3k9J08gSw19UGlFrlCmHs4FTIuLV7OSwa5Upi9HAfwcOAjYCbpV0W0T8oergmqxMWRwKLAQOBLYHbpB0S0Ssqjq4NtTwcbOdEkOZrjJGSncapfZT0h7AD4D3RMTTTYqtmcqUwxTgkpQUtgQOl7Q6In7WnBCbpuz346mIeAF4QdLNwJ5AtyWGMmVxDHBWZJXsD0p6BNgZ+E1zQmwrDR8326kqqUxXGVcBR6dW9rcBz0ZEb7MDbYJBy0LSJOBy4JNdeEbYZ9ByiIhtI6InInqAnwIndmFSgHLfjyuBAySNlrQxsD9wX5PjbIYyZbGE7MoJSVsBOwEPNzXK9tHwcbNtrhiiTlcZkk5I888hu+vkcOBB4EWys4KuU7IsTgfeCPxbOlteHV3Wo2TJchgRypRFRNwn6TrgLuA14AcR0XXd1Zf8XPwDcJ6ku8mqUk6JiK7silvSxcBUYEtJjwFfBcbA0I+b7hLDzMwK2qkqyczM2oATg5mZFTgxmJlZgRODmZkVODGYmVmBE4MZkHpmXSjpHkk/kbSxpJ56PVaadTMnBrPMSxGxV0TsDrwMnNDqgMxaxYnBbG23AJPT8ChJ3099+v9C0kYAko6XdEfq7/+y9KQxkj6arjruTF1SIGmUpG+k5e+S9JnW7JZZOU4MZjmSRpP9vsXdadIOwPciYjdgJfDhNP3yiNg3IvYk63bi2DT9dODQNP0DadqxZN0Q7AvsCxwvadvq98ZsaJwYzDIbSVoIzCfrZ+eHafojEbEwDS8AetLw7pJuSV0uHAXslqb/iqwrhuPJumsAeDdZXzULgdvJujLZocqdMVsXbdNXklmLvRQRe+UnpD6o/pyb9CpZd9YA5wFHRMSdkj5F1lcNEXGCpP2B9wILJe1F1lfP5yPi+ip3wGy4+IrBbGjGAr2SxpBdMQAgafuIuD0iTgeeIuvu+Hrgs2lZJO0o6Q2tCNqsDF8xmA3NV8iqhR4la48Ym6Z/Q9IOZFcJc8h+j/gusiqo36afmlwBHNHsgM3Kcu+qZmZW4KokMzMrcGIwM7MCJwYzMytwYjAzswInBjMzK3BiMDOzAicGMzMr+P9AxE8y9/wSuAAAAABJRU5ErkJggg==\n",
"text/plain": [
"