{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Working with a full instrument response" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the last tutorial we saw how to use an RSP file, which contains the response for a single detector and a single coordinate. If you don't provide `bc-rsp` with the `--direction` flag, it will generate an HDF5 file containning the detector response matrices (DRMs) for all detectors and multiple directions covering the full sky (in a HEALPix grid). If you want to use the GBM software, a convenient way to get an RSP file from a full `DRM` is to use the script `bc-rsp-extractor`.\n", "\n", "You can also get information directy from the `DRM`. The class to access this file is `bctools.io.DRM`. For example, let's get the effective area at 100 keV as a function of the angle with respect to to the normal of the first Single Quarter Detector:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/imartin5/burstcube/software/gbm-data-tools/gbm/plot/lal_post_subs.py:184: UserWarning: Basemap not installed. Some functionality not available.\n", " warnings.warn('Basemap not installed. Some functionality not available.')\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Detectors: ['SQD0', 'SQD1', 'SQD2', 'SQD3']\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAGwCAYAAACkfh/eAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABxu0lEQVR4nO3dd1hUV/oH8O8dysDQ+9BUmhU7SsCaosZUY6JJTNH0rDHRlDXrmk1M08TduNldsyluNNVNTIz5mU2TGHvvKCqCVBGkdxhg5vz+GGYEEQVm4E75fp5nnkfuDHdergPzzjnveY8khBAgIiIisiEKuQMgIiIiMjcmOERERGRzmOAQERGRzWGCQ0RERDaHCQ4RERHZHCY4REREZHOY4BAREZHNcZQ7ALnodDqcP38eHh4ekCRJ7nCIiIioA4QQqKqqQkhICBSK9sdp7DbBOX/+PMLDw+UOg4iIiLogNzcXYWFh7d5vtwmOh4cHAP0F8vT0lDkaIiIi6ojKykqEh4cb38fbY7cJjmFaytPTkwkOERGRlblaeQmLjImIiMjmMMEhIiIim8MEh4iIiGyORSY4ffr0gSRJbW5PPfUUAP0SsSVLliAkJASurq6YOHEiUlJSZI6aiIiILIVFJjgHDhxAfn6+8ZaUlAQAmDFjBgBg+fLlWLFiBVauXIkDBw5ArVZj0qRJqKqqkjNsIiIishAWmeAEBARArVYbb//73/8QFRWFCRMmQAiBd999F4sXL8b06dMRGxuLTz/9FLW1tVi7dq3coRMREZEFsMgEp6WGhgZ88cUXePjhhyFJEjIzM1FQUIDJkycbH6NUKjFhwgTs3r273fNoNBpUVla2uhEREZFtsvgE5/vvv0d5eTnmzJkDACgoKAAABAUFtXpcUFCQ8b7LWbZsGby8vIw3djEmIiKyXRaf4Hz88ceYOnUqQkJCWh2/tMGPEOKKTX8WLVqEiooK4y03N7db4iUiIiL5WXQn4+zsbPz222/47rvvjMfUajUA/UhOcHCw8XhhYWGbUZ2WlEollEpl9wVLREREFsOiR3DWrFmDwMBA3HzzzcZjERERUKvVxpVVgL5OZ9u2bUhMTJQjTCIiIrIwFjuCo9PpsGbNGsyePRuOjhfDlCQJCxYswNKlSxETE4OYmBgsXboUKpUKs2bNkjFiIiIishQWm+D89ttvyMnJwcMPP9zmvoULF6Kurg5z585FWVkZ4uPjsWnTpqvuLHolQgjoBOCguPLmXURERGT5JCGEkDsIOVRWVsLLywsFRaX4PqUUH27PQH2jFvdf0xtPTYyGl8pJ7hCJiIjoEob374qKCnh6erb7OLtPcIYt/h5lTa0HstSeLvjbjKEYG+MvU3RERER0OYUlZQjy971qgmPRRcY9oaSmAeG+rvjbjKFYPScOEf5uKKisx/0f78OSjSmob9TKHSIREZHd0zRp8fmeLNz0j+0derzF1uD0lFduHYgHJwyAk4M+17sm0g9LfzqFL/bm4JPdWdiZXoy/zxyGwWFeMkdKRERkfxq1Oqw/dA7/+j0deeV10GkaOvR9dj9F1d4Q15bUQiz8NhlFVRo4KiQsuCEGT06IgqOD3Q96ERERdbsmrQ7/d/Q8/rE5DTmltQCAQA8lHokPwpOThrAGpz0dKVIqrWnA4g3H8fMJ/RYQw3t54+8zh6GPv1tPhkpERGQ3dDqBH5L1iU1GUQ0AwN/dGU9OiML91/RGQ10Ni4yvpKNV2EIIbDiSh1f+LwVVmia4OjngpVsGYNboXlfcGoKIiIg6TqcT+DWlAH//7QzOXKgGAPionPDEhCg8mNAbKmd9VQ1XUV1FRy+QwbmyWrzwzTHszSgFAFzbLwBv3zUEgR4u3R0qERGRzRJC4LdThfh70hmczK8EAHi6OOKxcZGYM6YPPFxat21hgnMVnU1wAH12uXpXJpb/kooGrQ4+Kicsmz4YN8YGX/2biYiIyEgIgW1nivD3pDM4dq4CAOCudMTDY/rgkXGR8HK9fD86JjhX0ZUExyC1oAoLvj6KU82Z5p0jwvDKbQPh6cLmgERERFezO70Y7ySdwaHsMgCAq5MD5ozpg8fHRcLHzfmK38sE5ypMSXAA/Xr8d39LwwfbzkIIINTbFe/MHIprIv26IVoiIiLrtz+zFCuSUo3lHkpHBR64pjeenBgFf3dlh87BBOcqTE1wDA5kleK5dUeRW1oHSQIeGxeJ5yf3hdLRwYzREhERWa8jOWVYkXQGO9KKAQDODgrcOzocc6+NRpBn52pZmeBchbkSHACo1jTh9R9O4uuDuQCA/moP/P3uYRgQbNp5iYiIrNmJvAqsSDqD308XAgAcFRJmxIXj6euiEeLt2qVzMsG5CnMmOAZJJy/gT+uTUVLTACcHCc9P7ofHxkVyh3IiIrIrpwsq8fekM/g15QIAwEEhYfrwUDxzfQzCfVUmnZsJzlV0R4IDAMXVGvxp/XH8dkr/nzq6jy/emTnU5P9QIiIiS5deWIW//5aGH5PzAQCSBNw+NATzb+iLCDM1yWWCcxXdleAA+qVv3xw8h1d/SEFNgxbuSke8cutA3DUyjM0BiYjI5mQV1+Afm9Pwf0fzoGvOKm4eHIwFN8QgJsjDrM/FBOcqujPBMcgpqcXz3xzFgSz9MrjJA4OwbPpg+HWwUpyIiMiS5ZbW4l+/p2H94TxomzObyQOD8Oykvt1Wh8oE5yp6IsEBAK1O4KPtGViRlIpGrYC/uzPevnMIrh8Q1G3PSURE1N0yi2tw6792olrTBEDf4f+5Sf0wOMyrW5+XCc5V9FSCY5ByvgLPfX0MqReqAAD3jg7HSzcPhJvSsdufm4iIyJyEELh31V7szSjFoBBPvHZ7LEb29umR5+7o+7eiR6IhDArxwv/NG4PHxkVAkoD/7s/F1H/swKHsUrlDIyIi6pT1h/OwN6MUrk4O+OD+kT2W3HQGE5we5OLkgMU3D8SXj8YjxMsFOaW1mPHBHvz119No1OrkDo+IiOiqdDqBf29NBwA8fX20xa4SZoIjg8Qof/zy7HhMHxEKnQDe23IWy385LXdYREREV/X76UJkFNXAw8URDyb0kTucdjHBkYmnixNWzByGd2YMBQB8vDMTx3LL5Q2KiIjoKj7akQEAmBXfC+4WXEfKBEdmd44Mw+3DQqATwEvfn4Cd1nwTEZEVOJZbjv2ZpXBUSHgoMULucK6ICY4FePmWgXBxUuB4XgX2ZbLomIiILNOq5tGb24aGQO3VuU0yexoTHAvg567EnSPCAACrd2bKHA0REVFbuaW1+Om4fguGR8dFyhzN1THBsRAPjekDAEg6dQHZJTXyBkNERHSJNbuyoBPA2Gh/DAzp/v5xpmKCYyGiAz0woW8AhAA+2Z0ldzhERERGFbWN+OpADgDgsfGWP3oDMMGxKI+M1RdsfXPwHKrqG2WOhoiISG/t/hzUNmjRL8gD42P85Q6nQ5jgWJBxMf6IDnRHtaYJ6w6ekzscIiIiNDTp8MlufX3oo+MiIEmSzBF1DBMcCyJJEh4eox/F+WR3pnFnViIiIrn8cOw8LlRqEOihxG3DQuQOp8OY4FiYO4aHwlvlhNzSOiSdvCB3OEREZMeEEMal4XPG9IHS0UHmiDqOCY6FcXV2wKzRvQAAq3dxyTgREclnZ3oxThdUQeXsgPtG95Y7nE5hgmOBHkzoA0eFhP2ZpTiRVyF3OEREZKc+2q4fvZkZFw4vlZPM0XQOExwLpPZywc1DggFwFIeIiORxKr8SO9KKoZAurvK1JkxwLNRDzcXGPxw7j8LKepmjISIie/OfHfoP2FNjgxHuq5I5ms5jgmOhhoV7Y2RvHzRqBb7Ymy13OEREZEcuVNZj47E8APql4daICY4FMwwJfrEvB/WNWpmjISIie/HJ7iw0agVG9/HF8F4+cofTJUxwLNjkgUEI9XZFaU0DNh49L3c4RERkB2o0TfiyeebAWkdvACY4Fs3RQYHZifpleat3ZUIINv4jIqLute5gLirrmxDh74YbBgTJHU6XMcGxcHeP6gWVswNOF1Rh99kSucMhIiIb1qTV4eOd+uLiR8ZGQKGwjm0ZLocJjoXzcnXCjJFhAIDVO7lknIiIus8vKQU4V1YHXzdn3DkiTO5wTMIExwrMaV4yvvl0ITKLa2SOhoiIbJEQAquaG/s9cE1vuDpbz7YMl8MExwpE+Lvh+v6BAIA1bPxHRETd4EBWGY6dq4CzowIPJFjXtgyXwwTHSjzcvGT8m4PnUFHbKHM0RERkawzbMtw5Igz+7kqZozEdExwrkRjlh/5qD9Q1avH1wRy5wyEiIhtytqgam09fAGDdS8NbssgEJy8vD/fffz/8/PygUqkwbNgwHDp0yHi/EAJLlixBSEgIXF1dMXHiRKSkpMgYcfeTJAkPN9fifLo7G01ancwRERGRrfh4ZyaEAG4YEIioAHe5wzELi0twysrKMGbMGDg5OeHnn3/GyZMn8c4778Db29v4mOXLl2PFihVYuXIlDhw4ALVajUmTJqGqqkq+wHvAbcNC4OvmjLzyOvyackHucIiIyAaUVGuw/tA5AMBj4yJljsZ8LC7BefvttxEeHo41a9Zg9OjR6NOnD66//npERUUB0I/evPvuu1i8eDGmT5+O2NhYfPrpp6itrcXatWtljr57uTg54P74XgC4yzgREZnH53uzoWnSYWiYF0ZH+ModjtlYXIKzceNGxMXFYcaMGQgMDMTw4cOxatUq4/2ZmZkoKCjA5MmTjceUSiUmTJiA3bt3t3tejUaDysrKVjdrdP81veHkIOFQdhmO5pbLHQ4REVmx+kYtPttj2JYhEpJkvY39LmVxCU5GRgbef/99xMTE4Ndff8WTTz6JZ555Bp999hkAoKCgAAAQFNS6fXRQUJDxvstZtmwZvLy8jLfw8PDu+yG6UaCnC24dGgKAS8aJiMg03x3OQ2lNA0K9XTE1Vi13OGZlcQmOTqfDiBEjsHTpUgwfPhxPPPEEHnvsMbz//vutHndplimEuGLmuWjRIlRUVBhvubm53RJ/TzAUG/+YnI+CinqZoyEiImuk0wn8Z4d+afjDYyPg6GBxKYFJLO6nCQ4OxsCBA1sdGzBgAHJy9Euj1Wp9hnnpaE1hYWGbUZ2WlEolPD09W92sVWyofp60SSfw2Z4sucMhIiIrtPl0ITKKa+Dh4oi7R1nnrMaVWFyCM2bMGKSmprY6dubMGfTure+qGBERAbVajaSkJOP9DQ0N2LZtGxITE3s0VjkZRnHW7s9BXYNW5miIiMjarGoevbkvvjfclY4yR2N+FvcTPfvss0hMTMTSpUsxc+ZM7N+/Hx999BE++ugjAPqpqQULFmDp0qWIiYlBTEwMli5dCpVKhVmzZskcfc+ZNDAI4b6uyC2tw4YjeZjVvLqKyJrpdAJNOoEmnQ5NOgGtVv+11nDskq+1OoFGbeuvL36fzvhYXzdnjInyt+qdkYnM6WhuOfZnlsJRIWFOYh+5w+kWFpfgjBo1Chs2bMCiRYvw2muvISIiAu+++y7uu+8+42MWLlyIuro6zJ07F2VlZYiPj8emTZvg4eEhY+Q9y0EhYU5iBF7/30ms3pWJe0eH21T1O9mG/Io6fLwjEzvTi1HXqG1OUFonIo0tvhai+2KZndAbr94e231PQGRFDKM3tw0LgdrLReZouockRHf+SbFclZWV8PLyQkVFhdXW41TVNyJh2e+o1jTh04dHY0LfALlDIjISQmDS37cjvbDa5HM5OUhwUEhwVCjgoJDafO2okODoIMFBoYCjQmp1TIKEPRklAIC3pg/GPaM52kn2Lbe0FhP+ugU6Afw8fxwGBFvXe2BH378tbgSHOs7DxQkz4sKwZlcWVu/MZIJDFmV/ZinSC6vhrnTE8ruGIMjT5WLy4dCcgBgSlBYJi+E+w9cKqe2qyc761+Y0vJN0Bu9tTWeCQ3Zv9a5M6AQwLsbf6pKbzmCCY+UeSozAJ7uzsO1MEdILqxAdaD/TdGTZ1h3Ut36/dWgwbhocLGssD42NwD82pyG3tA65pbUI91XJGg+RXCpqG/H1AX2bFFvaluFyLG4VFXVOLz8VJg3QL49fvStL3mCImlXVN+Kn4/kAgBlx8i8/dVc6Yli4NwBgV3qxvMEQyWjt/hzUNmjRX+2BcTH+cofTrZjg2ICHx+qXjH93+BzKahpkjoZI34SyrlGL6EB3DG9OLOQ2Jlr/x3wnExyyUw1NOnyyW98B39a2ZbgcJjg2ID7CFwODPVHfqMN/D+TIHQ4R1h3UD4HPjAuzmD+iY5s/re45WwKdzi7XVpCd++HYeVyo1CDIU4nbmrf8sWVMcGyAJEnGUZzPdmejUauTOSKyZ+mFVTicUw4HhYQ7hofJHY7R0DBvqJwdUFLTgNMFVXKHQ9SjhBDGpeFzEiPg7Gj7b/+2/xPaiVuHBsPfXYmCynpj7QORHL5pLi6+rn8gAjyUMkdzkbOjAvERvgCA3Wc5TUX2ZUdaMU4XVEHl7IBZdrKSkAmOjVA6OuCBa/TbWazemQk7bW9EMmvU6rD+cB4AYKYFFBdfinU4ZK8Mozd3jwqHl8pJ5mh6BhMcG3LfNb3g7KDAsXMVOJxTLnc4ZIe2phahuFoDf3clJvazvL5MhgRnf2YpGpo4lUv24VR+JXakFUMhXdzH0B4wwbEh/u5K3D5MXzi2elemzNGQPTIUF985IhRODpb356VfkAf83JxR26DF0dxyucMh6hH/2aF/P5g6ONiuekBZ3l8gMomh2PiXEwXIK6+TORqyJ4VV9fj9dCEAYEac5RQXt6RQSEjkNBXZkYKKemw8pp82ftzGG/tdigmOjRkQ7InEKD9odQKf7c6SOxyyI98fyYNWJzCil7dFd9QeE+UHANjNBIfswCe7s9CoFRjdxxdDLaQnVU9hgmODDHOs/92fgxpNk8zRkD0QQhi3ZrDE4uKWDHU4R3PLUc3fD7Jh1ZomrN2XDQB4bLx9jd4ATHBs0nX9A9HHT4XK+iZ8d/ic3OGQHTiSW470wmq4Ojng5iHy7jt1NeG+KvT2U6FJJ7A/s0TucIi6zboDuaisb0Kkvxuu7x8odzg9jgmODVIoJDzUPIqzelcWu7ZSt/umubj4psHB8HCx/CWoiVHNdThpTHDINjVpdfh4p764+JFxEVAoLKOjeE9igmOj7hoZBg8XR2QW12DrmUK5wyEbVtvQhB+O6ZtLzrTQ4uJLjW2epmLDP7JVv6ToF5r4ujnjzhHW8XtpbkxwbJSb0hH3jNLXQqzemSVvMGTTfj5egGpNE/r4qTC6uVOwpUtoLjQ+XVCFoiqNzNEQmZcQAqu26xv7PZjQGy5ODjJHJA8mODbswYQ+UEj65bCp3HuHuomh982MuHCL2VjzanzdnDEoxBMAR3HI9uzPLMWxcxVQOiqMHe7tERMcGxbuq8KUQWoA+u0biMwtq7gG+zJLoZCA6SNC5Q6nUwyrqXZxuTjZmFXNjf2mjwiDn7vl7AfX05jg2LhHmhv/bTiah5JqDsWTeX17SL9Kb3zfAAR7ucocTedcTHBKuHcb2YyMompsPn0BAPDoOPvZluFymODYuJG9fTAkzAsNTTqs3ZcjdzhkQ7Q6YUxwLL33zeWM6uMDJwcJeeV1yC6plTscIrP4eGcmhABuGBCIqAB3ucORFRMcGydJkrHx32d7s7nBIJnNjrQiFFTWw0flhOsHWF+PDZWzI0b08gHAbRvINpRUa4wfOh61s20ZLocJjh24aXAwAj2UKKrS4H/J5+UOh2zEN82di6cND4XS0TpXaXC5ONmSL/bmQNOkw+BQL8RbyYrG7sQExw44OyowO7EPAMPwJesNyDSlNQ3YdLIAADBjpPVNTxkkGhOcEjbEJKtW36jF53uzAOhrb6xlRWN3YoJjJ+4d3QtKRwVSzlfiQFaZ3OGQlfv+SB4atQKDQ70wsHm5tTUaGuYFd6UjymsbcTK/Uu5wiLrs+yN5KK5uQIiXC24abNnbpfQUJjh2wtfN2biMl0vGyRT6jTX1vW+spXNxexwdFLgmUj+UzzocslY6ncB/mv+uPzQmAk4OfGsHmODYFcP+VJtOFuBcGVeNUNecyKvE6YIqODsqcNtQ6+p9cznsh0PWbuuZQqQXVsND6Yh7RlvvlLG5McGxI32DPHBNpC904mL/EqLOMoze3DhIDS+V5W+seTWGBOdAVinqG7UyR0PUeau260dv7hkdbhWb3fYUJjh25u7m/am+OXiORZXUafWNWvzf0TwA1tn75nJiAt0R4KFEfaMOh3NYn0bW5UReBfZklMBBIWHOGPtu7HcpJjh2ZmpsMDxcHJFXXoddXBpLnfRrSgEq65sQ6u2KxOYNK62dJEkY0/yz7E4vkTkaos75zw79ppo3Dw5GqLd1dRPvbkxw7IyLkwNuHxYCAPj6QK7M0ZC1MfS+uWtkGBQK21mGapimYqExWZPz5XX4X3I+AOAxNvZrgwmOHbo7rhcAYFPKBZTVNMgcDVmL3NJa46jfXSOte/XUpQwJTvK5clTWN8ocDVHHfLI7C006gfgIXwwO85I7HIvDBMcOxYZ6YkCwJxq0OnzfXE9BdDXrD5+DEMCYaD+E+6rkDsesQrxdEenvBp0A9p7lNBVZvqr6Rvy3eX9Bjt5cHhMcOyRJEu5u7l/y9YFcdjamq9LphHF6ylaKiy81pkVXYyJL9/3R86jSNCHS3w3X9be+veB6AhMcOzVteCicHRU4XVCF43kVcodDFm5PRgnyyuvg4eKIKYPUcofTLcZE6wuNWYdD1mBdcw3lrPheNlUPZ05McOyUt8rZ+EbFYmO6GkPvm9uHhcDFyTo31ryahEh/SBKQXliNgop6ucMhatfJ85U4nlcBJwcJ00fYVj2cOTHBsWN3N081bDx6HnUNbHBGl1dR24ifT+g31rTV6SkA8FI5YXCovlCTu4uTJTN84Jg0MAi+bs4yR2O5mODYscQoP4T5uKJK04SfT+TLHQ5ZqI3J59HQpEN/tYcxAbBVXC5Olq6+UYsNR2yr2WZ3YYJjxxQKCTNG6n9BOE1F7fmm+dPijLhwSJJtz/WPiWouNE4vYfE9WaRfUwpQUdeIEC8XjIsJkDsci8YEx87dFRcGSQL2ZZYiq7hG7nDIwpzKr0TyOf1c/7TmBpG2LK6PD5wdFSiorMfZIv4+kOUxTE/dNTIMDiwuviImOHYu1NvV+Cngm0McxaHWDEvDbxgQBD93pczRdD8XJwfE9fYBwDocsjy5pbXY1bydyAxOT10VExwyFht/e+gcmrQ6maMhS9HQpMOGI7bd++ZyjHU4aUxwyLIYpottsdlmd2CCQ7hhYCB8VE64UKnB9rQiucMhC7H51AWU1TYiyFOJcTH+cofTY8Y2Jzh7Mkqg1bEOhyyDVifwzSH7+8BhCiY4BKWjA+4YfrGzMRFwca7/zhFhcHSwnz8VsaFe8HRxRFV9E5tgksXYkVaE/Ip6eLk62WyzTXNz7Mo3bdy4sdPfM2nSJLi6cit3S3X3qHCs3pWJzacKUVSlQYCH7ddbUPsKKuqx7Yx+NM/e5vodFBISovzwa8oF7EovxrBwb7lDIjJ+4Jhmw802za1LCc60adM69XhJkpCWlobIyKtvCLZkyRK8+uqrrY4FBQWhoEDfaEwIgVdffRUfffQRysrKEB8fj/feew+DBg3qVEzUWj+1B4aGe+NYbjk2HDmHx8dHyR0SyWj94XPQCWB0H19E+LvJHU6PGxPtb0xwnro2Wu5wyM6VVGuQdPICAGDmKPv6wGGKLo87FxQUQKfTdeimUnWuGGrQoEHIz8833o4fP268b/ny5VixYgVWrlyJAwcOQK1WY9KkSaiqqurqj0LNDMXG3IDTvgkhWvS+sc828IZC44PZZahvZJdvkteGI3lo1ArEhnpiUIhtN9s0py4lOLNnz+7UdNP9998PT0/PDj/e0dERarXaeAsI0C9jFkLg3XffxeLFizF9+nTExsbi008/RW1tLdauXdvpn4Nau3VoMFydHHC2qAaHc8rkDodkciCrDFkltXBzdsBNg4PlDkcWkf5uUHu6oKFJh4NZ/F0g+QghjNNTd9vZdLGpupTgrFmzBh4eHh1+/Pvvvw9//46vwkhLS0NISAgiIiJwzz33ICMjAwCQmZmJgoICTJ482fhYpVKJCRMmYPfu3Vc8p0ajQWVlZasbtebh4mR8Q2Oxsf0y/DG9ZUgI3JRdmsW2epIkcdsGsghHc8tx5kI1lI4K3DYsVO5wrIrFLY2Ij4/HZ599hl9//RWrVq1CQUEBEhMTUVJSYqzDCQoKavU9LWt02rNs2TJ4eXkZb+HhzIQv5+7m+d3/JeejWtMkczTU06o1TfgxWb8v2cxR9jk9ZTAm2g8AG/6RvAwfOKbGquHl6iRzNNal0wlOXV0d8vLy2hxPSUkxS0BTp07FnXfeicGDB+OGG27Ajz/+CAD49NNPjY+5dD8cIcRV98hZtGgRKioqjLfcXI5QXM6oPj6I9HdDbYMWPyaflzsc6mE/Jp9HXaMWkQFuGNHLR+5wZGUYwTmeV4Hy2gaZoyF7VNvQhB+OGT5w8EN5Z3Uqwfn222/Rt29f3HTTTRgyZAj27dtnvO+BBx4we3AA4ObmhsGDByMtLQ1qtX7t/6WjNYWFhW1GdS6lVCrh6enZ6kZtSZJkXBbMaSr7s+7gxUZitr6x5tUEebogJtAdQgB7M0rkDofs0I/NI+m9fFW4JsJP7nCsTqcSnDfeeAOHDx/GsWPHsHr1ajz88MPG4t7uWnWj0Whw6tQpBAcHIyIiAmq1GklJScb7GxoasG3bNiQmJnbL89ujO0eGwkEh4XBOOdILuTrNXqQXVuNQdhkcFBKmD+dcPwDW4ZCsDNNTM+PCoODGmp3WqQSnsbHRuKIpLi4O27dvx4cffojXXnvNbJ/2XnjhBWzbtg2ZmZnYt28f7rrrLlRWVmL27NmQJAkLFizA0qVLsWHDBpw4cQJz5syBSqXCrFmzzPL8BAR6uODafoEAOIpjTwybrV7bLwCBni4yR2MZDAmOYYNDop5ytqgaB7LKoJCAu0ZyeqorOpXgBAYGIjk52fi1n58fkpKScOrUqVbHTXHu3Dnce++96NevH6ZPnw5nZ2fs3bsXvXv3BgAsXLgQCxYswNy5cxEXF4e8vDxs2rSpU6u66OoMxcbfHc5DQxM34LR1jVod1h/S19bZW+fiK4mP9IVCAjKLa5BXXid3OGRHDKM3E/oGQO3FDxxd0ak1oJ9//jkcHVt/i7OzM/773/9i3rx5Zgnoq6++uuL9kiRhyZIlWLJkiVmejy7v2n4BCPBQoqhKg99PX8CNsfbZD8VebEstQnG1Bv7uzriuf6Dc4VgMTxcnDA33xpGccuxKL+Ymh9QjWn7guHtUL5mjsV6dGsEJCwszFvpeasyYMWYJiCyDo4MCd47gBpz2wvBp8Y7hoXCyo401O2JMlH6aajfrcKiHbDldaPzAcf0AfuDoKrN08aqvr0dycjIKCwuh07WezrjtttvM8RQkg5lxYfhg21lsO1OEgop6DpPaKP0oXSEATk9dzphof6zcko5dZ0s61JKCyFSGDxzTR4TxA4cJTE5wfvnlFzz44IMoLm776UaSJGi13MfFWkUGuGN0H1/szyrFt4dyMe+6GLlDom7w/ZE8NOkEhoV7o28Qa9kuNaK3N1ycFCiq0iCtsJrXiLpVYWU9tqQWAQCnRE1kcmo4b948zJgxA/n5+W022WRyY/0MzaXWHTwHnY4bcNqalvvc8I/p5SkdHTCqjy8AYGcap6moe317+By0OoGRvX0QHegudzhWzeQEp7CwEM8999xVG+2RdbppsBruSkfklNZibyaXytqao7nlSCushouTArcMZSF5e8Y2Lxfntg3UnYQQ+Ka52SY31jSdyQnOXXfdha1bt5ohFLJEKmdH3Do0BACwjsXGNsfQufim2GB4unCfm/YY+uHszShFo5ZtE6h77M8sRWZxDdycHXDzEH7gMJXJNTgrV67EjBkzsGPHDgwePBhOTq3/SD7zzDOmPgXJ7O5R4fjv/hz8fKIAr9Y1csM3G1HXoMUPx/T7jbG4+MoGBnvCW+WE8tpGJJ8rx8jevnKHRDbo6+bp4luGhMBNaZY1QHbN5Cu4du1a/Prrr3B1dcXWrVtbrTCQJIkJjg0YGuaFfkEeSL1QhY1H8/BAQh+5QyIz+PnExX1u4iP4hn0lCoWExCg//HS8ALvSS5jgkNlV1jfip+PcWNOcTJ6ieumll/Daa6+hoqICWVlZyMzMNN4yMjLMESPJTJIk4y+c4RMGWT9DcfGMkdznpiO4LxV1px+OnUd9ow7Rge4Y0ctb7nBsgskJTkNDA+6++24oFFyrb8v0DeAknMirRMr5CrnDIRNll9Rgb0YpJAm4c2SY3OFYBUPDvyM5ZahtaJI5GrI1hhrHu+PC2WvJTEzOSmbPno2vv/7aHLGQBfN1c8bkgfou1iw2tn7fHtIXF4+LCUCIt6vM0ViH3n4qhHq7olErsD+zVO5wyIacLqjEsXMVcFRIuGNEqNzh2AyTa3C0Wi2WL1+OX3/9FUOGDGlTZLxixQpTn4IsxMxR4fjxeD6+P3oei24aABcnB7lDoi7Q6oQxwZkZx9GbjpIkCWOi/bDu4DnsPluCif3YQp/Mw7Adzg0DguDvrpQ5GtthcoJz/PhxDB8+HABw4sQJkwMiyzU22h8hXi44X1GPX1MKcPswftKwRjvTi5FfUQ9vlRMmDWT/qs4YE+2PdQfPseEfmY2mSYsNRwwba7K42JxMTnC2bNlijjjICjgoJNwVF45/bk7DuoO5THCslKG4eNqwUCgdOQrXGYnNdTgn8ytRWtMAXzdnmSMia5d08gLKaxuh9nTB+L4BcodjU0yuwVm2bBlWr17d5vjq1avx9ttvm3p6sjAzmgtSd6WXILe0VuZoqLPKahqQlHIBADCD01OdFuChRH+1fi8qdjUmczBMT901MgwOXM1oViYnOB9++CH69+/f5vigQYPwwQcfmHp6sjDhviqMifYDAHzDJeNW5/+O5qFBq8OgEE8MCvGSOxyrZFguvovLxclE58pqjW0HuBec+Zmc4BQUFCA4uG1L6YCAAOTn55t6erJAhl/Ebw7pN4Uj62HYmoF/TLvOkODvSufebGSabw+dgxBAQqQfevmp5A7H5pic4ISHh2PXrl1tju/atQshISGmnp4s0JRBani5OiG/oh470orkDoc66EReBU7mV8LZQYHbh/F3s6tGR/jBUSEhp7SW07TUZTpdi401WVzcLUxOcB599FEsWLAAa9asQXZ2NrKzs7F69Wo8++yzeOyxx8wRI1kYFycHTGt+g1zHaSqrYZhSnDwoCN4qFsd2lbvSEcPCvQFwmoq6btfZYuSV18HDxRE3xqrlDscmmbyKauHChSgtLcXcuXPR0NAAAHBxccGLL76IRYsWmRwgWaaZo8Lx6Z5sJJ28wNUkVqC+UYvvj+o31uT0lOnGRPvjYHYZdqYX457RveQOh6yQobh42rBQ9hTrJiaP4EiShLfffhtFRUXYu3cvjh07htLSUrz88svmiI8s1KAQL8SGeqJRK4w9HMhyJZ28gIq6RoR4uRiLZKnrDNdwz9kS6FiHRp1UVtOATc2rGTk91X3MtoGUu7s7Ro0ahdjYWCiV7MRoD+5uHglYdyAXQvCPvCUzTCVyKap5DAv3hsrZASU1DThdUCV3OGRlvm9ezTgw2BOxoVzN2F26lOAkJydDp9N1+PEpKSloauLmdLbmtmGhUDoqkHqhCsfOcQNOS5VXXmdcinrXSH5aNAdnRwVGR/gCYD8c6hwhhHF6iqM33atLCc7w4cNRUtLxJZIJCQnIycnpylORBfNydcLU5uK4r7kBp8XSj7BxKaq5jW2eptrJQmPqhON5FThdUAVnRwWmsRt8t+pSkbEQAn/5y1+gUnXsj6Wh+Jhsz8xR4fj+6Hn8cOw8/nLLAKicTa5bJzNq0uqMyee98SyGNSdDHc6+jFI0NOng7Gi2GX+yYYbfxxsHqeGlcrrKo8kUXXo3Gj9+PFJTUzv8+ISEBLi6unblqcjCXRPhh16+KuSU1uKn4wW4ayTb/1uSralFKKish6+bM6YM4saa5tQvyAN+bs4oqWnA0dxy45QVUXvqGrTY2LyakdNT3a9LCc7WrVvNHAZZK4VCwsy4MPxt0xmsO5DLBMfC/He/fmr4rpFh3FjTzBQKCYnR/vjh2HnsTC9mgkNX9dPxfFRpmhDu64qESD+5w7F5HFMlk901MhwKCdifVYqMomq5w6Fm58vrsCW1EABwDz8tdosxUfo3qd2sw6EO+Lp5NePMkeFQcDVjt2OCQyZTe7lgQt8AABf3OiL5fX0gF7rm4uLIAHe5w7FJhjqco7nlqNZwpSi1L7O4BvszS6GQgLviONLdE5jgkFkY5pPXHz6HJm3HWwhQ92jS6oy9b1hc3H3CfVXo5atCk05gfyY336T2GX4fx/cNQLAXa1J7AhMcMovr+gfBz80ZRVUabEnlBpxy23amCPkV9fBRObG4uJsZRnF2pjHBoctr0uqw/lDzxprcKqXHMMEhs3B2VGD6CH1PB/bEkd/afSwu7iljopvrcNjwj9qxNbUIhVUa+Lk54/oB/MDRU8zWtOTkyZPIyclp0/PmtttuM9dTkIW7e1Q4Vu3IxJbUQhRW1iPQ00XukOxSy+Lie7kRZLdLjNKP4JwuqEJhVT0CPfi6p9YMxcV3DA9lv6QeZHKCk5GRgTvuuAPHjx+HJEnGPYkkSV8hrtVqTX0KshLRgR4Y0csbh3PKsf5wHv4wMUrukOzSuoP64uJrIn1ZXNwDfN2cMSjEEynnK7HnbAluZ3daaqGwqh6/n9Z/4GDvm55lcio5f/58RERE4MKFC1CpVEhJScH27dsRFxfHfjl2yPAL/M1BbsAph1adizl602MMdTi7uFycLvHd4TxodQLDe3kjJshD7nDsiskJzp49e/Daa68hICAACoUCCoUCY8eOxbJly/DMM8+YI0ayIjcPCYHK2QEZxTU4kFUmdzh2p2Vx8Y3N+4RR97uY4JQwsScjIQTWGTbWZHFxjzM5wdFqtXB31w+D+/v74/x5fRvq3r17d2o7B7IN7kpH3DIkGACLjeXAzsXyGNXHB04OEvLK65BdUit3OGQhDmaXIaO4BipnB9wyNETucOyOyQlObGwskpOTAQDx8fFYvnw5du3ahddeew2RkZEmB0jWxzBN9dPxfFTVN8ocjf04X15nnOu/h9NTPUrl7IgRvXwAcHdxusjwIe/mwcFwV3Ij4p5mcoLz0ksvQafTN3Z7/fXXkZ2djXHjxuGnn37CP//5T5MDJOszopcPogPdUdeoxS8nCuQOx260LC6OYnFxjzNMU3G5OAFAVX0jfkzOB8DiYrmYnOBMmTIF06dPBwBERUXh5MmTKC4uRmFhIa677jqTAyTrI0kSbmqu/9jKpn89QqsTLC6W2cUEpwQ6Hetw7N3/kvNR16hFZIAbRvb2kTscu2SWBfk7duzA/fffj4SEBOTl5cHX1xdffPEFdu7caY7TkxWa0C8QALAjrYhbN/SAbWcKW3QuZnGxHIaGecFd6Yjy2kaczK+UOxyS2dctiosNbVOoZ5mc4Kxfvx5TpkyBq6srjhw5Ao1GAwCoqqrC0qVLTQ6QrNOwcG94q5xQWd+Eo7nlcodj8wydi+8cEQYXJxYXy8HRQYFrIn0BsA7H3p25UIWjueVwVEiYPoIba8rF5ATnjTfewAcffIBVq1bBycnJeDwxMRGHDx829fRkpRwUEsbF6HcY5zRV98qvuFhczI015cV+OARcHL25rn8gAjyUMkdjv0xOcFJTUzF+/Pg2xz09PVFeXm7q6cmKTeyrT3C2pzHB6U7rDpyDTgDxESwulpshwTmQVYr6RnZxt0cNTTpsOJIHgMXFcjM5wQkODkZ6enqb4zt37uQycTtn+GN/Iq8CFXVcLt4d9MXF+umpWRy9kV1MoDsCPJSob9ThcA4bXdqj305dQGlNAwI9lJjQ/CGP5GFygvPEE09g/vz52LdvHyRJwvnz5/Hll1/ihRdewNy5c00OcNmyZZAkCQsWLDAeE0JgyZIlCAkJgaurKyZOnIiUlBSTn4vMS+3lgkh/N+gEsD+zVO5wbNK2M4U4z+JiiyFJEsZENe8unl4iczQkB8P01F0jw+DowI015WTy1V+4cCGmTZuGa6+9FtXV1Rg/fjweffRRPPHEE5g3b55J5z5w4AA++ugjDBkypNXx5cuXY8WKFVi5ciUOHDgAtVqNSZMmoaqqyqTnI/NLMPyxZ2+QbrF2n/6PKYuLLYdh5JKFxvbnfHmdcUp+JrdmkJ1Z0ss333wTxcXF2L9/P/bu3YuioiK8/vrrJp2zuroa9913H1atWgUfn4s9BIQQePfdd7F48WJMnz4dsbGx+PTTT1FbW4u1a9ea+qOQmSVG6f/Y7znLT7Pmpi8uvgCAnYstiSHBST5Xjkp28rYr3x46B9FcD9fH303ucOyeSQlOY2Mjrr32Wpw5cwYqlQpxcXEYPXq0cW8qUzz11FO4+eabccMNN7Q6npmZiYKCAkyePNl4TKlUYsKECdi9e3e759NoNKisrGx1o+5nWDZ7uqAKxdUamaOxLS2Li6MDWVxsKUK8XY1Ts3uZ2NsNnU5g3cHm3jcsLrYIJiU4Tk5OOHHihNmbGH311Vc4fPgwli1b1ua+ggJ96/+goKBWx4OCgoz3Xc6yZcvg5eVlvIWH8wXYE/zcleiv9gAA7M3gH3tzYXGxZUuM1k/Ncrm4/diTUYJzZXXwUDpiamyw3OEQzDBF9eCDD+Ljjz82RywAgNzcXMyfPx9ffPEFXFxc2n3cpUmVEOKKidaiRYtQUVFhvOXmcqfrnmKYptrNT7Nms/1MEc5X1MObxcUWaayhHw5f83bDUFx827AQuDqzHs4SmLy9aUNDA/7zn/8gKSkJcXFxcHNrPe+4YsWKTp3v0KFDKCwsxMiRI43HtFottm/fjpUrVyI1NRWAfiQnOPhillxYWNhmVKclpVIJpZINl+SQEOWH1bsyOVxvRl+yc7FFuybSD5IEpBdWo6CiHmqv9j+skfWrqG3ELyn6GQROT1kOkxOcEydOYMSIEQCAM2fOtLqvK1NX119/PY4fP97q2EMPPYT+/fvjxRdfRGRkJNRqNZKSkjB8+HAA+iRr27ZtePvtt7v4U1B3Gh3hC4UEZBTXIL+iDsFernKHZNUKKuqNxcXcWNMyeaucMTjUC8nnKrD7bDHb9du474/moaFJh/5qDwwO9ZI7HGpmcoKzZcsWc8Rh5OHhgdjY2FbH3Nzc4OfnZzy+YMECLF26FDExMYiJicHSpUuhUqkwa9Yss8ZC5uHl6oTBoV44dq4Ce86W8I+9idYdzIVO6BNHFhdbrjHR/kg+V4Gd6UxwbJ1xY81R3FjTkpic4BicPHkSOTk5aGhoMB6TJAm33nqruZ7CaOHChairq8PcuXNRVlaG+Ph4bNq0CR4eHmZ/LjKPhCh/HDtXgd1McEyi1Ql8tb+5uJijNxZtTJQ/3t96FrvTS65aI0jW60ReBU7mV8LZQYFpw0LlDodaMDnBycjIwB133IHjx49DkiQIIQBcnJ7Sak3fj2Xr1q2tvpYkCUuWLMGSJUtMPjf1jMQoP3yw7Sz2nOUfe1O0LC6+MZbFxZYsro8PnB0VKKisx9miGo622SjD6M3kQUHwcXOWORpqyeRVVPPnz0dERAQuXLgAlUqFlJQUbN++HXFxcW0SE7JfcX184OQgIa+8DjmltXKHY7XW7mdxsbVwcXJAXG99k1J28rZN9Y1afH9Uv7HmPaM4omppTE5w9uzZg9deew0BAQFQKBRQKBQYO3Ysli1bhmeeecYcMZINUDk7Yni4/o89uxp3jb64uBAAcO9ortSwBsZtG9KY4NiiX04UoKq+CWE+rkhs3paGLIfJCY5WqzV2Lvb398f58+cBAL179zYu6SYCgGuM+1IxwemKdQdzodUJjO7ji+hA1ptZA0OCsyejBE1anczRkLkZpqdmjAyHQsFpd0tjcoITGxuL5ORkAEB8fDyWL1+OXbt24bXXXkNkZKTJAZLtSGyR4Bhqtahj9J2L9X9M2bnYegwO9YKHiyOq6puQeoGbAduS8+V12JNRAkkC7orjwglLZHKC89JLL0Gn038yeeONN5CdnY1x48bhp59+wj//+U+TAyTbMbyXN5SOChRXa5BeWC13OFZle1oR8srr4OXK4mJr4qCQMEDtCQB8zduY307pe1GN7OWDUG/29rJEJq+imjJlivHfkZGROHnyJEpLS+Hj48OVMtSK0tEBo/r4Ymd6MXafLUFMEKdZOmotOxdbrahAd+zPKkXaBSY4tmRTij7BmTyo/Q76JC+TR3Aux9fXl8kNXVaCcZqKRZcd1bK4eFY8i4utTUzz8nCO4NiOirpG4+bBkwZyRNVSdUuCQ9QeQx3O3oxS6HSsw+mIb1hcbNUM/W/SClmDYyu2phaiSScQE+iOCH+3q38DyYIJDvWowaFecFc6oqKuESfzK+UOx+JpdQJfNRcX38vRG6sUE6RPcLJLatHQxJVUtmDTSf301KSBnJ6yZExwqEc5OigwOsIXAPvhdETL4uKpscFyh0NdoPZ0gbvSEU06geySGrnDIRNpmrTYlloEgAmOpWOCQz0ukXU4HfZfFhdbPUmSEGWcpmIdjrXbc7YE1ZomBHooMTTMW+5w6ArMkuDs2LED999/PxISEpCXp29b/fnnn2Pnzp3mOD3ZGEOh8f7MUjSy+Vm7LlTWYzM7F9sEFhrbjqTm6akbBgaxuZ+FMznBWb9+PaZMmQJXV1ccOXIEGo0GAFBVVYWlS5eaHCDZngFqT3irnFDToEXyuQq5w7FY6w7oi4tH9fHhknorF80RHJug0wljgjOZ01MWz+QE54033sAHH3yAVatWwcnJyXg8MTERhw8fNvX0ZIMUCgkJkYbVVKzDuZyWxcXsXGz9OIJjG5LzKlBYpYG70tE4Ek2Wy+QEJzU1FePHj29z3NPTE+Xl5aaenmwU++Fc2Q4WF9sUwwjO2aJqaNkewWolnSwAAEzoFwClI2viLJ3JCU5wcDDS09PbHN+5cyf3oqJ2GQqND2aVob5RK3M0lsfQuXj6iFAWF9uAMB8VlI4KNDTpcK6sVu5wqIuM3Ys5PWUVTE5wnnjiCcyfPx/79u2DJEk4f/48vvzyS7zwwguYO3euOWIkGxQV4I4ADyU0TTocySmXOxyL0rK4eNZoTk/ZAgeFhMiA5jocbtlglTKLa5BWWA1HhYSJ/QLlDoc6wOS9qBYuXIiKigpce+21qK+vx/jx46FUKvHCCy9g3rx55oiRbJAkSUiM8sP/HT2PPWeLOZ/dgqFzMYuLbUtMoDtO5VcivagaN4AjANbGMD11TaQfvFydrvJosgRmWSb+5ptvori4GPv378fevXtRVFSE119/3RynJht2sR8OC40NtDqB/+5v7lzM0RubYig05giOdUpi92KrY3KC89BDD2Hz5s1wdXVFXFwcRo8eDXd3d3PERjYuMcofAHA0txw1miaZo7EMLYuLbxrM4mJbEm1cScU9qaxNcbUGh7LLADDBsSYmJzglJSW4+eabERYWhueffx5Hjx41Q1hkD8J9VQjzcUWTTuBg8x8Pe/ff/SwutlWGPanSC6shBFdSWZPfTxVCJ4DYUE+EeLvKHQ51kMkJzsaNG1FQUIBXXnkFhw4dwsiRIzFw4EAsXboUWVlZZgiRbJmhHw6Xi+uLi387ZehczOkpW9Pbzw2OCgk1DVrkV9TLHQ51gnFzzQFqmSOhzjBLDY63tzcef/xxbN26FdnZ2XjooYfw+eefIzo62hynJxuWGK1PcLjx5sXi4rjePujL4mKb4+SgQB9/NwBs+GdNahuasCNNv7nm5EGcnrImZt1ss7GxEQcPHsS+ffuQlZWFoCC+GOjKEiL1dTgn8ipQUdsoczTy0bUoLmbnYtsVHcAtG6zNjrRiaJp0CPNxRX81P3hYE7MkOFu2bMFjjz2GoKAgzJ49Gx4eHvjhhx+Qm5trjtOTDVN7uSAywA06AezLtN9RnB3pxcgrr4OniyOLi21Yyzocsg4Xm/upIUncXNOamNwHJywsDCUlJZgyZQo+/PBD3HrrrXBxcTFHbGQnEqP8kFFUg91nSzB5kH3Oca/dlw0AmD4ijMXFNowrqaxLk1aH309zebi1MjnBefnllzFjxgz4+PiYIx6yQ4lR/vhib47dbrxZ2KK4mNNTtq3lruJCCI4IWLhD2WUoq22Et8oJo/rwPc7amDxF9fjjjzO5IZNc07yS6nRBFYqrNTJH0/O+OXSOxcV2IirAHZIElNc2oqSmQe5w6CoMq6eu6x8IRwezlqxSD+jSCM5zzz2H119/HW5ubnjuueeu+NgVK1Z0KTCyH75uzuiv9sDpgirszSjBLUNC5A6px+iLi/W9b7g03Pa5ODkg3EeFnNJapF2ohr+7Uu6QqB1CCGP3Ym6uaZ26lOAcOXIEjY2Nxn+3h8Ov1FGJUf44XVCF3WftK8HZkV6Mc2X64uKbh7C42B7EBLojp7QW6UXV3IPNgqVeqEJOaS2UjgqM7xsgdzjUBV1KcLZs2XLZfxN1VWKUH1bvyrS7fjj/3WfoXMziYnsRHeiOzacLkX6BhcaWLKl59dTYaH+onE0uVyUZmDypmJOT027b8ZycHFNPT3ZidKQvFBKQWVyD/Io6ucPpEYWV9Ug6pf8jyukp+2FcSVXEpeKWzPC7yeZ+1svkBCciIgJFRUVtjpeUlCAiIsLU05Od8HRxwuAwbwD209XYUFw8srcP+rGBmN2I5q7iFi+/og7J5yogScB1/ZngWCuTE5z2ljpWV1ezHw51SmKUYV8q209wWhYXz+LojV0xJDiFVRpU1Nlv925L9ltzcfGIXj4I8GAhuLXq8sSiYfWUJEn4y1/+ApVKZbxPq9Vi3759GDZsmMkBkv1IiPTD+1vPYs/ZEpvvEbKTxcV2y8PFCcFeLsivqEd6YTVG9mabDUuziaunbEKXExzD6ikhBI4fPw5nZ2fjfc7Ozhg6dCheeOEF0yMkuxHXxwdODhLyyuuQU1qL3n5ucofUbdayuNiuRQe6Nyc4VUxwLExFXaNxmpzdi61blxMcw+qphx56CP/85z/h4cEaAjKNytkRw8N9sD+rFLvPlthsgqPvXMziYnsWHeiOHWnF3JPKAm1NLUSTTiA60B2RzZujknUyuQYnJiYG33zzTZvjq1evxttvv23q6cnOJNhBHc43h86hicXFdq3llg1kWQzN/Th6Y/1MTnA++ugj9O/fv83xQYMG4YMPPjD19GRnDIXGhjocW6PTCXx1gJ2L7V1MoD6x5QiOZdE0abE1Vb8qmAmO9TM5wSkoKEBwcNsiyYCAAOTn55t6erIzw3p5w8VJgeJqjU3+8d+ZXozc0jp4uDji5sEsLrZXhhGcc2V1qG1okjkaMtibUYpqTRMCPJQY1ty2gqyXyQlOeHg4du3a1eb4rl27EBJiPy33yTyUjg6I6+0LwDanqQxLw+8cEQZXZxYX2ytfN2f4uekXZmQU1cgcDRkknSwAANwwIAgKhe2u4rQXJic4jz76KBYsWIA1a9YgOzsb2dnZWL16NZ599lk89thj5oiR7MzFOpximSMxr8KqeuP8/j2jw2WOhuQWZazD4ZYNlkCna7G5JrsX2wSTN9hYuHAhSktLMXfuXDQ0NAAAXFxc8OKLL2LRokUmB0j2x1CHszejFFqdgIONfJL65qC+uHhEL2/0V3vKHQ7JLCbQHfszS9nR2EIcz6vAhUoN3JwdjH+DyLqZnOBIkoS3334bf/nLX3Dq1Cm4uroiJiYGSiW7P1LXDA71grvSERV1jTiVX4nYUC+5QzJZy+LiWfG9ZY6GLEGMYU8qG6w1s0aG0ZuJ/QKhdOT0sS0weYrKwN3dHaNGjUJsbCyTGzKJo4MC8RGGOhzbmKbadZbFxdRaNFdSWZRNzfU3XD1lO8yS4OzYsQP3338/EhMTkZeXBwD4/PPPsXPnzk6f6/3338eQIUPg6ekJT09PJCQk4OeffzbeL4TAkiVLEBISAldXV0ycOBEpKSnm+DHIgiS0WC5uC4ydi4eHsriYAAAxQfoRnOzSWmiatDJHY9+yimtw5kI1HBQSru0XKHc4ZCYmJzjr16/HlClT4OrqisOHD0Oj0QAAqqqqsHTp0k6fLywsDG+99RYOHjyIgwcP4rrrrsPtt99uTGKWL1+OFStWYOXKlThw4ADUajUmTZqEqioW6tmSxCh/AMD+zFI0anUyR2OawsqLxcX3xrP3DekFeijhoXSEVieQVVwrdzh2zfD7eU2kL7xUTjJHQ+ZicoLzxhtv4IMPPsCqVavg5HTxhZGYmIjDhw93+ny33norbrrpJvTt2xd9+/bFm2++CXd3d+zduxdCCLz77rtYvHgxpk+fjtjYWHz66aeora3F2rVrTf1RyIL0V3vAR+WEmgYtks9VyB2OSb46kGvsXMziYjKQJAnRQazDsQTG7sUDOD1lS0xOcFJTUzF+/Pg2xz09PVFeXm7SubVaLb766ivU1NQgISEBmZmZKCgowOTJk42PUSqVmDBhAnbv3n3Fc2k0GlRWVra6keVSKCRcE2mYprLeOpwmrc7Y++aBa1hcTK1FB3CpuNxKqjU4mF0KAJg0SC1zNGROJic4wcHBSE9Pb3N8586diIyM7NI5jx8/Dnd3dyiVSjz55JPYsGEDBg4ciIICfRFYUFDrLDsoKMh4X3uWLVsGLy8v4y08nH1ILF2iDexLtfl0IfIr6uHr5oypg/nHk1qL4QiO7DafLoROAINCPBHq7Sp3OGRGJic4TzzxBObPn499+/ZBkiScP38eX375JV544QXMnTu3S+fs168fjh49ir179+IPf/gDZs+ejZMnTxrvl6TWfVGEEG2OXWrRokWoqKgw3nJzc7sUG/WchOY6nIPZZahvtM4izC/2ZgMAZsaFc+kptcE9qeTHzTVtl1ka/VVUVODaa69FfX09xo8fD6VSiRdeeAHz5s3r0jmdnZ0RHR0NAIiLi8OBAwfwj3/8Ay+++CKAtvtfFRYWthnVuZRSqeTydSsTFeCGQA8lCqs0OJJTblxZZS0yi2uwI60YkgTcx+JiugzDnlQZRTVo0urg6GC2zh3UAXUNWuxI02+uOXkgR1htTZd+m5KTk6HTXVzZ8uabb6K4uBj79+/H3r17UVRUhNdff91sQQohoNFoEBERAbVajaSkJON9DQ0N2LZtGxITE832fGQZJElqsbu49dXhGEZvJvYNQLivSuZoyBKFervCxUmBBq0OuWV1codjd3akFaG+UYdQb1cMCPaQOxwysy4lOMOHD0dxsf4NJzIyEiUlJVCpVIiLi8Po0aPh7u7e5YD+/Oc/Y8eOHcjKysLx48exePFibN26Fffddx8kScKCBQuwdOlSbNiwASdOnMCcOXOgUqkwa9asLj8nWa4EK63DOZVfic/2ZAEAHkzoI2ssZLkUCglRhkLjCyw07mmbWuw9dbUyB7I+XZqi8vb2RmZmJgIDA5GVldVqNMdUFy5cwAMPPID8/Hx4eXlhyJAh+OWXXzBp0iQA+imxuro6zJ07F2VlZYiPj8emTZvg4cHs2xYZ+uEczS1HjaYJbkqTZ1W7XaNWh+fXHUOjVuCGAUGY2C9A7pDIgsUEuiPlfCXSi6ox+eoPJzPR6gR+P10IgPU3tqpL7xZ33nknJkyYgODgYEiShLi4ODg4XL6AMiMjo1Pn/vjjj694vyRJWLJkCZYsWdKp85J1CvdVIczHFefK6nAgqxQTraDL6Mrf03EyvxLeKicsnR7LT4Z0RYY6nHRuutmjDmWXobSmAV6uThjdx1fucKgbdCnB+eijjzB9+nSkp6fjmWeewWOPPcYRFOo2iVF+WHfwHPacLbH4BOdEXgXe26Jvm/Da7bEI9HCROSKydMY9qYqY4PSkTSn61iLX9w9kcbeN6lKCk5ycjMmTJ+PGG2/EoUOHMH/+fCY41G0So/yx7uA5i6/D0TRp8fy6Y2jSCUyNVePWIdxUk64uusWu4jqdgELBEb/uJoRA0ikuD7d1JhcZb9u2DQ0NDWYNiqglQ6FxyvkKVNQ2yhxN+/65OQ2pF6rg5+aMN6Zxaoo6prefCk4OEmobtDhfwZVUPeHMhWpkl9TC2VGB8X1ZI2erupTgGIqMAZi9yJjoUkGeLogMcINOAPsyLXMU52huOd7fehYA8Ma0WPi5s+cSdYyTgwIR/m4A2PCvpySd1E9PjY32t4qFC9Q1FldkTHQ5iVF+yCiqwe6zJZhsYfvF1Ddq8fy6o9AJ4LahIZg6mFNT1DnRge44c6Ea6YXVFl9nZgsM3Ysnc3rKprHImKxCYpQ/vtibgz0WWIezIukMzhbVIMBDiVdvGyR3OGSF9IXGBRzB6QEFFfU4dq4CkgRcz93DbVqXx+ZuvPFGAGCRMfUIw87iqReqUFSlQYCHZUwBHcouxaod+lHKZXcMho+bs8wRkTUyFBqnMcHpdobi4uHh3hbzd4S6R5fXxt10002oqKjAmjVr4OHhgTfffBPl5eXG+0tKSjBw4EBzxEgEXzdnDAj2BADszbCMUZy6Bi1e+CYZQgB3jgjDDRzupi6KabGSSgghczS2zbA83NKmusn8upzg/Prrr9BoNMav3377bZSWlhq/bmpqQmpqqmnREbVg3JfKQhKc5b+eRmZxDdSeLnj5Vibz1HUR/m5QSEBFXSOKqjVX/wbqksr6RuMHJC4Pt31dTnAu/ZTBTx3U3S5uvCl/grM3owRrdmUBAN66czC8XJ3kDYismouTA3o1b8jKjsbdZ2tqERq1AlEBbsY9wMh2sX0jWY1REb5QSEBmcQ3Ol8vXL6RG04Q/fnsMAHDPqHCueiGzMDb8Y0fjbmNYPTVpIKen7EGXExxJkto0MmNjM+pOni5OGBzmDUDeUZxlP59CbmkdQr1dsfjmAbLFQbbFsGVDGkdwukVDkw5bubmmXenyKiohBObMmQOlUl+FXl9fjyeffBJubvqGVS3rc4jMJTHKD8dyy7H7bAnuHBnW48+/M60YX+zNAQAsv2sIPFw4NUXm0bLQmMxvb0YJqjRN8HdXYni4t9zhUA/ocoIze/bsVl/ff//9bR7z4IMPdvX0RJeVGOWH97eexZ6zxRBC9OioYVV9I15cnwwAeOCa3hgT7d9jz022j0vFu9fF6alA7vdlJ7qc4KxZs8accRB1SFxvXzg5SDhfUY+c0lr09nPrsed+88dTyCuvQy9fFf40tX+PPS/Zh6jmBKe4WoPy2gZ4q9hTyVyEEC26F7P+xl6wyJisiquzA4b38gGAHt1dfEtqIb46kAsA+OtdQ7h/DZmdu9IRIV4uADhNZW7H8ypQUFkPlbODcfNesn1McMjqJDR3Ne6pBKeithF/ap6aemhMH8RH8g8kdY/oIH2hMRMc8zKM3kzsFwAXp8vvm0i2hwkOWZ2L/XCKe6T/0qv/S8GFSg0i/N2wcAqnpqj7RAewDqc7bEox1N9w9ZQ9YYJDVmdYL2+4OClQXN3Q7W8ESScv4LvDeVBIwN9mDIGrMz/9UfeJCWKCY27ZJTVIvVAFB4WEa9mzyq4wwSGro3R0wKg+vgCA3enF3fY8ZTUNWPTdcQDAY+MiMbK3b7c9FxFwcSXVWSY4ZmOYnoqP8GXhtp1hgkNWyVAo2J11OK9sTEFxtQbRge54dlLfbnseIgPDFFVeeR1qNE0yR2MbNp3k9JS9YoJDVikxSt+DZl9mKbQ689fh/Hw8HxuPnYeDQsI7M4ayMJF6hI+bM/zd9c1Tz3LLBpOV1jTgYJZ+E2gmOPaHCQ5ZpdgQT3goHVFR14hT+ZVmPXdxtQaLvz8BAHhyQiSGsusp9aDoQH1vJ27ZYLrNpy5AJ4CBwZ4I81HJHQ71MCY4ZJUcHRQYHdFch3PWfHU4Qgj85fsTKK1pQH+1B565PsZs5ybqiJjmPam46abpkjg9ZdeY4JDV6o46nB+S8/HziQI4KiT8bcZQKB05NUU9y7hlA0dwTFLXoMX2tCIAwORBTHDsERMcslqGOpz9maVo1OpMPl9hVT1e/j/91NS866IRG+pl8jmJOuvipptVMkdi3XamF6O+UYdQb1cMDPaUOxySARMcslr91R7wUTmhtkGL5HPlJp1LCIE/f3cc5bWNGBTiiaeujTZPkESdZBjBySmtRX2jVuZorNemlAIA+umpntyUlywHExyyWgqFZJym2mPiNNV3h/Pw26lCODlIeGfmUDg58FeD5BHgoYSniyN0AsgsrpE7HKuk1Qn8froQADCZ9Td2i3/FyaolNE9TmVKHU1BRjyU/pAAAFtzQF/3VHM4m+UiShBjuSWWSwzllKKlpgKeLI0ZFsEGnvWKCQ1bNsPHmweyyLg3nCyHw4vpkVNU3YWi4N54YH2nuEIk6jXtSmcYwPXX9gCCOxtox/s+TVYsKcEOghxINTToczinr9PevO5iLbWeK4OyowDszhsCRfwzJAhj2pOKWDZ0nhGD3YgLABIesnCRJLXYX79w01bmyWrz+v1MAgBcm90V0c/8RIrlFGZaKcyVVp6UVViO7pBbOjgqM7xsgdzgkIyY4ZPUSu1CHY5iaqtY0YWRvHzwyllNTZDkMS8Uzi2vQZIYWCPbE0NxvTJQf3JWOMkdDcmKCQ1bPsJLqWG55hzco/GJfDnall8DFSYG/zRgKBwWXkZLlCPFyhauTAxq1AtmltXKHY1UM01OTB6lljoTkxgSHrF64rwrhvq5o0gkcaN5Y70pySmqx7Cf91NSLN/ZHhL9bd4dI1CkKhcSOxl1wobIex3LLIUnA9QMC5Q6HZMYEh2yCYTXV1epwdDqBP357DLUNWsRH+GJ2Qp8eiI6o8wwJDncV7zjD9NSwcG8EerjIHA3JjQkO2YSO1uF8uicL+zJLoXJ2wF/vGgoFp6bIQl0cwWGhcUcZp6cGcnqKmOCQjTDU4Zw4X4GK2sbLPiajqBpv/3IaALDopgHo5afqsfiIOsu4JxVHcDqkqr4Re84WA+DycNJjgkM2IcjTBVEBbhAC2JvZdhRHqxP447fJqG/UYWy0P+6P7yVDlEQdF23cdLMaOp2QORrLt+1MERq1ApEBbsZrR/aNCQ7ZDMM01eXqcD7emYFD2WVwVzri7buGcPM9sni9fFVwdlCgvlGHvPI6ucOxeJtS2NyPWmOCQzajvYZ/6YVV+NumMwCAv9wyAKHerj0eG1FnOToojCv8uCfVlTU06bAllZtrUmtMcMhmXNO8kir1QhWKqjQAgCatDs+vO4aGJh0m9gvAzLhwOUMk6pToIHY07oh9mSWoqm+Cv7sSw8J95A6HLAQTHLIZPm7OGBCs3wl8b4Z+FOfD7Rk4dq4Cni6OeGs6p6bIuhg23eQIzpUZloffMCCQTTvJiAkO2RTDNNXus8U4XVCJd3/TT00tuW0Q1F7si0HWJSaIu4pfjRDCmOBMHsTpKbqICQ7ZFMPmehuO5OHptUfQqBW4YUAQ7hgeKnNkRJ3XciWVEFxJdTkn8iqRX1EPlbODcaEBEWCBCc6yZcswatQoeHh4IDAwENOmTUNqamqrxwghsGTJEoSEhMDV1RUTJ05ESkqKTBGTJRkf44+ESD/UN+qQVlgNb5UTlk6P5dQUWaUIfzcoJKCqvgmFzXVl1FrSyQIAwPiYALg4OcgcDVkSi0twtm3bhqeeegp79+5FUlISmpqaMHnyZNTU1Bgfs3z5cqxYsQIrV67EgQMHoFarMWnSJFRVsRDP3kmShKXTB8PZUf/Sfu32WLZsJ6uldHRAHz+upLqSTZyeonZY3F7yv/zyS6uv16xZg8DAQBw6dAjjx4+HEALvvvsuFi9ejOnTpwMAPv30UwQFBWHt2rV44okn5AibLEiEvxu+eCQehVX1uHlwsNzhEJkkKtAdGcU1SLtQhTHRnIJpKaekFqcLquCgkHBdf26uSa1Z3AjOpSoqKgAAvr6+AIDMzEwUFBRg8uTJxscolUpMmDABu3fvbvc8Go0GlZWVrW5ku0ZH+OKWISGcmiKrxy0b2repeXpqdB9feKucZY6GLI1FJzhCCDz33HMYO3YsYmNjAQAFBfoXdFBQ6+HIoKAg432Xs2zZMnh5eRlv4eHsh0JElu/ipptMcC5lWD3F7sV0ORad4MybNw/Jycn473//2+a+Sz+ZCyGu+Gl90aJFqKioMN5yc3PNHi8RkbnFBHoAYA3OpUqqNTiQVQqACQ5dnsXV4Bg8/fTT2LhxI7Zv346wsDDjcbVaDUA/khMcfLG+orCwsM2oTktKpRJKpbL7AiYi6gZRgfoi45KaBpTWNMDXjVMxALD5VCF0AhgU4olwX5Xc4ZAFsrgRHCEE5s2bh++++w6///47IiIiWt0fEREBtVqNpKQk47GGhgZs27YNiYmJPR0uEVG3Ujk7GvdP4yjORYb6m8kD1TJHQpbK4hKcp556Cl988QXWrl0LDw8PFBQUoKCgAHV1+t10JUnCggULsHTpUmzYsAEnTpzAnDlzoFKpMGvWLJmjJyIyP0NHYyY4ejWaJmxPKwYATInl9BRdnsVNUb3//vsAgIkTJ7Y6vmbNGsyZMwcAsHDhQtTV1WHu3LkoKytDfHw8Nm3aBA8Pjx6Oloio+0UHuGNrahE33Wy27UwRGpp06O2nQr8g/t2ny7O4BKcj7cglScKSJUuwZMmS7g+IiEhmHMFp7dcU/fTUlEFqtoKgdlncFBUREbXWck8qe9fQpMPvpwsBAFPYvZiugAkOEZGFiw7QT8PkV9Sjqr5R5mjktSejBFX1TQjwUGJ4uI/c4ZAFY4JDRGThvFROCPDQt7k4W1RzlUfbLp1OYPXOTAD63jcKBaenqH1McIiIrECMsaOx/RYav7s5DdvOFEHpqMCDCb3lDocsHBMcIiIrEG3ne1JtSinAPzenAQCWTR+M/mpPmSMiS8cEh4jIChg33bTDPanSC6vx3LpjAIA5iX0wfUTYVb6DiAkOEZFViDbsSWVnIziV9Y14/LODqNY0IT7CF4tvHiB3SGQlmOAQEVkBwxRVTmkt6hu1MkfTM3Q6gWe/OoqM4hqEeLngvftGwMmBb1vUMXylEBFZAX93Z3irnCAEcNZORnH+sTkNm08XwtlRgQ8eGAl/d26YTB3HBIeIyApIkoToAPtp+LcppQD/MBQV3zEYQ8K85Q2IrA4THCIiK2EvWzZcWlR850gWFVPnMcEhIrISUXYwglNZ34jHP2dRMZmOCQ4RkZWIad45O81GExydTuC5r48io4hFxWQ6vnKIiKyEYSVVVnENGrU6maMxv39sTsNvp1hUTObBBIeIyEqEeLnAzdkBTTqB7BLb2pOKRcVkbkxwiIishCRJF7dssKFpKhYVU3dggkNEZEWijJtu2kaCw6Ji6i5McIiIrEhMoO0UGrcsKg5mUTGZGV9JRERWxJamqP75+8Wi4g9ZVExmxgSHiMiKGHYVP1tUDa1OyBxN1yWdvIB3f9MXFS9lUTF1AyY4RERWJNxXBWdHBTRNOuSV1ckdTpekF1bj2a+PAtAXFd/FomLqBkxwiIisiINCQqS/GwAgrbBK5mg6r2VR8WgWFVM3YoJDRGRlrLUOR19UfMxYVPxvFhVTN+Iri4jIyljrSip9UfEFFhVTj2CCQ0RkZQy7iltTgtOyqPjNabEsKqZuxwSHiMjKGKaozhZWQwjLX0nVsqh4dkJvzIgLlzcgsgtMcIiIrEwfPzc4KCRUa5pQUFkvdzhXVFXfiCdaFBW/dMtAuUMiO8EEh4jIyjg7KtDbTwXAsguNdTqBZ78+hrOGTsWzWFRMPYevNCIiKxRjBXtS/ev3dGNR8Qf3j0SAB4uKqecwwSEiskLGpeJFlpng/HbyAv7+2xkA+qLioeHe8gZEdocJDhGRFTIsFU+3wBGcs0UsKib5McEhIrJCljqCU1XfiMc/O4gqTRNG92FRMcmHCQ4RkRWKCnCHJAGlNQ0oqdbIHQ6A5k7F61oUFbNTMcmIrzwiIivk6uyAMB9XAJbT8O9fv6cj6SSLiskyMMEhIrJS0QGWsycVi4rJ0jDBISKyUjFBzYXGMic4LYuKH2RRMVkIJjhERFbKEkZwLi0q/guLislCMMEhIrJS0cZNN6tkef6WRcVqTxYVk2XhK5GIyEoZlopfqNSgsr6xx5+/VVHxAywqJsvCBIeIyEp5ujghyFOfVPT0NNXmUxeLit+YFothLComC8MEh4jIivV0R2OdTuCzPVmY++VhAPqi4pksKiYL5Ch3AERE1HXRge7YmV7cIx2ND2WX4bX/ncSx3HIAwA0DAllUTBaLCQ4RkRWLNu4q3n2FxufKavH2L6n44dh5AICbswP+OKUfZif2gSRJ3fa8RKZggkNEZMW6c0+qGk0TPth2Fh9tz4CmSQdJAmaODMfzU/oi0MPF7M9HZE5McIiIrFhMc4JzrqwOdQ1auDo7mHxOnU5g/eFz+OuvqSis0u9zFR+h73ETG+pl8vmJeoJFFhlv374dt956K0JCQiBJEr7//vtW9wshsGTJEoSEhMDV1RUTJ05ESkqKPMESEcnIz10JH5UThNB3FDbV/sxS3P7eLvzx22QUVmnQy1eFD+4fia8ev4bJDVkVi0xwampqMHToUKxcufKy9y9fvhwrVqzAypUrceDAAajVakyaNAlVVfI0uyIikpNxJZUJS8VzS2vx1JeHMfPDPTieVwF3pSMWTe2PpOfG48ZYNWttyOpY5BTV1KlTMXXq1MveJ4TAu+++i8WLF2P69OkAgE8//RRBQUFYu3YtnnjiiZ4MlYhIdlGB7tifVdqljsZV9Y3499az+HhnJhqadFBIwD2je+G5SX3h787GfWS9LDLBuZLMzEwUFBRg8uTJxmNKpRITJkzA7t27201wNBoNNBqN8evKyspuj5WIqCcY6nA6M4Kj1Ql8czAXf9t0BsXV+r+NY6L98NLNAzEg2LNb4iTqSVaX4BQUFAAAgoKCWh0PCgpCdnZ2u9+3bNkyvPrqq90aGxGRHGKMe1J1LMHZc7YEr/3vJE7l6z/oRfi74c83DcANAwI5FUU2w+oSHINLfwmFEFf8xVy0aBGee+4549eVlZUID2f3TSKyfoal4tkltWho0sHZ8fLlldklNVj60yn8mnIBAODh4oj518fgwYQ+7X4PkbWyugRHrVYD0I/kBAcHG48XFha2GdVpSalUQqnkfDIR2R61pwvclY6o1jQhq6QGfYM8Wt1fWd+Ilb+nY82uTDRqBRwUEmaN7oVnJ/WFr5uzTFETdS+rS9kjIiKgVquRlJRkPNbQ0IBt27YhMTFRxsiIiOQhSRKijB2NL05TNWl1+GJvNq7961Z8tD0DjVqBcTH++Hn+OLw+LZbJDdk0ixzBqa6uRnp6uvHrzMxMHD16FL6+vujVqxcWLFiApUuXIiYmBjExMVi6dClUKhVmzZolY9RERPIZEuqFY7nlWH/4HG4eEoydacV4/X8nkdq8hUNkgBv+cvNATOwXwDobsguSEELIHcSltm7dimuvvbbN8dmzZ+OTTz6BEAKvvvoqPvzwQ5SVlSE+Ph7vvfceYmNjO/wclZWV8PLyQkVFBTw9uWKAiKxbZnENblixDVqdwDWRvtibUQoA8HJ1wrM3xOC+a3rDycHqBu2J2ujo+7dFJjg9gQkOEdmaP35zDN8cOgcAcFBIeOCa3lhwQwy8VZyKItvR0fdvi5yiIiKizlswqS9OnK9EqLcr/jS1H6IDPa7+TUQ2igkOEZGNCPV2xc/zx8kdBpFF4IQsERER2RwmOERERGRzmOAQERGRzWGCQ0RERDaHCQ4RERHZHCY4REREZHOY4BAREZHNYYJDRERENocJDhEREdkcJjhERERkc5jgEBERkc1hgkNEREQ2hwkOERER2RwmOERERGRzHOUOQC5CCABAZWWlzJEQERFRRxnetw3v4+2x2wSnpKQEABAeHi5zJERERNRZVVVV8PLyavd+u01wfH19AQA5OTlXvED2pLKyEuHh4cjNzYWnp6fc4VgEXpO2eE1a4/Voi9ekLV6Ttrp6TYQQqKqqQkhIyBUfZ7cJjkKhLz/y8vLii+0Snp6evCaX4DVpi9ekNV6PtnhN2uI1aasr16QjAxMsMiYiIiKbwwSHiIiIbI7dJjhKpRKvvPIKlEql3KFYDF6TtnhN2uI1aY3Xoy1ek7Z4Tdrq7msiiautsyIiIiKyMnY7gkNERES2iwkOERER2RwmOERERGRzmOAQERGRzbHLBOff//43IiIi4OLigpEjR2LHjh1yh9RjlixZAkmSWt3UarXxfiEElixZgpCQELi6umLixIlISUmRMWLz2759O2699VaEhIRAkiR8//33re7vyDXQaDR4+umn4e/vDzc3N9x22204d+5cD/4U5nW1azJnzpw2r5trrrmm1WNs6ZosW7YMo0aNgoeHBwIDAzFt2jSkpqa2eoy9vU46ck3s7XXy/vvvY8iQIcZGdQkJCfj555+N99vbawS4+jXpydeI3SU4X3/9NRYsWIDFixfjyJEjGDduHKZOnYqcnBy5Q+sxgwYNQn5+vvF2/Phx433Lly/HihUrsHLlShw4cABqtRqTJk1CVVWVjBGbV01NDYYOHYqVK1de9v6OXIMFCxZgw4YN+Oqrr7Bz505UV1fjlltugVar7akfw6yudk0A4MYbb2z1uvnpp59a3W9L12Tbtm146qmnsHfvXiQlJaGpqQmTJ09GTU2N8TH29jrpyDUB7Ot1EhYWhrfeegsHDx7EwYMHcd111+H22283JjH29hoBrn5NgB58jQg7M3r0aPHkk0+2Ota/f3/xpz/9SaaIetYrr7wihg4detn7dDqdUKvV4q233jIeq6+vF15eXuKDDz7ooQh7FgCxYcMG49cduQbl5eXCyclJfPXVV8bH5OXlCYVCIX755Zcei727XHpNhBBi9uzZ4vbbb2/3e2z9mhQWFgoAYtu2bUIIvk6EaHtNhODrRAghfHx8xH/+8x++RlowXBMhevY1YlcjOA0NDTh06BAmT57c6vjkyZOxe/dumaLqeWlpaQgJCUFERATuueceZGRkAAAyMzNRUFDQ6voolUpMmDDBbq5PR67BoUOH0NjY2OoxISEhiI2NtenrtHXrVgQGBqJv37547LHHUFhYaLzP1q9JRUUFgIub9PJ10vaaGNjr60Sr1eKrr75CTU0NEhIS+BpB22ti0FOvEbvabLO4uBharRZBQUGtjgcFBaGgoECmqHpWfHw8PvvsM/Tt2xcXLlzAG2+8gcTERKSkpBivweWuT3Z2thzh9riOXIOCggI4OzvDx8enzWNs9XU0depUzJgxA71790ZmZib+8pe/4LrrrsOhQ4egVCpt+poIIfDcc89h7NixiI2NBcDXyeWuCWCfr5Pjx48jISEB9fX1cHd3x4YNGzBw4EDjm7E9vkbauyZAz75G7CrBMZAkqdXXQog2x2zV1KlTjf8ePHgwEhISEBUVhU8//dRY6GXP18egK9fAlq/T3Xffbfx3bGws4uLi0Lt3b/z444+YPn16u99nC9dk3rx5SE5Oxs6dO9vcZ6+vk/auiT2+Tvr164ejR4+ivLwc69evx+zZs7Ft2zbj/fb4GmnvmgwcOLBHXyN2NUXl7+8PBweHNllgYWFhmyzbXri5uWHw4MFIS0szrqay5+vTkWugVqvR0NCAsrKydh9j64KDg9G7d2+kpaUBsN1r8vTTT2Pjxo3YsmULwsLCjMft+XXS3jW5HHt4nTg7OyM6OhpxcXFYtmwZhg4din/84x92/Rpp75pcTne+RuwqwXF2dsbIkSORlJTU6nhSUhISExNlikpeGo0Gp06dQnBwMCIiIqBWq1tdn4aGBmzbts1urk9HrsHIkSPh5OTU6jH5+fk4ceKE3VynkpIS5ObmIjg4GIDtXRMhBObNm4fvvvsOv//+OyIiIlrdb4+vk6tdk8ux9dfJ5QghoNFo7PI10h7DNbmcbn2NdKok2QZ89dVXwsnJSXz88cfi5MmTYsGCBcLNzU1kZWXJHVqPeP7558XWrVtFRkaG2Lt3r7jllluEh4eH8ed/6623hJeXl/juu+/E8ePHxb333iuCg4NFZWWlzJGbT1VVlThy5Ig4cuSIACBWrFghjhw5IrKzs4UQHbsGTz75pAgLCxO//fabOHz4sLjuuuvE0KFDRVNTk1w/lkmudE2qqqrE888/L3bv3i0yMzPFli1bREJCgggNDbXZa/KHP/xBeHl5ia1bt4r8/Hzjrba21vgYe3udXO2a2OPrZNGiRWL79u0iMzNTJCcniz//+c9CoVCITZs2CSHs7zUixJWvSU+/RuwuwRFCiPfee0/07t1bODs7ixEjRrRa5mjr7r77bhEcHCycnJxESEiImD59ukhJSTHer9PpxCuvvCLUarVQKpVi/Pjx4vjx4zJGbH5btmwRANrcZs+eLYTo2DWoq6sT8+bNE76+vsLV1VXccsstIicnR4afxjyudE1qa2vF5MmTRUBAgHBychK9evUSs2fPbvPz2tI1udy1ACDWrFljfIy9vU6udk3s8XXy8MMPG99LAgICxPXXX29MboSwv9eIEFe+Jj39GpGEEKJzYz5EREREls2uanCIiIjIPjDBISIiIpvDBIeIiIhsDhMcIiIisjlMcIiIiMjmMMEhIiIim8MEh4iIiGwOExwiIiKyOUxwiMjiLFmyBMOGDbviY7KysiBJEo4ePdojMXXWxIkTsWDBApPOIUkSJEmCt7f3FR/XkevVWYbrK0mS2c9N1BOY4BBZgMLCQjzxxBPo1asXlEol1Go1pkyZgj179rR63O7du3HTTTfBx8cHLi4uGDx4MN555x1otdpWjzO8MUmSBDc3N8TExGDOnDk4dOhQm+c+fvw4JkyYAFdXV4SGhuK1116D3A3OX3jhBWzevNn49Zw5czBt2jT5ApLRmjVrcObMmR5/3vDwcOTn5+P555/v8ecmMgcmOEQW4M4778SxY8fw6aef4syZM9i4cSMmTpyI0tJS42M2bNiACRMmICwsDFu2bMHp06cxf/58vPnmm7jnnnvaJCVr1qxBfn4+UlJS8N5776G6uhrx8fH47LPPjI+prKzEpEmTEBISggMHDuBf//oX/va3v2HFihU99rNfjru7O/z8/GSNwVJ4e3sjMDCwx5/XwcEBarUa7u7uPf7cRGZh8s5aRGSSsrIyAUBs3bq13cdUV1cLPz8/MX369Db3bdy4UQAQX331lfEYALFhw4Y2j33wwQeFh4eHKC0tFUII8e9//1t4eXmJ+vp642OWLVsmQkJChE6nu2wsr7zyyhU3XdTpdOLtt98WERERwsXFRQwZMkR88803xu83bOz522+/iZEjRwpXV1eRkJAgTp8+3eo5hg4d2u7zbdmyRWRmZgoAYv369WLixInC1dVVDBkyROzevbvd6yiEEO+8846IjY0VKpVKhIWFiT/84Q+iqqrKeP+aNWuEl5eX+OWXX0T//v2Fm5ubmDJlijh//rzxMY2NjeLpp58WXl5ewtfXVyxcuFA8+OCD4vbbbzc+ZsKECWL+/PnGrzUajfjjH/8oQkJChEqlEqNHjxZbtmy5Yqzt/T8uW7ZMBAYGCnd3d/Hwww+LF1980Xi9DFavXi369+8vlEql6Nevn3jvvfda3b9r1y4xdOhQoVQqxciRI8WGDRsEAHHkyJFWj2v5f0FkTZjgEMmssbFRuLu7iwULFrRKNFr67rvvBIB237z79u3b6s21vTfGI0eOCADi66+/FkII8cADD4jbbrut1WMOHz4sAIiMjIzLPldVVZXIz8833v72t78JlUpl3CX5z3/+s+jfv7/45ZdfxNmzZ8WaNWuEUqk0JnCGBCc+Pl5s3bpVpKSkiHHjxonExETjc7R8U62qqhIzZ84UN954o/E5NRqNMcHp37+/+N///idSU1PFXXfdJXr37i0aGxsvG7sQQvz9738Xv//+u8jIyBCbN28W/fr1E3/4wx+M969Zs0Y4OTmJG264QRw4cEAcOnRIDBgwQMyaNcv4mDfeeEP4+vqK7777Tpw6dUo8+eSTwtPT84oJzqxZs0RiYqLYvn27SE9PF3/961+FUqkUZ86caTfWy/0/fv3118LZ2VmsWrVKnD59WixevFh4eHi0SkI++ugjERwcLNavXy8yMjLE+vXrha+vr/jkk0+EEEJUVlYKX19fcf/994uUlBTx008/ib59+zLBIZvCBIfIAnz77bfCx8dHuLi4iMTERLFo0SJx7Ngx4/1vvfWWACDKysou+/233XabGDBggPHr9hKcuro6AUC8/fbbQgghJk2aJB577LFWj8nLy7tiMtXSnj17hIuLizFhqq6uFi4uLm2+95FHHhH33nuvEKL1CI7Bjz/+KACIuro6IUTbN9XZs2e3Sh6EEMYE5z//+Y/xWEpKigAgTp06ddXYDdatWyf8/PyMX69Zs0YAEOnp6cZj7733nggKCjJ+HRQUJP76178av25qahK9evVqN8FJT08XkiSJvLy8Vs99/fXXi0WLFrUb2+X+HxMSEsSTTz7Z6lh8fHyr6xUeHi7Wrl3b6jGvv/66SEhIEEII8f777ws/Pz/j9RZCiFWrVjHBIZvCGhwiC3DnnXfi/Pnz2LhxI6ZMmYKtW7dixIgR+OSTT1o9TrRT/CuEgLOz81Wfx/D9kiQZj7X8d3uPuZycnBxMmzYNL7zwAmbOnAkAOHnyJOrr6zFp0iS4u7sbb5999hnOnj3b6vuHDBli/HdwcDAAfbF1Z3X2PFu2bMGkSZMQGhoKDw8PPPjggygpKUFNTY3xMSqVClFRUa3OazhnRUUFLly4gNGjRxvvd3BwwMiRI9t9zsOHD0MIgb59+7a6Ltu2bWtzXa7m1KlTSEhIaHWs5ddFRUXIzc3FI4880uq53njjDeNzpaamYsiQIXBxcTF+X8ufh8gWOModABHpubi4YNKkSZg0aRJefvllPProo3jllVcwZ84cxMTEANC/uSUmJrb53tOnT3doKe+pU6cAABEREQAAtVqNgoKCVo8xvJEHBQW1e56amhrcdtttSEhIwGuvvWY8rtPpAAA//vgjQkNDW32PUqls9bWTk5Px34ZkyvD9ndGZ82RnZ+Omm27Ck08+iddffx2+vr7YuXMnHnnkETQ2Nl72nIbzXppctpcYXo5Op4ODgwMOHToEBweHVveZu4jX8LOvWrUK8fHxre4zPLcQolPxE1kjjuAQWaiBAwcaRxWmTJkCX19fvPPOO20et3HjRqSlpWHOnDlXPee7774LT09P3HDDDQD0n/y3b9+OhoYG42M2bdqEkJAQ9OnT57LnEELg/vvvh06nw+eff97qjXLgwIFQKpXIyclBdHR0q1t4eHgnfvrWnJ2d2yyF74qDBw+iqakJ77zzDq655hr07dsX58+f79Q5vLy8EBQUhP379xuPabVaHDlypN3vGT58OLRaLQoLC9tcF7Va3annHzBgAPbu3dvqWMuvg4KCEBoaioyMjDbPZUhs+/fvj+TkZGg0GuP3HTx4sFNxEFk6juAQyaykpAQzZszAww8/jCFDhsDDwwMHDx7E8uXLcfvttwMA3Nzc8OGHH+Kee+7B448/jnnz5sHT0xObN2/GH//4Rzz66KO46aabWp23vLwcBQUF0Gg0OHPmDD788EN8//33+Oyzz4yN42bNmoVXX30Vc+bMwZ///GekpaVh6dKlePnll9udolqyZAl+++03bNq0CdXV1aiurgagf+P38PDACy+8gGeffRY6nQ5jx45FZWUldu/eDXd3d8yePbtL16hPnz749ddfkZqaCj8/P3h5eXXpPFFRUWhqasK//vUv3Hrrrdi1axc++OCDTp/n6aefxrJlyxAdHY3+/fvjX//6F8rKytq9Zn379sV9992HBx98EO+88w6GDx+O4uJi/P777xg8eHCb/7srmT9/PmbPno24uDiMHTsWX375JVJSUhAZGWl8zJIlS/DMM8/A09MTU6dOhUajwcGDB1FWVobnnnsOs2bNwuLFi/H444/jT3/6E3JycvC3v/0NwNWnJomshlzFP0SkV19fL/70pz+JESNGCC8vL6FSqUS/fv3ESy+9JGpra1s9dvv27WLKlCnC09PTuGT6rbfeanNOw30AhIuLi4iKihKzZ88Whw4davPY5ORkMW7cOKFUKoVarRZLlixpd4m4EPri2ZbnN9xaLhP/xz/+Ifr16yecnJxEQECAmDJliti2bZsQ4mKRccuCacPqrszMTCFE28LWwsJCMWnSJOHu7t5mmXjLoljDkvsrLb9esWKFCA4OFq6urmLKlCnis88+axWPYZl4S4Yl1AaNjY1i3rx5wtPTU/j4+IgXX3xRzJgxQ9xzzz2trlPLVVQNDQ3i5ZdfFn369BFOTk5CrVaLO+64QyQnJ7cbK9opFn/zzTeFv7+/cHd3F7NnzxYLFy5sUwj85ZdfimHDhglnZ2fh4+Mjxo8fL7777jvj/bt27RJDhgwRzs7OYuTIkWLt2rUCQKvl+kKwyJislyQEJ16JrFF9fT1uv/125ObmYtu2bQgICJA7JLul0+kwYMAAzJw5E6+//rrZzitJEjZs2NAjXZy//PJLPPTQQ6ioqICrq6vx+JIlS/D9999b7JYYRO1hgkNkxerr6/Huu+8iJiYGd955p9zh2I3s7Gxs2rQJEyZMgEajwcqVK7FmzRocO3YMAwYMMNvzSJIEFxcX+Pn54dy5c2Y7LwB89tlniIyMRGhoKI4dO4Z58+Zh4sSJ+OKLLwDoV8kNHDgQDQ0NGDhwIBMcsjpMcIiIOik3Nxf33HMPTpw4ASEEYmNj8dZbb2H8+PFmfZ709HQA+tVPhgJhc1m+fDn+/e9/o6CgAMHBwZg2bRrefPNNqFQqAEBTUxOysrIA6FfAmVIkTiQHJjhERERkc7hMnIiIiGwOExwiIiKyOUxwiIiIyOYwwSEiIiKbwwSHiIiIbA4THCIiIrI5THCIiIjI5jDBISIiIpvz/9mVXflWdgQMAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import bctools.data as data\n", "from bctools.io import InstrumentResponse\n", "from bctools.util.coords import SpacecraftCoords\n", "import bctools.util.units as u\n", "\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "theta = np.linspace(0,360, 360)*u.deg\n", "energy = 100 * u.keV\n", "\n", "# The best way to work with them is to use them as context manager\n", "# The alternative to using context manager is to remember to call close() when you are done with the file.\n", "with InstrumentResponse(data.path/'sims/drm.h5') as irf: \n", "\n", " # You can get a list of the detector IDs\n", " print(\"Detectors: {}\".format(irf.detectors))\n", " \n", " # SQD0 zenith is pointing at 45d,45d in the spacecraft coordinates\n", " eff_area = []\n", " \n", " for angle in theta:\n", " \n", " coords = SpacecraftCoords(45*u.deg + angle, 45*u.deg)\n", " \n", " drm = irf.get_drm(\"SQD0\", coords = coords)\n", " \n", " \n", " eff_area += [drm.effective_area(energy)]\n", " \n", " \n", "fig,ax = plt.subplots()\n", "\n", "ax.plot(theta/u.deg, eff_area);\n", "\n", "ax.set_xlim(0,360)\n", "\n", "ax.set_xlabel(\"SQD0 zenith angle [deg]\")\n", "ax.set_ylabel(\"Effective area [cm$^2$]\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The effective area is not symmetric around 180$^\\circ$ because the *shadowing* effect from the other simulated components is taken into account.\n", "\n", "This plot is not very smooth because the underlying DRM used for this example has low resolution. You can see this with" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'all'" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "with InstrumentResponse(data.path/'sims/drm.h5') as irf: \n", "\n", " # You can get a list of the detector IDs\n", " print(\"Detectors: {}\".format(irf.detectors))\n", " \n", " # SQD0 zenith is pointing at 45d,45d in the spacecraft coordinates\n", " eff_area = []\n", " \n", " for angle in theta:\n", " \n", " coords = SpacecraftCoords(45*u.deg + angle, 45*u.deg)\n", " \n", " drm = irf.get_drm(\"SQD0\", coords = coords)\n", " \n", " \n", " eff_area += [drm.effective_area(energy)]" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['__array__',\n", " '__class__',\n", " '__delattr__',\n", " '__dict__',\n", " '__dir__',\n", " '__doc__',\n", " '__eq__',\n", " '__format__',\n", " '__ge__',\n", " '__getattribute__',\n", " '__gt__',\n", " '__hash__',\n", " '__init__',\n", " '__init_subclass__',\n", " '__le__',\n", " '__lt__',\n", " '__module__',\n", " '__ne__',\n", " '__new__',\n", " '__reduce__',\n", " '__reduce_ex__',\n", " '__repr__',\n", " '__setattr__',\n", " '__sizeof__',\n", " '__str__',\n", " '__subclasshook__',\n", " '__weakref__',\n", " '_drm',\n", " '_eff_area',\n", " 'axes',\n", " 'channel_effective_area',\n", " 'effective_area',\n", " 'energy_channels',\n", " 'fold_spectrum',\n", " 'from_gbm_rsp',\n", " 'from_gbm_rsp_file',\n", " 'load',\n", " 'photon_energies',\n", " 'to_gbm_rsp',\n", " 'to_gbm_rsp_file',\n", " 'write']" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dir(drm)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/imartin5/.local/lib/python3.7/site-packages/astropy/io/fits/column.py:649: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.\n", " array = np.array(array)\n" ] } ], "source": [ "drma.to_gbm_rsp().write(\"./\", \"coso.rsp\")" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Help on method write in module gbm.data.drm:\n", "\n", "write(directory, filename=None, **kwargs) method of gbm.data.drm.RSP instance\n", " Writes a single-DRM RSP object to a FITS file.\n", " \n", " Args:\n", " directory (str): Where to write the file\n", " filename (str, optional): The filename to write to\n", "\n" ] } ], "source": [ "help(drm.to_gbm_rsp().write)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Underlying HEALPix NSIDE = 2\n", "This means only 48 directions were simulated\n" ] } ], "source": [ "with InstrumentResponse(data.path/'sims/drm.h5') as irf:\n", " \n", " nside = irf.nside\n", " \n", "print(\"Underlying HEALPix NSIDE = {}\".format(nside))\n", "print(\"This means only {} directions were simulated\".format(12*nside*nside))" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.12" } }, "nbformat": 4, "nbformat_minor": 4 }