{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Preparing an extended source for LeakageLib fitting\n",
"\n",
"The steps to process extended source data before fitting are the same as for point sources, with the extra step of extracting spatial weights. This should be done with the `leakagelib_cxo` script, which generates spatial weights from a Chandra image. Note that `leakagelib_cxo` is NOT a replica of IXPEobssim; IXPEobssim incorporates IXPE's PSF and generates a simulated event file. `leakagelib_cxo` does not (and should not).\n",
"\n",
"This extended source section analyzes Crab OBS ID 02006001. Please download the data from Heasarc. (We will skip particle weights, so you do not need to download the level 1 data). After unzipping the files, the next step is to run `leakagelib_cxo`. Run the help script to see what the script requires."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\u001b[93m>>> PyXSPEC is not installed, you will no be able to use it.\u001b[0m\n",
"usage: clean [-h] --cxo-evt CXO_EVT [CXO_EVT ...] --cxo-arf CXO_ARF\n",
" [CXO_ARF ...] --ixpe-evt IXPE_EVT --ixpe-arf IXPE_ARF\n",
" [--expmap EXPMAP] --output OUTPUT [--width WIDTH] [--elow ELOW]\n",
" [--ehigh EHIGH] [--centerx CENTERX] [--centery CENTERY]\n",
" [--reg-src REG_SRC] [--reg-bkg REG_BKG]\n",
" [--clobber | --no-clobber]\n",
"\n",
"Creates a CXO image, adjusted to the IXPE band. Before running this command,\n",
"you should use the ciao tools merge_obs and mkwarf to create a merged image\n",
"(including an exposure map) and an arf for each observaton. These are required\n",
"as inputs for this script.\n",
"\n",
"options:\n",
" -h, --help show this help message and exit\n",
" --cxo-evt CXO_EVT [CXO_EVT ...]\n",
" List of Chandra event files to use\n",
" --cxo-arf CXO_ARF [CXO_ARF ...]\n",
" List of Chandra ARFs to use\n",
" --ixpe-evt IXPE_EVT IXPE event file (only one is necessary. Try DU1)\n",
" --ixpe-arf IXPE_ARF IXPE ARF (only one is necessary, for simplicity. Try\n",
" DU1)\n",
" --expmap EXPMAP Chandra merged observation ARF (Optional)\n",
" --output OUTPUT Name of the output file\n",
" --width WIDTH Width of the image, in arcseconds. Default: as big as\n",
" the CXO image (Optional)\n",
" --elow ELOW Low end of the energy range (keV). (Default: 2)\n",
" --ehigh EHIGH High end of the energy range (keV). (Default: 8)\n",
" --centerx CENTERX IXPE pixel on which the fit will be centered in the x\n",
" direction. (Default: 300)\n",
" --centery CENTERY IXPE pixel on which the fit will be centered in the y\n",
" direction. (Default: 300)\n",
" --reg-src REG_SRC Source region. Should be saved in CIAO format in fk5\n",
" coordinates (Optional)\n",
" --reg-bkg REG_BKG Background region. Should be saved in CIAO format in\n",
" fk5 coordinates (Optional)\n",
" --clobber, --no-clobber\n",
" Overwrite files\n"
]
}
],
"source": [
"!python -m leakagelib_cxo -h"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As you see, it requires\n",
"* CXO event file(s) to create the Source object from\n",
"* IXPE ARF and CXO ARF(s) to find the expected fluxes, taking the different effective areas into account. Provide multiple ARFs if providing multiple CXO event files\n",
"* IXPE event file to orient the CXO image relative to the IXPE pointing.\n",
"* IXPE energy range used for fitting (2-8 keV by default)\n",
"* IXPE pixel used for centering (300,300 by default)\n",
"\n",
"and optionally\n",
"* The width of the output spatial weight image\n",
"* A CXO source region to extract an image for\n",
"* A CXO background region to estimate background flux\n",
"* CXO Exposure map to exposure correct the event file\n",
"\n",
"`ciao` provides tutorials on how to extract ARFs and exposure maps for Chandra data, so we do not cover that here. If you are using more than one event file, I suggest using the full exposure map generated by the ciao tool merge_obs."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To make this example, I downloaded a short Crab observation from the Chandra archive (obsid 13146) and placed inside a `src.reg` region file containing\n",
"```\n",
"circle(5:34:31.8424,+22:00:52.502,0.5')\n",
"```\n",
"(This region was made in DS9 and saved with the ciao format with fk5 coordinates)\n",
"\n",
"I then ran the following script\n",
"```\n",
"chandra_repro . repro\n",
"asphist \"repro/pcadf13146_000N001_asol1.fits\" asp.hist evtfile=\"repro/acisf13146_repro_evt2.fits[sky=region(src.reg)]\" clobber=yes\n",
"sky2tdet \"repro/acisf13146_repro_evt2.fits[sky=region(src.reg)]\" asp.hist tdet[wmap] clobber=yes\n",
"mkwarf tdet arf.arf weightfile=none feffile=none spectrum=none egrid=0.3:10:0.1 pbkfile=none clobber=yes\n",
"```\n",
"\n",
"I also downloaded the IXPE data for the Crab from OBSID 02006001 and ran this script to generate an IXPE ARF:\n",
"```\n",
"ixpecalcarf evtfile=\"event_l2/ixpe02006001_det1_evt2_v01.fits\" attfile=\"hk/ixpe02006001_det1_att_v01.fits\" arfout=\"arf.arf\" specfile=None clobber=yes\n",
"```\n",
"Then the following code will create a source image. To keep things simple I did not apply background subtraction or exposure correction."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\u001b[93m>>> PyXSPEC is not installed, you will no be able to use it.\u001b[0m\n",
"Exposure map corrections will not be made because you did not pass an exposure map\n",
"You did not provide a background region. The Chandra image will not be background subtracted\n",
"Image saved to data/13146/crab.fits. To use it, load the source with LeakageLib using the following code:\n",
" \n",
"import leakagelib\n",
"source = leakagelib.Source.load_file(\"data/13146/crab.fits\")\n",
"\n",
"This gives a leakagelib.Source object, which you can use as a source in your LeakageLib fit.\n"
]
}
],
"source": [
"!python -m leakagelib_cxo\\\n",
" --output \"data/13146/crab.fits\"\\\n",
" --cxo-evt \"data/13146/repro/acisf13146_repro_evt2.fits\"\\\n",
" --cxo-arf \"data/13146/arf.arf\"\\\n",
" --ixpe-evt \"data/02006001/event_l2/ixpe02006001_det1_evt2_v01.fits\"\\\n",
" --ixpe-arf \"data/02006001/arf.arf\"\\\n",
" --reg-src \"data/13146/src.reg\"\\\n",
" --width 150\\\n",
" --clobber"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"
\n",
"\n",
"**Usage Notes**\n",
"\n",
"- This example doesn't subtract background because Crab is so bright, but most analyses should use the reg-bkg argument. Otherwise LeakageLib will include the CXO background as part of the source flux and try to fit for its polarization, which is incorrect.\n",
"- Since all LeakageLib fit sources need to be the same size, your \"width\" argument should be large enough to contain the entire data set, not just this one component.\n",
"- In order to make sure all sizes have the same width and center, make sure the width, centerx, and centery arguments are the same for all sources.\n",
"- If you later find that the spatial weight image is misaligned with the data, change centerx and centery to shift the weight image.\n",
"
\n",
"\n",
"A source object for LeakageLib is then loaded simply using"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"\u001b[93m>>> PyXSPEC is not installed, you will no be able to use it.\u001b[0m\n"
]
}
],
"source": [
"import leakagelib\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"source = leakagelib.Source.load_file(\"data/13146/crab.fits\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Displaying it,"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAj4AAAItCAYAAAA5cHCZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA5CUlEQVR4nO3de3RU9b3+8ScmJpNAAkaBJOWiUTEKKALFA4sjFTTFekMRELVUEThaC0gPlYqgRwT0LJSKxYXghfag7VFRutQDXqAHUbEqFX7FgpAKwRCCNwwQyIUk398fLqaOE/CzPZMZmO/7tRZ/MPNkz3f2XHjYM/nsFOecEwAAgAeOS/QCAAAA4oXiAwAAvEHxAQAA3qD4AAAAb1B8AACANyg+AADAGxQfAADgjbREL6C5tG7dWrW1tcrPz0/0UgAAQDOrqKhQRkaGKisrj5hLSdYBhpmZmaqpqVFKohcCAACanZMUCoVUXV19xFzSHvHJz89X6bZtykr0QgAAQLM7IJk+5eE7PgAAwBsUHwAA4A2KDwAA8AbFBwAAeIPiAwAAvEHxAQAA3qD4AAAAb1B8AACANyg+AADAGxQfAADgDYoPAADwBsUHAAB4g+IDAAC8QfEBAADeoPgAAABvUHwAAIA3KD4AAMAbFB8AAOANig8AAPAGxQcAAHiD4gMAALxB8QEAAN6g+AAAAG9QfAAAgDcoPgAAwBsUHwAA4A2KDwAA8AbFBwAAeIPiAwAAvEHxAQAA3qD4AAAAb1B8AACANyg+AADAGxQfAADgDYoPAADwBsUHAAB4g+IDAAC8QfEBAADeoPgAAABvUHwAAIA3KD4AAMAbFB8AAOCNuBafDz74QJdffrlyc3OVlZWlrl276uGHH47IrFmzRv369VNWVpby8vI0fvx4VVVVxXOZAAAgSaXF64Zee+01XXbZZTr33HM1bdo0tWzZUh9//LF27NgRzqxfv14DBw7UmWeeqTlz5mjHjh164IEHVFJSouXLl8drqQAAIEnFpfjs3btXI0eO1CWXXKIlS5bouOOaPtA0ZcoUnXDCCVq1apVycnIkSSeffLLGjBmj1157TcXFxfFYLgAASFJx+ajrD3/4gz799FPNnDlTxx13nPbv36/GxsaIzN69e/X666/r+uuvD5ceSRo5cqRatmypZ599Nh5LBQAASSwuxWfFihXKyclReXm5zjjjDLVs2VI5OTm65ZZbVFNTI0nasGGD6uvr1atXr4ifTU9PV/fu3bVu3boj3kZhYWHEn7Kysma7PwAA4NgUl+JTUlKi+vp6XXHFFfrxj3+s559/XqNGjdKjjz6qG2+8UZJUUVEhScrPz4/6+fz8fO3cuTMeSwUAAEksLt/xqaqq0oEDB3TzzTeHf4vrqquuUl1dnRYsWKDp06erurpakpSRkRH186FQKHz94WzdujXi74WFhSrdti1G9wAAACSDuBzxyczMlCSNGDEi4vJrr71WkvTOO++EM7W1tVE/X1NTE74eAADg+4pL8SkoKJAktWvXLuLytm3bSpK++uqr8Edchz7y+qaKiorwNgAAAL6vuBSfnj17SpLKy8sjLj/0vZ02bdqoa9euSktL09q1ayMydXV1Wr9+vbp37x6PpQIAgCQWl+IzbNgwSdITTzwRcfnjjz+utLQ0/ehHP1KrVq104YUX6qmnntK+ffvCmcWLF6uqqkpDhw6Nx1IBAEASi8uXm88991yNGjVKTz75pOrr69W/f3+tWrVKzz33nO64447wx1gzZ85U37591b9/f40dO1Y7duzQgw8+qOLiYg0aNCgeSwUAAEksxTnn4nFDBw8e1KxZs7Ro0SLt3LlTnTp10q233qrbbrstIvfWW29p8uTJ+uCDD5Sdna1hw4bpvvvuU3Z2dqDbO/RbXVkxvA9Asoj1od7UGG8vkRpivL3G744AiIEDkk4+5ZSo3/L+trgVn3ij+ACHR/E5PIoPcGyyFp+4np0dAAAgkSg+AADAGxQfAADgDYoPAADwBsUHAAB4g+IDAAC8QfEBAADeoPgAAABvUHwAAIA34nKuLgDRrP/rON6YSw9w2yFjLiPGt53I/2lZJyjXGXO1Md6eJB2Mcc46hZrp0vAJR3wAAIA3KD4AAMAbFB8AAOANig8AAPAGxQcAAHiD4gMAALxB8QEAAN6g+AAAAG9QfAAAgDeY3AwYWScoZxlzrROUk6RWxly2MWedBJ1qzAVhnWJsnaB8wJjbF+NckOxeY856X2pinJOYBo2jF0d8AACANyg+AADAGxQfAADgDYoPAADwBsUHAAB4g+IDAAC8QfEBAADeoPgAAABvUHwAAIA3mNyMpGRt9NbJxJKUZ8ydbsydZcx1MObaGHOSfbq0dVq1dUpvgzEXZEKwNWudYrzfmNtjzFUac5J9crN1m9bc58bcbmMuyG1bHz8mQSNWOOIDAAC8QfEBAADeoPgAAABvUHwAAIA3KD4AAMAbFB8AAOANig8AAPAGxQcAAHiD4gMAALzB5GYcU6yThHONOev0ZEn6V2Ou2Jjrbh2ffLIxFzLmJPsIZeu4Y+vo37222P4AY3pjPe041hOZrTkpcVOorROZy4w5SdpqzO005iqNuYPGHPzFER8AAOANig8AAPAGxQcAAHiD4gMAALxB8QEAAN6g+AAAAG9QfAAAgDcoPgAAwBsUHwAA4A0mN+OoYJ3IXGDM/dCYu9yYk6SrrIu8wpg7z5izjqEOMrLWOiI41iN9jdtrYZzwLEktjFOj8z615RqM27MOq6405iT7w2J9qK056/Rr632WpBJjboMxt9mY22XM1RpzSD4c8QEAAN6g+AAAAG9QfAAAgDcoPgAAwBsUHwAA4I2EFZ+ZM2cqJSVFXbt2jbpuzZo16tevn7KyspSXl6fx48erqqoqAasEAADJJCG/zr5jxw7NmjVLLVq0iLpu/fr1GjhwoM4880zNmTNHO3bs0AMPPKCSkhItX748AasFAADJIiHFZ9KkSfqXf/kXNTQ06Isvvoi4bsqUKTrhhBO0atUq5eTkSJJOPvlkjRkzRq+99pqKi4sTsWQAAJAE4v5R1+rVq7VkyRI99NBDUdft3btXr7/+uq6//vpw6ZGkkSNHqmXLlnr22WfjuFIAAJBs4nrEp6GhQePGjdPo0aPVrVu3qOs3bNig+vp69erVK+Ly9PR0de/eXevWrTvstgsLCyP+XlZmHSULAAB8Edfi8+ijj2r79u1asWJFk9dXVFRIkvLz86Ouy8/P15tvvtms60PsWc/y0MaYs56KYoQx9+MOxqAkXWnM9TfmTjfmQsZcgzEXxH5jznpqC+u5ESqNOcl+uoyttljqdlsub6cx96UtJ8l+zooDtpj14bM+LNacJBV+d0SS/TQ0Od8dkSStN+aMD58k8+7GMSJuxefLL7/UXXfdpWnTpqlNm6b/mauurpYkZWRkRF0XCoXC1zdl69bId7XCwkKVbtv2f1gxAABINnH7js/UqVOVm5urcePGHTaTmZkpSaqtjT59XE1NTfh6AACA7yMuR3xKSkq0cOFCPfTQQ9q5858HGGtqanTw4EGVlpYqJycn/BHXoY+8vqmiokIFBdaDogAAANHicsSnvLxcjY2NGj9+vE455ZTwn3fffVdbtmzRKaecounTp6tr165KS0vT2rVrI36+rq5O69evV/fu3eOxXAAAkKTicsSna9euWrp0adTlU6dO1b59+zR37lydeuqpatWqlS688EI99dRTmjZtmrKzsyVJixcvVlVVlYYOHRqP5QIAgCQVl+Jz0kknafDgwVGXH5rl883rZs6cqb59+6p///4aO3asduzYoQcffFDFxcUaNGhQPJYLAACS1FF3ktIePXpoxYoVyszM1MSJE7Vw4ULddNNNWrJkSaKXBgAAjnEJOWXFIatWrWry8n79+untt9+O72IAAEDSO+qO+AAAADSXhB7xwbEpSFtubcydZcxZT1H742xr0JiTpPOMOevI2lxjznpfWhhzkpQaPSS0ada3iHpbrC56RleTrJOgJanUmCsx5v5hzFknRge5L3uNOeMI5RbWXKUtlxdgdHMH47jjPOP2Whtz1mnx7xpzklRuzBmf3UgwjvgAAABvUHwAAIA3KD4AAMAbFB8AAOANig8AAPAGxQcAAHiD4gMAALxB8QEAAN6g+AAAAG8wuRmBBRkQ3N6Ysw5FHmC94VhPWZbsd7zGmDNOtjVPbk5tZQxK0g+MOetc3Za2WHqV8Wb/brxdSVmf2nIh4/asuSxj7nNjTrI/J/Ybc9ZJy9ZcgPsSMmYLK225LOO+abDFzLtQsr+krbvnYIDbRuxxxAcAAHiD4gMAALxB8QEAAN6g+AAAAG9QfAAAgDcoPgAAwBsUHwAA4A2KDwAA8AbFBwAAeIPJzQg73pgLMh/4ZGPuDGMuzzrFuMCYs07fleyjXo2DhNVozFnXmL3HGJSUcqIxeJIxZ30ErTNwreOTJeUYpzyfsd2Ws74QUo25IKPOrbvHmqs05qwjh637RrI/hMbXdJ7xdfUvlbZckIHau425WA/Utr5FIBiO+AAAAG9QfAAAgDcoPgAAwBsUHwAA4A2KDwAA8AbFBwAAeIPiAwAAvEHxAQAA3qD4AAAAbzC5GWHWQavWmb+SlBfrbVqnGFun6h603rDsU3Ct27ROwQ0yXdqqYKstl3KycYNdjLn2xlyQt6aWtlimcYRyl822XIsGW876IpDsE5mtI4J3GnPW6dLpxlyQbQYY0m1RaNw3Pwzw2jfO/NYuY876MNcacwiGIz4AAMAbFB8AAOANig8AAPAGxQcAAHiD4gMAALxB8QEAAN6g+AAAAG9QfAAAgDcoPgAAwBsUHwAA4A1OWeEBa7u1nkEhyIR56zaNw//tM9ytM+H3WW9YARZpZD0VxW5jLsjpN6z3peP/MwZbG3NnGHPG01BIkk415qxvdyfZYifvMOY+Nd6upAbj+Ra+NG7PeGYS83PR+oKWpAxjzvpctL6m29hihdbTecj+rC0x5mL9km405vA1jvgAAABvUHwAAIA3KD4AAMAbFB8AAOANig8AAPAGxQcAAHiD4gMAALxB8QEAAN6g+AAAAG8wudkDqTHOBWEdtmqdZFpTacuFrBsMMrk5yGRki2xjzrpG686W7PclZBwR3PZd4wati+xizElSe2POOJFZnYy5L4y5AJObUyttubYfGzdYbosdMG4uyOul0piL9RuUcbp0C+PmJOkHxlwHY267MWfd3UxuDoYjPgAAwBtxKT7vv/++fvGLX6hLly5q0aKFOnbsqGHDhmnLli1R2U2bNmnQoEFq2bKlcnNz9dOf/lSff/55PJYJAACSXFw+6vrP//xPvf322xo6dKjOPvts7dq1S/PmzVOPHj30l7/8RV27dpUk7dixQ+eff75atWqlWbNmqaqqSg888IA2bNig9957T+np6fFYLgAASFJxKT6//OUv9Yc//CGiuAwfPlzdunXT/fffr6eeekqSNGvWLO3fv19//etf1bFjR0lS7969ddFFF+l3v/udxo4dG4/lAgCAJBWXj7r69u0bdbTm9NNPV5cuXbRp06bwZc8//7wuvfTScOmRpAsvvFCdO3fWs88+G4+lAgCAJJawLzc75/Tpp5/qpJO+/k2L8vJyffbZZ+rVq1dUtnfv3lq3bl28lwgAAJJMwn6d/emnn1Z5ebmmT58uSaqoqJAk5efnR2Xz8/O1e/du1dbWKiMjo8ntFRYWRvy9rKwsxisGAADHuoQc8fnoo4906623qk+fPvrZz34mSaqurpakJotNKBSKyAAAAHwfcT/is2vXLl1yySVq1aqVlixZotTUr6dSZWZmSpJqa2ujfqampiYi05StW7dG/L2wsFCl27bFatkAACAJxLX47NmzRxdffLEqKyv15ptvqqCgIHzdoY+4Dn3k9U0VFRXKzc097MdciA3r4b+6ANusNOasH0yWGnNF1g3mGnOSfdSrdQeFAtx2rFkHKGcZc9nGCcGZ1p1onbIs2ac8Fxlz1p1jnZ5casxJ9pm+lbZYyPi4WB/nIBNFrOOErbvb+royTiUP8nGH9W3COrnZuj3rBLtYD5VPdnErPjU1Nbrsssu0ZcsWrVixQmeddVbE9T/4wQ/Upk0brV27Nupn33vvPXXv3j1OKwUAAMkqLt/xaWho0PDhw/XOO+/oueeeU58+fZrMDRkyRC+//HLEF5NXrlypLVu2aOjQofFYKgAASGJxOeLz7//+73rxxRd12WWXaffu3eGBhYdcf/31kqQpU6boueee0wUXXKAJEyaoqqpKs2fPVrdu3XTjjTfGY6kAACCJxaX4rF+/XpL00ksv6aWXXoq6/lDx6dChg9544w398pe/1K9//Wulp6frkksu0YMPPsj3ewAAwP9ZXIrPqlWrzNkuXbro1Vdfbb7FAAAAbyVscjMAAEC8UXwAAIA3KD4AAMAbFB8AAOCNhJ2kFPET63YbfVKRw9ttzO0y5rZ+d0SSVPSRMZhtzElSG2POOonWOtn2gDEX5L5Yb7vSmNtnzGV+agx+YcxJUoMxd6ox1zrGuXpjTpJ2GHNVtpj1uWN9/KzPbSnYiHeL1NjmggxOt760rBOZrTnrGvcbc5L9pZ/MOOIDAAC8QfEBAADeoPgAAABvUHwAAIA3KD4AAMAbFB8AAOANig8AAPAGxQcAAHiD4gMAALzB5OZjmLW1WgeeWid6Bhneah0I+7kxZ53cvMF4Z7ptMG5Qkk435qx3utKY+9KYyzPmJKmDMXfQmDM/KYwTh82zvCWp1JizTo0+yZhrZ8xlGHOStMeYM96XSuPmrM9Z6/NBktKNOetYZOtzLMsWSw3w3/7jje8n1knLOcac9ZljfY+XmNwsccQHAAB4hOIDAAC8QfEBAADeoPgAAABvUHwAAIA3KD4AAMAbFB8AAOANig8AAPAGxQcAAHiD4gMAALzBKSsQWF2A7H5jrtKY227M/T9jLnTAGJR0+mZjsMCYO96Ys470D3BfzLP1rae2sJ7y4LMGW67tx8YNSpL1gVlvzLUw5mqNOeurQDKfqqPauE3j7rae5sF8Ng/J/kZhzVmf38bXVUOAczdYz5ZhPXWE9SVtfZkiGI74AAAAb1B8AACANyg+AADAGxQfAADgDYoPAADwBsUHAAB4g+IDAAC8QfEBAADeoPgAAABvMLn5GGadEmrNWQUYeGqeeLrHmDPOtVUrY66NMSdJrY13ps1O4watk5utuYPGnCS1NuY6GXMZxtxuY878jJDUdocx+KkxZ50PXm/MVRlzkvktOd24OesQ6h/EeHuSfY3W5635uRPbmw2StR5JsA7Ktu5u61uEFOx+JyuO+AAAAG9QfAAAgDcoPgAAwBsUHwAA4A2KDwAA8AbFBwAAeIPiAwAAvEHxAQAA3qD4AAAAbzC5GYE1BMhap4TuN+YqjTnrjF7rkGVJKjDmQsYJz8dbx1rH+HYlSbnGnHVys3V0rHX6bsiYk6S25cagdYLyF8acdXJzEO1tsVTjqOV2xn1jfT5Yc5J9ZPznxpz1OWa83SDTjq1Z66TlWE9kjvV0/mTHER8AAOANig8AAPAGxQcAAHiD4gMAALxB8QEAAN6g+AAAAG9QfAAAgDcoPgAAwBsUHwAA4I2jcnJzbW2t7rrrLi1evFhfffWVzj77bM2YMUMXXXRRopeGgKxTnhuNOetw4kpjbpcxJ0llxpx16HC2MWedynrcAWNQUrp1tPU+Y846itY6ytu6syXpjO22XM5m4wbzjDnrnW5tzEn2UdldbLFM41t8pvGVlVNpy0lSTa0tZx1PbH3uGKUG+G9/yPgGFesJykxkbh5H5RGfG264QXPmzNF1112nuXPnKjU1VT/5yU/01ltvJXppAADgGHbUHfF577339N///d+aPXu2Jk2aJEkaOXKkunbtqttvv11r1qxJ8AoBAMCx6qg74rNkyRKlpqZq7Nix4ctCoZBuuukmvfPOOyorC3IMHAAA4J+OuiM+69atU+fOnZWTkxNxee/evSVJ69evV4cOHaJ+rrCwMOLvFCQAAPBtR90Rn4qKCuXn50ddfuiynTt3xntJAAAgSRx1R3yqq6uVkZERdXkoFApf35StW7dG/L2wsFCl27bFfoEAAOCYddQd8cnMzFRtbfSvQdbU1ISvBwAA+D6OuuKTn5+vioqKqMsPXVZQUBDvJQEAgCRx1BWf7t27a8uWLdq7d2/E5e+++274egAAgO/jqPuOz9VXX60HHnhACxcuDM/xqa2t1aJFi3Teeec1+RtdvrK21ljnrINWJSk9xrdtZR3yWhlgm9av1cd6InOstyfJvoOso7etw47rjLmt3x0J+9CY6/uqMdjKmOtuzAV5m21pzFknPNcbc1XGnHXfSGqxxZazPhetU8SNuTrruHjZl2idLG+9K9btWV+m+NpRV3zOO+88DR06VHfccYc+++wznXbaafr973+v0tJSPfHEE4leHgAAOIYddcVHkv7rv/5L06ZNizhX18svv6zzzz8/0UsDAADHsKOy+IRCIc2ePVuzZ89O9FIAAEASOeq+3AwAANBcKD4AAMAbFB8AAOANig8AAPAGxQcAAHiD4gMAALxB8QEAAN44Kuf4ILaspzKwnooiyKkRYt2sA52WwcA6Ol6SPjfmWhtzIWMuw5gLtG+sM+53G3NfGnPWU1u0MOYkqcyY+8i4yKJlxg1aH8GTjbkgWhtzZxhz+425UmNO0nHGU1ZYn4vW80YY78oB4+ak2J+KYu93RyTZHxVOWREMR3wAAIA3KD4AAMAbFB8AAOANig8AAPAGxQcAAHiD4gMAALxB8QEAAN6g+AAAAG9QfAAAgDeY3HwMazTmrFM9a40564TnILdtZZ2MGusJz0G2ad0/1ly2MReI9YEpMeZWGHPnGHNBHkDrRN9SY65guy2X87Fxg9YJz5J9RrB1m+2NuSpj7gtjTol7wRgFeW+K8dBo8/uY9dlg/bcAX+OIDwAA8AbFBwAAeIPiAwAAvEHxAQAA3qD4AAAAb1B8AACANyg+AADAGxQfAADgDYoPAADwBpObj2HWyaPWqaOxvl3JPrzVusZYTygNcl+s/0uwztRtbcx1MOYCsY6E3WzMlRlz1knQRcacJLUz5qwjsD835nLeNQatU5El+7PHeqdjPbl5lzEn+3hi64vf+mZizAUZDm59n7C+rKy7xjpNP9YT8pMdR3wAAIA3KD4AAMAbFB8AAOANig8AAPAGxQcAAHiD4gMAALxB8QEAAN6g+AAAAG9QfAAAgDeY3OwB61TPIJNME8U6GdXa6IM0f+v+sc7erTTmrINtrftGkkLG0bENB2w56yTa1p8ag9bpyZLU3ZizPoBZxtxB453pYL3TklqcaAyeYcztN+a+sMXcB8btyT7N2/rkSeAbmXVivPU1aM3VGXMIhiM+AADAGxQfAADgDYoPAADwBsUHAAB4g+IDAAC8QfEBAADeoPgAAABvUHwAAIA3KD4AAMAbTG4+hlmniVpz1sGoQdqydZtW1tu2TkYNMuTVmjUOOzav0Tp71zoAV7JPZLau0Xrb+40b/MFa4waD3PiXxtxOY+4fxlyhMSdJpxkXefoaWy69lS1XvceW22iLSZJKjDnr/rZO8660xYK8XoybNG/T+h5hndpufY/H1zjiAwAAvEHxAQAA3qD4AAAAb1B8AACANyg+AADAGxQfAADgDYoPAADwRlyKz8qVKzVq1Ch17txZWVlZKiws1OjRo1VRUdFkfs2aNerXr5+ysrKUl5en8ePHq6qqKh5LBQAASSwuAwwnT56s3bt3a+jQoTr99NO1detWzZs3Ty+//LLWr1+vvLy8cHb9+vUaOHCgzjzzTM2ZM0c7duzQAw88oJKSEi1fvjweywUAAEkqLsVnzpw56tevn4477p8HmAYNGqT+/ftr3rx5mjFjRvjyKVOm6IQTTtCqVauUk5MjSTr55JM1ZswYvfbaayouLo7HkgEAQBKKS/E5//zzm7wsNzdXmzZtCl+2d+9evf7665o4cWK49EjSyJEjNXHiRD377LMUn+8h1qeiCHKah1h/lmodzW4d9W49JYNkv9/WU0xYx9bvNua2GnOSdHyArEWsT01iPIGCJKndZlvuROupETpYb9iYO9mYk6TTjbnTjLk2xj1pfTJaTxshSWXG3K7Y3naD8bwRlcabDZK1vlat7xGxfl3hawk7V1dVVZWqqqp00kknhS/bsGGD6uvr1atXr4hsenq6unfvrnXr1h12e4WFkSfEKSuzvuoAAIAvEvZbXQ899JDq6uo0fPjw8GWHvuycn58flc/Pz9fOndb/sgEAAEQLfMSnsbFRdXV1pmxGRoZSUlKiLl+9erXuueceDRs2TAMGDAhfXl1dHf65bwuFQuHrm7J1a+SB/sLCQpVu22ZaJwAA8EPgIz6rV69WZmam6c/mzdEfvn/00Ue68sor1bVrVz3++OMR12VmZkqSamtro36upqYmfD0AAMD3EfiIT1FRkRYtWmTKfvsjq7KyMhUXF6tVq1ZatmyZsrOzm8w3Nd+noqJCBQUFQZcLAAAQFrj45OXl6YYbbgh8Q19++aWKi4tVW1urlStXNvk9nq5duyotLU1r167VsGHDwpfX1dVp/fr1EZcBAAAEFZcvN+/fv18/+clPVF5ermXLlun005v+nc1WrVrpwgsv1FNPPaV9+/75O4mLFy9WVVWVhg4dGo/lAgCAJBWXX2e/7rrr9N5772nUqFHatGlTxOyeli1bavDgweG/z5w5U3379lX//v01duxY7dixQw8++KCKi4s1aNCgeCwXAAAkqbgUn/Xr10uSnnzyST355JMR13Xq1Cmi+PTo0UMrVqzQ5MmTNXHiRGVnZ+umm27SfffdF4+lAgCAJBaX4lNaWhoo369fP7399tvNsxgclnUqcjJNE7VOeJYk2xAH+zYrjblSYy7IUF2rkDHXyphrYcxZBwlL9im4+4wTfQs22nLp1onDnxpzkn3asTVn/X2Q7O+OSLK/SQTJWseIG5+M1gnrQV77xqeOOWd9fgfZ3bBL2ABDAACAeKP4AAAAb1B8AACANyg+AADAGxQfAADgDYoPAADwBsUHAAB4g+IDAAC8QfEBAADeiMvkZiRWrKd/JlNbTuR9qTHmdhpzQSZqW58T1oG+HYy5dsZckMfFOvjXOi3Xmku3jv4NMiLY+qSI9ZjuLGMuPcA2rSO1Y/wGZb3ZLwNs0zqk2/qwWNeYTFPyjybJ9G8YAADAEVF8AACANyg+AADAGxQfAADgDYoPAADwBsUHAAB4g+IDAAC8QfEBAADeoPgAAABvMLkZgQUZtJpqzFkbuHV71px16q8khYw56xBc6/asw3z3GXPNsc1Y/w/Kuj5JqjXmrM8J61Td1taJzEHG71oXaWW9M9YHui7Abe825spssYbNttxG481ac5L0D2POOuHZ+vyO9dR9fI0jPgAAwBsUHwAA4A2KDwAA8AbFBwAAeIPiAwAAvEHxAQAA3qD4AAAAb1B8AACANyg+AADAG0xuRph1Smish8sG2aZ10rI1l27MSfZJy9ac9T5bh+Vah/RK0gFjzjqc2Po/KOuUZev0a0lqYcy1M+bMU3WNwYKtxg1KSrXucOsiYz0e/HNjTjJPZN5qvM8fGW/WOOBZG4w5SbI+hNZh1daHGc2DIz4AAMAbFB8AAOANig8AAPAGxQcAAHiD4gMAALxB8QEAAN6g+AAAAG9QfAAAgDcoPgAAwBsUHwAA4A1OWYHAGgJkrc3aevqGWJ82wnpqiyDZWJ/Sw3qaB+tpKCT76S1iPVrfusYgj4v1sbaebcF62oFyYy7PmJOkAuNpHvKMOeupP6zPMet9luynebDmthtz1jVatydJu4w56+vKenogNA+O+AAAAG9QfAAAgDcoPgAAwBsUHwAA4A2KDwAA8AbFBwAAeIPiAwAAvEHxAQAA3qD4AAAAbzC5GYEFmTpqzVqnQVsn+lon1lqn/kpStjFnXWONMbfXmNtnzAW5bevk5lg/fkFYJ2V/GeOccXiyWhtzkpQb421a97d1orZ1+rVkn6Bs3WasJ28HuS/W1xYTmY8NHPEBAADeoPgAAABvUHwAAIA3ElJ8xowZo5SUFF166aVNXv/iiy+qR48eCoVC6tixo+6++27V19fHeZUAACDZxL34rF27Vr/73e8UCjX9tdLly5dr8ODBat26tX77299q8ODBmjFjhsaNGxfnlQIAgGQT19/qcs5p/PjxGjlypFauXNlkZtKkSTr77LP12muvKS3t6+Xl5ORo1qxZmjBhgoqKiuK5ZAAAkETiesRn8eLF+vDDDzVz5swmr9+4caM2btyosWPHhkuPJP385z+Xc05LliyJ11IBAEASilvx2bdvnyZPnqwpU6YoLy+vycy6deskSb169Yq4vKCgQO3btw9fDwAA8H3E7aOu6dOnKzMzUxMnTjxspqKiQpKUn58fdV1+fr527tx52J8tLCyM+HtZmXW8GAAA8EXg4tPY2Ki6ujpTNiMjQykpKdqyZYvmzp2rP/7xj8rIyDhsvrq6Ovxz3xYKhbR3r3WGLY4W1om+1knC1lxrY66FMRcka50kvN+Ys06irTTmJPtEZuvkX+v2rPsmCOtha+v03T3GnHXyb5DnWKwnjsf6uRhkOnilMWfd39Z3f+v2rO8lkv35jWND4OKzevVqXXDBBabspk2bVFRUpAkTJqhv374aMmTIEfOZmZmSpNra2qjrampqwtc3ZevWrRF/LywsVOm2baZ1AgAAPwQuPkVFRVq0aJEpm5+frz//+c965ZVX9MILL6i0tDR8XX19vaqrq1VaWqrc3Fzl5OSEP+KqqKhQhw4dIrZVUVGh3r17B10uAABAWODik5eXpxtuuMGc/+STTyRJV111VdR15eXlOuWUU/Sb3/xGt912m7p37y7p61k/3yw5O3fu1I4dOzR27NigywUAAAhr9i83DxgwQEuXLo26fOzYserUqZPuvPNOdevWTZLUpUsXFRUVaeHChfq3f/s3paZ+/Qn1/PnzlZKSoquvvrq5lwsAAJJYsxefjh07qmPHjlGX33bbbWrXrp0GDx4ccfns2bN1+eWXq7i4WNdcc40+/PBDzZs3T6NHj9aZZ57Z3MsFAABJ7Kg7Semll16qF154Qbt379a4ceP0wgsvaMqUKXrkkUcSvTQAAHCMi+spK77pm190/rbBgwdHHQkCAAD4vzrqjvgAAAA0F4oPAADwRsI+6oIfGo056xRV6+RY62Rbay5I1jrl1Tr515oLMlXXyjp523qfE/k/LesUY+u0auu040pjTrI/x9IDbNPCNotfOhBgm9bnozVn3d/W9xLrexOSD0d8AACANyg+AADAGxQfAADgDYoPAADwBsUHAAB4g+IDAAC8QfEBAADeoPgAAABvUHwAAIA3mNyMo4J1iqp1cqx12nGQ6a3WybFWe4w562Rb6/RkKfb/44n1FFzrxOggrJObrfvGOiHYeruSfWp0kG1aWPd3kMnNsZ60bH1+M5EZ34UjPgAAwBsUHwAA4A2KDwAA8AbFBwAAeIPiAwAAvEHxAQAA3qD4AAAAb1B8AACANyg+AADAGxQfAADgDU5ZgWOKdWx9c5zmodKYs55OoM6YC7JGq1ifEqI5TjFhZT1FQaz3o/V/jUFOLxHrU1FYWR8/6+klJPv+bo7nN3AkHPEBAADeoPgAAABvUHwAAIA3KD4AAMAbFB8AAOANig8AAPAGxQcAAHiD4gMAALxB8QEAAN5gcjOSknWa74EA27ROrU3U9N0gkul/PImcGn20s74OrILs61jfNhAryfT+BwAAcEQUHwAA4A2KDwAA8AbFBwAAeIPiAwAAvEHxAQAA3qD4AAAAb1B8AACANyg+AADAG0xuBoysk2iPhYm11v/xHDTmEjmtOpkmNx8Lzx3gWMcRHwAA4A2KDwAA8AbFBwAAeIPiAwAAvEHxAQAA3qD4AAAAb1B8AACANyg+AADAGxQfAADgDSY3Ax6K9YRgJg4DOFZwxAcAAHgjrsVnxYoVGjBggFq1aqXs7Gz17NlTzzzzTFTuxRdfVI8ePRQKhdSxY0fdfffdqq+vj+dSAQBAEorbR12LFi3STTfdpIsuukizZs1SamqqNm/erLKysojc8uXLNXjwYP3oRz/Sb3/7W23YsEEzZszQZ599pvnz58druQAAIAmlOOdcc99IaWmpzjrrLI0ZM0Zz5849YrZLly46/vjjtXbtWqWlfd3Lpk6dqlmzZmnjxo0qKioy3WZhYaFKt21T1v959QAA4Gh3QNLJp5yirVu3HjEXl4+6Hn30UTU0NGj69OmSpKqqKjXVtzZu3KiNGzdq7Nix4dIjST//+c/lnNOSJUvisVwAAJCk4lJ8VqxYoaKiIi1btkzt27dXdna2TjzxRE2bNk2Njf/8fZB169ZJknr16hXx8wUFBWrfvn34+qYUFhZG/Pn2R2gAAABx+Y5PSUmJUlNTdeONN+r222/XOeecoxdeeEEzZsxQfX297rvvPklSRUWFJCk/Pz9qG/n5+dq5c2c8lgsAAJJU4OLT2Niouro6UzYjI0MpKSmqqqpSY2Oj7r//fk2ePFmSNGTIEO3evVtz587VlClTlJ2drerq6vDPfVsoFNLevXsPe1vf/kzv0Hd8AAAADgn8Udfq1auVmZlp+rN582ZJUmZmpiRpxIgREdsaMWKEqqurwx9hHcrV1tZG3W5NTU34egAAgO8j8BGfoqIiLVq0yJQ99JFVQUGBSkpK1K5du4jr27ZtK0n66quvIvIVFRXq0KFDRLaiokK9e/cOulwAAICwwMUnLy9PN9xwQ6Cf6dmzp0pKSlReXq7CwsLw5Ye+s9OmTRtJUvfu3SVJa9eujSg5O3fu1I4dOzR27NigywUAAAiLy291DR8+XJL0xBNPhC9rbGzUokWLlJubq549e0r6eoZPUVGRFi5cqIaGhnB2/vz5SklJ0dVXXx2P5QIAgCQVl9/quuKKKzRw4EDdd999+uKLL3TOOefoT3/6k9566y0tWLAg4svMs2fP1uWXX67i4mJdc801+vDDDzVv3jyNHj1aZ555ZjyWCwAAklRcJjdLXw8tnDp1qp555hnt3r1bZ5xxhiZPnqzrrrsuKvunP/1J99xzjzZt2qQ2bdrohhtu0F133aXjjz/efHtMbgYAwB/Wyc1xKz7xRvEBAMAfR9UpKwAAAI4GFB8AAOANig8AAPAGxQcAAHiD4gMAALxB8QEAAN6g+AAAAG9QfAAAgDcoPgAAwBsUHwAA4A2KDwAA8AbFBwAAeIPiAwAAvEHxAQAA3qD4AAAAb1B8AACANyg+AADAGxQfAADgDYoPAADwBsUHAAB4g+IDAAC8QfEBAADeoPgAAABvUHwAAIA3KD4AAMAbFB8AAOANig8AAPAGxQcAAHiD4gMAALxB8QEAAN6g+AAAAG9QfAAAgDcoPgAAwBsUHwAA4A2KDwAA8AbFBwAAeIPiAwAAvEHxAQAA3qD4AAAAb6Q451yiF9EcMjMzVVNTo5RELwQAADQ7JykUCqm6uvqIubT4LCf+MjIyJEn5+fkJXsn3V1ZWJknq0KFDgldydGB/RGOfRGOfRGJ/RGOfREuGfVJRURH+t/9IkvaITzIoLCyUJG3dujXBKzk6sD+isU+isU8isT+isU+i+bRP+I4PAADwBsUHAAB4g+IDAAC8wXd8AACANzjiAwAAvEHxAQAA3qD4AAAAb1B8AACANyg+AADAGxQfAADgDYpPgn3wwQe6/PLLlZubq6ysLHXt2lUPP/xwRGbNmjXq16+fsrKylJeXp/Hjx6uqqipBK46vmTNnKiUlRV27do26zof98v777+sXv/iFunTpohYtWqhjx44aNmyYtmzZEpXdtGmTBg0apJYtWyo3N1c//elP9fnnnydg1YlRW1uryZMnq6CgQJmZmTrvvPP0+uuvJ3pZCbVy5UqNGjVKnTt3VlZWlgoLCzV69GhVVFQ0mffhNfVNY8aMUUpKii699NImr3/xxRfVo0cPhUIhdezYUXfffbfq6+vjvMr4WLFihQYMGKBWrVopOztbPXv21DPPPBOVS4p94pAwr776qktPT3fnnXeemzNnjlu4cKGbPHmy+9WvfhXOrFu3zoVCIXfuuee6+fPnuzvvvNNlZGS4QYMGJXDl8VFWVuaysrJcixYtXJcuXSKu82W/DBkyxOXl5blx48a5xx57zN17772uXbt2rkWLFm7Dhg3hXFlZmTvppJPcqaee6ubOnetmzpzpTjjhBHfOOee42traBN6D+LnmmmtcWlqamzRpkluwYIHr06ePS0tLc2+++Wail5YwPXv2dKeccoq7/fbb3WOPPebuuOMOl52d7dq1a+cqKioisr68pg55//33XVpamguFQu6SSy6Jun7ZsmUuJSXFXXDBBW7hwoVu3Lhx7rjjjnM333xzAlbbvJ588kmXkpLiiouL3bx589z8+fPdbbfd5mbPnh2RS5Z9QvFJkD179rh27dq5K6+80jU0NBw2d/HFF7v8/Hy3Z8+e8GWPPfaYk+ReffXVeCw1YYYPH+4GDBjg+vfvH1V8fNkvb7/9dlRx2bJli8vIyHDXXXdd+LJbbrnFZWZmuu3bt4cve/31150kt2DBgritN1HeffddJynijbq6utqdeuqprk+fPglcWWK98cYbUe8vb7zxhpPk7rzzzojLfXlNOedcY2Oj69Onjxs1apTr1KlTk8XnrLPOcuecc447ePBg+LI777zTpaSkuE2bNsVzuc1q27ZtLjMz040fP/47s8myTyg+CTJ//nwnyW3cuNE551xVVVXUG9SePXtcWlpaxBEg55yrra11LVu2dDfddFPc1htvb7zxhktNTXV/+9vfooqPz/vlkB49ergePXqE/962bVs3dOjQqFznzp3dwIED47m0hPjVr37lUlNTI/7Rds65WbNmOUnuk08+SdDKjk65ubnuqquuCv/dt9fU73//e5edne0qKiqaLD5///vfnST3yCOPRFxeXl7uJLl77703nsttVpMnT3bp6emusrLSOefcvn37XGNjY1QumfYJ3/FJkBUrVignJ0fl5eU644wz1LJlS+Xk5OiWW25RTU2NJGnDhg2qr69Xr169In42PT1d3bt317p16xKx9GbX0NCgcePGafTo0erWrVvU9b7ul0Occ/r000910kknSZLKy8v12WefRe0PSerdu3fS7w9JWrdunTp37qycnJyIy3v37i1JWr9+fQJWdXSqqqpSVVVV+Pkj+fWa2rdvnyZPnqwpU6YoLy+vycyh+/vt/VFQUKD27dsn1f5YsWKFioqKtGzZMrVv317Z2dk68cQTNW3aNDU2NoZzybRPKD4JUlJSovr6el1xxRX68Y9/rOeff16jRo3So48+qhtvvFGSwl9AzM/Pj/r5/Px87dy5M65rjpdHH31U27dv17333tvk9b7ul0OefvpplZeXa/jw4ZK+e3/s3r1btbW1cV1jvFVUVBz2/ktK+udEEA899JDq6urCzx/Jr9fU9OnTlZmZqYkTJx4249P+KCkpUVlZmW688UaNGjVKS5Ys0cUXX6wZM2bozjvvDOeSaZ+kJXoBvqqqqtKBAwd08803h3+L66qrrlJdXZ0WLFig6dOnq7q6WpKUkZER9fOhUCh8fTL58ssvddddd2natGlq06ZNkxkf98shH330kW699Vb16dNHP/vZzyR99/44lGnq+mRxuPv3zft/rGtsbFRdXZ0pm5GRoZSUlKjLV69erXvuuUfDhg3TgAEDwpcfi6+p77M/tmzZorlz5+qPf/zjEV8P37U/9u7d+/0W3cy+zz6pqqpSY2Oj7r//fk2ePFmSNGTIEO3evVtz587VlClTlJ2dfczuk6ZwxCdBMjMzJUkjRoyIuPzaa6+VJL3zzjvhTFP/W6+pqQlfn0ymTp2q3NxcjRs37rAZH/eLJO3atUuXXHKJWrVqpSVLlig1NVXSd++Pb2aSVWZmZtLf/9WrVyszM9P0Z/PmzVE//9FHH+nKK69U165d9fjjj0dcdyy+pr7P/pgwYYL69u2rIUOGHHHbx+L+kL7fPjncv0UjRoxQdXV1+COsY3WfNIUjPglSUFCgv//972rXrl3E5W3btpUkffXVVzr11FMlqcmZGxUVFSooKGj+hcZRSUmJFi5cqIceeijisGlNTY0OHjyo0tJS5eTkhA+1+rJfJGnPnj26+OKLVVlZqTfffDPiPn7X/sjNzU3qoz3S1/ugvLw86vJD+yQZnhNFRUVatGiRKfvtjyPKyspUXFysVq1aadmyZcrOzm4yfyy9poLujz//+c965ZVX9MILL6i0tDR8XX19vaqrq1VaWqrc3Nyo95gOHTpEbKuioiL83bGjzfd5jhQUFKikpOSI/xZ9M3+s7ZMmJfrb1b769a9/7SS5lStXRly+cuVKJ8k9/fTTrrKy8oi/aTFq1Kh4LrnZ/e///q+TdMQ/EyZM8G6/VFdXu3/91391WVlZbs2aNU1m2rRpc9jf6howYEBzLzHhJk2a1ORvdc2cOdP73+r64osvXFFRkWvbtq3bsmVLkxkfXlOLFi36zveX3/zmN8455z788MMj/gbT9OnTE3APmsc111zjJLmPP/444vInnnjCSXJvv/22cy659gnFJ0E++OADJ8lde+21EZePGDHCpaWlufLycuecc4MGDXL5+flu79694czjjz/uJLnly5fHdc3N7fPPP3dLly6N+tOlSxfXsWNHt3TpUve3v/3NOefPfqmvr3eXX365S0tLc//zP/9z2NzNN9/sMjMzI/6BX7FihZPk5s+fH4+lJtRf/vKXqDk+NTU17rTTTnPnnXdeAleWWFVVVa53794uOzvbrV279ojZZH9Nbd++vcn3lzZt2rhevXq5pUuXun/84x/hfFFRkTvnnHNcfX19+LKpU6e6lJSU8BiSZLB06VInyU2ZMiV8WUNDg+vXr5/Lzc11NTU14cuTZZ9QfBJo1KhRTpIbNmyYe+SRR9zQoUOdJHfHHXeEM3/9619dRkZGxDTVUCjkiouLE7jy+GpqgKEv+2XChAlOkrvsssvc4sWLo/4c8sknn7gTTzzRnXrqqe7hhx92s2bNcieccILr1q1bxBtXMhs6dGj4qMWCBQtc3759XVpamnvjjTcSvbSEueKKK5wkN2rUqKjnztKlSyOyvrymvu1wAwxfeukll5KS4gYMGOAWLlzoxo8f74477jg3ZsyYBKyy+TQ2NrqBAwe6lJQUN3bsWPfII4+4iy66qMnhp8myTyg+CVRXV+f+4z/+w3Xq1Mkdf/zx7rTTTgsfav2mN9980/Xt29eFQiHXpk0bd+utt0b8ryzZNVV8nPNjv/Tv3/+Ih+a/6cMPP3TFxcUuKyvLtW7d2l133XVu165dCVp5/FVXV7tJkya5vLw8l5GR4X74wx+6V155JdHLSqhOnTod9rnTqVOnqLwPr6lvO1zxce7royHdu3d3GRkZrn379m7q1Kmurq4uzitsfvv27XMTJkxweXl5Lj093XXr1s099dRTTWaTYZ+kOOdcrL83BAAAcDTi19kBAIA3KD4AAMAbFB8AAOANig8AAPAGxQcAAHiD4gMAALxB8QEAAN6g+AAAAG9QfAAAgDcoPgAAwBsUHwAA4A2KDwAA8Mb/B6Sil9ZLhRwgAAAAAElFTkSuQmCC",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots()\n",
"ax.pcolormesh(source.pixel_centers, source.pixel_centers, np.flip(source.source, axis=1))\n",
"ax.set_aspect(\"equal\")\n",
"ax.set_xlim(source.pixel_centers[-1], source.pixel_centers[0])\n",
"ax.set_ylim(source.pixel_centers[0], source.pixel_centers[-1]);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This image is sharper than the normal IXPE data because this is a model for where the source photons arrive from, not where IXPE detects them. LeakageLib can uses this model to separate sources which overlap.\n",
"\n",
"The map doesn't include the pulsar because of pileup, so it is a nebula-only spatial weight map (which is what we want to measure the nebula polarization). The pulsar's readout streak is barely visible and should be removed in a more careful analysis."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "base",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.2"
}
},
"nbformat": 4,
"nbformat_minor": 2
}