{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Source Injector" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The source injector can produce mock simulated data independent of the MEGAlib software.\n", "\n", "Standard data simulation requires the users to install and use MEGAlib to convolve the source model with the detector effects to generate data. The source injector utilizes the response generated by intensive simulation, which contains the statistical detector effects. With the source injector, you can convolve response, source model, and orientation to gain the mock data quickly. \n", "\n", "The advantages of using the source injector include:\n", "\n", "- No need to install and use MEGAlib\n", "- Get the data much faster than the standard simulation pipeline\n", "- The generated data are in the format that can be used for spectral fitting, localization, imaging, etc.\n", "\n", "The disadvantages are:\n", "\n", "- The data are binned based on the binning of the response, which means that you lost the unbinned event distribution as you will get from the MEGAlib pipeline.\n", "- If the response is coarse, the data you generated might not be precise." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%%capture\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from scipy.interpolate import interp1d\n", "import astropy.units as u\n", "from pathlib import Path\n", "from astropy.coordinates import SkyCoord\n", "from astromodels.functions.function import Function1D, FunctionMeta\n", "from cosipy import SpacecraftFile, SourceInjector\n", "from histpy import Histogram\n", "from threeML import Powerlaw, Band, Model, PointSource\n", "from cosipy.threeml.custom_functions import SpecFromDat\n", "from cosipy.util import fetch_wasabi_file\n", "import shutil\n", "import os\n", "import h5py as h5\n", "\n", "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "data_dir = Path(\"\") # Current directory by default. Modify if you want a different path" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Get the data\n", "The data can be downloaded by running the cells below. Each respective cell also gives the wasabi file path and file size. " ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "%%capture\n", "zipped_response_path = data_dir/\"SMEXv12.Continuum.HEALPixO3_10bins_log_flat.binnedimaging.imagingresponse.nonsparse_nside8.area.good_chunks_unzip.earthocc.zip\"\n", "response_path = data_dir/\"SMEXv12.Continuum.HEALPixO3_10bins_log_flat.binnedimaging.imagingresponse.nonsparse_nside8.area.good_chunks_unzip.earthocc.h5\"\n", "\n", "# download response file ~839.62 MB\n", "if not response_path.exists():\n", " \n", " fetch_wasabi_file(\"COSI-SMEX/DC2/Responses/SMEXv12.Continuum.HEALPixO3_10bins_log_flat.binnedimaging.imagingresponse.nonsparse_nside8.area.good_chunks_unzip.earthocc.zip\", zipped_response_path)\n", "\n", " # unzip the response file\n", " shutil.unpack_archive(zipped_response_path)\n", " \n", " # delete the zipped response to save space\n", " os.remove(zipped_response_path)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "%%capture\n", "orientation_path = data_dir/\"20280301_3_month_with_orbital_info.ori\"\n", "\n", "# download orientation file ~684.38 MB\n", "if not orientation_path.exists():\n", " fetch_wasabi_file(\"COSI-SMEX/DC2/Data/Orientation/20280301_3_month_with_orbital_info.ori\", orientation_path)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Inject a source response" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Method 1 : Define the point source" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this method, we are setting up an analytical function (eg: a power law model) to simulate the spectral characteristics of a point source:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
01:12:02 WARNING The current value of the parameter beta (-2.0) was above the new maximum -2.15. parameter.py:794\n", "\n" ], "text/plain": [ "\u001b[38;5;46m01:12:02\u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m The current value of the parameter beta \u001b[0m\u001b[1;38;5;251m(\u001b[0m\u001b[1;37m-2.0\u001b[0m\u001b[1;38;5;251m)\u001b[0m\u001b[1;38;5;251m was above the new maximum \u001b[0m\u001b[1;37m-2.15\u001b[0m\u001b[1;38;5;251m.\u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=475164;file:///Users/shengyong/conda_envs/source_injector/lib/python3.10/site-packages/astromodels/core/parameter.py\u001b\\\u001b[2mparameter.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=190355;file:///Users/shengyong/conda_envs/source_injector/lib/python3.10/site-packages/astromodels/core/parameter.py#794\u001b\\\u001b[2m794\u001b[0m\u001b]8;;\u001b\\\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Defind the Crab spectrum\n", "\n", "alpha_inj = -1.99\n", "beta_inj = -2.32\n", "E0_inj = 531. * (alpha_inj - beta_inj) * u.keV\n", "xp_inj = E0_inj * (alpha_inj + 2) / (alpha_inj - beta_inj)\n", "piv_inj = 100. * u.keV\n", "K_inj = 7.56e-4 / u.cm / u.cm / u.s / u.keV\n", "\n", "spectrum_inj = Band()\n", "\n", "spectrum_inj.alpha.min_value = -2.14\n", "spectrum_inj.alpha.max_value = 3.0\n", "spectrum_inj.beta.min_value = -5.0\n", "spectrum_inj.beta.max_value = -2.15\n", "spectrum_inj.xp.min_value = 1.0\n", "\n", "spectrum_inj.alpha.value = alpha_inj\n", "spectrum_inj.beta.value = beta_inj\n", "spectrum_inj.xp.value = xp_inj.value\n", "spectrum_inj.K.value = K_inj.value\n", "spectrum_inj.piv.value = piv_inj.value\n", "\n", "spectrum_inj.xp.unit = xp_inj.unit\n", "spectrum_inj.K.unit = K_inj.unit\n", "spectrum_inj.piv.unit = piv_inj.unit" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "# Define the coordinate of the point source\n", "source_coord = SkyCoord(l = 184.5551, b = -05.7877, frame = \"galactic\", unit = \"deg\")\n", "\n", "# define the Crab point source\n", "point_source = PointSource('Crab', l = source_coord.l.deg, b = source_coord.b.deg, spectral_shape=spectrum_inj)\n", "\n", "# define the model. The model can contain multiple point sources.\n", "model = Model(point_source)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Read orientation file" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "# Read the 3-month orientation\n", "# It is the pointing of the spacecraft during the the mock simlulation\n", "ori = SpacecraftFile.parse_from_file(orientation_path)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Get the expected counts and save to a data file" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "# Define an injector by the response\n", "injector = SourceInjector(response_path = response_path)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 21.4 s, sys: 3.73 s, total: 25.1 s\n", "Wall time: 26.4 s\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkQAAAG7CAYAAAAizIoLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAqlElEQVR4nO3dfXRU9YH/8c8QkgGZQTCxIU0QhXiQVKtdC4isSnSgyUGyWuJDy+6Rg0YgW/EIFQo17rYJBVoOoujGhdWNa2slLOKhHB7cEZCH+rRbkIVABKHhoUkwkcAMYIYm9/cHv8waE+JkMpM7w/f9OsfjzL3fufMZ2ls+/d4nh2VZlgAAAAzWw+4AAAAAdqMQAQAA41GIAACA8ShEAADAeBQiAABgPAoRAAAwHoUIAAAYj0IUgi+//FKVlZX68ssv7Y4CAACigEIUgqqqKhUUFKiqqsruKAAAIAooRAAAwHg97Q4Qy7xer7xer/x+v91RAABAFFGIOuDxeOTxeFRZWamCggK74wAAgCjhkBkAADAehQgAABiPQgQAAIxHIQIAAMajEAEAAONRiAAAgPEoRAAAwHgUIgAAYDwKEQAAMB6FCAAAGI9CBAAAjEchAgAAxuPhrui0VftW6dmtz8rX6LM7Sqe4nW4VZxcrPyvf7igAgBhDIUKnPbv1WR2oO2B3jM7zSUVbiihEAIA2KETotJaZoR6OHkpzpdmcJjTV/mo1W81xN6sFAOgeFCKELc2VpuMzj9sdIyQZSzJ0wnfC7hgAgBjFSdUAAMB4FCIAAGA8ChEAADAe5xB1wOv1yuv1yu/32x0FAABEEYWoAx6PRx6PR5WVlSooKLA7DgAAiBIOmQEAAONRiAAAgPEoRAAAwHgUIgAAYDwKEQAAMB6FCAAAGI9CBAAAjEchAgAAxuPGjDBKtb9aGUsy7I4RMrfTreLsYuVn5dsdBQAuaxQiGMHtdEs+qdlq1gnfCbvjhM4nFW0pohABQJRRiGCE4uxiFW0pkq/RZ3eUkFX7q9VsNcdVZgCIVxQiGCE/Kz/uZlkylmTE12wWAMQxTqoGAADGoxABAADjUYgAAIDxKEQAAMB4FCIAAGA8ChEAADAehQgAABiPQgQAAIzHjRlttmrfKj279dm4uhtxtb/a7ggAAEQUhchmz259VgfqDtgdIyxup9vuCAAARASFyGYtM0M9HD2U5kqzOU3oWp7CDgDA5YBCFCPSXGk6PvO43TEAADASJ1UDAADjUYgAAIDxKEQAAMB4Rp1D9MYbb2j16tXy+/3KyMjQsmXLdMUVV9gdCwAA2MyYQvTWW2/pww8/1L/8y7/oW9/6lg4fPqyePY35+QAAoANGNIKmpia9/vrrevHFF5WamipJGjJkiM2pAABArIjJQnTu3Dm9+eabqqio0P79++Xz+TR37lzl5ua2GRsIBPTKK6/onXfekc/n05AhQ/TYY49p+PDhwTGff/65GhsbtXXrVpWXl8vlcunhhx/WhAkTuvNnAQCAGBWTJ1WfPn1aZWVlqqqqUmZmZodjFyxYoPLyco0dO1YzZsxQjx49NHv2bO3Zsyc45vPPP5ff79exY8dUXl6uX/7yl1q+fLk++eSTaP8UAAAQB2KyECUnJ2vNmjVatWqVpk+ffslxFRUVevfdd/X444+rsLBQeXl5Wrp0qQYMGKDS0tLgOKfTKUmaPHmynE6nhgwZonvuuUcffPBB1H8LAACIfTFZiJKSkpScnPyN49577z0lJCQoLy8vuMzpdGr8+PHat2+famtrJUkDBw5UYmKiHA5HcNxXXwMAALPFZCEK1cGDB5WRkaE+ffq0Wj5s2DBJ0qFDhyRJvXv31l133aX/+I//UCAQ0J///Gdt3rxZt912W7vbraurU2VlZfCfqqqq6P4QAABgq5g8qTpU9fX17c4ktSyrq6sLLnvqqae0aNEiTZgwQVdeeaUeffRR3Xzzze1ud+3atSorK4tKZgAAEHviuhA1NjYqMTGxzfKkpKTg+hZut1slJSUhbTcvL0+jR48Ovq+qqgr5swAAIP7EdSFyOp26cOFCm+WBQCC4PhwpKSlKSUnpUjYAABA/4vocouTkZNXX17dZ3rKMUgMAAEIR14UoMzNTx48f19mzZ1str6ioCK4HAAD4JnFdiMaMGaOmpiatXbs2uCwQCGj9+vXKysoKPqYDAACgIzF7DlHLU+lbDn/t3LlTJ0+elCRNnDhRLpdLWVlZys7O1vLly9XQ0KD09HRt3LhRNTU1mjNnTpczeL1eeb1e+f3+Lm8LAADErpgtRCtXrlRNTU3w/bZt27Rt2zZJ0rhx4+RyuSRJ8+bNU2pqqjZt2iS/36/Bgwdr0aJFuuWWW7qcwePxyOPxqLKyUgUFBV3eHgAAiE0xW4jKy8tDGud0OlVYWKjCwsIoJwIAAJeruD6HCAAAIBJidoYIwEXV/mplLMmwO0bI3E63irOLlZ+Vb3cUAAgZhQiIUW6nW/JJzVazTvhO2B0ndD6paEsRhQhAXKEQdYCrzGCn4uxiFW0pkq/RZ3eUkFX7q9VsNcdVZgCQKEQd4ioz2Ck/Kz/uZlkylmTE12wWAPx/nFQNAACMRyECAADGoxABAADjUYgAAIDxKEQAAMB4XGXWAS67BwDADBSiDnDZPQAAZuCQGQAAMB6FCAAAGI9CBAAAjEchAgAAxqMQAQAA41GIAACA8bjsvgPchwgAADNQiDrAfYgAADADh8wAAIDxKEQAAMB4FCIAAGA8ChEAADAehQgAABiPQgQAAIxHIQIAAMajEAEAAONxY8YOcKdqAADMQCHqAHeqBgDADBwyAwAAxqMQAQAA41GIAACA8ShEAADAeBQiAABgPAoRAAAwHoUIAAAYj0IEAACMRyECAADGoxABAADj8eiODvAsMwAAzEAh6gDPMgMAwAwUIgARV+2vVsaSDLtjhMztdKs4u1j5Wfl2RwFgEwoRgIhxO92ST2q2mnXCd8LuOKHzSUVbiihEgMEoRAAipji7WEVbiuRr9NkdJWTV/mo1W81xlRlA5FGIAERMflZ+3M2yZCzJiK/ZLABRwWX3AADAeBQiAABgPAoRAAAwHoUIAAAYj0IEAACMRyECAADGoxABAADjUYgAAIDxKEQAAMB43Km6A16vV16vV36/3+4oAAAgiihEHfB4PPJ4PKqsrFRBQYHdcQAAQJRwyAwAABiPQgQAAIxHIQIAAMajEAEAAONRiAAAgPEoRAAAwHgUIgAAYDwKEQAAMB6FCAAAGI9CBAAAjEchAgAAxqMQAQAA41GIAACA8ShEAADAeBQiAABgPAoRAAAwHoUIAAAYj0IEAACMRyECAADG62l3gFjm9Xrl9Xrl9/vtjgIAAKKIQtQBj8cjj8ejyspKFRQU2B0HAABECYfMAACA8ShEAADAeBQiAABgvIifQ7RhwwZ9+umncrvdGj9+vFJTUyP9FQAAABEVdiF6++239dprr8nhcOill15SWlqa/vmf/1lbt25tNWbFihWUIgAAENPCPmS2e/duffHFF0pKSlJaWppOnDihLVu2SJIsy5JlWTpz5ozefPPNiIUFAACIhrBniA4dOiSHw6GbbrpJkrRr1y5JUkJCgkaOHKkPPvhATU1N+p//+Z/IJAWAKKr2VytjSYbdMULmdrpVnF2s/Kx8u6MAl4WwC9GpU6ckSd/61rckSUeOHJEk/e3f/q1+8YtfaNGiRVq/fr1qa2sjEBMAosPtdEs+qdlq1gnfCbvjhM4nFW0pohABERJ2ITp//rwkqXfv3pKk48ePy+FwaPDgwZKk9PR0SVJTU1NXMwJA1BRnF6toS5F8jT67o4Ss2l+tZqs5rjIDsS7sQtSnTx/5fD7t3r1b48ePV0VFhSQpI+PilLPPd3FH7du3bwRiAkB05Gflx90sS8aSjPiazQLiQNgnVbfMBP33f/+37r//fp05c0aSNGzYMEnS559/LklKSUnpakYAAICoCrsQ5eTkyLIsSQr+OysrS9/+9rdlWZZ27dolh8OhG264ITJJAQAAoiTsQ2a5ubmqrq7W2rVrFQgE9N3vflczZ86UJB0+fFh9+vRRnz59NHz48IiFBQAAiIYu3al6ypQpmjJlSpvlQ4YM0W9/+9uubBoAAKDbhF2IysrKJEnDhw/Xd77znTbrGxoagpfcDx06NNyvAQAAiLqwC9G///u/y+FwqHfv3u0Woo0bN+rll1+Ww+EI3sEaAAAgFkX84a4tWh7fAQAAEOvCvsqsI42Njdq7d280Ng0AABBxnZohGjNmTKv3lmWptLRUpaWll/zMFVdcEVYwAACA7tKpQmRZlhwOR6tDYZc6LOZwOORwONo9vwgAACCWdPqQWajnBVmWpauuukrTpk3rdCgAAIDu1KkZop/97GfB1wsXLpTD4dDdd9/d5uaLPXv21NVXX63vfOc7SkxMjExSAACAKOlUIcrNzQ2+XrhwoSzL0tChQ1stBwAAiDdhX3a/cuVKSTzNHgAAxL+wC9GAAQOCr8+fPy+fz3fJ84tSU1PD/RoAAICo69KNGTds2KA33nhDx44du+QY7lQNAABiXdiFaN26dVq8eLGk0K88AwAAiEVdOoeo5b5E8WDGjBmqqKhQQkKCJOm73/2ufvOb39icCgAAxIKwC1F1dbUcDoecTqcKCgp0zTXXKDExMaYL0uzZszVu3Di7YwAAgBgTdiFyu906deqU/u7v/k75+fmRzAQAANCtwn646+233y7LstTY2BjJPJKkc+fO6dVXX9VPf/pTjR8/Xnfeeac2bNjQ7thAIKDS0lLdf//98ng8mjp1qj7++ON2xy5btkwTJkzQzJkz9dlnn0U8NwAAiE9hF6IpU6YoOTlZGzZs0B//+MdIZtLp06dVVlamqqoqZWZmdjh2wYIFKi8v19ixYzVjxgz16NFDs2fP1p49e1qNmzZtmlauXKn//M//1Pe//309/fTTOnfuXERzAwCA+BT2IbNf/vKXcjqdqq+v17x585ScnKy0tLTgScstHA6Hli5d2qltJycna82aNUpOTtaBAwf0+OOPtzuuoqJC7777rqZPn64f/ehHkqQf/OAHmjx5skpLS1VaWhocm5WVFXz94x//WOvXr9e+ffvaPHYEAACYJ+xCtHv37uAT7S3LUl1dnerr61uNCfcqtKSkJCUnJ3/juPfee08JCQnKy8sLLnM6nRo/fryWL1+u2traS94UsiU3AABAl27M+PVC0d0F4+DBg8rIyFCfPn1aLR82bJgk6dChQ0pNTZXP59OBAwd08803y+FwaM2aNfL5fK1mjb7q6+Wuqqoqej8CAADYLuxClJOTE8kcYamvr293JqllWV1dnSSpqalJy5cv19GjR9WzZ09lZmZq0aJFcrlc7W537dq1Kisri1puAAAQW8IuRHPnzo1kjrA0NjYqMTGxzfKkpKTgeknq16+fVqxYEfJ28/LyNHr06OD7qqoqlZSUdDEtAACIVV06ZGY3p9OpCxcutFkeCASC68ORkpKilJSULmUDAADxI+xCVFtbG/LYaD3tPjk5WZ9//nmb5S3n/1BqAABAKMIuRA8++GBIV5BF82n3mZmZ2rVrl86ePdvqxOqKiorgegAAgG8S9o0ZW1iW9Y3/RMuYMWPU1NSktWvXBpcFAgGtX79eWVlZXZ6Z8nq9+tnPfqZly5Z1NSoAAIhhEb3svkXLzFFXytDq1avl9/uDh7927typkydPSpImTpwol8ulrKwsZWdna/ny5WpoaFB6ero2btyompoazZkzJ+zvbuHxeOTxeFRZWamCgoIubw8AAMSmsAvR888/32ZZIBDQ8ePHtWbNGh07dkyjRo3Sww8/HNb2V65cqZqamuD7bdu2adu2bZKkcePGBS+ZnzdvnlJTU7Vp0yb5/X4NHjxYixYt0i233BLW9wIAAPOEXYguVThGjBihcePG6ZFHHtEHH3ygCRMmhLX98vLykMY5nU4VFhaqsLAwrO8BAADo8jlE7XG5XLrppptkWZZ+97vfReMrAAAAIiYqhejs2bPav3+/pIuPzwAAAIhlYR8ye/LJJ9td3tjYqKNHj+rcuXOS/u+u0QAAALGqy0+7b0/LU+4dDoeGDx8edji7eb1eeb1e+f1+u6MAAIAoispl9y3rMjIy4vpkZy67BxDLqv3VyliSYXeMkLmdbhVnFys/K9/uKEAbEX/avcPhkMvl0rBhw3THHXdwyAwAIsztdEs+qdlq1gnfCbvjhM4nFW0pohAhJsX10+4BwETF2cUq2lIkX6PP7ighq/ZXq9lqjqvMMEtcP+0eAEyUn5Ufd7MsGUsy4ms2C8bpciH65JNPVF5ern379snv98vlcunGG2/UAw88oJtvvjkSGQEAAKKqS4WovLxcpaWlrR7ieurUKe3YsUM7d+7U9OnT9eCDD0YkKAAAQLSEXYj279+v0tJSNTc3t3v5fXNzs0pLS3XTTTdp2LBhXQppFy67BwDADGHfqXr16tXBMtSrVy+NGTNGP/zhDzVmzBj16tVL0sVL7996662Ihe1uHo9HCxcu1BNPPGF3FAAAEEVhzxDt2bNHkpSamqp//dd/Vb9+/YLrTp06palTp6q2tlaffPJJl0MCAABEU9gzRF988YUcDoc8Hk+rMiRJ/fv319ixY4PjAAAAYlnYhSghIUHSxQe5tqdlecs4AACAWBV2IUpLS5NlWdqwYYM++uijVus++ugjrV+/Xg6HQ2lpaV0OCQAAEE1hn0M0YsQIHTlyRI2NjZo9e7b69eun/v3769SpU2poaAg+4HXkyJGRzAsAABBxYc8QPfjgg3K73ZIuXk126tQpHTlyRKdOnQrek8jlcumBBx6ITFIAAIAoCXuGKCUlRSUlJSoqKtKZM2ckKViEJKlv374qKSlRSkpK11PahPsQAQBghi7dqfqWW27R73//e23YsEH79u3TmTNn1LdvX914443KycmRy+WKVE5beDweeTweVVZWqqCgwO44AAAgSrr8LLOWw2IcGgMAAPGq04WotrZWkuR0Otvcf6hFQ0ODGhsbJV28cSMAAEAs69RJ1R9++KEeeughPfTQQ9qxY8clx+3YsUMPPfSQHn74Ye5UDQAAYl6nCtHWrVtlWZb69++v3NzcS47LyclRv379ZFmWNm/e3OWQAAAA0dSpQlRRUSGHw6ERI0Z0eAfqnj17asSIEbIsS7t37+5qRgAAgKjqVCGqqamRJKWnp3/j2JYxLZ8BAACIVZ0qRH/9618ltb7fUKifAQAAiFWdusrM7XaroaFBn3766TeObRnTcjfreMSNGQEAMEOnZoiuueYaWZal999/X5999tklx3322Wd6//335XA4NHDgwC6HtIvH49HChQv1xBNP2B0FAABEUacK0a233ipJampq0qxZs7R9+/Y2Y3bs2KGnn35aTU1NkqTvf//7EYgJAAAQPZ06ZJaXl6ff/va3unDhgk6dOqWioiK53W5lZGRIko4fPy6fzxc8xygxMVETJkyIfGoAAIAI6tQMUf/+/VVYWCjLsuRwOGRZls6cOaP9+/dr//79OnPmTHCdw+HQtGnTdNVVV0UrOwAAQER0qhBJ0v3336+pU6fK4XBIUvDfX33tcDj02GOPaeLEiRGKCQAAED1hPdz1xz/+sW6//XatXr1af/rTn1RXVydJSklJ0a233qof/vCHuvbaayOZEwAAIGrCftr9tddeq1mzZkUyCwAAgC06fcgMAADgckMhAgAAxqMQAQAA41GIAACA8cI+qdoEPMsMAAAzUIg64PF45PF4VFlZqYKCArvjAACAKOGQGQAAMB6FCAAAGI9CBAAAjEchAgAAxqMQAQAA41GIAACA8ShEAADAeBQiAABgPAoRAAAwHoUIAAAYj0d3AAC6TbW/WhlLMuyO0Slup1vF2cXKz8q3OwqiiEIEAIg6t9Mt+aRmq1knfCfsjtM5PqloSxGF6DJHIQIARF1xdrGKthTJ1+izO0qnVPur1Ww1x11udB6FqANer1der1d+v9/uKAAQ1/Kz8uNyhiVjSUb8zWghLBSiDng8Hnk8HlVWVqqgoMDuOAAAIEq4ygwAABiPQgQAAIxHIQIAAMajEAEAAONRiAAAgPEoRAAAwHgUIgAAYDwKEQAAMB6FCAAAGI9CBAAAjEchAgAAxqMQAQAA41GIAACA8ShEAADAeBQiAABgPAoRAAAwHoUIAAAYj0IEAACMRyECAADGoxABAADj9bQ7QCzzer3yer3y+/12RwEAAFFEIeqAx+ORx+NRZWWlCgoK7I4DAACihENmAADAeBQiAABgPAoRAAAwHoUIAAAYj0IEAACMRyECAADGoxABAADjUYgAAIDxKEQAAMB4FCIAAGA8ChEAADAehQgAABiPQgQAAIxHIQIAAMajEAEAAONRiAAAgPEoRAAAwHgUIgAAYDwKEQAAMF5PuwMAABDrqv3VyliSYXeMkLmdbhVnFys/K9/uKHGDQgQAwCW4nW7JJzVbzTrhO2F3nND5pKItRRSiTqAQAQBwCcXZxSraUiRfo8/uKCGr9ler2WqOq8yxgEIEAMAl5Gflx90sS8aSjPiazYoRnFQNAACMRyECAADGoxABAADjGVeI9u7dq7vuukuvvfaa3VEAAECMMKoQNTc368UXX9QNN9xgdxQAABBDjLrK7A9/+IOGDRums2fP2h0FAADEkJicITp37pxeffVV/fSnP9X48eN15513asOGDe2ODQQCKi0t1f333y+Px6OpU6fq448/bjPu9OnTWrVqlaZMmRLt+AAAIM7EZCE6ffq0ysrKVFVVpczMzA7HLliwQOXl5Ro7dqxmzJihHj16aPbs2dqzZ0+rcStWrNADDzwgt9sdzegAACAOxWQhSk5O1po1a7Rq1SpNnz79kuMqKir07rvv6vHHH1dhYaHy8vK0dOlSDRgwQKWlpcFxn376qQ4cOKB77723O+IDAIA4E5PnECUlJSk5Ofkbx7333ntKSEhQXl5ecJnT6dT48eO1fPly1dbWKjU1Vbt379axY8c0ceJESZLf71dCQoL+8pe/aO7cuVH7HQAAID7EZCEK1cGDB5WRkaE+ffq0Wj5s2DBJ0qFDh5Samqq8vDzdc889wfUvvPCC0tLSNGnSpHa3W1dXp/r6+uD7qqqqKKQHAACxIq4LUX19fbszSS3L6urqJEm9evVSr169guudTqd69+59yfOJ1q5dq7KyssgHBgAAMSmuC1FjY6MSExPbLE9KSgqub8+8efM63G5eXp5Gjx4dfF9VVaWSkpIuJAUAALEsrguR0+nUhQsX2iwPBALB9eFISUlRSkpKl7IBAID4EZNXmYUqOTm51bk+LVqWUWoAAEAo4roQZWZm6vjx423uPF1RURFcDwAA8E3iuhCNGTNGTU1NWrt2bXBZIBDQ+vXrlZWVpdTUVBvTAQCAeBGz5xCtXr1afr8/ePhr586dOnnypCRp4sSJcrlcysrKUnZ2tpYvX66Ghgalp6dr48aNqqmp0Zw5c7qcwev1yuv1yu/3d3lbAAAgdsVsIVq5cqVqamqC77dt26Zt27ZJksaNGyeXyyXp4hVjqamp2rRpk/x+vwYPHqxFixbplltu6XIGj8cjj8ejyspKFRQUdHl7AAAgNsVsISovLw9pnNPpVGFhoQoLC6OcCAAAXK7i+hwiAACASKAQAQAA41GIAACA8WL2HKJYwFVmAACYgULUAa4yAwDADBwyAwAAxqMQAQAA41GIAACA8ShEAADAeBQiAABgPK4y6wCX3QMAYAYKUQe47B4AADNwyAwAABiPQgQAAIxHIQIAAMajEAEAAONRiAAAgPEoRAAAwHhcdt8B7kMEAIAZKEQd4D5EAACYgUNmAADAeBQiAABgPAoRAAAwHoUIAAAYj0IEAACMRyECAADGoxABAADjUYgAAIDxuDFjB7hTNQAAZqAQdYA7VQMAYAYOmQEAAONRiAAAgPEoRAAAwHgUIgAAYDwKEQAAMB6FCAAAGI9CBAAAjEchAgAAxqMQAQAA41GIAACA8Xh0Rwd4lhkAIF5V+6uVsSTD7hghczvdKs4uVn5Wvi3fTyHqAM8yAwDEG7fTLfmkZqtZJ3wn7I4TOp9UtKWIQgQAALquOLtYRVuK5Gv02R0lZNX+ajVbzbZmphABAHAZyc/Kt22WJVwZSzJsn83ipGoAAGA8ChEAADAehQgAABiPQgQAAIxHIQIAAMajEAEAAONRiAAAgPEoRAAAwHgUIgAAYDwKEQAAMB6FCAAAGI9nmXXA6/XK6/XK7/fbHQUAAEQRhagDHo9HHo9HlZWVKigosDsOAACIEgpRCBobGyVJVVVVEd92UkOS+pzto6SmJFVWVkZ8+wAAxLpo/104aNAg9erVq8MxDsuyrIh/82XmnXfeUUlJid0xAABAGFasWKGhQ4d2OIZCFIKGhgZ99NFHevvtt/Xkk0+G/Llly5bpiSee6HBMVVWVSkpK9Mwzz2jQoEFdjXpZCOXPzS52ZIvGd0Zqm13ZTjif7exn2AfDE8v7oNT9+aL1fSbsh6GOjfZ+GMoMEYfMQtCvXz+NGzdOmzdv/saG+VUulyvk8YMGDerUti9nnflz6252ZIvGd0Zqm13ZTjif7exn2AfDE8v7oNT9+aL1fSbsh53dvp37IZfdd4LH44nqeFwUy39udmSLxndGaptd2U44n2Uf7B6x/ufW3fmi9X0m7Iex/t+lr+KQmc1armAL5fgmgMhjHwTsFwv7ITNENktOTtbkyZOVnJxsdxTASOyDgP1iYT9khggAABiPGSIAAGA8ChEAADAehSjGBQIBLVy4UPn5+crJydG0adO0d+9eu2MBRvnNb36j++67Tzk5OXrkkUe0c+dOuyMBxtq7d6/uuusuvfbaaxHdLucQxbjz589r5cqVys3N1dVXX60tW7Zo6dKlWrlypa644gq74wFGqKqqUlpampKSkrR//37NnDlTb775pq688kq7owFGaW5uVmFhoSzL0u23365HHnkkYttmhijG9e7dW5MnT1Zqaqp69Oihe+65Rz179tSxY8fsjgYYY9CgQUpKSpIkORwOXbhwQXV1dTanAszzhz/8QcOGDYvK3ay5U3WEnTt3Tm+++aYqKiq0f/9++Xw+zZ07V7m5uW3GBgIBvfLKK3rnnXfk8/k0ZMgQPfbYYxo+fPglt3/s2DH5fD6lp6dH82cAcSta++CSJUu0fv16BQIB3XbbbRo8eHB3/BwgLkVjPzx9+rRWrVql0tJSLVu2LOKZmSGKsNOnT6usrExVVVXKzMzscOyCBQtUXl6usWPHasaMGerRo4dmz56tPXv2tDu+sbFRJSUlmjRpklwuVzTiA3EvWvvgzJkztWnTJj333HMaPny4HA5HtH4CEPeisR+uWLFCDzzwgNxud3RCW4ioxsZGq66uzrIsy9q/f791xx13WOvXr28zbt++fdYdd9xhvfHGG8FlX375pfXwww9b06ZNazP+woUL1uzZs61f/OIXVnNzc/R+ABDnorUPftWcOXOsP/7xj5ENDlxGIr0fVlZWWo8++qj117/+1bIsy5o/f75VVlYW0czMEEVYUlJSSHfafO+995SQkKC8vLzgMqfTqfHjx2vfvn2qra0NLm9ublZJSYkcDofmzZvH/zMFOhCNffDrmpqadOLEiYjkBS5Hkd4Pd+/erWPHjmnixIm67777tHnzZr3xxhtasGBBxDJzDpFNDh48qIyMDPXp06fV8mHDhkmSDh06pNTUVEnS4sWLVV9fr8WLF6tnT/4jAyIh1H3Q7/fr/fff1+jRo5WUlKTt27dr165devzxx+2IDVxWQt0P8/LydM899wTXv/DCC0pLS9OkSZMiloW/XW1SX1/fbntuWdZyBUtNTY3WrVunpKSkVg3617/+tW6++ebuCQtchkLdBx0Oh9atW6fnnntOlmUpPT1dRUVFuv7667s1L3A5CnU/7NWrl3r16hVc73Q61bt374ieT0QhskljY6MSExPbLG+5tLexsVGSNGDAAG3btq1bswEmCHUf7NOnj55//vluzQaYItT98OvmzZsX8SycQ2QTp9OpCxcutFkeCASC6wFED/sgYL9Y2g8pRDZJTk5WfX19m+Uty1JSUro7EmAU9kHAfrG0H1KIbJKZmanjx4/r7NmzrZZXVFQE1wOIHvZBwH6xtB9SiGwyZswYNTU1ae3atcFlgUBA69evV1ZWVvAKMwDRwT4I2C+W9kNOqo6C1atXy+/3B6f8du7cqZMnT0qSJk6cKJfLpaysLGVnZ2v58uVqaGhQenq6Nm7cqJqaGs2ZM8fO+EDcYx8E7Bdv+yFPu4+CBx98UDU1Ne2uW7lypdLS0iRdPHu+5fktfr9fgwcP1mOPPaYRI0Z0Z1zgssM+CNgv3vZDChEAADAe5xABAADjUYgAAIDxKEQAAMB4FCIAAGA8ChEAADAehQgAABiPQgQAAIxHIQIAAMajEAEAAONRiAAAgPF4uCuAy151dbUeeuihNsszMzP16quvSpJ27dqlJ598Mrjuq89aihXz5s3Tjh072ix//vnn9b3vfc+GRMDlg0IEoMu+XiYuJScnR/PmzeuGRLHp7//+73X06FFJ0qBBg/T666+3O+78+fO67777dP78eUnSHXfcofnz53dbTsBEFCIAxvmHf/gHuVwu9e/fv1u/NycnR8uXL5ckVVVVqbKyUkOHDm0zbvv27cEy1PI5ScrNzdVNN90kSSotLe2GxIA5KEQAIu7uu+9u9y/6wYMH25CmrXvvvdeWw2E/+MEP9G//9m9qbm6WJG3cuLHdP6eNGzcGX/fr10+jRo2SdHGmqAWFCIgsChGAiBs5cqRyc3M7HPP183qee+45HT58WG+//bZqa2uVnp6uSZMmady4cTp//rxeeeUVbd68WWfOnNGgQYM0efLkVgUhmsrLy/Xiiy8G30+ePFlTpkyRJJ09e1Zr1qzR9u3bdfToUTU2Nuqqq67S3/zN3+hHP/qRrrvuuuDnrr76ag0fPlwffvihJGnz5s36x3/8R/Xs+X//U1xXV6c//elPwffjxo1rtR5AdLCXAYgJL7/8siorK4Pvjxw5opKSEp07d04bNmzQ/v37g+sOHjyoZ555RkuWLNGtt94a1Vxr1qxpVYamTp2qSZMmSZKOHTumWbNmqaamptVnTp48qY0bN2rz5s36+c9/ruzs7OC6nJycYCE6deqUPv744+AMkCT913/9V3AGSdI3FksAkUEhAhBxH374oRoaGtosv/vuu5WamtruZyorKzVy5EjdcMMNWrdunerr6yVJS5YskSSNHj1a1113nVavXq3z58/Lsiz9/ve/j2ohWrdunZYuXRp8P2PGDOXn50uSmpqa9MwzzwTLUL9+/eTxeNS3b1999NFH2rt3rwKBgObPn6+hQ4fq29/+tqSLh71cLpf8fr+ki4fHvlqI3nnnneDr66+/XkOGDIna7wPwfyhEACJu8+bN2rx5c5vlN9xwwyUL0fDhw/XrX/9aDodDV199tRYvXhxcN2rUKC1YsECSZFmWfve730mSDhw4EIX0F23atEllZWWyLEsOh0OzZs1SXl5ecP3777+vI0eOSJISEhL00ksvaeDAgZIunrT96KOP6vDhwwoEAnrrrbf0k5/8RJKUlJQkj8ejt99+W5K0c+dO+f1+uVwuHTp0SJ999lnwO5gdAroPhQhATPB4PHI4HJKkAQMGtFr31UNO6enpwdc+ny9qeVruT9SjRw/NmTOnTTn53//93+Drpqam4GG09uzdu7fV+5ycnGAhCgQC2rp1q+69915t2rQpOCYxMVFjx47t6s8AECIKEYCImzt3bqdnN1JSUoKvExMTL7kuISEh+NqyrDAThq5nz57tzmp1pox9/fBhVlaWrr32Wv35z3+WdHE2Kjc3V16vNzhm1KhRuvLKK8PKDKDzKEQAYkJHV1J9tQR1l2uuuUZHjx5VIBDQ3Llz9dxzzykrKyu43u12B18nJSXp0UcfveS2XC5Xm2U5OTl6+eWXJUl79uxpdd6UxOEyoLtRiACgHf/0T/+k+fPn6/Dhwzp//ryefvppvfDCC8GTnG+88cbg2EAgoOuuu0633XZbm+1UVFS0mfGSLl5Ov2LFCjU1NcmyLL300kvBdVdddZVGjhwZhV8F4FIoRAAi7lJXmblcLk2YMKH7A4XB5XJp8eLFmj59umpra+Xz+TRr1iwtW7ZMAwcO1KhRozRo0CBVVVVJkn7+85/rzjvv1LXXXqvm5mb95S9/0SeffKKamhrNnTtX119/favtp6SkaPjw4frggw8kSV9++WVw3dixY7n3ENDN2OMARNylrjIbMGBA3BQi6WJpWbx4sX7yk5/o9OnT+uKLL/TUU0/ppZdeUmpqqn71q18F70N04cIFvfvuu53afm5ubrAQfX05gO7Vw+4AABDLBg0apEWLFql3796SLt508amnnlJ9fb0GDhyosrIyTZ8+XTfeeKPcbrcSEhJ0xRVXaMiQIbr33ns1f/58eTyedrc9evRo9e3bt9WyoUOHxswjTgCTOKzuuEwDAGz09ceEfPXhri0PTo0H27dv1/HjxyW1fpbZ888/r+9973t2xQIuCxwyA2Cc119/XZKUmZkZV4Vow4YN2rFjh90xgMsSh8wAAIDxOGQGAACMxwwRAAAwHoUIAAAYj0IEAACMRyECAADGoxABAADjUYgAAIDxKEQAAMB4FCIAAGA8ChEAADDe/wMaF2asJqJ5VAAAAABJRU5ErkJggg==", "text/plain": [ "