diff --git a/source_injector_Poisson.ipynb b/source_injector_Poisson.ipynb
new file mode 100644
index 0000000..9e47974
--- /dev/null
+++ b/source_injector_Poisson.ipynb
@@ -0,0 +1,1385 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "id": "a0226ddc",
+ "metadata": {},
+ "source": [
+ "# Source injector"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "c3cf1f63-9d4c-4922-9555-c6a33bdc6fc1",
+ "metadata": {},
+ "source": [
+ "The \"source injector\" is a cosipy module that will generate mocked binned data based on the detector response and a source hypothesis (provided by the users). This should result in the same output as simulating the source using MEGAlib, but be much quicker. MEGAlib is only needed when the event-by-event data is required, or to create the detector response itself. \n",
+ "\n",
+ "The goal of this notebook is to get an idea how the source injector will work in practice. We need to take it from here to something that is user friendly and compatible with the rest of the modules."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "12989d27-f8eb-4764-947f-e5197b13b3b5",
+ "metadata": {},
+ "source": [
+ "First, let's load all dependecies:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "id": "1863fe19-1d2b-4d9d-aacc-6e1eba99b882",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "
16:14:44 WARNING The naima package is not available. Models that depend on it will not be functions.py:48\n",
+ "available \n",
+ "
\n"
+ ],
+ "text/plain": [
+ "\u001b[38;5;46m16:14:44\u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m The naima package is not available. Models that depend on it will not be \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=928783;file:///Users/thomassiegert/.virtualenvs/cosipy/lib/python3.9/site-packages/astromodels/functions/functions_1D/functions.py\u001b\\\u001b[2mfunctions.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=720320;file:///Users/thomassiegert/.virtualenvs/cosipy/lib/python3.9/site-packages/astromodels/functions/functions_1D/functions.py#48\u001b\\\u001b[2m48\u001b[0m\u001b]8;;\u001b\\\n",
+ "\u001b[38;5;46m \u001b[0m \u001b[1;38;5;251mavailable \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[2m \u001b[0m\n"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "text/html": [
+ "
WARNING The GSL library or the pygsl wrapper cannot be loaded. Models that depend on it functions.py:69\n",
+ "will not be available. \n",
+ "
\n"
+ ],
+ "text/plain": [
+ "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m The GSL library or the pygsl wrapper cannot be loaded. Models that depend on it \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=367225;file:///Users/thomassiegert/.virtualenvs/cosipy/lib/python3.9/site-packages/astromodels/functions/functions_1D/functions.py\u001b\\\u001b[2mfunctions.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=685736;file:///Users/thomassiegert/.virtualenvs/cosipy/lib/python3.9/site-packages/astromodels/functions/functions_1D/functions.py#69\u001b\\\u001b[2m69\u001b[0m\u001b]8;;\u001b\\\n",
+ "\u001b[38;5;46m \u001b[0m \u001b[1;38;5;251mwill not be available. \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[2m \u001b[0m\n"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "/Users/thomassiegert/.virtualenvs/cosipy/lib/python3.9/site-packages/numba/core/decorators.py:262: NumbaDeprecationWarning: \u001b[1mnumba.generated_jit is deprecated. Please see the documentation at: https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-generated-jit for more information and advice on a suitable replacement.\u001b[0m\n",
+ " warnings.warn(msg, NumbaDeprecationWarning)\n"
+ ]
+ },
+ {
+ "data": {
+ "text/html": [
+ "
16:14:45 WARNING The ebltable package is not available. Models that depend on it will not be absorption.py:36\n",
+ "available \n",
+ "
\n"
+ ],
+ "text/plain": [
+ "\u001b[38;5;46m16:14:45\u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m The ebltable package is not available. Models that depend on it will not be \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=270317;file:///Users/thomassiegert/.virtualenvs/cosipy/lib/python3.9/site-packages/astromodels/functions/functions_1D/absorption.py\u001b\\\u001b[2mabsorption.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=135787;file:///Users/thomassiegert/.virtualenvs/cosipy/lib/python3.9/site-packages/astromodels/functions/functions_1D/absorption.py#36\u001b\\\u001b[2m36\u001b[0m\u001b]8;;\u001b\\\n",
+ "\u001b[38;5;46m \u001b[0m \u001b[1;38;5;251mavailable \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[2m \u001b[0m\n"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "text/html": [
+ "
WARNING We have set the min_value of K to 1e-99 because there was a postive transform parameter.py:704\n",
+ "
\n"
+ ],
+ "text/plain": [
+ "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m We have set the min_value of K to \u001b[0m\u001b[1;37m1e-99\u001b[0m\u001b[1;38;5;251m because there was a postive transform \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=798320;file:///Users/thomassiegert/.virtualenvs/cosipy/lib/python3.9/site-packages/astromodels/core/parameter.py\u001b\\\u001b[2mparameter.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=575920;file:///Users/thomassiegert/.virtualenvs/cosipy/lib/python3.9/site-packages/astromodels/core/parameter.py#704\u001b\\\u001b[2m704\u001b[0m\u001b]8;;\u001b\\\n"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "text/html": [
+ "
WARNING We have set the min_value of K to 1e-99 because there was a postive transform parameter.py:704\n",
+ "
\n"
+ ],
+ "text/plain": [
+ "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m We have set the min_value of K to \u001b[0m\u001b[1;37m1e-99\u001b[0m\u001b[1;38;5;251m because there was a postive transform \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=487302;file:///Users/thomassiegert/.virtualenvs/cosipy/lib/python3.9/site-packages/astromodels/core/parameter.py\u001b\\\u001b[2mparameter.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=443077;file:///Users/thomassiegert/.virtualenvs/cosipy/lib/python3.9/site-packages/astromodels/core/parameter.py#704\u001b\\\u001b[2m704\u001b[0m\u001b]8;;\u001b\\\n"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "text/html": [
+ "
WARNING We have set the min_value of K to 1e-99 because there was a postive transform parameter.py:704\n",
+ "
\n"
+ ],
+ "text/plain": [
+ "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m We have set the min_value of K to \u001b[0m\u001b[1;37m1e-99\u001b[0m\u001b[1;38;5;251m because there was a postive transform \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=180693;file:///Users/thomassiegert/.virtualenvs/cosipy/lib/python3.9/site-packages/astromodels/core/parameter.py\u001b\\\u001b[2mparameter.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=200087;file:///Users/thomassiegert/.virtualenvs/cosipy/lib/python3.9/site-packages/astromodels/core/parameter.py#704\u001b\\\u001b[2m704\u001b[0m\u001b]8;;\u001b\\\n"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "text/html": [
+ "
WARNING We have set the min_value of K to 1e-99 because there was a postive transform parameter.py:704\n",
+ "
\n"
+ ],
+ "text/plain": [
+ "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m We have set the min_value of K to \u001b[0m\u001b[1;37m1e-99\u001b[0m\u001b[1;38;5;251m because there was a postive transform \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=842666;file:///Users/thomassiegert/.virtualenvs/cosipy/lib/python3.9/site-packages/astromodels/core/parameter.py\u001b\\\u001b[2mparameter.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=634849;file:///Users/thomassiegert/.virtualenvs/cosipy/lib/python3.9/site-packages/astromodels/core/parameter.py#704\u001b\\\u001b[2m704\u001b[0m\u001b]8;;\u001b\\\n"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "text/html": [
+ "
WARNING We have set the min_value of F to 1e-99 because there was a postive transform parameter.py:704\n",
+ "
\n"
+ ],
+ "text/plain": [
+ "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m We have set the min_value of F to \u001b[0m\u001b[1;37m1e-99\u001b[0m\u001b[1;38;5;251m because there was a postive transform \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=687689;file:///Users/thomassiegert/.virtualenvs/cosipy/lib/python3.9/site-packages/astromodels/core/parameter.py\u001b\\\u001b[2mparameter.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=139954;file:///Users/thomassiegert/.virtualenvs/cosipy/lib/python3.9/site-packages/astromodels/core/parameter.py#704\u001b\\\u001b[2m704\u001b[0m\u001b]8;;\u001b\\\n"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "text/html": [
+ "
WARNING We have set the min_value of K to 1e-99 because there was a postive transform parameter.py:704\n",
+ "
\n"
+ ],
+ "text/plain": [
+ "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m We have set the min_value of K to \u001b[0m\u001b[1;37m1e-99\u001b[0m\u001b[1;38;5;251m because there was a postive transform \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=351491;file:///Users/thomassiegert/.virtualenvs/cosipy/lib/python3.9/site-packages/astromodels/core/parameter.py\u001b\\\u001b[2mparameter.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=963954;file:///Users/thomassiegert/.virtualenvs/cosipy/lib/python3.9/site-packages/astromodels/core/parameter.py#704\u001b\\\u001b[2m704\u001b[0m\u001b]8;;\u001b\\\n"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "/Users/thomassiegert/.virtualenvs/cosipy/lib/python3.9/site-packages/numba/core/decorators.py:262: NumbaDeprecationWarning: \u001b[1mnumba.generated_jit is deprecated. Please see the documentation at: https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-generated-jit for more information and advice on a suitable replacement.\u001b[0m\n",
+ " warnings.warn(msg, NumbaDeprecationWarning)\n"
+ ]
+ },
+ {
+ "data": {
+ "text/html": [
+ "
WARNING The cthreeML package is not installed. You will not be able to use plugins which __init__.py:94\n",
+ "require the C/C++ interface (currently HAWC) \n",
+ "
\n"
+ ],
+ "text/plain": [
+ "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m The cthreeML package is not installed. You will not be able to use plugins which \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=895505;file:///Users/thomassiegert/.virtualenvs/cosipy/lib/python3.9/site-packages/threeML/__init__.py\u001b\\\u001b[2m__init__.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=501476;file:///Users/thomassiegert/.virtualenvs/cosipy/lib/python3.9/site-packages/threeML/__init__.py#94\u001b\\\u001b[2m94\u001b[0m\u001b]8;;\u001b\\\n",
+ "\u001b[38;5;46m \u001b[0m \u001b[1;38;5;251mrequire the C/C++ interface \u001b[0m\u001b[1;38;5;251m(\u001b[0m\u001b[1;38;5;251mcurrently HAWC\u001b[0m\u001b[1;38;5;251m)\u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[2m \u001b[0m\n"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "text/html": [
+ "
WARNING Could not import plugin HAWCLike.py. Do you have the relative instrument __init__.py:144\n",
+ "software installed and configured? \n",
+ "
\n"
+ ],
+ "text/plain": [
+ "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m Could not import plugin HAWCLike.py. Do you have the relative instrument \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=555226;file:///Users/thomassiegert/.virtualenvs/cosipy/lib/python3.9/site-packages/threeML/__init__.py\u001b\\\u001b[2m__init__.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=712719;file:///Users/thomassiegert/.virtualenvs/cosipy/lib/python3.9/site-packages/threeML/__init__.py#144\u001b\\\u001b[2m144\u001b[0m\u001b]8;;\u001b\\\n",
+ "\u001b[38;5;46m \u001b[0m \u001b[1;38;5;251msoftware installed and configured? \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[2m \u001b[0m\n"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "text/html": [
+ "
WARNING Could not import plugin FermiLATLike.py. Do you have the relative instrument __init__.py:144\n",
+ "software installed and configured? \n",
+ "
\n"
+ ],
+ "text/plain": [
+ "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m Could not import plugin FermiLATLike.py. Do you have the relative instrument \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=484813;file:///Users/thomassiegert/.virtualenvs/cosipy/lib/python3.9/site-packages/threeML/__init__.py\u001b\\\u001b[2m__init__.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=984403;file:///Users/thomassiegert/.virtualenvs/cosipy/lib/python3.9/site-packages/threeML/__init__.py#144\u001b\\\u001b[2m144\u001b[0m\u001b]8;;\u001b\\\n",
+ "\u001b[38;5;46m \u001b[0m \u001b[1;38;5;251msoftware installed and configured? \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[2m \u001b[0m\n"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "text/html": [
+ "
16:14:46 WARNING Env. variable OMP_NUM_THREADS is not set. Please set it to 1 for optimal __init__.py:387\n",
+ "performances in 3ML \n",
+ "
\n"
+ ],
+ "text/plain": [
+ "\u001b[38;5;46m16:14:46\u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m Env. variable OMP_NUM_THREADS is not set. Please set it to \u001b[0m\u001b[1;37m1\u001b[0m\u001b[1;38;5;251m for optimal \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=274895;file:///Users/thomassiegert/.virtualenvs/cosipy/lib/python3.9/site-packages/threeML/__init__.py\u001b\\\u001b[2m__init__.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=103135;file:///Users/thomassiegert/.virtualenvs/cosipy/lib/python3.9/site-packages/threeML/__init__.py#387\u001b\\\u001b[2m387\u001b[0m\u001b]8;;\u001b\\\n",
+ "\u001b[38;5;46m \u001b[0m \u001b[1;38;5;251mperformances in 3ML \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[2m \u001b[0m\n"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "text/html": [
+ "
WARNING Env. variable MKL_NUM_THREADS is not set. Please set it to 1 for optimal __init__.py:387\n",
+ "performances in 3ML \n",
+ "
\n"
+ ],
+ "text/plain": [
+ "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m Env. variable MKL_NUM_THREADS is not set. Please set it to \u001b[0m\u001b[1;37m1\u001b[0m\u001b[1;38;5;251m for optimal \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=808420;file:///Users/thomassiegert/.virtualenvs/cosipy/lib/python3.9/site-packages/threeML/__init__.py\u001b\\\u001b[2m__init__.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=183276;file:///Users/thomassiegert/.virtualenvs/cosipy/lib/python3.9/site-packages/threeML/__init__.py#387\u001b\\\u001b[2m387\u001b[0m\u001b]8;;\u001b\\\n",
+ "\u001b[38;5;46m \u001b[0m \u001b[1;38;5;251mperformances in 3ML \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[2m \u001b[0m\n"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "text/html": [
+ "
WARNING Env. variable NUMEXPR_NUM_THREADS is not set. Please set it to 1 for optimal __init__.py:387\n",
+ "performances in 3ML \n",
+ "
\n"
+ ],
+ "text/plain": [
+ "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m Env. variable NUMEXPR_NUM_THREADS is not set. Please set it to \u001b[0m\u001b[1;37m1\u001b[0m\u001b[1;38;5;251m for optimal \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=208938;file:///Users/thomassiegert/.virtualenvs/cosipy/lib/python3.9/site-packages/threeML/__init__.py\u001b\\\u001b[2m__init__.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=247579;file:///Users/thomassiegert/.virtualenvs/cosipy/lib/python3.9/site-packages/threeML/__init__.py#387\u001b\\\u001b[2m387\u001b[0m\u001b]8;;\u001b\\\n",
+ "\u001b[38;5;46m \u001b[0m \u001b[1;38;5;251mperformances in 3ML \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[2m \u001b[0m\n"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "# We'll use histpy's Histogram and mhealpy's HealpixMap as the basis\n",
+ "# develop cosipy. These object (or a derivative) will be passed around by \n",
+ "# the different modules.\n",
+ "from histpy import Histogram\n",
+ "from mhealpy import HealpixMap\n",
+ "\n",
+ "# Needed for coordinate conversion.\n",
+ "# cosipy uses astropy coordinates, with a custom\n",
+ "# SpacecraftFrame (coordinate frame attached to COSI)\n",
+ "from astropy.coordinates import SkyCoord\n",
+ "from cosipy.coordinates import SpacecraftFrame, Attitude\n",
+ "\n",
+ "# cosipy uses astropy units\n",
+ "import astropy.units as u\n",
+ "\n",
+ "#Other standard libraries\n",
+ "import numpy as np\n",
+ "\n",
+ "import matplotlib.pyplot as plt"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "e713e7ed-3956-403e-919d-5f0e1330d1ce",
+ "metadata": {},
+ "source": [
+ "## Signal from a source"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "ebe2b03c-e736-46e2-bcbb-bfb8d3fc3424",
+ "metadata": {},
+ "source": [
+ "Here we obtain the expected number of counts for a given source hypothesis.\n",
+ "\n",
+ "We need the detector response, which can convert from physical values --e.g. source position, spectrum-- to detected counts. For this tutorial we'll use the detector response for an idealized detector, called the \"Compton sphere. The goals of this tutorial is not to get realistics simulations, but to understand the mechanics of of the code.. The response can be found in the sftp server at `/uploads/cosipy/test_data/FlatContinuumIsotropic.LowRes.binnedimaging.imagingresponse.area.nside8.cosipy.h5.zip` (remember to unzip it). "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "id": "a5ff96a4",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "response_path = \"/Users/thomassiegert/COSItools/cosipy/cosipy/test_data/FlatContinuumIsotropic.LowRes.binnedimaging.imagingresponse.area.nside8.cosipy.h5\"\n",
+ "\n",
+ "from cosipy.response import FullDetectorResponse\n",
+ "\n",
+ "response = FullDetectorResponse.open(response_path)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "49f2b36b-b919-4f8b-91ee-fb761590ad71",
+ "metadata": {},
+ "source": [
+ "These is the source hypothesis. Let's inject a GRB-like event, lasting only 1s"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "id": "7c583e3b-6779-4712-94ef-d42b54ecf69f",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "duration = 1*u.s\n",
+ "\n",
+ "# The attitude defines the rotation of the spacecraft with respect\n",
+ "# to the inertial ICRS frame. During flight this information will come from the\n",
+ "# spacecraft telemetry data, but for now let's assume the spacecraft is aligned with the ICRS.\n",
+ "coord = SkyCoord(ra = 20*u.deg, dec = 40*u.deg, \n",
+ " frame = 'icrs', \n",
+ " attitude = Attitude.identity())"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "d27769bb-e830-43ab-a372-52500d71985f",
+ "metadata": {},
+ "source": [
+ "The detector response contains the response as a function of the incoming direction in the spacecraft coordinates. As the spacecraft (SC) rotates, the fixed sky location (RA, Dec) of a given source will move in SC coordinate system, and we need to integrate all locations weighted by the time spent in that direction. These is encoded in the \"dwell time map\", which will be produced by the \"spacecraft orientation\" cosipy module, based on the telemetry information. For our case though, we can assume the GRB occurred at a fixed location because the duration is very short. \n",
+ "\n",
+ "Side note: at some point I was incorrectly calling the \"dwell time map\" a \"exposure map\", which might be confusing."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "id": "28548161",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# The dwell time map has the same pixelation (base) as the detector response.\n",
+ "# We start with an empty map\n",
+ "dwell_time_map = HealpixMap(base = response, \n",
+ " unit = u.s, \n",
+ " coordsys = SpacecraftFrame())\n",
+ "\n",
+ "# Here we add duration of the GRB, at the location where it happend (in spacecraft coordinates)\n",
+ "# We use interpolation to get a better result\n",
+ "pixels, weights = dwell_time_map.get_interp_weights(coord)\n",
+ "\n",
+ "for p,w in zip(pixels, weights):\n",
+ " dwell_time_map[p] += w*duration\n",
+ " \n",
+ "# Without interpolating, it would look like these. Give it a try.\n",
+ "# pixel = dwell_time_map.ang2pix(coord)\n",
+ "# dwell_time_map[pixel] += duration"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "8b32b999-486b-4c05-b608-2df0028a58f9",
+ "metadata": {},
+ "source": [
+ "This is how the dwell time map looks. The red dot shows the exact location of the GRB."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "id": "718ff8f4",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ ""
+ ]
+ },
+ "execution_count": 5,
+ "metadata": {},
+ "output_type": "execute_result"
+ },
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnUAAAGaCAYAAABpDbWeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAABcSAAAXEgFnn9JSAABRJ0lEQVR4nO3deXxU1cH/8c+9d2aSyUIgCRD2XVBBUVlU3IoWKWprtXZRazcrtVWr1i7Wn/ap1lpr+7Q+TzeXWrVuT11q61bQgtYFQUVlkVW2AAkQErJPZube8/tjMpMEAmSZZJLJ9/16zWtm7nLuCUOSb8655xzLGGMQERERkV7NTnUFRERERKTzFOpERERE0oBCnYiIiEgaUKgTERERSQMKdSIiIiJpQKFOREREJA0o1ImIiIikAYU6ERERkTSgUCciIiKSBhTqRERERNKAQp2IiIhIGlCoExEREUkDvlRXQET6JmMM4XCYuro6QqEQ4XCYcDhMQ0ND4nVr2yKRCK7rJh6e57V43n8bgGVZWJaFbduJ9/HXtm0n9juOg9/vx+fz4ff7Ew+fz0cgEEhsj7/PzMwkMzOTYDDY4tnn049WEel++skjIh0SiUSorq6mqqoq8dz8dW1tLXV1dYd8uK6b6i+jS/j9/haBLxgMkp2dTU5OTuKRm5vb4n3zR15eHpmZman+MkSkl7GMMSbVlRCR1ItEIlRUVFBRUUF5eTnl5eWJ1/HtzcNbfX190q4dCAQI10bAs8CzsDzAWC3fe/H3jfsAyzS+NsQeNHsdf59gWr5t/jpemNVYhG1i721aPBsLsBuPsw04BmMbsvODhEKhpIbUjIwM+vfv3+KRl5d3wLb8/HwKCgrIyMhI2rVFpHdSqBNJc67rUlFRwZ49exKP3bt3s2fPHsrKyhKhrbq6ut1lW5aFCQNRGytqxZ4jNkQtLNcG18JyLXCtxm0W7L/dtbAOSFi9j4mHRqcp7OEYsA3GMeDzYtt9sdfxbfH3OIb+w/pRXV1NNBpt9/Vzc3MpLCykoKCAwsLCFq+bP/v9/uR/8SLSIyjUifRy9fX1lJaWUlJSws6dOyktLW0R4MrKytreguQBkcZgFo49H/C6WYAjmh6BrCcxNIZBfyzwGb8Xe93Ks/F7BPr5CIfDbSrbsiwKCwsZPHgwRUVFrT4Hg8Eu/gpFpKso1In0cNFolN27d1NSUpIIbvHXJSUlVFRUHLYM27bx6sFqcLDCNlaDEwtqDU6L0KaQ1vvEQ6AJeBBwMRle0+tA7PXgIwrYu3cvkUjksOXl5eUxePBghgwZwrBhwxg+fDjDhg1j2LBhFBYWJgaYiEjPo1An0gMYY6ioqKC4uLjFY9u2bezcufPw3XFRCyvkxB4NDjTYBwY4hbU+zdDY+pfhYTJcTKYLGW7ifXZRJjU1NYcsIxAIHBD0hg8fzvDhwxk0aBCWpf9jIqmkUCfSjaLRKDt27GDz5s1s2bKlRYCrra096HmBQIDIPhcr5EBjcEuEuJATu39NpJOM0xj4MrxY6AtGMZkuQ48ZRElJySG78YPBIKNGjWLUqFGMHj068Xro0KE4jtONX4VI36VQJ9IFPM+jtLSUTZs2sWXLlsTz1q1bD9oFZlkWpt7Gqnew6x2seh9W4zMNammT1DKYWMtesPGR6WKCUYZPLWLHjh0HDXyBQIDhw4cngt7YsWMZN24cQ4cOVVeuSJIp1Il0Ul1dHRs3bmTDhg1s2LAhEeBCoVDrJ7hg1flij2bBzap3sIyCm/Q+xjKxkJcVjT2Csdf+fPuggziCwSBjx45l/PjxjB8/nnHjxjF27FiysrK6ufYi6UOhTqQd9u3bx/r16xMBbv369ezYsYNWv428ZuGt8WHXOtDgqNVN+oR4656X1SzwZR867A0bNoxx48YxYcIEjjjiCCZNmsSAAQO6ueYivZNCnchBVFVVsWbNGtasWcPatWvZsGEDe/bsaf3gBhu71odV48eqjbfCKbyJtMZgYq152RFMdhQvO0rBEf0oKytr9fjBgwczceJEJk2axKRJk5g4cSK5ubndXGuRnk+hTgQIh8Ns3LiRNWvW8NFHH7FmzRq2b99+wHGWZUGtHQtuNf6mIBfVvUEinWV8HiY7gpcdxeREGT5tEMXFxa22hA8bNiwR8I466igmTpyoVTWkz1Ookz5p165drFixgo8++oiPPvqIjRs3tjqAwap3sKr9WNWNAa7Wp5GmIt3IOF6sNS83gsmJUHRsATt37jzgOJ/PxxFHHMHkyZMTj8LCwhTUWCR1FOok7Xmex5YtW1ixYkXisXv37gMPjFjY8QBX7VcLnEgPZXweJieClxPF5EboPyGb8vLyA44rKipKBLwpU6YwZswYfD5fCmos0j0U6iTthMNh1q1bx4oVK1i5ciUrV648YF1Tx3Hw9rUMcYR0D5xIb5QYkNEvgtcvwthTh7Fp0yY8z2txXFZWFscccwzHHXccU6dOZcKECQp5klYU6qTXi0ajrF+/nuXLl/Pee++xcuXKA0fWuWBVB7Cr/NiVsSBneWqFE0lXxmlszWsMelnD/QdM8K2QJ+lGoU56HWMMmzdvToS4Dz/88MDljSIWdmUAq8qPXRWI3QunOeBE+iyDid2blxfG5EXIGuU/4OdGVlYWU6ZM4fjjj2fatGmMGzdOEyRLr6JQJ73C7t27WbZsGcuXL2f58uUH3j8TjYU4e18AqzKAVaeuVBE5uLaEvPz8fKZNm8b06dOZNm0aBQUFKaqtSNso1EmPFIlEWLFiBUuXLmXp0qVs3ry55QEuWFWxEGdXBrBqfApxItJhLUJe/zAZwxzq6+tbHDN+/HimT5/O9OnTmTJliqZQkR5HoU56jF27dvH222+zdOlS3nvvvRY/UG3bxuxzYiFuXyB2T5y6U0WkixjLYPpF8Po3MP6sEaxbt67F/oyMDKZOncpJJ53EySefTFFRUYpqKtJEoU5SxnVdVq5cyZtvvsnbb7/N1q1bWx4QtrErAtgVGbEgp+lFRCRFjM/D6x/GG9BA/tE5B6x+MXbsWGbNmsXJJ5/MpEmTcBwnRTWVvkyhTrpVKBTinXfe4Y033uCtt96isrIysc+2bUyFkwhyVq26VEWk5zEYTJaLN6CBoz8zllWrVrWYPqV///6JFrzp06eTlZWVwtpKX6JQJ12uoqKCt956izfeeIN33nmn5XQjESvWEre3sTVOqzWISC9jfB7egAa8/AayxrQccOH3+zn++OM5/fTTOeWUU+jfv3/qKippT6FOusSuXbt49dVXef3111m5cmXLtRtDNs7eTOzyDKxKv1rjRCRtJO7Fy2+gaMYAduzYkdhn2zbHHnssZ5xxBqeeeqqWMZOkU6iTpNmzZw+vvvoqixcvZtWqVS32WTW+WGtcubpVRaRvMBhM0MUraGDcnKGsX78+sc+yLCZPnszpp5/O6aefzuDBg1NYU0kXCnXSKWVlZbz22mssXry4RYucZVmwz4dT1tgi16CbhkWkbzMZUdzCBiadO4rVq1e32Ddp0iTOOussZs+erRY86TCFOmm3yspKFi9ezKJFi/jwww9bdK1alX7sskycvRlYYQU5EZHWmICLW9DA5AvGsmLFihZ/EB9//PGcddZZnH766eTk5KS4ptKbKNRJm4TDYZYsWcLChQtZsmQJ0Wg0sc+q8mOXZeCUZSrIiYi0k/G7eIUNTPrMqBa3rgQCAU488UTOOussTjrpJE12LIelUCcHZYxh1apVLFy4kEWLFlFdXZ3YZ9X4sPdkxoKculZFRJLCZLi4A0OMOKOQLVu2JLZnZ2dz+umnM3fuXI499tjYLS4i+1GokwPs3LmTBQsWsHDhwhYjt2iwcfZkYu8OYtf5UldBEZE0l1i2bGCI/OOy2bNnT2LfsGHDmDt3LnPnztUAC2lBoU4AaGho4D//+Q/PP/8877//ftMO14p1re7OxKoMaNSqiEg3M8SmSXEHhcgca1FXVwfE7r+bPn068+bN45RTTiEQCKS4ppJqCnV93ObNm3nuuedYuHAhVVVVQOPI1XI/zu5M7L2ZWJ6CnIhIT2Bsg1cYYvLnx/LBBx8ktvfr149PfvKTzJs3jwkTJqSugpJSCnV9UH19Pa+++irPPfdcy/nkQjbOriDOrqAGPIiI9HAmM4o7KET+CVktumcnTZrE+eefz+zZs8nMzExhDaW7KdT1IZs3b+bvf/87L7/8MrW1tQA4joPZ5cPeFcSuUPeqiEhvYzCY/mFOveZ4Xn/9dSKRCAA5OTnMmzePz3zmM4wYMSLFtZTuoFCX5lzX5c033+Tpp59uea9cvdPYKpeJFVGrnIhIOjA+D3dwPYNOzqWkpCSxfdq0aZx//vmcfPLJ+Hwa6JauFOrSVGVlJS+88ALPPvsspaWlQGzdQXb7cUqCGvQgIpLGDAZvQJgZ3zyKJUuWJCY3HjhwIJ/+9Kf59Kc/zYABA1JcS0k2hbo0s3HjRp5++mlefvllwuFwbGPEwikN4pRmaU45EZE+xmS4uEV15E4JsG/fPiA2sfGcOXO46KKLGDNmTGorKEmjUJcGjDG8/fbbPP744y1GQ1k1PpydWdhlGsEqItLXGSs2cnbC+cNYs2ZNYvvMmTP5/Oc/z7Rp0zSpcS+nUNeLRSIRXnnlFZ544gk2b94MNA58KPXhlGRhVfnVxSoiIi0YDCY3wklXT+H1119PdM2OGTOGL3zhC5x11lma866XUqjrhWpra3nuued48sknm4axRxu7WHdmaToSERFpE5MZxR1aR2BcbLorgPz8fC644AI++9nPkpubm+IaSnso1PUiZWVlPPXUU/zzn/+kpqYmtjFs4+zIwikNYrl2aisoIiK9knE83KJ68mcEE40F2dnZfPazn+Wiiy7SoIpeQqGuF9i1axePPvooL7zwQmL+IavOwdmRjb07E8uoi1VERDovft/diLkFidt6AoEA5513Hl/84he11mwPp1DXg5WUlPDII4/w0ksvEY1GAbCq/Djbs7DLM3S/nIiIdAmDwctvYPwFQ1i7di0Qu2f77LPP5tJLL2X48OEprqG0RqGuB9qxYwd//etfWbBgAa7rAmDt8+MrzsGu1M2rIiLSPQwGkxdm8mVjEhPYx8PdZZddxtChQ1NcQ2lOoa4HKSkp4cEHH2ThwoVNYa4igK84G7tKYU5ERFLHyw1z/BVHsHTpUiAW7ubNm8dll12mbtkeQqGuB9i7dy+PPPII//jHPxLdrHZ5AKc4G7taYU5ERHoOLzfMsd8Yx7vvvguA3+/n3HPP5dJLL2XgwIEprl3fplCXQtXV1Tz++OM89dRThEIhoLFlbmsOdo0/xbUTERE5OK9fmMlfGZ3olg0EAlx44YVceumlmgolRRTqUiAUCvHUU0/x2GOPJaYmsar8sTCne+ZERKQX8fLCHHnJCFasWAFAbm4ul156KRdccAEZGRkprl3folDXjTzPY8GCBdx///2JeYCsWh/O1myNZhURkV7LYPAGhBlx7oDEVCiDBg3iG9/4BnPmzMFxNCl+d1Co6ybLly/n97//PRs2bIhtCNmxlrk9mQpzIiKSFgwGb1CIAbMyEo0XY8eO5dvf/jYzZsxIce3Sn0JdF9u6dSt//OMfeeutt2IbohZOcXZsOS9NGiwiImnI2AZ3SB3BKSRuMzr55JP5zne+w4gRI1Jcu/SlUNdFqqureeCBB3j22WdxXRfHcTDFAXzbcrCiWs5LRETSn/F5fObXZ/DMM8/gui4+n4+LLrqIyy67jOzs7FRXL+0o1CVZ/L65P/3pT1RUVABg7w3gbMnFrveluHYiIiLdzwtGOf7b4xNz3A0YMIBvfvObfOpTn9L9dkmkUJdE69at47e//S2rV68GYuuz+jblYu/T6B8RERF3QAND5uZSXFwMwKRJk/je977HxIkTU1yz9KBQlwRVVVXcd999/POf/8QYA66Fs033zYmIiOzPWLH77TKnGGpra7Ftm/PPP5/LL7+cnJycVFevV1Oo6wRjDK+88gr/+7//y759+wCwd2fi25KDFVZzsoiIyMEYv8tpN03llVdeASA/P5+rr76a2bNnY1lqEOkIhboOKikp4b//+78T9wdYtQ6+j/tpjVaRZnxDitgzZwwDHlrSsQJOPIZT713G68dkJrdiItJjeHkNDP5UDtu3bwdg2rRpXH/99QwfPjzFNet9FOraKRqN8tRTT/HAAw/ElvbyiE1Rsj1bXa0ijXxDitj+hbFgQUaFYcCD7Q915qRjeebJewhaAe7Ye5SCnUgaM5bBHV6LMyFCOBwmIyODyy+/nM997nMaSNEOCnXtsHHjRu68807WrVsHgFXpx7exn0a1ijTjG1LE9i+OxTTO3JO5t32hzpx0LE8++SccLHLsWJD7bcVoXjq6fxfUVkR6EpMZZcr80bz33nsAHHXUUfzoRz9i9OjRqa1YL6FQ1wbRaJRHH32UBx98ENd1IWrh25yDvSuo1SBEmtk/0EH7Ql28dS4e5uIU6kT6DoPBG1xPxlSP2tpa/H4/X/3qV/nSl76Ez6dGlEPRLLiHsWXLFr797W/z5z//Gdd1sfdmEHivAGdXlgKdSDOtBTqAUL5FxVdOOuz5Bwt0ItK3WFg4u7KILA5y0kknEYlEuO+++5g/fz4ff/xxqqvXo6ml7iBc1+XJJ5/k/vvvJxwOx1rnPs7VWq0i+2l+/9z+gS4uo9yQ/5fWW+vMycfy5N9adrfuTy11In2TweANDJE1PTZ9mN/v54orruCiiy7CttUutT/9i7SipKSE7373u/zhD38gHA5jlwcILC/A2aPuVpHmEq1zzsED3aHEA12eHVQLnYgcwMLC2RMk9HKAk08+mUgkwu9//3uuv/56du/enerq9TgKdfv597//zTe+8Q1WrFgBroVvQy6+j/pr3jmR/Rysu7Wtmgc6EZFDsSIO7/5yI74NuWRmZrJ8+XK++tWv8u9//zvVVetRdMdho/r6eu6++25efPFFAKwqP/51eVgNCnMi+4t3ubYn0PU3IeaxmWPYQzA3gwnOczj/mw0X+2GgfhSJyKHF77VzK6Mc+eUxrFmzhp/+9KcsWbKE66+/nqysrFRXMeV0Tx2wfv16br31VrZt24ZlWdhbs3C2ZaurVaQViUDXxr93MiJhfrToWc5f/Q5+zz1gv/EDX+iHua0QMltPibqnTkSaM5bh4vvm8fDDD+N5HiNGjODWW29l3Lhxqa5aSvXpUGeM4ZlnnuEPf/gDkUgEGmz86/OwK7UqhEhrOhLo/vT0fUzffvgRa+akTMyjQyF4YLBTqBOR1ni5Yfqf5WPPnj0EAgG++93vcu655/bZZcb67D119fX13Hbbbdx9991EIpHYVCXvFyjQiRyEb9jQdgU6gB8tfpbp2z/mcH85GsBaEsK6pazV/f2dOpzxY9p+YRHpE+zqAJXPwYknnkg4HOauu+7itttuo66uLtVVS4k+2VK3bds2br75ZjZv3gwGnM05ODs175zI4fhGjaD4cyPadGxBbRWv3HNbq12uB2P8YN4b3eIeu4erCnns0rnYlXW4Gza1t8oi0gcYDO6wOhhfj+u6jBgxgttvv73PrUTR51rqXn31Va644opYoGuw8a8cgG+n7p8TaZNo2wPahSuXtivQAVgR4LGqxPt4oLNCUQj41VonIq2ysPDtyMZe3o+BAwdSXFzM/Pnzef3111NdtW7VZ0Kd67r88Y9/5JZbbqGurg6r0k/gg3zsKnW3irRVdMdORjxV3KZjpxV3bOZ36636xGsPOxboAGNZ0EfvkxGRtol3xx533HHU19dz00038cADD+B5Xqqr1i36RKirq6vjpptu4vHHHwfA2Z6Ff+UArIimKxFptza21mVFGjpWfm3sh++DVYN44stnt9yn1joROQwrarP6dzu48MILAXjwwQe56aabqK2tTXHNul7ah7qSkhKuvPJK3nrrLfDAt7Yfvi256m4V6aC2tNYZB2qCGR27QLadCHRWfaRluZYFGQEFOxE5JMtYPH/dG/jW9yMQCPDmm28yf/58duzYkeqqdam0DnUrVqxg/vz5sfvnwjb+Ffk4ZZq9XqTTDtFaZxwI5xqWjenYfFHvnTCw1UCXKF/dsCLSRs7uIGZZDgMHDmTbtm1861vfYtWqVamuVpdJ21C3cOFCrr32Wvbt24dV44vdP1fjT3W1RNLCwVrr4oEO4MkTZhCx23eLg/HD2guLDhroEtQNKyJtZNf4qXzeMHHiRCorK7n22mtZvHhxqqvVJdIu1BljePTRR/nZz35GNBrFLsvAvyJfa7eKJNt+rXXNAx1AWW4/nj5+emzfYYqK7193wWAevv6cw15a3bAi0h5WxGHzXyqYNWsW4XCYn/zkJzz66KOk26xuaTVPneu6/O53v+Ppp58GYgMinC05un9OpIs0n7du/1AHkBGJ8OeH7mPmlsOPhC2dkcs3vO8Sbcf4CqshrLnrRKTNDIZP//Y0nnrqKQA+/elPc9111+E46dHwkzahrqGhgdtvv51XX30VAGdTDr6d2amtlEia8w0povjisa0GuriMSIT/98KzXLi89bVfw47DxxcO5IZNl7cr0AFYxkBDGHfj5o5UX0T6qOiQOsyEWjzP4xOf+AQ33XQTgUDvn+IsLUJdbW0tN954Ix988EFshOv6PJyyzFRXS6RP8I0dzcavDTnscQU11Vz07lJmbv6Y7HADtYEMlo4Zx7NHnkhwKYe/j+4gFOxEpCPcghDWMXVEIhGmT5/Oz372M4LB3j2YsteHuurqam644QbWrFkDUQv/mv5av1WkG/lGjWDj5cM7dG5mucWQ372L5ffB2JEdKsMKR2JdsL37R5mIpICX14BvRgP19fUcffTR3HnnnfTr1y/V1eqwXj1QYt++fVx77bWxQBex8K8coEAn0s2i27Yz7qHSjp1swETCeHV1sHFLu09P3FOnQCciHWBXZhBdEiQ3N5fVq1dz9dVXs3fv3lRXq8N6bajbu3cv11xzDRs2bIjNQbcyH7tWU5aIdDtjcD/e0u5gl1luMeT37zYrpn3BzIp3uSrQiUgn2DV+Qq8FKCgoYPPmzVx77bW9Ntj1ylC3e/durrnmGrZs2QINdqyFrs6X6mqJ9F3GYEWibT483u1qIuGmIhoaYNO2dl9XRKSz7Dof1f+CQYMGsXXrVq699lrKy8tTXa1263X31JWXl3P11VdTXFwMIZvAygFYDQp0IilnWTjjRvPxV4oOeVjm3lgLXfNA15ydlXXY++vUSiciXcFkRsn7lM2ePXsYNWoUd999N/n5+amuVpv1qpa6yspKrrvuumaBLl+BTqSnaENr3eECHYBxD74EGSjQiUjXsUI+Kl/yGDhwYK9sses1oa6mpoYbbrghto5rg01g1QCshvSYLFAkXRxq0ERbAh0cuhtWgU5EulrzYLdlyxZuuOEGqqurU12tNukVoa6uro4f/OAHrFu3LjYoYtUArJBa6ER6nIMMmsjcazHkj+8dNtDFeXV1BwQ7BToR6S5WyMe+lzzy8/PZuHEjN954Iw0N7ZwdPQV6fKiLRqPccsstrFq1KjYP3er+2PUKdCI9ljFYDU0TCScCXTt/IDbvhlWgE5HuZod8VP8bsrOzWbFiBT/5yU+IRts+ICwVenSoM8bwq1/9imXLloFr4V81QNOWiPQC0e07Eq11lke7Ax200g2rQCci3cyu9RNeEiAQCPDWW2/xy1/+Es/zUl2tg+rRo18feOABHnzwQTDg+6g/TkVGqqskIm1lWdgZGRhjOhTqmpcDKNSJSMq4+SHMlBpc1+WSSy5h/vz5qa5Sq3psS93zzz8fC3SAb2OuAp1Ib2MMXijUuUDXWI4CnYikklOeibUmG4BHH32Ul156KcU1al2PDHXvvPMOv/71rwFwtmXj7MpKcY1ERESkL3N2B7nssssAuOuuu/jggw9SW6FW9Lju1x07dnDFFVdQXV2NvTsT3/p+WFiprpaIiIj0cQbDrJuPZvHixfTr148//elPDB8+PNXVSuhRoa6uro4rr7ySzZs3Y1X58a8cgGUU6ERERKRnMLZh7NfzWbt2LSNHjuRPf/oTOTk5qa4W0IO6X40x/PznP09MLuxfm6dAJyIiIj2K5VlseqiMgQMHsm3bNu644w56SvtYjwl1Dz/8MP/5z3/AIxbowlotQkRERHoeK+Kw72UXv9/P66+/zhNPPJHqKgE9JNQtX76cBx54AADfx/2wqwMprpGIiIjIwdk1fq6++moA7r333h4xcCLl99RVVFTw9a9/nb1792KXZuLfmJfK6oiIiIi0icFw+k+P5eWXXyY/P58///nPFBQUpKw+KW2p8zyPn//85+zduxerzsG3qV8qqyMiIiLSZhYWr/30Q8aMGUN5eTm33357SlecSGmo+9vf/sbSpUvBA9/aPCxPAyNERESk97A8i+1PVZGRkcG7777L008/nbK6pCzUrV+/nnvuuQcA36Zc7Dqt6SoiIiK9j13v49vf/jYA99xzT2wmj1TUIxUXjUQi3HHHHbiui12WgV0aTEU1RERERJLi9xc+ysyZMwmHw/zsZz8jHA53ex1SEur++te/8vHHH0PEwvexVowQERGR3s3C4v3/+Zi8vDw2bNjAX/7yl26vQ7eHug0bNvDXv/4ViE1fYkV6xKwqIiIiIp1iRRxql8ReP/HEE2zYsKFbr9+tiSoajfKLX/yiqdu1LKM7Ly8iIiLSpZzyTM444wxc1+Wuu+7Cdd1uu3a3hrpnnnkmllojFr6Pc9XtKiIiImnnrV+uIjs7m7Vr1/Lss89223W7LdSVlZU1rRqxJRcromXAREREJP1YYYf58+cDsdUmdu/e3S3X7bZQd88991BXV4dV7cPeldldlxURERHpdr+74BGOPvpo6uvr+cMf/tAt1+yWULdixQoWLFgABo12FRERkbRnYbH+kRIsy2LRokWsWrWqy6/Z5aHO8zzuvvvu2MV2BbFrNMmwiIiIpD+71s+8efMA+P3vf48xpmuv16WlA4sXL44Njoha+LbkdPXlRERERHqMl29ZSjAYZPXq1SxatKhLr9WloS4ajXL//fcD4OzIwopqTjoRERHpO6yIw8UXXwzExhd05UoTXZqyXnjhBXbs2AFhG2dHVldeSkRERKRH+us3nqOwsJDS0lJefPHFLrtOl4W6hoYGHnroIQCc4mwsT610IiIi0vdYnsUll1wCwCOPPEIkEumS63RZ0nr++ecpKyuDkI1TGuyqy4iIiIj0eH/8/OMUFBSwe/duXnrppS65RpeEumg0yv/93/8B4NuejWU0hYmIiIj0XZaxEvfWPfLII0Sj0aRfo0tC3WuvvUZpaSlELOzdaqUTERERuedL/0d+fj6lpaW8/PLLSS8/6aHOGMNjjz0GgLMzC8tTK52IiIiI5VlcdNFFADz11FNJn7cu6aFu+fLlsXnpXAunRCNeRUREROL+8rVnyMjIYMOGDXz44YdJLTvpoe7vf/97rOBdmZqXTkRERKQZK2ozZ84cINZal0xJTV1lZWW8+eabADilaqUTERER2d+/fvoWAG+88UZsppAkSWqoe/HFF3FdF6vKj13nS2bRIiIiImnBrvMxZcoUPM/jX//6V/LKTVZBxhheeOEFAM1LJyIiInIIa57ZCsBLL72UtAETSQt1q1evpqSkBFwLuywzWcWKiIiIpB27LINgMEhxcTGrVq1KTplJKQX497//HStwb4amMRERERE5BMu1Of3004GmDNVZSQl10WiURYsWxQrco1Y6ERERkcN55ddvA7FFGzzP63R5SQl1K1eupKKiIraCxL5AMooUERERSWv2vgA5OTns3bs3KV2wSQl1b70VG5prl2donVcRERGRNrCMxaxZswB49dVXO11eUkLdkiVLYoWVZySjOBEREZE+4d+/jXXBLl26tNNldTrUFRcXs23bNvBQ16uIiIhIO9j7AjiOQ3FxMTt37uxcWZ2tzHvvvQeAVeXHcrUsmIiIiEhbWa7N0UcfDcA777zTqbI6ncLef//9WEGVaqUTERERaa8ZM2YAKQ51xhg+/PDDWEEKdSIiIiLt9tA1TwOwatWqTq0u0alQV1xcTHl5OXhgVfs7U5SIiIhIn2TV+PH5fJSXl8dW5+qgToW6tWvXJiqjqUxERERE2s/yLI444giATs1X16lQt379+lghNb7OFCMiIiLSp8UHS6xZs6bDZSQl1Fk16noVERER6ai/374AgE2bNnW4jA6HOmMMGzduBMCqVUudiIiISEfFs9SmTZs6PFiiw6GuoqKCmpoaMGDVKdSJiIiIdJRV78O2bSorK9m7d2+HyuhwqCsuLo69aLA1SEJERESkEyzPYsiQIQBs3769Q2V0ONTFL2jVq5VOREREpLPioa6j05p0ONTF1yezQk5HixARERGRRkVFRQCUlpZ26PwOh7p4f6/VoPVeRURERDrrpd8tBlIQ6srLywGwImqpExEREemseKaKZ6z26nSoI6yWOhEREZFOi8QGnlZVVXXo9A4nsvgFrYhCnYiIiEhnxTNVZWVlh87vcCKrr6+PvXA1nYmIiIhIp0VjsazbW+rioc7yFOpEREREOiueqcLhcIfO71Coc1236YJqqRMRERHpPC/2FIlEOrRUWIdDXULHlicTERERkVYYY5puc2uHDoU629bgCBEREZGu0pEu2M6nM/W+ioiIiHSeaQpVPl/7l2FVS52IiIhIDxMMBtt9TodDXSAQiL1xdFOdiIiISKc5sZESGRkZOE77V+zqcJNbdnY2AKaxAiIiIiLSccYXayjLysrq0PkdDnWJC/rUUiciIiLSaU6KQl1OTg7QlCpFREREpOOML9b7mZub26HzOxzqCgoKYi/87qEPFBEREZHDC8RCXWFhYYdO73CoGzhwIAAmQ/fUiYiIiHSWaQx1iYazdupwqIunSBNQS52IiIhIZ539rdOBFIS6wYMHA2AyFepEREREOmvnzp0ADBkypEPndzjUjRw5EgATVKgTERER6axt27YBMGLEiA6d3+FQl7hghqe56kREREQ6wTge5eXlQFPDWXt1ONTl5uaSn58fq4ha60REREQ6LJ6l8vPzE9PGtVenFnEdO3ZsrCLZkc4UIyIiItKnmZxYlho/fnyHy+hUqJs0aRIAXq5CnYiIiEhHnf3d0wCYOHFih8voVKg78sgjATC50c4UIyIiItKnrVu3DmhqMOuIpLTUmawoxtZyYSIiIiLtZRyPzZs3AylsqRs4cCBFRUVggekX7kxRIiIiIn2S1y+C67oMHTqUQYMGdbicToU6gOOPPz5Wof4KdSIiIiLtddHPPgXA1KlTO1VO8kJdnkKdiIiISHu9//77QFOm6qikhTqTE8X4NAmxiIiISFsZv5sYJJHyUFdYWMiECRPAAi+/obPFiYiIiPQZXn4YYwwTJ06ksLCwU2V1OtQBnHLKKbGKFSjUiYiIiLTViV87BoBZs2Z1uqzkhrr+DZraRERERKQNjG149913gR4U6saPHx+b2sRRF6yIiIhIW3j5IUKhEEOGDOnU8mBxSQl1lmXxyU9+EgBvUH0yihQRERFJa9O/djQAn/zkJ7Esq9PlJSXUAU2hrn9Yo2BFREREDsH4PJYtWwY0ZajOSlqoGz16dGxpCxu8gaFkFSsiIiKSdtyBIVzX5YgjjmDUqFFJKTNpoQ5g7ty5ALhFdRg0YEJERERkfwbD8Nn5AJx77rlJKzepoe7ss88mMzMTk+1i+kWSWbSIiIhIWjD9ImzdupVgMJi0rldIcqjLyclJVM4dUpfMokVERETSwqnXHgfAWWedRXZ2dtLKTWqoAzj//POB2ETEJuAmu3gRERGRXstkuLz22mtAU2ZKlqSHugkTJnDccceBDe4wtdaJiIiIxJ33i9NxXZdp06bFlllNoqSHOoBLLrkEALeoXtObiIiIiBCbxuT5558H4OKLL056+V0S6qZPnx6b3sQxuEPVWiciIiJy8Z/OoaGhgYkTJ3LCCSckvfwuCXWWZTW11g2twzhqrRMREZG+yzgeTz/9NBDr0UzGChL765JQB3DaaacxcuRI8BndWyciIiJ92hf+MJfq6mpGjhzJqaee2iXX6LJQZ9s2l19+ORAbMGH8GgkrIiIifY/xuzz55JMAXHHFFTiO0yXX6bJQB3D66adz5JFHgmOIjqztykuJiIiI9Ejz7jiFUCjE0Ucf3WWtdNDFoc6yLK688koAvKJ6vMxoV15OREREpEcxmVGee+45AObPn98l99LFdWmoA5g6dSonnXQSWOCOqenqy4mIiIj0GNO/OwnXdZk5cyZTp07t0mt1eagD+Na3voXjOHgFDbgDGrrjkiIiIiIp5Q5o4M0338RxHL7zne90+fW6JdSNGTOGiy66CIDo2GqMbbrjsiIiIiIpYWzDwLOCAHz+859n9OjRXX7Nbgl1AF/96lcZOHAgBF3c4Ro0ISIiIunr4nvnUVJSwsCBA/nKV77SLdfstlCXlZXF1VdfDYA7vFaDJkRERCQteZlRHnvsMQCuuuoqsrKyuuW63RbqIDbFyYwZM8CG6IQqDOqGFRERkfRhMBzxlSLC4TDTp0/njDPO6LZrW8aYbk1WJSUlfPWrX6W+vh5ncw6+HdndeXkRERGRLhMdVos7pobs7GwefPBBBg8e3G3X7taWOoAhQ4Y0dcOOqsELqhtWREREej8vGMWZGAZi3a7dGeggBaEO4JxzzmHmzJmxbtgjKtUNKyIiIr2awTD+0kGEw2FmzpzJvHnzur0OKQl1lmXxwx/+kNzcXExuFFdLiImIiEgv9oV7zmbNmjXk5OTwgx/8oEtXjjiYlIQ6gMLCQq6//noA3BG1eHmalFhERER6H69/A4888ggAN9xwQ2wKtxRIWagDOPPMMznvvPPAgsjEKozfTWV1RERERNrF+F1yTrUwxvCZz3yG2bNnp6wuKQ11ANdccw1jx46FgEdkou6vExERkd7BYJh85UgqKioYN24cV111VUrrk/JQl5GRwa233kowGMT0j+j+OhEREekVvnjvXJYvX04wGOS//uu/yMjISGl9Uh7qAEaOHMn3vvc9ANyRtbgFoRTXSEREROTg3IIQDz/8MADf+973GDVqVIpr1ENCHcCcOXO46KKLAIgeUYWXHUlxjUREREQO5GVH8B0XG+D5hS98gTlz5qS4RjE9JtQBXHnllUybNg0cQ+TIfRi/l+oqiYiIiCQYn0fBpwKEQiGmT5/Ot771rVRXKaHblwk7nOrqaubPn8/27duxKv34Vw3AMt0/14uIiIhIc8YyHPntIXz44YcMGzaMe++9l9zc3FRXK6FHtdQB5Obmcscdd5CTk4PJixCdUKURsSIiIpJSBsOpt0zhww8/JCsrizvuuKNHBTrogS11ccuWLeOHP/whruvibM/Ct6Vn/cOJiIhI3xEdXY07vA7Hcbjrrrtit4v1MD2upS5uxowZ/OhHPwLAHV5HdKimOhEREZHuFx1aizu8DoAbb7yxRwY66MGhDuDss89m/vz5ALhja3ALNdWJiIiIdB+3MIQ7tgaAb33rWz1mpGtrenSoA7j44ou58MILAYgeUYk7QGvEioiISNdzBzRgjo4FugsvvJAvfelLKa7RofX4UGdZFldddRWf+MQnwIbokfvw8hTsREREpOt4/Ruwptbiui5nnXUWV111FZbVs2fj6PGhDsBxHG6++WZmzZoFNkSOqsTrF051tURERCQNef3C2CfUE4lEOPXUU/nxj3+M4ziprtZh9djRr60Jh8PcdNNNLF26FFwL/6r+2NWBVFdLRERE0oSXG8Z3Yoj6+npOPPFEbr/9dvx+f6qr1Sa9KtQBNDQ0cOONN/Luu+9C1MK/WsFOREREOs/LDRM4OUxtbS0nnHACv/jFL8jIyEh1tdqs14U6gFAoxA9+8AM++OCDWIvdR/2xKxXsREREpGO8vAac6SFCoRBTp07lzjvvJBgMprpa7dIr7qnbX2ZmJnfeeWfTOrFHV2hUrIiIiHSIO6AB64Q6QqEQM2bM4Je//GWvC3TQS0MdQDAY5I477uCUU06JjYo9ah9ugeaxExERkbZzC0KYY6oJh8Oceuqp/PznPyczMzPV1eqQXtn92lw0GuXnP/85r7zyChjwbeyHs6v3pWsRERHpXu7gOryJtXiex1lnncWPf/xjfD5fqqvVYb0+1AG4rsuvf/1rnn/+eQCcrdk4xdlY9Oz5ZERERKT7GQzuyFrckbElSM855xxuuOGGXjFtyaGkRagDMMZw33338cgjjwBgl2bi+7gfllGwExERkRhjGaLjq/AGx27Z+spXvsLXv/71Hj+xcFukTaiLe/bZZ/ntb3+L53lY5QH86/Kw3F5766CIiIgkibE9jrtuLEuXLsW2ba6//no+/elPp7paSZN2oQ7gjTfe4Kc//SkNDQ1YNT78H/XHCvfuJlURERHpOBNwGX1pPuvXrycjI4P/+q//iq1UlUbSMtQBfPTRR/zoRz9i37590GDjX9Mfu6Z3zAgtIiIiyePlRMidDeXl5eTl5XHnnXdy1FFHpbpaSZe2oQ5gx44d3HjjjWzZsgU88G3oh7NHI2NFRET6CndgPfaUEOFwmDFjxnDHHXcwdOjQVFerS6R1qAOora3ltttu46233gLAKc7C2ZqjkbEiIiJpzGBwR9XgjqgDYNasWdx8881kZWWluGZdJ+1DHcSmPPnzn//cNDK2PIBPAyhERETSknE8Trh+PEuWLAHg0ksv5fLLL8e20/v3fp8IdXGvvPIKv/jFLwiHw1Dv4F+bh12r++xERETShZcdoXBeBiUlJQQCAX7wgx8wZ86cVFerW/SpUAewbt06br75ZkpLS2P32X2ci70rqO5YERGRXsxg8AbXYx/dQDgcZujQodx6660cccQRqa5at+lzoQ6gqqqKn//854n77OxdjRMVewp2IiIivY2xDdFxTRMKn3zyydx0003k5uamuGbdq0+GOgDP83jssce4//77YxMV1/rwrc3Dru+9a76JiIj0NV4wyoiL8ti0aRO2bXP55Zdz8cUXp/39c63ps6Eu7v333+enP/0p5eXl4Fr4Nudgl6o7VkREpCczGLyieuyjYt2tAwYM4Cc/+QnHH398qquWMn0vxu7nuOOO4/7772fatGngGKLjq4keWYnxeamumoiIiLTC+Dxm3jiR6PhqwuEwM2fO5C9/+UufDnSglroEz/N48sknueeee4hGo7FVKNb3w67MSHXVREREpJGXFybvEw5lZWX4/X7mz5/P5z73uT7Z3bo/hbr9rF+/nltvvZVt27YB4GzPwtmWo0EUIiIiKWTs2GTC3vB6jDGMHDmSW265pU+Nbj0chbpWhEIhfve73/HPf/4TAKvOwbehH3Z1IMU1ExER6Xu83DBF5+ZQXFwMwHnnncdVV11FMKilP5tTqDuEN998k7vuuis2iMKAs6NxiTGjVjsREZGuZiyDO7IGMyqE53kUFhbygx/8gBNPPDHVVeuRFOoOo6qqiv/5n/9h4cKFQGOr3fo87BqtRCEiItJVvJwIw87vx5YtWwA4++yzueaaa/rc3HPtoVDXRq+//jq//vWvW7ba6V47ERGRpDK2hzuqFjMi1jqXn5/P9773PU499dRUV63HU6hrh8rKSu6++25eeeWV2IaQje/jfjgVGiErIiLSWe6ABgpmZ7Br1y4AzjrrLL773e+Sl5eX4pr1Dgp1HbBkyRJ+85vfxNaPBew9Gfg25WJFnBTXTEREpPcxfpfo2Gq8gQ0AFBUV8b3vfY+ZM2emuGa9i0JdB9XX1/PAAw/w1FNP4bouRC18W7QahYiISFvFV4XInGqoqanBcRwuuugivva1r2lkawco1HXS+vXr+dWvfsXatWsBsKp9+DblavoTERGRQ/Byw4z+XCEbNmwAYOLEiXz/+9/XvHOdoFCXBK7r8uyzz3L//fdTW1sLgL0rE9+WHHXJioiINGP8LtHRNXiDQwDk5OTwjW98g/PPPx/H0e/MzlCoS6Ly8nLuvfdeXnzxxdiGqIVTnI2zM0tz24mISJ9mLIM7pI6MyR51dXUAnHPOOVxxxRUMGDAgxbVLDwp1XeCjjz7it7/9bVOXbJ2DszkXuyKg++1ERKRPMRi8/DBD5/ZLLME5adIkrr32Wo466qgU1y69KNR1Ec/zeOmll7jnnnvYt28fANY+P77Nudi1mrhYRETSn5cdYfLXR/H+++8DkJeXx/z585k3bx62bae4dulHoa6LVVdX88gjj/D0008TDocBsHdn4tuag9WgewdERCT9mAyX6KgavEGx++YCgQAXXXQRl1xyCTk5OSmuXfpSqOsmpaWl3Hfffbz88suxDR44O7NwirOxXP21IiIivZ9xPNzhtTjjoomGjDlz5nD55ZdTVFSU4tqlP4W6brZ27Vr++Mc/JpqiiVg4OxoHU2jJMRER6YWMbXCH1hGcDDU1NQAce+yxfOc732HSpEkprl3foVCXAsYYlixZwp/+9KfEQsWEbZzt2TglQY2UFRGRXsFYBreonn7TfFRUVAAwZswYLr/8ck455RQsS7/PupNCXQq5rssrr7zCX/7yF3bu3BnbGLLxFWdj71a4ExGRnslg8AaHyJ+Vye7duwEYNmwYX/va1zjzzDM131yKKNT1ANFolBdffJGHHnqIPXv2xDbWO7FwtydT4U5ERHoEYxm8QfUMOj030RgxcOBAvvKVrzBv3jx8Pl+Ka9i3KdT1IA0NDfzjH//gkUceSUyDQsjGtz0be5da7kREJDWMZfAG11MwK4tdu3YBselJvvzlL/OZz3yGjIyMFNdQQKGuR6qrq+Mf//gH//d//0d5eXlsY0PjPXe7ghpQISIi3cLYBreojv7TM9i7dy8ABQUFfOlLX+K8884jGAymuIbSnEJdD9bQ0MBzzz3H448/3tQtG7ZxdmThlAY1FYqIiHQJ43i4Q+rJOdahsrISgEGDBnHJJZcwb948tcz1UAp1vUA4HOall17i0UcfpbS0NLYxauGUBmNToYR1Q6qIiHSeCbi4Q+vInGhRW1sLwNChQ7n00ks5++yz8fu1IlJPplDXi0SjUV5++WUef/zxpqlQPLDLMnF2ZGn5MRER6RAvO4I7rA5rSATXdQEYOXIkX/7ylznzzDM1AKKXUKjrhTzPY+nSpTzxxBNNkxgDVkUA344srH0BLHTfnYiIHJzBYPqHiQ6vw/QPJ7ZPnTqVL37xi5x44olan7WXUajr5datW8cTTzzBq6++mvjryqr14ezMik2HokEVIiLSjLEN3sAQI+YUsHnzZgAcx+GMM87gC1/4glaA6MUU6tJESUkJTz31FM8//zz19fWxjRELZ1cwtkpFg5rORUT6MpPh4g6pI/soH1VVVQAEg0HOO+88Pve5z2lt1jSgUJdmqquref7553n22WcpKSmJbTRgVwRigyrUNSsi0mcYDN6AMN6QOkxBhPiv/KKiIj772c9y7rnnkpubm+JaSrIo1KUp13VZtmwZzzzzDEuXLk1st+od7JJgbL47TYkiIpKWjM/DHVxP0Sl57NixI7F9+vTpfPazn+Wkk07SUl5pSKGuDyguLubZZ5/lpZdeoqamJrbRbRw1uyuIVeVX652ISC9nMJjcCG5RPf6RhoaGBgBycnL41Kc+xfnnn8+IESNSXEvpSgp1fUh9fT0vv/wyzzzzDJs2bUpst+oc7F1BnN2ZWBH95SYi0psYv4s7KMSw0/PZtm1bYvv48eP57Gc/y1lnnaWVH/oIhbo+yBjDRx99xPPPP8+iRYuaBlZ4YJdnYO8KYlfo3jsRkZ4qca9cUT3WoGhi9oPMzEw+8YlPcO655zJ58mQsSz/H+xKFuj6urq6ORYsW8cILL7B69eqmHQ02zq4g9p5M7HqNnBUR6Qm8zCje4HoGHJdFWVlZYvtRRx3FOeecw+zZs8nOzk5hDSWVFOokYdOmTbzwwgssWLAgMdwdwKr2Ye8O4pRlYkU0uEJEpDsZv4dbGMIbVI/JjSa25+XlcfbZZ3POOecwZsyYFNZQegqFOjlAOBzm9ddfZ8GCBbzzzjuJZv341Cj27kzsck1sLCLSVYxt8ApCeANDWAPdxM9hx3GYNm0a55xzDrNmzdJarNKCQp0cUkVFBYsWLWLhwoWsWbOmaYdrYZdl4OzJ1Nx3IiJJEF+2yx0UImOU1XS/MzBp0iTmzJnD7Nmzyc/PT2EtpSdTqJM2Ky4uZuHChSxcuLBpYmOAiIW9NxOnLEMBT0SkHQwGkxfGLWyg35EZ7Nu3L7Fv6NChzJkzh09+8pOaikTaRKFO2s0Yw6pVq1i4cCGvvvoqlZWVTTubB7zKAJZRwBMRac5YTUEud2Kgxc/QvLw8Zs+ezZw5czjqqKM0elXaRaFOOiUajfLBBx+wePFi/vOf/7QS8DJiAywU8ESkDzOWwesfxisMkT3eT3V1dWJfXl4ep512GqeffjrHH388Pp9mHJCOUaiTpIlGo3z44YeJgNe8G4GoFZsDrzwjNgeeligTkTRnHC82l1x+A1ljfU0r+gADBgzgtNNO44wzzuDYY49VkJOkUKiTLhEPeK+++ir/+c9/qKioaNrpgVUVwNkbC3lWg1axEJH0YDJc3PwGvIIG7IKmUasA+fn5nHbaaXziE5/gmGOO0dqrknQKddLlXNdlzZo1vPHGG7z55pts3bq1xX6r1ocdD3g1Pg20EJFeI77eqpcfa5Ez2dEW+0eOHMmsWbOYNWsWRx99tIKcdCmFOul2xcXFvPnmm7z55pusXLkSz/OadoYt7H2xLlp7X4YmOxaRHsf4Pbz+DXgDwuRNDLboiXAchylTpjBr1ixOPvlkjVqVbqVQJylVWVnJ22+/zRtvvMGyZctazMsEjatZVGRg7wtgVfs12EJEup2xGlvjBsSCnMlp2RqXlZXFzJkzmTVrFieeeCL9+vVLUU2lr1Ookx4jEomwevVqli5dyrJly9iwYUPLA6IW9r5ArBWvMgAhR121IpJ0BgOZbmy06oAwGcPtA/7gnDBhAtOnT2fGjBlMmTJFKztIj6BQJz3W3r17eeedd1i2bBnvvPNOy+lSAEI2dmUg8dCACxHpKJPh4uWFEw8yvRb78/LyEiFu+vTpFBQUpKimIgenUCe9guu6rF+/nmXLlvHuu++yevVqotGWXSDUO7GAt88fC3kRhTwRaZ0JNAtx/SOQ6bbY7zgORx99NDNmzGDmzJlMmDAB29Y9vtKzKdRJrxQKhVi5ciXvv/8+y5cvZ926dS2mDgCw6hysqgB2lR+7yq/uWpE+KtGd2i+ClxvB5IUxWQeGuIkTJ3Lcccdx3HHHMWXKFILBYIpqLNIxCnWSFmpra1mxYgXLly9n+fLlbNy4kQP+a4dt7Co/VpUfuyqAVevTwAuRNGQsg8mO4vULY/pF8PpFINCyO9WyLI444ohEiDvmmGPIzs5OUY1FkkOhTtJSVVUVH374IatWrWLlypWsXbv2wO5aF6zqxoBX7ceu9mNF1b0i0tsYnxdrgesXwesXJjDIpqGhocUxPp+PiRMnMnnyZI455himTp1Kbm5uimos0jUU6qRPaGhoYN26daxYsSIR9JqvvZgQsmPhrqYx5NX4sDwFPZGewtgeJicaC3E5Ebzc6AH3wwH069ePyZMnM3nyZKZMmcKkSZPIyMhIQY1Fuo9CnfRJnuexbdu2RMhbs2YN27ZtO7DL1oBV54vNl9cY9qw6dduKdId4N2osvMVCnJVrWk5Y3mjEiBGJADd58mRGjhypgQ3S5yjUiTSqqalh3bp1rFmzhrVr17JmzRr27Nlz4IEeWPW+2D15NT7sWn/stbpuRTrM+LzYfXDZkcYgF8XJ48DbJoCBAwcyadIkjjzySI488kiOOOIIdaWKoFAnckhlZWUtQt769eupqqpq/eCQ3RTw4g+NuBVpITESNTsaC2/ZEbzs6AHzwsXl5uYyadKkRIibNGkShYWF3Vxrkd5BoU6kHYwx7N69m40bN7Jx40bWr1/Pxo0bKSkpaf0Et7FVry7+cLDrfJpeRdJeIrxlRTFZLiYrGnsEo3CQKSSHDBnC+PHjmTBhAuPGjWPChAkMHjwYy9L3ikhbKNSJJEF1dTUff/wxGzZsSAS+rVu3Eg6HWz9h/7BX78QeIR+Wp19g0nsY22Ayo5igiwm2DHAc5I6EQCDAmDFjEsFt/PjxjBs3jpycnO6tvEiaUagT6SKu61JSUsLmzZvZsmVL4nHIsAcQsmPhrnnQq3diXbkaoCEpYCyDyWwMbUG3xWsyWu82hVh4GzVqFKNHj27xGDJkCD6frxu/ApG+QaFOpJs1D3tbt25ly5YtbN++neLi4tanWYkzQENjuAvZWPHXjc+EbXXpSocYDAS8WFjLaAxtGV7sdTCKnUWrI07jcnJyGD58+AEBrqioCMfRcn0i3UWhTqQHqaysTAS85s/bt2+nvr7+0Cd7YDU4sfv1GhyshsbgF7ZjYTBsg2sp+PUxBgOOwQQ8CDQLaxkuJjP2mgz3oF2lccFgkOHDhzNixAiGDx/e4pGXl6f73kR6AIU6kV7AGMPevXvZsWMHpaWllJaWUlJSQklJCaWlpezevfuAtW9b5VpYDXasVS/cGPzCja18YRsiNlZE4a83iIc1/B7G72ECsQcZbuNrN9b6luHFjjsMx3EYNGgQRUVFDBkyhMGDB1NUVMTQoUMZPnw4+fn5Cm4iPZxCnUgaiEaj7N27t0XQ27VrF3v27KGsrIw9e/ZQU1PT9gI9EgGvRdiLP0dtrIgVe442Put+v04xlgGfh/E1e24MbIng5vdiQc0Xez5c61pzOTk5FBYWUlhYyKBBgxgyZAhFRUWJR0FBge5zE+nlFOpE+oj6+nrKysoSIa954Nu7dy/79u2joqKCurq6jl3ApWXIa3zGtbBcC9zG154FUQu8+PamZ7zGh6HHtxSaWCXBbuzebHzGMRin6TV28/f7hbZm7w82zcfhBINB8vLyKCgoYODAgYngVlhY2OJ9MBhM6tcvIj2PQp2ItBAKhRIBL/5cXl6eeF1RUUFNTQ3V1dVUV1dTU1NzyJvoO8yjKeR5xFoCG1/Hgx+NrYNW/KeYadwOLfZjIJER4wcn3jceED/VajzGJhbYGl+bVrYlO3dalkVOTg65ubnk5OTQv39/BgwYQP/+/ROvm7/v37+/wpqIJCjUiUineJ5HbW1ti5AXf11VVUV9fT319fXU1dUd9nVv/XHkOA7BYDDxyMzMbPE+/sjKymoR2po/5+bmkp2drfVKRaTDFOpEpEcwxhCJRAiHw4TD4Ravm7+PRCI0NDQQjUbxPA/P8zDG4LouxpjENs/zWmyzbRvLslo8xx/7b/f5fAQCAfx+P36//5Cv4w8NIhCRVFOoExEREUkDaucXERERSQMKdSIiIiJpQKFOREREJA0o1ImIiIikAYU6ERERkTSgUCciIiKSBhTqRERERNKAQp2IiIhIGlCoExEREUkDCnUiIiIiaUChTkRERCQNKNSJiIiIpAGFOhEREZE0oFAnIiIikgYU6kRERETSgEKdiIiISBpQqBMRERFJAwp1IiIiImlAoU5EREQkDSjUiYiIiKQBhToRERGRNKBQJyIiIpIGFOpERERE0oBCnYiIiEgaUKgTERERSQO+VFdAJM4YQygUSnU1RETaJTMzE8uyUl0NEYU66TlCoRBnn312qqshItIuCxYsIBgMproaIup+FREREUkHaqmTHimwbBCWafybw7KxbAssG2wLLAvLju9r3G5ZYFtY8WMS+6zEOYkHNNtmt9wfOzGxzVhW058+zcpIbLeartV8m7FixST22bFyY9utxL74OaZxW2I/NJVhNx4f30/La7Q4p7H6xm5lX4vjaVHHpm3WAfsOOIfm9dhvPwfZfpDyDlaPA845VLmJ7ebA85udk9jfrCzTuJ1m58X2mWb1ie23mu9LHBvfZxJlWs2Pt0xiX+K/WHx7vLjGY2L/FUziffwcu/F9bF/sffy8xD7LYNF0nt24LfHAJM6zLVpsj53vNZ1H/HgPJ35O4/umsrxEeU6z8h1i2514eYljPZx4mcTr4TUdT1PZsTI9bGLXj+2Llec0brPwcOLnNzvHgdh5xK4T//eIv49dyzS+pnGfwW78d3GwsAGn8cO2sRq3WTiWhY2N1fjJRcIOF36zCJGeRKFOeibXavzxSizU0RjAGn9bNu2zwG5KMFYsITUWEv/tbnPAb+2mxNQyScTLPOC3PPtta34NWtm2/3k0hblmoe6Abc1CWPP3+1ex5fGtnGMfYt/BvowD6nGQL/tQ+w72T9XR8pqV2Vrg69JQ19p+9n9vEmU3r0fza7a2LxECaXZM8+MPOMe0ci3T4tE81DUFxcbHwfYRD36xIpsHwHj4g3g4IxGKmu+LhTqvKRRZzUNR7LVtWbHA1fhM4rWVOC9WDo1lxs+l8bzG7a3ta3aO0xhInUQ946HOHDbUNS/Pif970HKbTfM6NvsMRXoIdb+KiIiIpAGFOhEREZE0oFAnIiIikgYU6kRERETSgEKdiIiISBpQqBMRERFJAwp1IiIiImlA89RJz+QYjIlNOBqbd81q9mztNyFw/LnZa5pvM81et2Ffs0nLmqaIPdj2pmfT4jUtzjMAJr69qUyDBYbEuc33J8poMbla87q08t60qNJ+/x4Heex/bFvmojvUvjZfq437ml/ykOeZw5RpDlLHg08+3HJuuWb7Esd2fPLhpno0m6eOjs9TZ2g6z1im5YPYc2wfLbZ7lgHLayqT+LW8ZvPpNR7TuN9YXqI8WpTf+By/VuN7u/GY+DNwwDav2bd1/LVngUfTPHVe4zaLg81TZyUmDHZo+szi7+3Gc/af+67tkw9bNE0+3Nr3pUhqKdRJjxSesTvVVega8d+ZHbR/JhGJa/5fy0tlRTqseapWJ5JIR+g7R3qMUCiU6iqIiIj0Wmqpkx4jIyMj8fof//gHmZmZKayNJEsoFOIzn/kMoM81nehzbdKXv3bpWRTqpMewrKYOxczMTILBYAprI11Bn2t60ucq0jOo+1VEREQkDSjUiYiIiKQBhToRERGRNKBQJyIiIpIGLGNMJ2bNEhEREZGeQC11IiIiImlAoU5EREQkDSjUiYiIiKQBhToRERGRNKBQJyIiIpIGFOpERERE0oBCnYiIiEgaUKgTERERSQMKdSIiIiJpwJfqCkh6qqur44knnuC1116jtLQU27YZMWIEs2fP5sILL8Tv93e47PLych577DGWLFnCrl27yMjIYMyYMcydO5dzzjkHy7KS+JVIc13xue7Zs4c33niD999/nw0bNrBnzx4A8vPzOfroozn33HM54YQTkv2lSDNd+f26v1/96lf885//BKCoqIi//e1vSStbpK/TMmGSdKWlpVxzzTWUlpYCkJmZied5hMNhACZMmMBvf/tbcnNz2132unXruOGGG6isrAQgGAwSDodxXReAGTNmcMcddyT1l5DEdMXnumvXLj7/+c/T/MdQZmYmxhgaGhoS2+bNm8f3v/99HMdJ0lcjcV35/bq/5cuXc9111yU+b4U6keRS96skVTQa5Uc/+hGlpaUUFBTw3//93yxcuJCFCxfyk5/8hKysLDZs2MBtt93W7rJramr44Q9/SGVlJSNHjuTee+9lwYIFLFy4kGuvvRafz8eyZcv43//93y74yvq2rvpcPc/DGMMJJ5zAj3/8Y5555hkWLlzIggULePjhhznllFMAePHFF/nLX/7SFV9an9aV36/7C4VC3HXXXTiOw6RJk5JQexHZn0KdJNW//vUvNm3aBMBtt93GtGnTALBtmzPPPJMbbrgBgLfffpv33nuvXWU/8cQTlJeXk5GRwS9/+cvELwa/388FF1zA17/+dQCee+45iouLk/UlCV33uebm5nLffffxm9/8hrlz51JYWJgod/To0dx+++3MnDkTgCeffLJF6510Xld+v+7vvvvuY8eOHXzpS19i9OjRnSpLRFqnUCdJ9a9//QuA4447jsmTJx+w/8wzz2TIkCEtjm2rBQsWJMoYOnToAfsvuOACgsEgruvy8ssvt7fqcghd9bnm5OQwceLEg+63LIt58+YBUF9fz9atW9tTbTmMrvx+bW716tU8/fTTjBgxgssuu6zD5YjIoSnUSdKEQiFWrVoFwIknntjqMZZlJVpe3nnnnTaXvW3bNnbt2gWQOH9/WVlZHHPMMe0uWw6tKz/XtggEAonXnucltey+rLs+13A4zC9+8QuMMXz/+98nIyOjYxUWkcNSqJOk2bp1a+KX7pgxYw56XHxfeXk5VVVVbSo73kV0uLLHjh0LwJYtW9pUrhxeV36ubfH+++8DsW72ESNGJK3cvq67PteHHnqIrVu3cs455zB16tQO1VVE2kahTpKmrKws8XrgwIEHPS5+39T+5xzK3r1721V2bW0tdXV1bSpbDq0rP9fD2blzZ2L6i9mzZ5OdnZ2UcqV7Ptf169fz2GOPkZ+fz5VXXtn+SopIuyjUSdI0D1GH6mLJzMxs9ZxUlS2Hlqp/+4aGBn7yk58QCoXIy8tj/vz5nS5TmnT15xqNRrnzzjtxXZdrrrkmKVOiiMihKdSJSI8TjUa59dZbWbduHT6fj5tvvrlFi5H0fI8++igbNmzg5JNPZvbs2amujkifoFAnSZOVlZV4faipJ0KhUKvnpKpsObTu/rd3XZfbbruN119/HcdxuOWWW5gxY0aHy5PWdeXnumXLFh5++GGCwSDXXXddxyspIu2iUCdJ07wlJb7UU2ua35fT1taXgoKCdpWdnZ2tUJckXfm57s91XX72s5+xePFiHMfh//2//8cZZ5zRobLk0Lryc/3Nb35DJBLhy1/+Mrm5udTV1bV4xFeAMcYktkWj0Q5+JSISp7VfJWlGjRqFbdt4nsfmzZsPOk3C5s2bgdjanv369WtT2fFRrfHzDzZ5aXyUrCY3TZ6u/Fybi7fQLVq0KBHozjzzzE7VXQ6uKz/XkpISAO69917uvffegx63a9cu5s6dC8BVV13F5z//+fZ8CSKyH7XUSdJkZmYmJjBdunRpq8cYY1i2bBkA06dPb3PZI0aMYPDgwYcsu76+nhUrVrS7bDm0rvxc41zX5dZbb1Wg60bd8bmKSPdSqJOkiv/V/f777/PRRx8dsH/x4sXs3LmzxbFtYVkWZ599NgCLFi1KtAQ09/e//536+nocx+GTn/xkR6ovB9FVnys0tdA173JVoOseXfW5/u1vf+M///nPQR/xsoqKihLb1Eon0nkKdZJUc+fOZezYsRhjuPnmmxPrRXqex+LFi7nrrruA2KoQJ5xwQotzH3jgAU477TROO+20VkPbF7/4RfLz8wmFQvzwhz9k3bp1AEQiEZ599ln+/Oc/A3Deeedpktok66rPNX4PXbyF7pZbblGg60Zd+f0qIt1P99RJUvl8Pu644w6++93vUlpaynXXXUdmZiae5xEOhwGYMGECt9xyS7vLzsnJ4c477+SGG25gy5YtfPOb3yQrK4twOJy4yXr69OlcddVVSf2apOs+15UrV/Lvf/8biLXG3n333dx9990HPf7qq69W6Euirvx+FZHup1AnSTdkyBAefPBBnnjiCV577TVKS0vx+XyMGTOGM888kwsvvBC/39+hsidOnMhDDz3EY489xltvvcXu3bvJzMxk7NixzJ07l3nz5mHbaoDuCl3xuRpjEq+j0Sjl5eWHPD4eNCR5uvL7VUS6l2Wa/1QVERERkV5JTRoiIiIiaUChTkRERCQNKNSJiIiIpAGFOhEREZE0oFAnIiIikgYU6kRERETSgEKdiIiISBpQqBMRERFJAwp1IiIiImlAoU5EREQkDSjUiYiIiKQBhToRERGRNKBQJyIiIpIGFOpERERE0oBCnYiIiEgaUKgTERERSQMKdSIiIiJpQKFOREREJA38f4aMcFpQvzsUAAAAAElFTkSuQmCC",
+ "text/plain": [
+ "
"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "_,ax = dwell_time_map.plot(coord = SpacecraftFrame(attitude = Attitude.identity()));\n",
+ "\n",
+ "ax.scatter(coord.ra, coord.dec, transform=ax.get_transform('world'), color = 'red')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "d532c332-d566-4787-8830-d42271c22442",
+ "metadata": {},
+ "source": [
+ "The sum of all pixel in the dwell time map is simply the duration of the data that was integrated. In this case is just the duration of the GRB."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "id": "6bbffb4e-853a-4925-8987-7cfd57d33029",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/latex": [
+ "$1 \\; \\mathrm{s}$"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "execution_count": 6,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "sum(dwell_time_map)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "9c6baed9-14e7-48fc-84f5-af349de2b715",
+ "metadata": {},
+ "source": [
+ "The detector response is then convolved with the dwell time map to get the point source response:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "id": "b2d58300-deae-4c26-9469-e1ac07fa8313",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "psr = response.get_point_source_response(dwell_time_map)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "f817be35-5398-48b9-b088-d308195a302f",
+ "metadata": {},
+ "source": [
+ "The point source response is still quite generic. We obtained the response for a give location and duration, but we still convolved this with a given spectral assumption:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "id": "3d93229f-7c8e-4c9a-a582-54cc73bf9eec",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from threeML import *\n",
+ "from astropy.units import Quantity\n",
+ "\n",
+ "# This is not a very realistic spectrum, but help to illustrate the method.\n",
+ "# You can play with other parameters. \n",
+ "spectrum = Powerlaw(K = 1, index = -1.7)\n",
+ "\"\"\"index = -1.7\n",
+ "piv = 1 * u.keV\n",
+ "K = 1e0 / u.cm / u.cm / u.s / u.keV\n",
+ "spectrum.index.value = index\n",
+ "spectrum.K.value = K.value\n",
+ "spectrum.piv.value = piv.value \n",
+ "spectrum.K.unit = K.unit\n",
+ "spectrum.piv.unit = piv.unit\n",
+ "\"\"\" \n",
+ "# We project into the only event parameters that we can measure in COSI\n",
+ "signal = psr.get_expectation(spectrum).project(['Em', 'Phi', 'PsiChi'])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "7f6e60b5-65e4-4171-a8a7-e8e4f4d72008",
+ "metadata": {},
+ "source": [
+ "The result `signal` is histogram that contains the expected counts in measured energy (`Em`) and the \"Compton Data Space\": Compton scatter angle (`Phi`) and direction (in SC coordinates) of the scattered photon in he (`PsiChi`). For reference, see the following figure from [this](https://arxiv.org/abs/2102.13158) paper. The only different is that here we are using spacecraft coordinate instead of galactic coordinates. "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "9ea7ffcb-d555-4da0-913a-87c7087b60c0",
+ "metadata": {},
+ "source": [
+ "The `signal` object is a 3D histogram. Note that the last axis, `PsiChi`, is actually a 2D axis encoding the coordinates in a sphere as pixels in a HEALPix grid. So, in a sense, it's really a 4D histogram."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "id": "8d821639-d2c7-42a8-9f3b-d6ef53268356",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "array(['Em', 'Phi', 'PsiChi'], dtype=',\n",
+ " )"
+ ]
+ },
+ "execution_count": 10,
+ "metadata": {},
+ "output_type": "execute_result"
+ },
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkUAAAG7CAYAAADNDuE1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA2ZElEQVR4nO3de1RVdf7/8dfhqlwU56BoQnjBEDWtGc28jGLgFaWc1GyGJmu85K3SmVFrqTORM82UdjGzzJwfrTJBbanMjKJmhqKpaZoKZpgjY17h4AUEDhfP7w8X5xsDGB44F+D5WMuV5/N5733em5bHl3t/9j4Gi8ViEQAAQCPn5uwGAAAAXAGhCAAAQIQiAAAASYQiAAAASYQiAAAASYQiAAAASYQiAAAASYSiO1JUVKSTJ0+qqKjI2a0AAIA6Rii6A1lZWZo0aZKysrKc3QoAAKhjhCIAAAARigAAACQRigAAACQRigAAACQRigAAACQRigAAACQRigAAACQRigAAACQRigAAACQRigAAACQRigAAACQRigAAACQRigAAACQRigAAACQRigAAACRJHs5uAKit0/su6ND6TJUUlTq7FfwEzyYe6jm2k9r3buPsVgCgEkIR6r1D6zN17fwNZ7eBGjHr4LpMQhEAl0QoQr1XfobIYJCatvB2cjeoTuEVsywWcUYPgMsiFKHBaNrCW79e9pCz20A1PpnxuQpyzc5uAwCqxUJrAAAAEYoAAAAkueDls//85z/6f//v/+nkyZPKzc1VkyZNFBoaqscff1z9+vWrUHvmzBktW7ZMx44dk4eHh/r06aMZM2YoICCgQt3NmzeVmJiojRs3Kjc3V8HBwYqLi1N0dLQDjwwAALgylwtFFy9eVEFBgYYNG6bAwEAVFRUpNTVVL7zwgv7whz8oNjZWknT58mXNnDlTfn5+mjRpkgoLC5WYmKjTp09rxYoV8vT0tO5z5cqVWr16tUaNGqXOnTsrLS1N8fHxMhgMioqKctahAgAAF+JyoahPnz7q06dPhbFf/epXmjRpktauXWsNRR9//LGKior0wQcfKCgoSJIUERGh2bNna8uWLda67OxsJSUlafTo0Zo1a5YkaeTIkZo5c6aWL1+uyMhIubu7O/AIAQCAK6oXa4rc3d3VqlUr5efnW8dSU1PVt29fayCSpJ49eyokJEQ7d+60jqWlpam0tFSjR4+2jhkMBj3yyCPKzs5Wenq6Yw4CAAC4NJc7U1SusLBQZrNZN27c0J49e7R//34NGjRI0q2zP1euXFF4eHil7SIiIrRv3z7r68zMTDVt2lShoaGV6srnu3fvXmUPOTk5MplM1tdZWVm1Pi4AAOCaXDYUvfPOO0pOTpYkubm5acCAAdbLX+VBxWg0VtrOaDTq+vXrKi4ulpeXl0wmk1q0aCGDwVCpTroVfKqTnJyshISEujgcAADg4lw2FI0dO1aRkZHKycnRzp07VVZWppKSEkmS2XzrAXA/XkxdzsvLy1rj5eUls9n8k3XViY2NrXDHW1ZWlhYtWmT7QQEAAJflsqEoNDTUeslr2LBhmj17tubNm6cVK1bI2/vWVzmUh6QfKy4uliRrjbe3d43qqhIYGKjAwMDaHQgAAKgX6sVCa0mKjIzUt99+q7Nnz1ovff14vU85k8mkZs2aWc8EGY1G5ebmymKxVKqTROgBAACS6lEoKr/MlZ+fr5YtWyogIEAnT56sVHfixAmFhYVZX4eFhamoqKjSIumMjAzrPAAAgMuFoitXrlQaKy0t1datW+Xt7a127dpJkgYOHKi9e/fq0qVL1rpDhw7p7Nmz1rvUJKl///7y8PDQhg0brGMWi0WbNm1Sy5Yt1a1bN/sdDAAAqDdcbk3R4sWLdePGDfXo0UMtW7aUyWTS9u3b9d///lfTp0+Xj4+PJCkuLk5ffPGFnn/+eY0ZM0aFhYVas2aNOnTooOHDh1v316pVK40dO1Zr1qxRaWmpIiIitHv3bh09elQLFizgwY0AAECSC4aihx56SP/+97+1adMmXbt2TT4+PgoPD9czzzyj/v37W+uCgoK0dOlSLVu2TCtWrLB+99n06dOt64nKTZkyRf7+/kpOTlZKSoqCg4M1f/58DR482NGHBwAAXJTLhaKoqKgafx9Z+/bttWTJkp+sc3NzU1xcnOLi4mrbHgAAaKBcbk0RAACAMxCKAAAARCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQJHk4u4H/deLECaWkpOjw4cO6ePGimjVrpq5du2rixIkKCQmx1v31r39VSkpKpe3vvvtuffzxxxXGbt68qcTERG3cuFG5ubkKDg5WXFycoqOj7X48AACgfnC5UPTJJ5/o2LFjGjRokDp27CiTyaQNGzZo4sSJevfdd9WhQwdrrZeXl+bMmVNhe19f30r7XLlypVavXq1Ro0apc+fOSktLU3x8vAwGg6Kioux+TAAAwPW5XCgaN26cFi5cKE9PT+vYQw89pKeeekqrV6/WggULrOPu7u4aMmTIbfeXnZ2tpKQkjR49WrNmzZIkjRw5UjNnztTy5csVGRkpd3d3+xwMAACoN1xuTdG9995bIRBJUkhIiNq1a6esrKxK9WVlZbpx40a1+0tLS1NpaalGjx5tHTMYDHrkkUeUnZ2t9PT0umseAADUWy53pqgqFotFV65cUbt27SqMFxUVafjw4SoqKpK/v7+ioqL0zDPPyMfHx1qTmZmppk2bKjQ0tMK2ERER1vnu3btX+b45OTkymUzW11WFMgAA0DDUi1C0fft2ZWdn6+mnn7aOGY1GPf7447rnnntksVi0f/9+bdy4Ud9//73eeusteXjcOjSTyaQWLVrIYDBU2KfRaJR0K/hUJzk5WQkJCXV/QAAAwOW4fCjKysrSG2+8oa5du2rYsGHW8SlTplSoi4qKUkhIiFauXKnU1FTrAmqz2Vzpcpx0a5F2+Xx1YmNj1a9fvwq9LFq0qFbHAwAAXJPLrSn6MZPJpLlz58rX11cvv/zyTy6IHjdunNzc3HTw4EHrmLe3t0pKSirVFhcXW+erExgYqPDwcOuv/70EBwAAGg6XDUX5+fmaM2eO8vPztXjxYgUGBv7kNt7e3mrWrJmuX79uHTMajcrNzZXFYqlQW75WqCb7BQAADZ9LhiKz2ax58+bp7Nmz+tvf/lZpgXV1CgoKdO3aNQUEBFjHwsLCVFRUVGmRdEZGhnUeAADA5UJRWVmZ/vznPys9PV0vvfSSunXrVqnGbDaroKCg0viHH34oi8Wi3r17W8f69+8vDw8PbdiwwTpmsVi0adMmtWzZssr9AwCAxsflFlq/88472rNnj/r27au8vDxt27atwvyQIUOUm5ur3/3ud4qOjtbdd98tSTpw4ID27dun3r17q3///tb6Vq1aaezYsVqzZo1KS0sVERGh3bt36+jRo1qwYAEPbgQAAJJcMBSdOnVKkrR3717t3bu30vyQIUPk5+envn376quvvlJKSopu3ryptm3bavLkyRo/frzc3CqeAJsyZYr8/f2VnJyslJQUBQcHa/78+Ro8eLBDjgkAALg+lwtFS5cu/ckaf39/zZ8/v8b7dHNzU1xcnOLi4mrTGgAAaMBcbk0RAACAMxCKAAAARCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQRCgCAACQJHk4u4H/deLECaWkpOjw4cO6ePGimjVrpq5du2rixIkKCQmpUHvmzBktW7ZMx44dk4eHh/r06aMZM2YoICCgQt3NmzeVmJiojRs3Kjc3V8HBwYqLi1N0dLQDjwwAALiyGoWixx57rNZvNHbsWI0ZM+Yn6z755BMdO3ZMgwYNUseOHWUymbRhwwZNnDhR7777rjp06CBJunz5smbOnCk/Pz9NmjRJhYWFSkxM1OnTp7VixQp5enpa97ly5UqtXr1ao0aNUufOnZWWlqb4+HgZDAZFRUXV+tgAAED9V6NQdPHiRfn6+srPz8+mN7l8+bLy8/NrVDtu3DgtXLiwQqh56KGH9NRTT2n16tVasGCBJOnjjz9WUVGRPvjgAwUFBUmSIiIiNHv2bG3ZskWxsbGSpOzsbCUlJWn06NGaNWuWJGnkyJGaOXOmli9frsjISLm7u9t0XAAAoOGo8eWzcePGacKECTa9ycCBA2tce++991YaCwkJUbt27ZSVlWUdS01NVd++fa2BSJJ69uypkJAQ7dy50xqK0tLSVFpaqtGjR1vrDAaDHnnkEcXHxys9PV3du3e35bAAAEAD4nJriqpisVh05coVtWvXTtKtsz9XrlxReHh4pdqIiAjt27fP+jozM1NNmzZVaGhopbry+epCUU5Ojkwmk/X1j0MZAABoWGoUij766CM1b97c5jep7fbbt29Xdna2nn76aUmyBhWj0Vip1mg06vr16youLpaXl5dMJpNatGghg8FQqU66FXyqk5ycrISEBJv7BgAA9UeNQtHdd99dqzepzfZZWVl644031LVrVw0bNkySZDabJanCuqNyXl5e1hovLy+ZzeafrKtObGys+vXrV6GXRYsW2XwsAADAddX55bPi4mK5ubnJw6P2uzaZTJo7d658fX318ssvWxdEe3t7S5JKSkqqfP8f13h7e9eoriqBgYEKDAys3UEAAIB6waaHNx45ckSrVq1SXl6edezatWv64x//qKFDh2r48OF67733atVYfn6+5syZo/z8fC1evLhCOCm/9PXj9T7lTCaTmjVrZj0TZDQalZubK4vFUqlOEqEHAABIsjEUJSYm6rPPPpO/v7917J133tGBAwfUpk0b+fn5KTExUZ9//rlNTZnNZs2bN09nz57V3/72N+sC63ItW7ZUQECATp48WWnbEydOKCwszPo6LCxMRUVFlRZJZ2RkWOcBAABsCkWZmZkVbp03m83auXOnevXqpU8++USrV69Wq1attGnTpjved1lZmf785z8rPT1dL730krp161Zl3cCBA7V3715dunTJOnbo0CGdPXtWgwYNso71799fHh4e2rBhg3XMYrFo06ZNatmyZbX7BwAAjYtNC3+uX7+uli1bWl+np6eruLhYw4cPlyT5+Piob9++Sk1NveN9v/POO9qzZ4/69u2rvLw8bdu2rcL8kCFDJElxcXH64osv9Pzzz2vMmDEqLCzUmjVr1KFDB2sfktSqVSuNHTtWa9asUWlpqSIiIrR7924dPXpUCxYs4MGNAABAko2hyNvbWwUFBdbXhw8flsFg0H333Wcda9q0aYU1RzV16tQpSdLevXu1d+/eSvPloSgoKEhLly7VsmXLtGLFCut3n02fPt26nqjclClT5O/vr+TkZKWkpCg4OFjz58/X4MGD77g/AADQMNkUitq2bav9+/eruLhYBoNBO3bsUGhoaIXnBl26dEktWrS4430vXbq0xrXt27fXkiVLfrLOzc1NcXFxiouLu+N+AABA42DTmqJRo0bp3Llzevzxx/XEE0/o/PnzGjFiRIWa7777rtICaQAAAFdlUyiKiYnR+PHjVVxcrBs3bujhhx/W2LFjrfPHjx/X2bNn9fOf/7zOGgUAALAnmy6fGQwGTZ06VVOnTq1yPjw8XP/+97/VpEmTWjUHAADgKHb5QlhPT88qv1oDAADAVdl0+QwAAKChIRQBAACIUAQAACCJUAQAACCJUAQAACDpDkJRenq6PfsAAABwqhrfkj9t2jS1a9dOMTExGjJkiAICAuzYFgAAgGPV+ExR165ddebMGS1fvlxjxozRwoULtX//flksFnv2BwAA4BA1PlO0fPlynT17Vv/617+0fft2paamateuXTIajRo+fLhGjBihu+66y569AgAA2M0dLbQOCQnR1KlTtX79er3yyivq16+frl69qo8++ki//vWv9fzzz2v79u0qLi62V78AAAB2YdPXfLi5ualv377q27evrl69qq1bt2rz5s06fPiwjhw5ojfffFPR0dEaMWKEwsPD67pnAACAOlfrW/IDAgL02GOP6cMPP9R7772n2NhYSdKmTZv0zDPP1LpBAAAAR6jT5xQFBwcrNDRUgYGBslgsLMIGAAD1hk2Xz/7XV199pc2bNystLU0lJSWyWCzq0aOHYmJi6mL3AAAAdmdzKLpw4YK2bNmilJQUXb58WRaLRUajUcOGDdOIESMUHBxcl30CAADY1R2FouLiYn3xxRfavHmzvvnmG928eVPu7u7q16+fYmJi9OCDD8rNjW8OAQAA9U+NQ9GSJUu0Y8cOFRQUyGKxKCQkRCNGjNCwYcP0s5/9zJ49AgAA2F2NQ1FycrKaNGmioUOHKiYmRt27d7dnXwAAAA5V41D0hz/8QVFRUfLx8bFnPwAAAE5R41A0atSoKscLCgp09uxZFRUVqUePHnXWGAAAgCPV6u6zpUuXat++fbJYLDIYDNq5c6ck6dixY3r11Vc1e/Zs3X///XXWLAAAgL3YdKvYpUuXNHXqVO3bt0/9+/dX165dKzyoMSIiQteuXdNnn31WZ40CAADYk02h6B//+Ify8vK0dOlSvfzyy+rZs2eFeQ8PD3Xv3l3Hjx+vkyYBAADszaZQdODAAf3yl7/UvffeW21N69atlZ2dbXNjAAAAjmRTKMrLy1Pr1q1vW2OxWFRSUmJTUwAAAI5mUyhq0aKFfvjhh9vWnD59WkFBQTY1BQAA4Gg2haKePXvqyy+/1Pfff1/l/DfffKOvv/5aDz74YK2aAwAAcBSbbsn/7W9/q9TUVM2cOVPjx4/XuXPnJEn79u3T8ePHtXbtWjVv3lzjx4+v02YBAADsxaZQ1KZNGy1evFh//vOftWrVKhkMBlksFs2bN08Wi0VBQUGKj49XYGBgXfcLAABgFzY/vLFLly765JNPtHfvXmVkZCgvL08+Pj7q0qWL+vfvL09Pz7rsEwAAwK5sCkXHjh3TvffeKw8PDw0YMEADBgyosm7t2rUaN25crRoEAABwBJsWWr/wwgs6c+bMbWvWrl2r5cuX27J7AAAAh7MpFHl7e+uPf/xjtQ9nXL9+vd555x117dq1Vs0BAAA4ik2h6LXXXtONGzf0+9//XtevX68w9+mnn+rtt99Wly5d9Nprr9VJkwAAAPZmUyjq0KGD/va3v+nixYuaO3euzGazpFuBaOnSpYqIiNCSJUvk4+NTp80CAADYi02hSJK6d++uhQsX6ttvv9WCBQu0bt06LV26VOHh4Vq8eDGBCAAA1Cs2hyJJ6t+/v37/+99r//79eueddxQeHq433nhDfn5+ddUfAACAQ9TolvwjR45UOxccHKx+/frp2LFjeuKJJ5SZmVlh/r777rujhgoKCpSYmKiMjAydOHFCeXl5euGFFzR8+PAKdX/961+VkpJSafu7775bH3/8cYWxmzdvKjExURs3blRubq6Cg4MVFxen6OjoO+oNQO0VXjHrkxmfO7uNesmziYd6ju2k9r3bOLsVoEGqUSh67rnnZDAYbltjsVi0YMGCSuNffPHFHTV07do1JSQkKCgoSGFhYTp8+HC1tV5eXpozZ06FMV9f30p1K1eu1OrVqzVq1Ch17txZaWlpio+Pl8FgUFRU1B31B8A2nk08JJllsUgFuWZnt1NPmXVwXSahCLCTGoWiJ5988idDUV0xGo3asGGDjEajvv32W02ePLnaWnd3dw0ZMuS2+8vOzlZSUpJGjx6tWbNmSZJGjhypmTNnavny5YqMjJS7u3udHgOAynqO7aSD6zJVUlTq7FbqpcIrtwIlPz/AfmoUip5++ml792Hl5eUlo9FY4/qysjIVFRVVeYZIktLS0lRaWqrRo0dbxwwGgx555BHFx8crPT1d3bt3r3XfAG6vfe82nOGohU9mfM4ZNsDObP7uM1dQVFSk4cOHq6ioSP7+/oqKitIzzzxT4c63zMxMNW3aVKGhoRW2jYiIsM5XF4pycnJkMpmsr7OysuxwFAAAwBXU21BkNBr1+OOP65577pHFYtH+/fu1ceNGff/993rrrbfk4XHr0Ewmk1q0aFHp8l/52aicnJxq3yM5OVkJCQl2OwYAAOA66m0omjJlSoXXUVFRCgkJ0cqVK5WammpdQG02m+Xp6Vlpey8vL+t8dWJjY9WvXz/r66ysLC1atKgu2gcAAC6mVs8pcjXjxo2Tm5ubDh48aB3z9vZWSUlJpdri4mLrfHUCAwMVHh5u/fW/l+AAAEDD0aBCkbe3t5o1a1bh+9iMRqNyc3NlsVgq1JavFQoMDHRojwAAwDU1qFBUUFCga9euKSAgwDoWFhamoqKiSoukMzIyrPMAAAD1MhSZzWYVFBRUGv/www9lsVjUu3dv61j//v3l4eGhDRs2WMcsFos2bdqkli1bqlu3bg7pGQAAuDaXXGj96aefKj8/33qJa8+ePbp8+bIk6dFHH1VeXp5+97vfKTo6Wnfffbck6cCBA9q3b5969+6t/v37W/fVqlUrjR07VmvWrFFpaakiIiK0e/duHT16VAsWLODBjQAAQJINochkMun48eNyd3dXjx495O/vX2XdkSNHdOTIEU2YMOGOm0pKStLFixetr3ft2qVdu3ZJkoYMGSI/Pz/17dtXX331lVJSUnTz5k21bdtWkydP1vjx4+XmVvEE2JQpU+Tv76/k5GSlpKQoODhY8+fP1+DBg++4NwAA0DDdUShKTEzUBx98oNLSW4+Z9/Ly0hNPPKG4uLhKzwE6fPiwPvzwQ5tC0dq1a3+yZv78+TXen5ubm+Li4hQXF3fHvQAAgMahxmuKDhw4oHfffVdeXl4aOXKkHnnkEfn4+GjVqlWaN2+e9RZ3AACA+qjGZ4rWrVunJk2aaMWKFQoJCZEkTZ48WYsXL9aOHTs0b948vfLKK7d97g8AAICrqvGZom+//VYDBgywBiJJ8vHx0cKFC/XrX/9ahw4d0rx58277hGgAAABXVeNQVFhYqFatWlU5N2XKFD3xxBP6+uuvNXfuXIIRAACod2p8+SwwMFDZ2dnVzk+cOFGS9NFHH2nOnDkKDw+vfXcAAAAOUuNQ1L59ex06dOi2NT8ORsePH69dZwAAAA5U48tnffr0UU5Ojr788svb1k2cOFG//e1vrbftAwAA1Ac1PlMUGRkpi8WiJk2a/GTt7373O911110VHsAIAADgymocipo1a6aHH364xjsePny4TQ0BAAA4Q738QlgAAIC6VqsvhM3MzNSpU6dkMpmqXENkMBj05JNP1uYtAAAAHMKmUHTlyhXFx8fr8OHDkiSLxVJlHaEIAADUFzaFojfeeENff/21HnzwQUVFRcloNMrd3b2uewMAAHAYm0LRgQMHdP/99+vvf/97XfcDAADgFDaFIg8PD55YXYdO77ugQ+szVVLEs51sUXiFr5UBANSeTaGoe/fuyszMrOteGq1D6zN17fwNZ7dR73k2qdV9AwCARs6mv0UmT56s6dOn69NPP9Wjjz5a1z01OuVniAwGqWkLbyd3Uz95NvFQz7GdnN0GAKAesykUtWvXTsuWLdOMGTP06aefqmPHjvL19a2ydt68ebVqsDFp2sJbv172kLPbAACgUbIpFJ0/f14vvvii8vPzlZ+fr3PnzlVZZzAYCEUAAKBesCkUvfXWWzp//rwefvhhRUdHc0s+AACo92wKRd9884369u2r2bNn13U/AAAATmHTd595enoqJCSkrnsBAABwGptCUa9evXT8+PG67gUAAMBpbApF06ZNk8lk0vLly2U28+A8AABQ/9m0pujll1+Wn5+f1q5dq3/+858KDg6Wj49PpTqDwaA333yztj0CAADYnU2h6MiRI9bfFxQU6LvvvquyzmAw2NQUAACAo9kUilJTU+u6DwAAAKeyaU0RAABAQ2NTKCorK9ONGzd08+bN286XlZXVqjkAAABHsSkUJSQk6OGHH9b169ernM/Ly9PDDz+sjz76qFbNAQAAOIpNoWjv3r36+c9/roCAgCrnAwIC1LNnT6WlpdWmNwAAAIexKRRduHBBd999921rQkJCdOHCBZuaAgAAcDSbQlFpaanc3G6/qcFgUHFxsU1NAQAAOJpNoaht27b6+uuvb1vz9ddfq02bNjY1BQAA4Gg2haIBAwbo1KlTWrVqVaU7zMrKyvTBBx/o1KlTioyMrIseAQAA7M6mhzc+9thj2rFjhz766CPt2LFD999/v1q2bKns7GwdPnxY58+fV2hoqMaPH1/X/QIAANiFTaHIx8dHy5Yt05IlS7R7926dO3fOOufm5qaBAwdq9uzZVX4fGgAAgCuyKRRJt267f/nll5Wbm6uTJ08qPz9ffn5+6ty5s1q0aFGXPQIAANidzaGo3M9+9jP16dOnLnoBAABwGr77DAAAQDU8UxQfH6+BAwdq4MCBNr3JnWxfUFCgxMREZWRk6MSJE8rLy9MLL7yg4cOHV6o9c+aMli1bpmPHjsnDw0N9+vTRjBkzKj1p++bNm0pMTNTGjRuVm5ur4OBgxcXFKTo62qbjAQAADU+NzhTt2LFD//nPf2x+kzvZ/tq1a0pISFBWVpbCwsKqrbt8+bJmzpypc+fOadKkSRo/fry+/PJLzZ49WyUlJRVqV65cqffee0+9evXSc889p6CgIMXHx2vHjh02HxMAAGhYarymKDMzUykpKfbsRZJkNBq1YcMGGY1Gffvtt5o8eXKVdR9//LGKior0wQcfKCgoSJIUERGh2bNna8uWLYqNjZUkZWdnKykpSaNHj9asWbMkSSNHjtTMmTO1fPlyRUZGyt3d3e7HBQAAXFuNQ1FaWpr27Nlzx29gsVjuqN7Ly0tGo/En61JTU9W3b19rIJKknj17KiQkRDt37rSGorS0NJWWlmr06NHWOoPBoEceeUTx8fFKT09X9+7d76hHAADQ8NQoFM2bN6/Wb9SpU6da76Ncdna2rly5ovDw8EpzERER2rdvn/V1ZmammjZtqtDQ0Ep15fPVhaKcnByZTCbr66ysrLpoHwAAuKAahaKqFjk7U3lQqeqMktFo1PXr11VcXCwvLy+ZTCa1aNFCBoOhUp10K/hUJzk5WQkJCXXXOAAAcFm1fk6RM5jNZkmSp6dnpTkvLy9rjZeXl8xm80/WVSc2Nlb9+vWzvs7KytKiRYtq1TsAAHBN9TIUeXt7S1Klu8wkqbi4uEKNt7d3jeqqEhgYqMDAwFr3CwAAXF+9fHhj+aWvH6/3KWcymdSsWTPrmSCj0ajc3NxKC77LtyX0AAAAqZ6GopYtWyogIEAnT56sNHfixIkKzzcKCwtTUVFRpUXSGRkZ1nkAAIB6GYokaeDAgdq7d68uXbpkHTt06JDOnj2rQYMGWcf69+8vDw8PbdiwwTpmsVi0adMmtWzZUt26dXNo3wAAwDW55JqiTz/9VPn5+dZLXHv27NHly5clSY8++qj8/PwUFxenL774Qs8//7zGjBmjwsJCrVmzRh06dKhwt1yrVq00duxYrVmzRqWlpYqIiNDu3bt19OhRLViwgAc3AgAASS4aipKSknTx4kXr6127dmnXrl2SpCFDhsjPz09BQUFaunSpli1bphUrVli/+2z69OnW9UTlpkyZIn9/fyUnJyslJUXBwcGaP3++Bg8e7NDjAgAArsslQ9HatWtrVNe+fXstWbLkJ+vc3NwUFxenuLi42rYGAAAaqHq7pggAAKAuEYoAAABEKAIAAJBEKAIAAJBEKAIAAJBEKAIAAJBEKAIAAJBEKAIAAJBEKAIAAJBEKAIAAJDkol/zAQCoWuEVsz6Z8bmz26iXPJt4qOfYTmrfu42zW4GLIhQBQD3g2cRDklkWi1SQa3Z2O/WUWQfXZRKKUC1CEQDUAz3HdtLBdZkqKSp1div1UuGVW4GSnx9uh1AEAPVA+95tOMNRC5/M+JwzbPhJLLQGAAAQoQgAAEASoQgAAEASoQgAAEASoQgAAEASoQgAAEASoQgAAEASoQgAAEASoQgAAEASoQgAAEASoQgAAEASoQgAAEASoQgAAEASoQgAAEASoQgAAEASoQgAAEASoQgAAEASoQgAAEASoQgAAEASoQgAAEASoQgAAEASoQgAAEASoQgAAEASoQgAAEASoQgAAEASoQgAAECS5OHsBmx1+PBhPffcc1XOvfvuu+ratav19bFjx/Tee+/pu+++k6+vrwYNGqRJkybJx8fHUe0CAAAXV29DUblHH31UERERFcbatm1r/X1mZqZmzZql0NBQzZgxQ5cvX1ZSUpJ++OEHvfbaa45uFwAAuKh6H4p69OihyMjIaufff/99+fv7a+nSpfL19ZUktWnTRq+++qoOHDigBx54wEGdAgAAV9Yg1hQVFBSotLS00viNGzd08OBBDRkyxBqIJGno0KFq2rSpdu7c6cg2AQCAC6v3Z4peeeUVFRYWyt3dXd27d9fUqVPVuXNnSdLp06dVVlam8PDwCtt4enqqU6dOyszMvO2+c3JyZDKZrK+zsrLq/gAAAIBLqLehyMPDQwMHDtSDDz6o5s2b68yZM0pKStKMGTO0fPly3XPPPdZAYzQaK21vNBr1zTff3PY9kpOTlZCQYI/2AQCAi6m3oejee+/Vvffea33dv39/RUZG6qmnntL777+vxYsXy2w2S7p1Zuh/eXl5qbi4+LbvERsbq379+llfZ2VladGiRXV0BAAAwJXU21BUleDgYPXv31+7du1SWVmZvL29JUklJSWVaouLi+Xl5XXb/QUGBiowMNAuvQIAANfSIBZa/1irVq1UUlKioqIi62WzH68LKmcymQg8AADAqsGFovPnz8vLy0tNmzZV+/bt5e7urpMnT1aoKSkpUWZmpsLCwpzUJQAAcDX1NhRdvXq10tipU6e0Z88e9erVS25ubvLz81PPnj21bds2FRQUWOu2bt2qwsJCDRo0yIEdAwAAV1Zv1xT96U9/kre3t7p166YWLVrozJkz+uc//6kmTZpoypQp1rqJEydq+vTpmjlzpmJjY61PtO7Vq5d69+7txCMAAACupN6Gol/+8pfavn271q5dqxs3biggIEADBgzQhAkTFBwcbK0LDw/X66+/rvfee09vv/22fHx8FBMTUyE4AQAA1NtQNGbMGI0ZM6ZGtd27d9fy5cvt3BEAAKjP6u2aIgAAgLpEKAIAABChCAAAQBKhCAAAQBKhCAAAQBKhCAAAQBKhCAAAQBKhCAAAQBKhCAAAQBKhCAAAQBKhCAAAQBKhCAAAQBKhCAAAQBKhCAAAQBKhCAAAQBKhCAAAQBKhCAAAQJLk4ewGAABwlMIrZn0y43Nnt4FqeDbxUM+xndS+dxunvD+hCADQ4Hk28ZBklsUiFeSand0OqmXWwXWZhCIAAOyl59hOOrguUyVFpc5uBdUovHIrtDrz/xGhCADQ4LXv3cZpZx9QM5/M+NzpZ/FYaA0AACBCEQAAgCRCEQAAgCRCEQAAgCRCEQAAgCRCEQAAgCRCEQAAgCRCEQAAgCRCEQAAgCRCEQAAgCRCEQAAgCRCEQAAgCRCEQAAgCRCEQAAgCRCEQAAgCRCEQAAgCRCEQAAgCRCEQAAgCRCEQAAgCRCEQAAgCTJw9kNOEpxcbFWrVqlbdu2KS8vTx07dtTEiRPVq1cvZ7cGAABcQKM5U/TKK69o7dq1Gjx4sJ599lm5ublpzpw5Onr0qLNbAwAALqBRhKKMjAzt2LFDkydP1rRp0xQbG6s333xTrVu31rvvvuvs9gAAgAtoFKEoNTVV7u7uio2NtY55e3srJiZG6enpunTpkhO7AwAArqBRrCnKzMxUcHCwfH19K4xHRERIkk6dOqWgoKBK2+Xk5MhkMllfZ2Vl2bdRAADgNI0iFJlMJhmNxkrj5WM5OTlVbpecnKyEhAR7tgYAAFxEowhFZrNZnp6elca9vLys81WJjY1Vv379rK+zsrK0aNGiOu/Pp7l3hf8CANDYuMLfhY0iFHl7e6ukpKTSeHFxsXW+KoGBgQoMDLRrb5L0yF/6/XQRAAANmCv8XdgoFlobjcYKa4PKlY85IvgAAADX1ihCUVhYmH744QfduHGjwnhGRoZ1HgAANG6NIhRFRkaqrKxMycnJ1rHi4mJt3rxZXbp0qfLOMwAA0Lg0ijVFXbp00aBBg/T+++/r6tWratu2rVJSUnTx4kXNnTvX2e0BAAAX0ChCkSS9+OKLCgoK0tatW5Wfn68OHTro73//u+677z5ntwYAAFxAowlF3t7emjZtmqZNm+bsVgAAgAtqFGuKAAAAfgqhCAAAQIQiAAAASYQiAAAASYQiAAAASYQiAAAASYQiAAAASY3oOUV1wWw2S5KysrKc3AkAALhToaGhatKkSbXzhKI7cPHiRUnSokWLnNwJAAC4UytXrlR4eHi18waLxWJxYD/12tWrV3XgwAG1adNGXl5et619++23NXPmzBrtNysrS4sWLdL8+fMVGhpaF602Onfy83Y1rtC7o3qwx/vU1T5rsx9bt+VzwrFc4c+aLVylb0f0Ye/PCM4U1aGAgAANGTKkRrV+fn63TaNVCQ0NveNtcIstP29X4Qq9O6oHe7xPXe2zNvuxdVs+JxzLFf6s2cJV+nZEH87+jGChtZ1ER0c7u4VGpT7/vF2hd0f1YI/3qat91mY/tm7rCv/vG5P6+vN2lb4d0YezPyO4fOYCTp48qUmTJv3ktU4AjRefE4D9cabIBRiNRk2YMEFGo9HZrQBwUXxOAPbHmSIAAABxpggAAEASoQgAAEASt+TXG8XFxXr99dd18OBB5efnq127dpoxY4a6devm7NYAuIjXXntNe/bsUVFRkYKCgjR58mT169fP2W0B9QZriuqJwsJCJSUlafjw4WrZsqV27typN998U0lJSfLx8XF2ewBcQFZWlvXhsidOnNDs2bOVmJio5s2bO7s1oF7g8lk90bRpU02YMEFBQUFyc3NTVFSUPDw8dPbsWWe3BsBFhIaGWp+2bzAYVFJSopycHCd3BdQfXD6zk4KCAiUmJiojI0MnTpxQXl6eXnjhBQ0fPrxSbXFxsVatWqVt27YpLy9PHTt21MSJE9WrV69q93/27Fnl5eWpbdu29jwMAHZir8+I119/XZs3b1ZxcbEefPBBdejQwRGHAzQInCmyk2vXrikhIUFZWVkKCwu7be0rr7yitWvXavDgwXr22Wfl5uamOXPm6OjRo1XWm81mLVq0SL/5zW/k5+dnj/YB2Jm9PiNmz56trVu36o033lCvXr1kMBjsdQhAg0MoshOj0agNGzZo3bp1mjp1arV1GRkZ2rFjhyZPnqxp06YpNjZWb775plq3bq133323Un1paakWLlyotm3basKECXY8AgD2ZK/PCElyd3fXL37xCx06dEhffvmlvQ4BaHAIRXbi5eVVoyfPpqamyt3dXbGxsdYxb29vxcTEKD09XZcuXbKO37x5U4sWLZLBYNCLL77IvwCBeswenxH/q6ysTOfOnauTfoHGgFDkZJmZmQoODpavr2+F8YiICEnSqVOnrGOLFy+WyWTSSy+9JA8PloMBjUFNPyPy8/O1fft2FRQUqLS0VDt37tThw4fVo0cPh/cM1Ff8zepkJpOpyn8tlo+V3zly8eJF/etf/5KXl1eFfzG++uqrfOgBDVhNPyMMBoP+9a9/6Y033pDFYlHbtm21YMECderUyaH9AvUZocjJzGazPD09K42X31ZrNpslSa1bt9auXbsc2hsA56vpZ4Svr6/eeusth/YGNDRcPnMyb29vlZSUVBovLi62zgNovPiMAByHUORkRqNRJpOp0nj5WGBgoKNbAuBC+IwAHIdQ5GRhYWH64YcfdOPGjQrjGRkZ1nkAjRefEYDjEIqcLDIyUmVlZUpOTraOFRcXa/PmzerSpYuCgoKc2B0AZ+MzAnAcFlrb0aeffqr8/Hzrae49e/bo8uXLkqRHH31Ufn5+6tKliwYNGqT3339fV69eVdu2bZWSkqKLFy9q7ty5zmwfgJ3xGQG4FoPFYrE4u4mGaty4cbp48WKVc0lJSWrTpo2kW3ePlH+vUX5+vjp06KCJEyfqgQcecGS7AByMzwjAtRCKAAAAxJoiAAAASYQiAAAASYQiAAAASYQiAAAASYQiAAAASYQiAAAASYQiAAAASYQiAAAASYQiAAAASYQiAAAASYQiAI3UgAEDKvwym83WuS1btmjAgAHasmWLEzv8P5s2barQ61//+ldntwQ0SB7ObgBAw3bhwgU99thjt61p3bq11q5d66COKr7vsGHDJEnu7u52fa8DBw7oD3/4g3r16qUlS5bctjY+Pl6fffaZFixYoMGDBys8PFwTJkxQfn6+1q9fb9c+gcaMUATAIdq2bavBgwdXOefn5+fgbm5p3bq1nn76aYe8V8+ePRUUFKRDhw7p0qVLCgoKqrIuPz9fu3fvlp+fnwYMGCBJ6ty5szp37qwLFy4QigA7IhQBcIi2bds6LIC4Ijc3Nw0fPlwJCQlKSUnRk08+WWXdZ599JrPZrBEjRsjb29vBXQKNG2uKALicAQMG6Nlnn1V2drbi4+M1atQoDR06VHPmzNH58+clSWfOnNGLL76omJgYDR06VAsWLFBubq5d+7p8+bKefPJJRUdH64svvrCOX7lyRW+//bYef/xxRUVFadSoUZo/f75Onz5dYfsRI0bIYDBoy5YtslgsVb7H5s2bJUkxMTF2Ow4AVSMUAXBJeXl5mj59ui5cuKChQ4fq/vvv1759+zR79mydPn1a06ZNU2FhoUaMGKHOnTsrNTVVL730kt36OXPmjKZNm6bLly/rtddeU2RkpCTp3LlzmjhxotatW6e77rpLv/rVr/Tggw/qwIEDmjp1qjIyMqz7aN26tX7xi1/o/PnzOnz4cKX3OH36tL799lt16tRJ99xzj92OBUDVuHwGwCHOnTunf/zjH1XOde3aVb17964w9v3332vcuHGaMWOGdez111/Xxo0bNWPGDD311FMaO3asJMlisWju3Lnat2+fTp48qfDw8DrtPT09XXPnzpWHh4fefvtthYWFWef+8pe/KDc3V4sXL9YDDzxgHf/tb3+rSZMm6dVXX1VCQoJ1PCYmRgcPHtTmzZv185//vML7cJYIcC7OFAFwiHPnzikhIaHKX/v3769U37RpU02cOLHCWFRUlCSpefPmGjNmjHXcYDBY577//vs67fvLL7/UrFmz5O/vr+XLl1cIRN99952OHz+uoUOHVghEkhQSEqKRI0fq9OnTFS6j/fKXv1Tz5s2VmpqqGzduWMdLS0u1bds2eXl5VbsgHYB9caYIgEM88MADWrx4cY3rg4OD1aRJkwpjRqNRktShQwcZDIYq53JycmrZ6f/ZuXOnvvrqK3Xs2FGvvfaaWrRoUWG+/NLYlStXqjwL9t///tf63w4dOkiSNfSsX79en332mR5++GFJ0p49e3T16lVFR0fL39+/zo4BQM0RigC4JF9f30pj5c8Sut1caWlpnfWQnp6usrIyde/evVIgkqTr169LunU26csvv6x2P4WFhRVex8TEaP369dq8ebM1FHHpDHA+QhEAVGPy5MlKS0vT+vXr5e7urunTp1eYLw9nzz33nB599NEa77djx47q3LmzTpw4of/85z/y9/fXgQMH1KZNm0rrjAA4DmuKAKAaXl5e+stf/qI+ffooKSlJy5YtqzAfEREh6dYZpTtVfkbo3//+t7Zu3aqysjLrLfsAnINQBAC34eXlpUWLFqlv375au3at3n77betcly5d1KVLF+3YsUM7duyotO3Nmzd15MiRKvcbHR2tJk2aaNu2bdq8ebPc3NysXzkCwDm4fAbAIW53S74k/eY3v3HZJzh7enrq5Zdf1sKFC7Vu3TpZLBY9++yzkqSFCxfq+eef10svvaT169erU6dO8vb21uXLl3X8+HFdu3ZNn332WaV9+vr6auDAgdq6dauuXr2q3r17V/vVHwAcg1AEwCHKb8mvztixY102FEn/F4z+9Kc/af369bJYLHruued01113adWqVUpKStLu3bu1ZcsWubm5yWg0qkePHtaHPFYlJiZGW7dulXTradcAnMtgqe5Z8wDQgA0YMED33Xefli5d6uxWauzChQt67LHHNGzYML344ovObgdocDhTBKDROnLkiPWb6Ldv3+6yZ6o2bdqkJUuWOLsNoMEjFAFolCZMmFDhdflzjlxReHh4hX47derkvGaABozLZwAAAOKWfAAAAEmEIgAAAEmEIgAAAEmEIgAAAEmEIgAAAEmEIgAAAEmEIgAAAEmEIgAAAEmEIgAAAEnS/wftqm95h4pTIwAAAABJRU5ErkJggg==",
+ "text/plain": [
+ "
"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "signal.project('Em').plot()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "05960a53-f5ee-41f0-a3d3-1cb41d141f64",
+ "metadata": {},
+ "source": [
+ "This shape is a combination of the spectrum, the energy resolution, and the effective area of the detector as a function of energy.\n",
+ "\n",
+ "We can get the total number of events we expect from the GRB by summing over all bins:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "id": "08eb4fae-d2e7-4336-82ca-cf8b301a1252",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "698.1483813523819"
+ ]
+ },
+ "execution_count": 11,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "np.sum(signal)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "aea9bc5f-3a6d-45ab-96ae-bd1d767f0df8",
+ "metadata": {},
+ "source": [
+ "Now let's explore the CDS. It's easier to visualize if we take slices in energy and scatter angle. For reference, these are the bin edges:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 12,
+ "id": "7b7bc781-6100-4123-b805-d4c3742aff04",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/latex": [
+ "$[100,~200,~500,~1000,~2000,~5000] \\; \\mathrm{keV}$"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "execution_count": 12,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "signal.axes['Em'].edges"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 13,
+ "id": "1339c2c6-48b4-4ce0-9c6f-c33263b6bc03",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/latex": [
+ "$[0,~10,~20,~30,~40,~50,~60,~70,~80,~90,~100,~110,~120,~130,~140,~150,~160,~170,~180] \\; \\mathrm{{}^{\\circ}}$"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "execution_count": 13,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "signal.axes['Phi'].edges"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "1d4b2946-6165-4d08-b045-d80f9292748a",
+ "metadata": {},
+ "source": [
+ "This is the plot of the distribution of events within the energy range 1-2 MeV (bin 3) and scattered angles between 40-50deg (bin 4):"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 19,
+ "id": "db0b2090-2ad3-4807-87bc-9269b214d68d",
+ "metadata": {},
+ "outputs": [
+ {
+ "ename": "TypeError",
+ "evalue": "matplotlib.transforms.Bbox.from_bounds() argument after * must be an iterable, not int",
+ "output_type": "error",
+ "traceback": [
+ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+ "\u001b[0;31mTypeError\u001b[0m Traceback (most recent call last)",
+ "Cell \u001b[0;32mIn[19], line 8\u001b[0m\n\u001b[1;32m 5\u001b[0m fig \u001b[38;5;241m=\u001b[39m plt\u001b[38;5;241m.\u001b[39mfigure(dpi \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m150\u001b[39m)\n\u001b[1;32m 7\u001b[0m \u001b[38;5;66;03m# Try also other projections, e.g. projection = 'mollview'\u001b[39;00m\n\u001b[0;32m----> 8\u001b[0m ax \u001b[38;5;241m=\u001b[39m \u001b[43mfig\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43madd_subplot\u001b[49m\u001b[43m(\u001b[49m\u001b[43mprojection\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43morthview\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m)\u001b[49m\n\u001b[1;32m 10\u001b[0m m_signal\u001b[38;5;241m.\u001b[39mplot(ax, coord \u001b[38;5;241m=\u001b[39m SpacecraftFrame(attitude \u001b[38;5;241m=\u001b[39m Attitude\u001b[38;5;241m.\u001b[39midentity()))\n\u001b[1;32m 12\u001b[0m \u001b[38;5;66;03m# Location of the source\u001b[39;00m\n",
+ "File \u001b[0;32m~/.virtualenvs/cosipy/lib/python3.9/site-packages/matplotlib/figure.py:757\u001b[0m, in \u001b[0;36mFigureBase.add_subplot\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m 754\u001b[0m args \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mtuple\u001b[39m(\u001b[38;5;28mmap\u001b[39m(\u001b[38;5;28mint\u001b[39m, \u001b[38;5;28mstr\u001b[39m(args[\u001b[38;5;241m0\u001b[39m])))\n\u001b[1;32m 755\u001b[0m projection_class, pkw \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_process_projection_requirements(\n\u001b[1;32m 756\u001b[0m \u001b[38;5;241m*\u001b[39margs, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs)\n\u001b[0;32m--> 757\u001b[0m ax \u001b[38;5;241m=\u001b[39m \u001b[43mprojection_class\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mpkw\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 758\u001b[0m key \u001b[38;5;241m=\u001b[39m (projection_class, pkw)\n\u001b[1;32m 759\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_add_axes_internal(ax, key)\n",
+ "File \u001b[0;32m~/.virtualenvs/cosipy/lib/python3.9/site-packages/mhealpy/plot/axes.py:259\u001b[0m, in \u001b[0;36mOrthview.__init__\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m 257\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m__init__\u001b[39m(\u001b[38;5;28mself\u001b[39m, \u001b[38;5;241m*\u001b[39margs, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs):\n\u001b[0;32m--> 259\u001b[0m \u001b[38;5;28;43msuper\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[38;5;21;43m__init__\u001b[39;49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 260\u001b[0m \u001b[43m \u001b[49m\u001b[43mframe_class\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m \u001b[49m\u001b[43mkwargs\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mpop\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mframe_class\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mEllipticalFrame\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 261\u001b[0m \u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n",
+ "File \u001b[0;32m~/.virtualenvs/cosipy/lib/python3.9/site-packages/mhealpy/plot/axes.py:41\u001b[0m, in \u001b[0;36mHealpyAxes.__init__\u001b[0;34m(self, fig, rect, coord, flip, rot, **kwargs)\u001b[0m\n\u001b[1;32m 24\u001b[0m \u001b[38;5;250m\u001b[39m\u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m 25\u001b[0m \u001b[38;5;124;03mBase class for WCSAxes that behave similar to healpy's projections.\u001b[39;00m\n\u001b[1;32m 26\u001b[0m \n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 37\u001b[0m \u001b[38;5;124;03m of angle psi around this direction is applied.\u001b[39;00m\n\u001b[1;32m 38\u001b[0m \u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m 40\u001b[0m \u001b[38;5;66;03m# Get equivalent WCS FITS header\u001b[39;00m\n\u001b[0;32m---> 41\u001b[0m naxis1,naxis2 \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_get_naxis\u001b[49m\u001b[43m(\u001b[49m\u001b[43mfig\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mrect\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 43\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_autoscale:\n\u001b[1;32m 44\u001b[0m \u001b[38;5;66;03m# Gnom specifies the reosolution and the limits are adjusted based\u001b[39;00m\n\u001b[1;32m 45\u001b[0m \u001b[38;5;66;03m# on the figure sizes\u001b[39;00m\n\u001b[1;32m 46\u001b[0m \u001b[38;5;66;03m# For the others the limits are fixed and the resolution is adjusted\u001b[39;00m\n\u001b[1;32m 47\u001b[0m \u001b[38;5;66;03m# by scaling a normalized appropiate resolution\u001b[39;00m\n\u001b[1;32m 48\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_cdelt \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m=\u001b[39m \u001b[38;5;241m360\u001b[39m\u001b[38;5;241m/\u001b[39mnaxis1\n",
+ "File \u001b[0;32m~/.virtualenvs/cosipy/lib/python3.9/site-packages/mhealpy/plot/axes.py:140\u001b[0m, in \u001b[0;36mHealpyAxes._get_naxis\u001b[0;34m(self, fig, rect)\u001b[0m\n\u001b[1;32m 133\u001b[0m \u001b[38;5;250m\u001b[39m\u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m 134\u001b[0m \u001b[38;5;124;03mReturn approapiate NAXIS1 and NAXIS2 for a given figure and Bbox (or bounds)\u001b[39;00m\n\u001b[1;32m 135\u001b[0m \u001b[38;5;124;03m \u001b[39;00m\n\u001b[1;32m 136\u001b[0m \u001b[38;5;124;03maspect = naxis2/naxis1\u001b[39;00m\n\u001b[1;32m 137\u001b[0m \u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m 139\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(rect, Bbox):\n\u001b[0;32m--> 140\u001b[0m rect \u001b[38;5;241m=\u001b[39m \u001b[43mBbox\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mfrom_bounds\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mrect\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 142\u001b[0m naxis1 \u001b[38;5;241m=\u001b[39m fig\u001b[38;5;241m.\u001b[39mget_figwidth() \u001b[38;5;241m*\u001b[39m fig\u001b[38;5;241m.\u001b[39mdpi \u001b[38;5;241m*\u001b[39m rect\u001b[38;5;241m.\u001b[39mwidth\n\u001b[1;32m 143\u001b[0m naxis2 \u001b[38;5;241m=\u001b[39m fig\u001b[38;5;241m.\u001b[39mget_figheight() \u001b[38;5;241m*\u001b[39m fig\u001b[38;5;241m.\u001b[39mdpi \u001b[38;5;241m*\u001b[39m rect\u001b[38;5;241m.\u001b[39mheight\n",
+ "\u001b[0;31mTypeError\u001b[0m: matplotlib.transforms.Bbox.from_bounds() argument after * must be an iterable, not int"
+ ]
+ },
+ {
+ "data": {
+ "text/plain": [
+ "
"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "#Since `PsiChi` is encoded as pixel in a HEALPix grid, we need mhealpy to render it back to a sphere\n",
+ "m_signal = HealpixMap(signal.slice[{'Em':3, 'Phi':4}].project('PsiChi').todense().contents,\n",
+ " coordsys = SpacecraftFrame())\n",
+ "\n",
+ "fig = plt.figure(dpi = 150)\n",
+ "\n",
+ "# Try also other projections, e.g. projection = 'mollview'\n",
+ "ax = fig.add_subplot(projection = 'orthview')\n",
+ "\n",
+ "m_signal.plot(ax, coord = SpacecraftFrame(attitude = Attitude.identity()))\n",
+ "\n",
+ "# Location of the source\n",
+ "ax.scatter(coord.ra, coord.dec, transform=ax.get_transform('world'), color = 'red')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "290a9873-c011-4bbd-a2e5-3605f0a56931",
+ "metadata": {},
+ "source": [
+ "This is a horizontal slice of the Compton cone shown in the figure above, spread by detector effects and the finite size of the `Em` and `Phi` bins. Try selecting different `Phi` bins to see how these circle grows or shrinks, and relate that to the CDS figure.\n",
+ "\n",
+ "You can also try selecting different energy bins. The opening of the cone in the CDS is geometrically constrained and does not depend on the energy. This circle becomes more blurry at different energies though, which is related to the energy resolution and the bin width."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "79b32eb9-b8d4-49c1-9c38-52def6410278",
+ "metadata": {},
+ "source": [
+ "## Getting a fake background"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "b90c008e-40b3-46b0-a20d-fd3d85c7c58c",
+ "metadata": {},
+ "source": [
+ "The background from Compton telescopes can be complex, and in general we need to either simulate all the different components with MEGAlib and/or use real data to constrain it. For the purpose of having a toy background that we can use to develop our algorithms, let's use the detector response to simulate an (unrealistic) isotropic gamma-ray background. The final source injector should use a background model as input instead.\n",
+ "\n",
+ "We'll repurpose the point source convolution by generating an effective dwell time map with the same value for all pixels. Since all pixels have the same area, this is simulating an isotropic distribution"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 15,
+ "id": "4700ced8-79da-4d5f-87e3-8771f050caa4",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "iso_map = HealpixMap(base = response, \n",
+ " unit = u.s, \n",
+ " coordsys = SpacecraftFrame())\n",
+ "\n",
+ "# Filling all pixels with a constant. The actual value doesn't\n",
+ "# since we will renormalize it\n",
+ "iso_map[:] = 1*u.s"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 16,
+ "id": "b0dcee61-d3e5-478c-a2a5-206b80812f49",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Non-realistic spectrum\n",
+ "bkg_spectrum = Powerlaw(K = 1, index=-2)\n",
+ "\"\"\"bkg_index = -2\n",
+ "bkg_piv = 1 * u.keV\n",
+ "bkg_K = 1 / u.cm / u.cm / u.s / u.keV\n",
+ "bkg_spectrum.index.value = bkg_index\n",
+ "bkg_spectrum.K.value = bkg_K.value\n",
+ "bkg_spectrum.piv.value = bkg_piv.value\n",
+ "bkg_spectrum.K.unit = bkg_K.unit\n",
+ "bkg_spectrum.piv.unit = bkg_piv.unit\n",
+ "\"\"\" \n",
+ "iso_response = response.get_point_source_response(iso_map)\n",
+ " \n",
+ "bkg = iso_response.get_expectation(bkg_spectrum).project(['Em', 'Phi', 'PsiChi'])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "c0de926a-f12e-4c79-b7df-3cd94edd0c88",
+ "metadata": {},
+ "source": [
+ "Now, let's renormalize the background to a total rate of 1k Hz. This is again not realistic, but was chosen such that the signal will show up clearly enough above the background and we can work on the algorihtms."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 20,
+ "id": "02ae055f-9549-4f38-a8e6-d0e57af60cbe",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "90238.31928090738"
+ ]
+ },
+ "execution_count": 20,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "np.sum(bkg)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 21,
+ "id": "5c0cc346-7bd6-4d67-b98b-e94bcf5569aa",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "bkg = bkg * 1e3 / np.sum(bkg)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 22,
+ "id": "3bf9becd-f9af-4259-b3af-22d0aa3e100a",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "1000.0000000000001"
+ ]
+ },
+ "execution_count": 22,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "np.sum(bkg)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "71024a7d-1742-4a57-9b23-cd47bc1ae603",
+ "metadata": {},
+ "source": [
+ "These are the same plots as we did for the signal, so you can compare:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 23,
+ "id": "aa8589a2-cff1-4c68-bae8-c78c1691a0fb",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "(,\n",
+ " )"
+ ]
+ },
+ "execution_count": 23,
+ "metadata": {},
+ "output_type": "execute_result"
+ },
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkUAAAG7CAYAAADNDuE1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAwbUlEQVR4nO3deXRV1d3/8XdCSBiCQhMMSChzQ9CitqAyFEFAgUjUClhbVGoBy2SV9ifoEtsi1lbBARGLQh/6ODG5FNoyKcUoKFIRlEkMpaYUGcNgAiEhcH9/sHIf0wQMlyT3Bt6vtVh6997nnO8Ji8uHffY5JyoQCASQJEk6z0WHuwBJkqRIYCiSJEnCUCRJkgQYiiRJkgBDkSRJEmAokiRJAgxFkiRJgKHojBw9epQtW7Zw9OjRcJciSZLKmaHoDGRlZTFkyBCysrLCXYokSSpnhiJJkiQMRZIkSYChSJIkCTAUSZIkAYYiSZIkwFAkSZIEGIokSZIAQ5EkSRJgKJIkSQIMRZIkSYChSJIkCTAUSZIkAYYiSZIkwFAkSZIEGIokSZIAiAl3AdLZ2rZqJ2vmZXLsaGG4S9E3qF4jhnb9W9HsqobhLkWSSjAUqcpbMy+TQ18eDncZKpN8PpqbaSiSFJEMRaryimaIoqKgZr24MFejU8k7kE8ggDN6kiKWoUjnjJr14vjxlGvDXYZO4dWRf+fI/vxwlyFJp+RCa0mSJAxFkiRJgKFIkiQJMBRJkiQBhiJJkiTAUCRJkgQYiiRJkgBDkSRJEmAokiRJAgxFkiRJgKFIkiQJMBRJkiQBhiJJkiTAUCRJkgQYiiRJkgBDkSRJEmAokiRJAgxFkiRJgKFIkiQJMBRJkiQBhiJJkiTAUCRJkgQYiiRJkgBDkSRJEmAokiRJAgxFkiRJgKFIkiQJMBRJkiQBhiJJkiTAUCRJkgQYiiRJkgBDkSRJEmAokiRJAgxFkiRJgKFIkiQJMBRJkiQBhiJJkiTAUCRJkgQYiiRJkgBDkSRJEmAokiRJAgxFkiRJgKFIkiQJgJhwF/BN/vd//5fp06fTrFkz/vznPxfrW79+PX/84x/5/PPPqV27Nt26dWPIkCHUqlWr2LiCggJmzJjB0qVLycnJoUWLFgwePJj27dtX5qlIkqQIFtEzRXv27OHll1+mZs2aJfoyMzO57777OHr0KCNHjiQtLY2//OUv/PrXvy4x9rHHHmPOnDn07NmTe+65h+joaO6//34+/fTTyjgNSZJUBUT0TNHUqVNp06YNJ06c4NChQ8X6XnjhBerUqcPkyZOpXbs2AA0bNuTxxx9n9erVXHnllQBs2rSJZcuWMWzYMG677TYArr/+egYNGsTzzz/P888/X7knJUmSIlLEzhStW7eOjIwMRo0aVaLv8OHDfPTRR1x33XXBQAQnw07NmjVZvnx5sC0jI4Nq1aqRnp4ebIuLiyMtLY2NGzeye/fuij0RSZJUJUTkTNHx48d55plnSEtLo0WLFiX6t23bxvHjx0lJSSnWXr16dVq1akVmZmawLTMzk+Tk5GLhCSA1NRWArVu3kpSUVGod+/btIzs7O/g5Kysr5HOSJEmRLSJD0fz589m9ezdPPfVUqf1FQSUhIaFEX0JCAp988kmxsacaByeDz6ksWLCAmTNnnknpkiSpioq4UHTo0CH+9Kc/cccdd1C3bt1Sx+Tn5wMnZ4b+W2xsLAUFBcXGnmrc1/dVmvT0dDp16hT8nJWVxYQJE8p0HpIkqWqJuFA0ffp06tSpwy233HLKMXFxcQAcO3asRF9BQUEw8BSNPdW4r++rNImJiSQmJpa5dkmSVHVFVCjavn07f/nLXxg1alSxy1oFBQUUFhayc+dOateuHbz09fX1PkWys7OLBZmEhAT27t1b6jjA0CNJkoAIC0X79u3jxIkTPPPMMzzzzDMl+m+99Vb69evHXXfdRbVq1diyZQvXXnttsP/YsWNkZmbSrVu3YFvLli1Zu3Ythw8fLrbYetOmTcF+SZKkiApFzZo149FHHy3RPn36dI4cOcI999zDxRdfTHx8PO3atWPp0qXceeedwSdYL1myhLy8vGKhqGvXrsyaNYsFCxYEn1NUUFDAwoULadOmzSnvPJMkSeeXiApFdevW5Qc/+EGJ9rlz5wIU6xs8eDAjRoxg1KhRpKens2fPHmbPnk379u256qqrguPatGlDt27deOGFFzh48CCNGjVi8eLF7Nq1izFjxlT8SUmSpCohokLRmUhJSeHJJ5/kj3/8I88++yy1atUiLS2Nu+++u8TYBx98kKSkJJYsWUJubi7NmzfnD3/4A5dffnnlFy5JkiJSlQhFkydPLrW9bdu2TJ069Ru3j4uLY/jw4QwfPry8S5MkSeeIiH3NhyRJUmUyFEmSJGEokiRJAgxFkiRJgKFIkiQJMBRJkiQBhiJJkiTAUCRJkgQYiiRJkgBDkSRJEmAokiRJAgxFkiRJgKFIkiQJMBRJkiQBhiJJkiTAUCRJkgQYiiRJkgBDkSRJEmAokiRJAgxFkiRJgKFIkiQJMBRJkiQBhiJJkiTAUCRJkgQYiiRJkgBDkSRJEmAokiRJAgxFkiRJgKFIkiQJMBRJkiQBhiJJkiTAUCRJkgQYiiRJkgBDkSRJEmAokiRJAgxFkiRJgKFIkiQJMBRJkiQBhiJJkiTAUCRJkgQYiiRJkgBDkSRJEmAokiRJAgxFkiRJgKFIkiQJMBRJkiQBhiJJkiTAUCRJkgQYiiRJkgBDkSRJEmAokiRJAgxFkiRJgKFIkiQJMBRJkiQBEFOWQbfeeutZH6h///7069fvrPcjSZJUEcoUinbt2kXt2rWJj48P6SB79uwhNzc3pG0lSZIqQ5lCEcCAAQMYNGhQSAe55pprQtpOkiSpspQ5FFWWf/3rX/zP//wPW7ZsYf/+/dSoUYMmTZpw22230alTp2Jjv/jiC6ZMmcL69euJiYmhQ4cOjBw5krp16xYbd+LECWbNmsWbb77J/v37SU5OZuDAgfTo0aMSz0ySJEWyMoWil156iQsvvDDkg5zJ9rt27eLIkSP06tWLxMREjh49SkZGBg888AC/+tWvSE9PB05ekhs1ahTx8fEMGTKEvLw8Zs2axbZt25g2bRrVq1cP7vPFF1/klVdeoW/fvrRu3ZoVK1Ywfvx4oqKi6N69e8jnJUmSzh1lCkXf/va3z+ogZ7J9hw4d6NChQ7G2H/7whwwZMoQ5c+YEQ9HLL7/M0aNHmT59OklJSQCkpqYyevRoFi1aFBy3d+9eZs+ezc0338x9990HwA033MCoUaOYOnUqXbt2pVq1amd1fpIkqeor91vyCwoKKCwsLNd9VqtWjYsuuqjYYu2MjAw6duwYDEQA7dq1o3HjxixfvjzYtmLFCgoLC7n55puDbVFRUdx0003s3buXjRs3lmutkiSpagppTdG6detYs2YNAwYMoE6dOgAcOnSICRMm8NFHHxETE8Mtt9zCz3/+85ALy8vLIz8/n8OHD7Ny5Uo+/PBDunXrBpyc/Tlw4AApKSkltktNTWXVqlXBz5mZmdSsWZMmTZqUGFfU37Zt21Jr2LdvH9nZ2cHPWVlZIZ+PJEmKbCGFolmzZpGVlcXPfvazYNtzzz3H6tWradSoUXB9z3e+8x2uvfbakAp77rnnWLBgAQDR0dF06dIlePmrKKgkJCSU2C4hIYGvvvqKgoICYmNjyc7Opl69ekRFRZUYByeDz6ksWLCAmTNnhlS/JEmqWkIKRZmZmXz/+98Pfs7Pz2f58uW0b9+eiRMncuTIEQYNGsT8+fNDDkX9+/ena9eu7Nu3j+XLl3P8+HGOHTsWPB5QbDF1kdjY2OCY2NhY8vPzv3HcqaSnpxe74y0rK4sJEyaEdD6SJCmyhRSKvvrqK+rXrx/8vHHjRgoKCujduzcAtWrVomPHjmRkZIRcWJMmTYKXvHr16sXo0aMZO3Ys06ZNIy4uDiAYkr6uoKAAIDgmLi6uTONKk5iYSGJiYsjnIEmSqo6QFlrHxcVx5MiR4Oe1a9cSFRXF5ZdfHmyrWbMmOTk5Z11gka5du/LZZ5+xffv24KWvr6/3KZKdnc0FF1wQnAlKSEhg//79BAKBEuMAQ48kSQJCDEWNGjXiww8/pKCggGPHjrFs2TKaNGlSbI3P7t27qVevXrkVWnSZKzc3l/r161O3bl22bNlSYtzmzZtp2bJl8HPLli05evRoiUXSmzZtCvZLkiSFFIr69u3Ljh07uO2227j99tv58ssv6dOnT7Exn3/+OU2bNj3jfR84cKBEW2FhIUuWLCEuLi64z2uuuYb333+f3bt3B8etWbOG7du3B+9SA+jcuTMxMTG88cYbwbZAIMD8+fOpX78+l1566RnXKEmSzj0hrSlKS0tj+/btLFy4kIKCAm688Ub69+8f7N+wYQPbt28nLS3tjPc9ceJEDh8+zGWXXUb9+vXJzs7mrbfe4t///jcjRoygVq1aAAwcOJB33nmHe++9l379+pGXl8drr71G8+bNg2ubAC666CL69+/Pa6+9RmFhIampqbz33nt8+umnjBs3zgc3SpIkIMRQFBUVxbBhwxg2bFip/SkpKfztb3+jRo0aZ7zva6+9lr/97W/Mnz+fQ4cOUatWLVJSUvj5z39O586dg+OSkpKYPHkyU6ZMYdq0acF3n40YMSK4nqjI3XffTZ06dViwYAGLFy8mOTmZhx56iJ49e55xfZIk6dxUIS+ErV69eqm3wZdF9+7dy/w+smbNmjFp0qRvHBcdHc3AgQMZOHBgSDVJkqRzX7m/5kOSJKkqMhRJkiRhKJIkSQIMRZIkSYChSJIkCTiDULRx48aKrEOSJCmsynxL/vDhw2natClpaWlcd9111K1btwLLkiRJqlxlnim65JJL+OKLL5g6dSr9+vXj4Ycf5sMPPyzxolVJkqSqqMwzRVOnTmX79u389a9/5a233iIjI4N3332XhIQEevfuTZ8+fbj44osrslZJkqQKc0YLrRs3bsywYcOYN28ejz32GJ06deLgwYO89NJL/PjHP+bee+/lrbfeoqCgoKLqlSRJqhAhveYjOjqajh070rFjRw4ePMiSJUtYuHAha9euZd26dTz99NP06NGDPn36kJKSUt41S5IklbuzviW/bt263Hrrrfz5z3/mj3/8I+np6QDMnz+fn//852ddoCRJUmUo1+cUJScn06RJExITEwkEAi7CliRJVUZIl8/+2z/+8Q8WLlzIihUrOHbsGIFAgMsuu4y0tLTy2L0kSVKFCzkU7dy5k0WLFrF48WL27NlDIBAgISGBXr160adPH5KTk8uzTkmSpAp1RqGooKCAd955h4ULF/LJJ59w4sQJqlWrRqdOnUhLS+Pqq68mOto3h0iSpKqnzKFo0qRJLFu2jCNHjhAIBGjcuDF9+vShV69efOtb36rIGiVJkipcmUPRggULqFGjBtdffz1paWm0bdu2IuuSJEmqVGUORb/61a/o3r07tWrVqsh6JEmSwqLMoahv376lth85coTt27dz9OhRLrvssnIrTJIkqTKd1d1nkydPZtWqVQQCAaKioli+fDkA69ev5/HHH2f06NFcccUV5VasJElSRQnpVrHdu3czbNgwVq1aRefOnbnkkkuKPagxNTWVQ4cO8fbbb5dboZIkSRUppFD0pz/9iZycHCZPnswjjzxCu3btivXHxMTQtm1bNmzYUC5FSpIkVbSQQtHq1av5wQ9+wHe/+91TjmnQoAF79+4NuTBJkqTKFFIoysnJoUGDBqcdEwgEOHbsWEhFSZIkVbaQQlG9evX4z3/+c9ox27ZtIykpKaSiJEmSKltIoahdu3Z88MEH/POf/yy1/5NPPuHjjz/m6quvPqviJEmSKktIt+TfcccdZGRkMGrUKH70ox+xY8cOAFatWsWGDRuYM2cOF154IT/60Y/KtVhJkqSKElIoatiwIRMnTuQ3v/kNM2bMICoqikAgwNixYwkEAiQlJTF+/HgSExPLu15JkqQKEfLDG9u0acOrr77K+++/z6ZNm8jJyaFWrVq0adOGzp07U7169fKsU5IkqUKFFIrWr1/Pd7/7XWJiYujSpQtdunQpddycOXMYMGDAWRV4Pti2aidr5mVy7GhhuEupkvIO5Ie7BEnSOSCkUPTAAw8wZcoUmjZtesoxc+bMYerUqYaiMlgzL5NDXx4OdxlVXvUaIU98SpIUWiiKi4vj//2//8fUqVOpX79+if558+bx3HPPcemll551geeDohmiqCioWS8uzNVUTdVrxNCuf6twlyFJqsJCCkVPPPEEI0eO5Je//CVTpkzhggsuCPa9/vrrPPvss7Rp04Ynnnii3Ao9H9SsF8ePp1wb7jIkSTovhfScoubNm/P73/+eXbt2MWbMGPLzT67peP3115k8eTKpqalMmjSJWrVqlWuxkiRJFSWkUATQtm1bHn74YT777DPGjRvH3LlzmTx5MikpKUycONFAJEmSqpSQQxFA586d+eUvf8mHH37Ic889R0pKCk899RTx8fHlVZ8kSVKlKNOaonXr1p2yLzk5mU6dOrF+/Xpuv/12MjMzi/VffvnlZ1OfJElSpShTKPrFL35BVFTUaccEAgHGjRtXov2dd94JqTBJkqTKVKZQdOedd35jKJIkSarKyhSK7rrrroquQ5IkKazOaqG1JEnSucJQJEmShKFIkiQJMBRJkiQBhiJJkiTAUCRJkgQYiiRJkgBDkSRJEhBCKMrOziYjI4MVK1aQk5NzynHr1q1j5syZZ1ObJElSpSnTE62LzJo1i+nTp1NYWAhAbGwst99+OwMHDizxGpC1a9fy5z//mUGDBpVbsZIkSRWlzDNFq1ev5vnnnyc2NpYbbriBm266iVq1ajFjxgzGjh1LQUFBRdYpSZJUoco8UzR37lxq1KjBtGnTaNy4MQBDhw5l4sSJLFu2jLFjx/LYY48RFxdXYcVKkiRVlDLPFH322Wd06dIlGIgAatWqxcMPP8yPf/xj1qxZw9ixY8nPz6+QQiVJkipSmUNRXl4eF110Ual9d999N7fffjsff/wxY8aMMRhJkqQqp8yXzxITE9m7d+8p+wcPHgzASy+9xP33309KSsrZVydJklRJyhyKmjVrxpo1a0475uvBaMOGDWdXmSRJUiUq8+WzDh06sG/fPj744IPTjhs8eDB33HFH8LZ9SZKkqqDMM0Vdu3YlEAhQo0aNbxz7s5/9jIsvvphdu3adcUGbN29m8eLFrF27ll27dnHBBRdwySWXMHjw4GKLvAG++OILpkyZwvr164mJiaFDhw6MHDmSunXrFht34sQJZs2axZtvvsn+/ftJTk5m4MCB9OjR44zrkyRJ56Yyh6ILLriAG2+8scw77t27d0gFvfrqq6xfv55u3brRokULsrOzeeONNxg8eDDPP/88zZs3B2DPnj2MGjWK+Ph4hgwZQl5eHrNmzWLbtm1MmzaN6tWrB/f54osv8sorr9C3b19at27NihUrGD9+PFFRUXTv3j2kOiVJ0rnljJ5oXRkGDBjAww8/XCzUXHvttfz0pz/llVdeYdy4cQC8/PLLHD16lOnTp5OUlARAamoqo0ePZtGiRaSnpwOwd+9eZs+ezc0338x9990HwA033MCoUaOYOnUqXbt2pVq1apV8lpIkKdKcVSjKzMxk69atZGdnl7qGKCoqijvvvPOM9vnd7363RFvjxo1p2rQpWVlZwbaMjAw6duwYDEQA7dq1o3HjxixfvjwYilasWEFhYSE333xzsbpuuukmxo8fz8aNG2nbtu0Z1ShJks49IYWiAwcOMH78eNauXQtAIBAodVwooag0gUCAAwcO0LRpU+Dk7M+BAwdKve0/NTWVVatWBT9nZmZSs2ZNmjRpUmJcUf+pQtG+ffvIzs4Ofv56KJMkSeeWkELRU089xccff8zVV19N9+7dSUhIqNBLUG+99RZ79+7lrrvuAggGlYSEhBJjExIS+OqrrygoKCA2Npbs7Gzq1atX4oW1Rdvu27fvlMddsGABM2fOLKezkCRJkSykULR69WquuOIK/vCHP5R3PSVkZWXx1FNPcckll9CrVy+A4BOzv77uqEhsbGxwTGxsLPn5+d847lTS09Pp1KlTsVomTJgQ+slIkqSIFVIoiomJqZQnVmdnZzNmzBhq167NI488EpyNKnrp7LFjx0psU1BQUGxMXFxcmcaVJjExkcTExLM7CUmSVCWU+eGNX9e2bVsyMzPLu5ZicnNzuf/++8nNzWXixInFwknRpa+vr/cpkp2dzQUXXBCcCUpISGD//v0l1j0VbWvokSRJEGIoGjp0KFu2bOH1118v73qAk5e0xo4dy/bt2/n9738fXGBdpH79+tStW5ctW7aU2Hbz5s20bNky+Llly5YcPXq0xCLpTZs2BfslSZJCunzWtGlTpkyZwsiRI3n99ddp0aIFtWvXLnXs2LFjz2jfx48f5ze/+Q0bN27kd7/7HZdeemmp46655hoWL17M7t27g7flr1mzhu3btzNgwIDguM6dOzNlyhTeeOON4HOKAoEA8+fPp379+qfcvyRJOr+EFIq+/PJLHnzwQXJzc8nNzWXHjh2ljouKijrjUPTcc8+xcuVKOnbsSE5ODkuXLi3Wf9111wEwcOBA3nnnHe6991769etHXl4er732Gs2bNy/2NO2LLrqI/v3789prr1FYWEhqairvvfcen376KePGjfPBjZIkCQgxFD3zzDN8+eWX3HjjjfTo0aNcb8nfunUrAO+//z7vv/9+if6iUJSUlMTkyZOZMmUK06ZNC777bMSIEcH1REXuvvtu6tSpw4IFC1i8eDHJyck89NBD9OzZs1xqliRJVV9IoeiTTz6hY8eOjB49urzrYfLkyWUe26xZMyZNmvSN46Kjoxk4cCADBw48m9IkSdI5LKSF1tWrVy/xxnpJkqSqLKSZovbt27Nhw4byrkXSeSDvQD6vjvx7uMuokqrXiKFd/1Y0u6phuEuRzkkhhaLhw4czYsQIpk6dys9+9rPTPgBRkuDkX+iQTyAAR/af+knyOp18PpqbaSiSKkhIoeiRRx4hPj6eOXPm8Je//IXk5GRq1apVYlxUVBRPP/302dYo6RzQrn8rPpqbybGjheEupUrKO3AyUPrzkypOSKFo3bp1wf8/cuQIn3/+eanj/vslrJLOX82uaugMx1l4deTfnWGTKlhIoSgjI6O865AkSQqrkO4+kyRJOteEFIqOHz/O4cOHOXHixGn7jx8/flbFSZIkVZaQQtHMmTO58cYb+eqrr0rtz8nJ4cYbb+Sll146q+IkSZIqS0ih6P333+d73/sedevWLbW/bt26tGvXjhUrVpxNbZIkSZUmpFC0c+dOvv3tb592TOPGjdm5c2dIRUmSJFW2kEJRYWEh0dGn3zQqKoqCgoKQipIkSapsIYWiRo0a8fHHH592zMcff0zDhj6TRJIkVQ0hhaIuXbqwdetWZsyYUeIOs+PHjzN9+nS2bt1K165dy6NGSZKkChfSwxtvvfVWli1bxksvvcSyZcu44oorqF+/Pnv37mXt2rV8+eWXNGnShB/96EflXa8kSVKFCCkU1apViylTpjBp0iTee+89duzYEeyLjo7mmmuuYfTo0aW+D02SJCkShRSK4ORt94888gj79+9ny5Yt5ObmEh8fT+vWralXr1551ihJklThQg5FRb71rW/RoUOH8qhFkiQpbHz3mSRJEmUMRePHjycjIyPkg5zt9pIkSRWtTKFo2bJl/Otf/wr5IGe7vSRJUkUr85qizMxMFi9eXJG1SJIkhU2ZQ9GKFStYuXLlGR8gEAic8TaSJEmVrUyhaOzYsWd9oFatWp31PiRJkipKmUJR7969K7oOSZKksPKWfEmSJAxFkiRJgKFIkiQJMBRJkiQBhiJJkiTAUCRJkgQYiiRJkgBDkSRJEmAokiRJAgxFkiRJgKFIkiQJMBRJkiQBhiJJkiTAUCRJkgQYiiRJkgBDkSRJEmAokiRJAgxFkiRJgKFIkiQJMBRJkiQBhiJJkiTAUCRJkgQYiiRJkgBDkSRJEmAokiRJAgxFkiRJgKFIkiQJMBRJkiQBhiJJkiTAUCRJkgQYiiRJkgBDkSRJEgAx4S7gvx05coRZs2axadMmNm/eTE5ODg888AC9e/cuMfaLL75gypQprF+/npiYGDp06MDIkSOpW7dusXEnTpxg1qxZvPnmm+zfv5/k5GQGDhxIjx49KumsJElSpIu4UHTo0CFmzpxJUlISLVu2ZO3ataWO27NnD6NGjSI+Pp4hQ4aQl5fHrFmz2LZtG9OmTaN69erBsS+++CKvvPIKffv2pXXr1qxYsYLx48cTFRVF9+7dK+vUJElSBIu4UJSQkMAbb7xBQkICn332GUOHDi113Msvv8zRo0eZPn06SUlJAKSmpjJ69GgWLVpEeno6AHv37mX27NncfPPN3HfffQDccMMNjBo1iqlTp9K1a1eqVatWOScnSZIiVsStKYqNjSUhIeEbx2VkZNCxY8dgIAJo164djRs3Zvny5cG2FStWUFhYyM033xxsi4qK4qabbmLv3r1s3LixfE9AkiRVSRE3U1QWe/fu5cCBA6SkpJToS01NZdWqVcHPmZmZ1KxZkyZNmpQYV9Tftm3bUo+zb98+srOzg5+zsrLKo3xJkhSBqmQoKgoqpc0oJSQk8NVXX1FQUEBsbCzZ2dnUq1ePqKioEuPgZPA5lQULFjBz5szyK1ySJEWsKhmK8vPzAYotpi4SGxsbHBMbG0t+fv43jjuV9PR0OnXqFPyclZXFhAkTzqp2SZIUmapkKIqLiwPg2LFjJfoKCgqKjYmLiyvTuNIkJiaSmJh41vVKkqTIF3ELrcui6NLX19f7FMnOzuaCCy4IzgQlJCSwf/9+AoFAiXGAoUeSJAFVNBTVr1+funXrsmXLlhJ9mzdvpmXLlsHPLVu25OjRoyUWSW/atCnYL0mSVCVDEcA111zD+++/z+7du4Nta9asYfv27XTr1i3Y1rlzZ2JiYnjjjTeCbYFAgPnz51O/fn0uvfTSSq1bkiRFpohcU/T666+Tm5sbvMS1cuVK9uzZA8Att9xCfHw8AwcO5J133uHee++lX79+5OXl8dprr9G8efNirwS56KKL6N+/P6+99hqFhYWkpqby3nvv8emnnzJu3Dgf3ChJkoAIDUWzZ89m165dwc/vvvsu7777LgDXXXcd8fHxJCUlMXnyZKZMmcK0adOC7z4bMWJEcD1Rkbvvvps6deqwYMECFi9eTHJyMg899BA9e/as1POSJEmRKyJD0Zw5c8o0rlmzZkyaNOkbx0VHRzNw4EAGDhx4tqVJkqRzVJVdUyRJklSeDEWSJEkYiiRJkgBDkSRJEhChC60lSaXLO5DPqyP/Hu4yqqTqNWJo178Vza5qGO5SFKEMRZJUBVSvEQPkEwjAkf2nfpG1Tiefj+ZmGop0SoYiSaoC2vVvxUdzMzl2tDDcpVRJeQdOBkp/fjodQ5EkVQHNrmroDMdZeHXk351h0zdyobUkSRKGIkmSJMBQJEmSBBiKJEmSAEORJEkSYCiSJEkCDEWSJEmAoUiSJAkwFEmSJAGGIkmSJMBQJEmSBBiKJEmSAEORJEkSYCiSJEkCDEWSJEmAoUiSJAkwFEmSJAGGIkmSJMBQJEmSBBiKJEmSAEORJEkSYCiSJEkCDEWSJEmAoUiSJAkwFEmSJAGGIkmSJMBQJEmSBBiKJEmSAEORJEkSYCiSJEkCDEWSJEmAoUiSJAkwFEmSJAGGIkmSJMBQJEmSBBiKJEmSAEORJEkSYCiSJEkCDEWSJEkAxIS7AEmSKkvegXxeHfn3cJehU6heI4Z2/VvR7KqGYTm+oUiSdM6rXiMGyCcQgCP788Ndjk4pn4/mZhqKJEmqKO36t+KjuZkcO1oY7lJ0CnkHTobWcP4eGYokSee8Zlc1DNvsg8rm1ZF/D/ssngutJUmSMBRJkiQBhiJJkiTAUCRJkgQYiiRJkgBDkSRJEnAe3ZJfUFDAjBkzWLp0KTk5ObRo0YLBgwfTvn37cJcmSZIiwHkzU/TYY48xZ84cevbsyT333EN0dDT3338/n376abhLkyRJEeC8CEWbNm1i2bJlDB06lOHDh5Oens7TTz9NgwYNeP7558NdniRJigDnRSjKyMigWrVqpKenB9vi4uJIS0tj48aN7N69O4zVSZKkSHBerCnKzMwkOTmZ2rVrF2tPTU0FYOvWrSQlJZXYbt++fWRnZwc/Z2VlVWyhkiQpbM6LUJSdnU1CQkKJ9qK2ffv2lbrdggULmDlzZkWWJkmSIsR5EYry8/OpXr16ifbY2Nhgf2nS09Pp1KlT8HNWVhYTJkwo9/pqXRhX7L+SJJ1vIuHvwvMiFMXFxXHs2LES7QUFBcH+0iQmJpKYmFihtQHc9Ginbx4kSdI5LBL+LjwvFlonJCQUWxtUpKitMoKPJEmKbOdFKGrZsiX/+c9/OHz4cLH2TZs2BfslSdL57bwIRV27duX48eMsWLAg2FZQUMDChQtp06ZNqXeeSZKk88t5saaoTZs2dOvWjRdeeIGDBw/SqFEjFi9ezK5duxgzZky4y5MkSRHgvAhFAA8++CBJSUksWbKE3Nxcmjdvzh/+8Acuv/zycJcmSZIiwHkTiuLi4hg+fDjDhw8PdymSJCkCnRdriiRJkr6JoUiSJAlDkSRJEmAokiRJAgxFkiRJgKFIkiQJMBRJkiQB59FzispDfn4+AFlZWWGuRJIknakmTZpQo0aNU/Ybis7Arl27AJgwYUKYK5EkSWfqxRdfJCUl5ZT9UYFAIFCJ9VRpBw8eZPXq1TRs2JDY2NjTjn322WcZNWpUmfablZXFhAkTeOihh2jSpEl5lHreOZOfd6SJhNorq4aKOE557fNs9hPqtn5PVK5I+LMWikipuzLqqOjvCGeKylHdunW57rrryjQ2Pj7+tGm0NE2aNDnjbXRSKD/vSBEJtVdWDRVxnPLa59nsJ9Rt/Z6oXJHwZy0UkVJ3ZdQR7u8IF1pXkB49eoS7hPNKVf55R0LtlVVDRRynvPZ5NvsJddtI+L0/n1TVn3ek1F0ZdYT7O8LLZxFgy5YtDBky5BuvdUo6f/k9IVU8Z4oiQEJCAoMGDSIhISHcpUiKUH5PSBXPmSJJkiScKZIkSQIMRZIkSYC35FcZBQUFPPnkk3z00Ufk5ubStGlTRo4cyaWXXhru0iRFiCeeeIKVK1dy9OhRkpKSGDp0KJ06dQp3WVKV4ZqiKiIvL4/Zs2fTu3dv6tevz/Lly3n66aeZPXs2tWrVCnd5kiJAVlZW8OGymzdvZvTo0cyaNYsLL7ww3KVJVYKXz6qImjVrMmjQIJKSkoiOjqZ79+7ExMSwffv2cJcmKUI0adIk+LT9qKgojh07xr59+8JclVR1ePmsghw5coRZs2axadMmNm/eTE5ODg888AC9e/cuMbagoIAZM2awdOlScnJyaNGiBYMHD6Z9+/an3P/27dvJycmhUaNGFXkakipIRX1HPPnkkyxcuJCCggKuvvpqmjdvXhmnI50TnCmqIIcOHWLmzJlkZWXRsmXL04597LHHmDNnDj179uSee+4hOjqa+++/n08//bTU8fn5+UyYMIGf/OQnxMfHV0T5kipYRX1HjB49miVLlvDUU0/Rvn17oqKiKuoUpHOOoaiCJCQk8MYbbzB37lyGDRt2ynGbNm1i2bJlDB06lOHDh5Oens7TTz9NgwYNeP7550uMLyws5OGHH6ZRo0YMGjSoAs9AUkWqqO8IgGrVqvH973+fNWvW8MEHH1TUKUjnHENRBYmNjS3Tk2czMjKoVq0a6enpwba4uDjS0tLYuHEju3fvDrafOHGCCRMmEBUVxYMPPui/AKUqrCK+I/7b8ePH2bFjR7nUK50PDEVhlpmZSXJyMrVr1y7WnpqaCsDWrVuDbRMnTiQ7O5vf/va3xMS4HEw6H5T1OyI3N5e33nqLI0eOUFhYyPLly1m7di2XXXZZpdcsVVX+zRpm2dnZpf5rsait6M6RXbt28de//pXY2Nhi/2J8/PHH/dKTzmFl/Y6Iiorir3/9K0899RSBQIBGjRoxbtw4WrVqVan1SlWZoSjM8vPzqV69eon2ottq8/PzAWjQoAHvvvtupdYmKfzK+h1Ru3ZtnnnmmUqtTTrXePkszOLi4jh27FiJ9oKCgmC/pPOX3xFS5TEUhVlCQgLZ2dkl2ovaEhMTK7skSRHE7wip8hiKwqxly5b85z//4fDhw8XaN23aFOyXdP7yO0KqPIaiMOvatSvHjx9nwYIFwbaCggIWLlxImzZtSEpKCmN1ksLN7wip8rjQugK9/vrr5ObmBqe5V65cyZ49ewC45ZZbiI+Pp02bNnTr1o0XXniBgwcP0qhRIxYvXsyuXbsYM2ZMOMuXVMH8jpAiS1QgEAiEu4hz1YABA9i1a1epfbNnz6Zhw4bAybtHit5rlJubS/PmzRk8eDBXXnllZZYrqZL5HSFFFkORJEkSrimSJEkCDEWSJEmAoUiSJAkwFEmSJAGGIkmSJMBQJEmSBBiKJEmSAEORJEkSYCiSJEkCDEWSJEmAoUjSeapLly7FfuXn5wf7Fi1aRJcuXVi0aFEYK/w/8+fPL1br7373u3CXJJ2TYsJdgKRz286dO7n11ltPO6ZBgwbMmTOnkioqftxevXoBUK1atQo91urVq/nVr35F+/btmTRp0mnHjh8/nrfffptx48bRs2dPUlJSGDRoELm5ucybN69C65TOZ4YiSZWiUaNG9OzZs9S++Pj4Sq7mpAYNGnDXXXdVyrHatWtHUlISa9asYffu3SQlJZU6Ljc3l/fee4/4+Hi6dOkCQOvWrWndujU7d+40FEkVyFAkqVI0atSo0gJIJIqOjqZ3797MnDmTxYsXc+edd5Y67u233yY/P58+ffoQFxdXyVVK5zfXFEmKOF26dOGee+5h7969jB8/nr59+3L99ddz//338+WXXwLwxRdf8OCDD5KWlsb111/PuHHj2L9/f4XWtWfPHu6880569OjBO++8E2w/cOAAzz77LLfddhvdu3enb9++PPTQQ2zbtq3Y9n369CEqKopFixYRCARKPcbChQsBSEtLq7DzkFQ6Q5GkiJSTk8OIESPYuXMn119/PVdccQWrVq1i9OjRbNu2jeHDh5OXl0efPn1o3bo1GRkZ/Pa3v62wer744guGDx/Onj17eOKJJ+jatSsAO3bsYPDgwcydO5eLL76YH/7wh1x99dWsXr2aYcOGsWnTpuA+GjRowPe//32+/PJL1q5dW+IY27Zt47PPPqNVq1Z85zvfqbBzkVQ6L59JqhQ7duzgT3/6U6l9l1xyCVdddVWxtn/+858MGDCAkSNHBtuefPJJ3nzzTUaOHMlPf/pT+vfvD0AgEGDMmDGsWrWKLVu2kJKSUq61b9y4kTFjxhATE8Ozzz5Ly5Ytg32PPvoo+/fvZ+LEiVx55ZXB9jvuuIMhQ4bw+OOPM3PmzGB7WloaH330EQsXLuR73/teseM4SySFlzNFkirFjh07mDlzZqm/PvzwwxLja9asyeDBg4u1de/eHYALL7yQfv36BdujoqKCff/85z/Lte4PPviA++67jzp16jB16tRigejzzz9nw4YNXH/99cUCEUDjxo254YYb2LZtW7HLaD/4wQ+48MILycjI4PDhw8H2wsJCli5dSmxs7CkXpEuqWM4USaoUV155JRMnTizz+OTkZGrUqFGsLSEhAYDmzZsTFRVVat++ffvOstL/s3z5cv7xj3/QokULnnjiCerVq1esv+jS2IEDB0qdBfv3v/8d/G/z5s0BgqFn3rx5vP3229x4440ArFy5koMHD9KjRw/q1KlTbucgqewMRZIiUu3atUu0FT1L6HR9hYWF5VbDxo0bOX78OG3bti0RiAC++uor4ORs0gcffHDK/eTl5RX7nJaWxrx581i4cGEwFHnpTAo/Q5EkncLQoUNZsWIF8+bNo1q1aowYMaJYf1E4+8UvfsEtt9xS5v22aNGC1q1bs3nzZv71r39Rp04dVq9eTcOGDUusM5JUeVxTJEmnEBsby6OPPkqHDh2YPXs2U6ZMKdafmpoKnJxROlNFM0J/+9vfWLJkCcePHw/esi8pPAxFknQasbGxTJgwgY4dOzJnzhyeffbZYF+bNm1o06YNy5YtY9myZSW2PXHiBOvWrSt1vz169KBGjRosXbqUhQsXEh0dHXzliKTw8PKZpEpxulvyAX7yk59E7BOcq1evziOPPMLDDz/M3LlzCQQC3HPPPQA8/PDD3Hvvvfz2t79l3rx5tGrViri4OPbs2cOGDRs4dOgQb7/9dol91q5dm2uuuYYlS5Zw8OBBrrrqqlO++kNS5TAUSaoURbfkn0r//v0jNhTB/wWjX//618ybN49AIMAvfvELLr74YmbMmMHs2bN57733WLRoEdHR0SQkJHDZZZcFH/JYmrS0NJYsWQKcfNq1pPCKCpzqWfOSdA7r0qULl19+OZMnTw53KWW2c+dObr31Vnr16sWDDz4Y7nKkc44zRZLOW+vWrQu+if6tt96K2Jmq+fPnM2nSpHCXIZ3zDEWSzkuDBg0q9rnoOUeRKCUlpVi9rVq1Cl8x0jnMy2eSJEl4S74kSRJgKJIkSQIMRZIkSYChSJIkCTAUSZIkAYYiSZIkwFAkSZIEGIokSZIAQ5EkSRIA/x/MbYXiqWGwxwAAAABJRU5ErkJggg==",
+ "text/plain": [
+ "
"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "bkg.project('Em').plot()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 24,
+ "id": "eee3d467-d795-43b3-857e-b6cc1c09c1f7",
+ "metadata": {},
+ "outputs": [
+ {
+ "ename": "TypeError",
+ "evalue": "matplotlib.transforms.Bbox.from_bounds() argument after * must be an iterable, not int",
+ "output_type": "error",
+ "traceback": [
+ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+ "\u001b[0;31mTypeError\u001b[0m Traceback (most recent call last)",
+ "Cell \u001b[0;32mIn[24], line 5\u001b[0m\n\u001b[1;32m 1\u001b[0m m_bkg \u001b[38;5;241m=\u001b[39m HealpixMap(bkg\u001b[38;5;241m.\u001b[39mslice[{\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mEm\u001b[39m\u001b[38;5;124m'\u001b[39m:\u001b[38;5;241m3\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mPhi\u001b[39m\u001b[38;5;124m'\u001b[39m:\u001b[38;5;241m4\u001b[39m}]\u001b[38;5;241m.\u001b[39mproject(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mPsiChi\u001b[39m\u001b[38;5;124m'\u001b[39m)\u001b[38;5;241m.\u001b[39mtodense()\u001b[38;5;241m.\u001b[39mcontents)\n\u001b[1;32m 3\u001b[0m fig \u001b[38;5;241m=\u001b[39m plt\u001b[38;5;241m.\u001b[39mfigure(dpi \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m150\u001b[39m)\n\u001b[0;32m----> 5\u001b[0m ax \u001b[38;5;241m=\u001b[39m \u001b[43mfig\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43madd_subplot\u001b[49m\u001b[43m(\u001b[49m\u001b[43mprojection\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43morthview\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m)\u001b[49m\n\u001b[1;32m 7\u001b[0m m_bkg\u001b[38;5;241m.\u001b[39mplot(ax)\n",
+ "File \u001b[0;32m~/.virtualenvs/cosipy/lib/python3.9/site-packages/matplotlib/figure.py:757\u001b[0m, in \u001b[0;36mFigureBase.add_subplot\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m 754\u001b[0m args \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mtuple\u001b[39m(\u001b[38;5;28mmap\u001b[39m(\u001b[38;5;28mint\u001b[39m, \u001b[38;5;28mstr\u001b[39m(args[\u001b[38;5;241m0\u001b[39m])))\n\u001b[1;32m 755\u001b[0m projection_class, pkw \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_process_projection_requirements(\n\u001b[1;32m 756\u001b[0m \u001b[38;5;241m*\u001b[39margs, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs)\n\u001b[0;32m--> 757\u001b[0m ax \u001b[38;5;241m=\u001b[39m \u001b[43mprojection_class\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mpkw\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 758\u001b[0m key \u001b[38;5;241m=\u001b[39m (projection_class, pkw)\n\u001b[1;32m 759\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_add_axes_internal(ax, key)\n",
+ "File \u001b[0;32m~/.virtualenvs/cosipy/lib/python3.9/site-packages/mhealpy/plot/axes.py:259\u001b[0m, in \u001b[0;36mOrthview.__init__\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m 257\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m__init__\u001b[39m(\u001b[38;5;28mself\u001b[39m, \u001b[38;5;241m*\u001b[39margs, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs):\n\u001b[0;32m--> 259\u001b[0m \u001b[38;5;28;43msuper\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[38;5;21;43m__init__\u001b[39;49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 260\u001b[0m \u001b[43m \u001b[49m\u001b[43mframe_class\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m \u001b[49m\u001b[43mkwargs\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mpop\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mframe_class\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mEllipticalFrame\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 261\u001b[0m \u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n",
+ "File \u001b[0;32m~/.virtualenvs/cosipy/lib/python3.9/site-packages/mhealpy/plot/axes.py:41\u001b[0m, in \u001b[0;36mHealpyAxes.__init__\u001b[0;34m(self, fig, rect, coord, flip, rot, **kwargs)\u001b[0m\n\u001b[1;32m 24\u001b[0m \u001b[38;5;250m\u001b[39m\u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m 25\u001b[0m \u001b[38;5;124;03mBase class for WCSAxes that behave similar to healpy's projections.\u001b[39;00m\n\u001b[1;32m 26\u001b[0m \n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 37\u001b[0m \u001b[38;5;124;03m of angle psi around this direction is applied.\u001b[39;00m\n\u001b[1;32m 38\u001b[0m \u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m 40\u001b[0m \u001b[38;5;66;03m# Get equivalent WCS FITS header\u001b[39;00m\n\u001b[0;32m---> 41\u001b[0m naxis1,naxis2 \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_get_naxis\u001b[49m\u001b[43m(\u001b[49m\u001b[43mfig\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mrect\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 43\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_autoscale:\n\u001b[1;32m 44\u001b[0m \u001b[38;5;66;03m# Gnom specifies the reosolution and the limits are adjusted based\u001b[39;00m\n\u001b[1;32m 45\u001b[0m \u001b[38;5;66;03m# on the figure sizes\u001b[39;00m\n\u001b[1;32m 46\u001b[0m \u001b[38;5;66;03m# For the others the limits are fixed and the resolution is adjusted\u001b[39;00m\n\u001b[1;32m 47\u001b[0m \u001b[38;5;66;03m# by scaling a normalized appropiate resolution\u001b[39;00m\n\u001b[1;32m 48\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_cdelt \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m=\u001b[39m \u001b[38;5;241m360\u001b[39m\u001b[38;5;241m/\u001b[39mnaxis1\n",
+ "File \u001b[0;32m~/.virtualenvs/cosipy/lib/python3.9/site-packages/mhealpy/plot/axes.py:140\u001b[0m, in \u001b[0;36mHealpyAxes._get_naxis\u001b[0;34m(self, fig, rect)\u001b[0m\n\u001b[1;32m 133\u001b[0m \u001b[38;5;250m\u001b[39m\u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m 134\u001b[0m \u001b[38;5;124;03mReturn approapiate NAXIS1 and NAXIS2 for a given figure and Bbox (or bounds)\u001b[39;00m\n\u001b[1;32m 135\u001b[0m \u001b[38;5;124;03m \u001b[39;00m\n\u001b[1;32m 136\u001b[0m \u001b[38;5;124;03maspect = naxis2/naxis1\u001b[39;00m\n\u001b[1;32m 137\u001b[0m \u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m 139\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(rect, Bbox):\n\u001b[0;32m--> 140\u001b[0m rect \u001b[38;5;241m=\u001b[39m \u001b[43mBbox\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mfrom_bounds\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mrect\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 142\u001b[0m naxis1 \u001b[38;5;241m=\u001b[39m fig\u001b[38;5;241m.\u001b[39mget_figwidth() \u001b[38;5;241m*\u001b[39m fig\u001b[38;5;241m.\u001b[39mdpi \u001b[38;5;241m*\u001b[39m rect\u001b[38;5;241m.\u001b[39mwidth\n\u001b[1;32m 143\u001b[0m naxis2 \u001b[38;5;241m=\u001b[39m fig\u001b[38;5;241m.\u001b[39mget_figheight() \u001b[38;5;241m*\u001b[39m fig\u001b[38;5;241m.\u001b[39mdpi \u001b[38;5;241m*\u001b[39m rect\u001b[38;5;241m.\u001b[39mheight\n",
+ "\u001b[0;31mTypeError\u001b[0m: matplotlib.transforms.Bbox.from_bounds() argument after * must be an iterable, not int"
+ ]
+ },
+ {
+ "data": {
+ "text/plain": [
+ "
"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "m_bkg = HealpixMap(bkg.slice[{'Em':3, 'Phi':4}].project('PsiChi').todense().contents)\n",
+ "\n",
+ "fig = plt.figure(dpi = 150)\n",
+ "\n",
+ "ax = fig.add_subplot(projection = 'orthview')\n",
+ "\n",
+ "m_bkg.plot(ax)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "3bcca3df-33c4-4843-aac1-27facf67cafe",
+ "metadata": {},
+ "source": [
+ "Note: I actually don't understand what causes the strip in the middle. Maybe it's a beating pattern caused by converting from FISBEL to HEALPix during the creation of the detector response. I plan to generate a detector response using HEALPix directly, and will revisit this then."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "61f927b2-1132-4afe-971a-cc649a832ebb",
+ "metadata": {},
+ "source": [
+ "## Putting it togetter"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "2b0ce845-5b17-48e9-8282-97e6d703bcbc",
+ "metadata": {},
+ "source": [
+ "Once we obtain the expected signal, it's easy to add it do the background to simulate how the observed data would look like"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "7d4f3132",
+ "metadata": {},
+ "source": [
+ "Here we draw a Poisson sample to get the actual counts for both sky (signal) and background"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 25,
+ "id": "2b549107",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def Poisson_Sample_Sparse_Histogram(signal):\n",
+ " from histpy import Histogram\n",
+ " from sparse import COO\n",
+ " pdata = np.random.poisson(signal.todense()[:].value)\n",
+ " psignal = Histogram(signal.axes,\n",
+ " COO(np.nonzero(pdata),\n",
+ " data=pdata[np.nonzero(pdata)],\n",
+ " shape=pdata.shape),\n",
+ " sparse=True)\n",
+ " return psignal"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 26,
+ "id": "d6e1970f",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "psignal = Poisson_Sample_Sparse_Histogram(signal)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 27,
+ "id": "27810ec2",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "