{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "title",
   "metadata": {},
   "source": [
    "# Pointing Cut Example\n",
    "\n",
    "This example reads a DC4 unbinned source FITS file, builds good time intervals for a source using a 60 degree field-of-view cut, and writes a binned `.hdf5` file containing events either inside or outside that FOV selection."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "overview",
   "metadata": {},
   "source": [
    "The FOV selection uses `GoodTimeInterval.from_pointing_cut`, with the default example value `max_offaxis = 60 * u.deg`, to keep events when the source is inside the FOV.\n",
    "\n",
    "**When making a pointing cut in an analysis (e.g., spectral analysis), the data, background, and orientation files need to all be cut self-consistently.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "imports",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "import numpy as np\n",
    "import astropy.units as u\n",
    "from astropy.coordinates import SkyCoord\n",
    "from astropy.time import Time\n",
    "import matplotlib.pyplot as plt\n",
    "from histpy import Histogram\n",
    "\n",
    "from cosipy import BinnedData\n",
    "from cosipy.event_selection import GoodTimeInterval\n",
    "from cosipy.spacecraftfile import SpacecraftHistory\n",
    "from cosipy.util import fetch_wasabi_file\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "download-input-files",
   "metadata": {},
   "outputs": [],
   "source": [
    "fetch_wasabi_file(\n",
    "    \"COSI-SMEX/DC4/Data/Mock_Dataset/dc4_mock_dataset_3months_unbinned_data_filtered_with_SAAcut_time_ordered.fits.gz\",\n",
    "    checksum=\"762b7ad849b501ddaa31799ec19d14b9\",\n",
    ")\n",
    "\n",
    "fetch_wasabi_file(\n",
    "    \"COSI-SMEX/DC4/Data/Backgrounds/Total_DC4_BG_3months_unbinned_data_filtered_with_SAAcut_withSAAbck.fits.gz\",\n",
    "    checksum=\"73c8b23e43684da3b35c25499254202b\",\n",
    ")\n",
    "\n",
    "fetch_wasabi_file(\n",
    "    \"COSI-SMEX/DC4/Data/Orientation/DC4_final_530km_3_month_with_slew_15sbins_GalacticEarth_SAA.fits\",\n",
    "    checksum=\"ca94ff1d7a73c1f41479aaf598807673\",\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "configuration",
   "metadata": {},
   "outputs": [],
   "source": [
    "### Input files to apply pointing cuts\n",
    "data_file = Path(\n",
    "    \"dc4_mock_dataset_3months_unbinned_data_filtered_with_SAAcut_time_ordered.fits.gz\"\n",
    ")\n",
    "\n",
    "bkg_file = Path(\n",
    "    \"Total_DC4_BG_3months_unbinned_data_filtered_with_SAAcut_withSAAbck.fits.gz\"\n",
    ")\n",
    "\n",
    "orientation_file = Path(\n",
    "    \"DC4_final_530km_3_month_with_slew_15sbins_GalacticEarth_SAA.fits\"\n",
    ")\n",
    "\n",
    "### Output names for the binned files after the pointing cut\n",
    "binned_output_file = Path(\n",
    "    \"dc4_mock_dataset_3months_binned_data_filtered_with_SAAcut_time_ordered.hdf5\"\n",
    ")\n",
    "binned_bkg_output_file = Path(\n",
    "    \"Total_DC4_BG_3months_binned_data_filtered_with_SAAcut_withSAAbck.hdf5\"\n",
    ")\n",
    "\n",
    "### This config is only used to initialize the DataIO object\n",
    "config_file = Path(\"inputs.yaml\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "922f56e4",
   "metadata": {},
   "source": [
    "Change details below for the source to be analyzed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "c6c6a1dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "source_name = \"NGC 4151\"\n",
    "source_coord = SkyCoord(l=155.07 * u.deg, b=75.06 * u.deg, frame=\"galactic\")\n",
    "\n",
    "### Default FOV cut used by the pointing-cut helper.\n",
    "max_offaxis = 60 * u.deg\n",
    "earth_occ = True"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "read-data",
   "metadata": {},
   "source": [
    "Read the unbinned FITS data and open the spacecraft orientation over the event time range."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "load-data-and-orientation",
   "metadata": {},
   "outputs": [],
   "source": [
    "analysis = BinnedData(str(config_file))\n",
    "\n",
    "analysis.cosi_dataset = analysis.get_dict(str(data_file))\n",
    "time_tags = analysis.cosi_dataset[\"TimeTags\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "344dd39a",
   "metadata": {
    "id": "cut-orientation-file"
   },
   "source": [
    "Load in the orientation file.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "7d5edebc",
   "metadata": {
    "id": "save-or-load-time-cut-orientation"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time-cut orientation pointings: 531,998\n"
     ]
    }
   ],
   "source": [
    "orientation = SpacecraftHistory.open(orientation_file)\n",
    "\n",
    "print(f\"Time-cut orientation pointings: {len(orientation.obstime):,}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "make-gti",
   "metadata": {},
   "source": [
    "Build the GTIs where the source is inside the FOV."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "build-gti",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Source: NGC 4151\n",
      "FOV cut: off-axis <= 60.0 deg\n",
      "GTI intervals: 675\n",
      "GTI livetime: 980,415.0 s\n",
      "Total livetime: 6,579,555.0 s\n",
      "Livetime fraction in FOV: 0.1490\n"
     ]
    }
   ],
   "source": [
    "source_gti = GoodTimeInterval.from_pointing_cut(\n",
    "    source_coord,\n",
    "    orientation,\n",
    "    max_offaxis,\n",
    "    earth_occ=earth_occ,\n",
    ")\n",
    "\n",
    "pointing_cut_orientation = orientation.apply_gti(source_gti)\n",
    "gti_livetime = pointing_cut_orientation.cumulative_livetime().to_value(u.s)\n",
    "total_livetime = orientation.cumulative_livetime().to_value(u.s)\n",
    "livetime_fraction = gti_livetime / total_livetime\n",
    "\n",
    "print(f\"Source: {source_name}\")\n",
    "print(f\"FOV cut: off-axis <= {max_offaxis.to_value(u.deg):.1f} deg\")\n",
    "print(f\"GTI intervals: {len(source_gti)}\")\n",
    "print(f\"GTI livetime: {gti_livetime:,.1f} s\")\n",
    "print(f\"Total livetime: {total_livetime:,.1f} s\")\n",
    "print(f\"Livetime fraction in FOV: {livetime_fraction:.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "mask-events",
   "metadata": {},
   "source": [
    "Apply the GTI to the mock-data file. This uses the same half-open convention as `select_data_time`: events are kept for `start <= TimeTags < stop`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "event-mask",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Events in FOV: 26,164,190\n",
      "Events outside FOV: 156,163,781\n"
     ]
    }
   ],
   "source": [
    "in_fov_mask = source_gti.contains(analysis.cosi_dataset[\"TimeTags\"])\n",
    "n_in_fov = np.count_nonzero(in_fov_mask)\n",
    "n_out_fov = len(in_fov_mask) - n_in_fov\n",
    "\n",
    "analysis.cosi_dataset = {\n",
    "    key: values[in_fov_mask]\n",
    "    for key, values in analysis.cosi_dataset.items()\n",
    "}\n",
    "selected_time_tags = analysis.cosi_dataset[\"TimeTags\"]\n",
    "analysis.tmin = float(np.min(selected_time_tags))\n",
    "analysis.tmax = float(np.max(selected_time_tags))\n",
    "\n",
    "print(f\"Events in FOV: {n_in_fov:,}\")\n",
    "print(f\"Events outside FOV: {n_out_fov:,}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "f82b36b3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wrote: dc4_mock_dataset_3months_binned_data_filtered_with_SAAcut_time_ordered.hdf5\n"
     ]
    }
   ],
   "source": [
    "analysis.get_binned_data(\n",
    "    make_binning_plots=False,\n",
    "    show_plots=False,\n",
    ")\n",
    "analysis.binned_data.write(binned_output_file, overwrite=True)\n",
    "\n",
    "print(f\"Wrote: {binned_output_file}\")\n",
    "\n",
    "del analysis.cosi_dataset, analysis.binned_data, in_fov_mask"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e5771856",
   "metadata": {
    "id": "apply-pointing-cut-background"
   },
   "source": [
    "Apply the same pointing-cut time intervals to the background file, following the same unbinned-selection and binning pattern used in the other tutorial notebooks.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "35469f50",
   "metadata": {
    "id": "cut-and-bin-background-file"
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "fill() discarded one or more values due to out-of-bounds coordinate in a dimension without under/overflow tracking\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Background events in FOV: 25,071,036\n",
      "Background events outside FOV: 143,577,508\n",
      "Wrote: Total_DC4_BG_3months_binned_data_filtered_with_SAAcut_withSAAbck.hdf5\n"
     ]
    }
   ],
   "source": [
    "background = BinnedData(str(config_file))\n",
    "\n",
    "background.cosi_dataset = background.get_dict(str(bkg_file))\n",
    "bkg_in_fov_mask = source_gti.contains(background.cosi_dataset[\"TimeTags\"])\n",
    "\n",
    "background.cosi_dataset = {\n",
    "    key: values[bkg_in_fov_mask]\n",
    "    for key, values in background.cosi_dataset.items()\n",
    "}\n",
    "selected_bkg_time_tags = background.cosi_dataset[\"TimeTags\"]\n",
    "background.tmin = float(np.min(selected_bkg_time_tags))\n",
    "background.tmax = float(np.max(selected_bkg_time_tags))\n",
    "\n",
    "background.get_binned_data(\n",
    "    make_binning_plots=False,\n",
    "    show_plots=False,\n",
    ")\n",
    "background.binned_data.write(binned_bkg_output_file, overwrite=True)\n",
    "\n",
    "print(f\"Background events in FOV: {np.count_nonzero(bkg_in_fov_mask):,}\")\n",
    "print(f\"Background events outside FOV: {len(bkg_in_fov_mask) - np.count_nonzero(bkg_in_fov_mask):,}\")\n",
    "print(f\"Wrote: {binned_bkg_output_file}\")\n",
    "\n",
    "del background.cosi_dataset, background.binned_data, bkg_in_fov_mask"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "plot-saved-counts-text",
   "metadata": {},
   "source": [
    "Load the saved pointing-cut binned files and compare the mock-data and background counts.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "plot-saved-counts",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxAAAAF/CAYAAADZxC9bAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAPelJREFUeJzt3Qd4VGXe/vFfAgklREBQQFgLZSkWUCm6KKICiihWUGxYiAUVXXURWPuLa30tgOKCKKuuUkRcVJoUBVFcRRGkBLCgSEeJJJEAIf/rfvyfvMkwDCeQZMr5fq7rXCEzZ84852Qcn/s8LamgoKDAAAAAAMCHZD87AQAAAAABAgAAAECJ0AIBAAAAwDcCBAAAAADfCBAAAAAAfCNAAAAAAPCNAAEAAADANwKED9u3b7fMzEz3EwAAAAgyAoQPq1evtoyMDPczWnJycqL23gAQDXzvAQianDip7xEg4sTu3bujXQQAKFd87wEImt1xUt8jQAAAAADwjQABAAAAwLeK/ncNnhkzZrgtOzs72kUBAAAAYgIBIoJOnTq5TTMwaRA1AACAJz8/33bu3MkFQanZsWNHmc/6mZKSYhUqVDigYxAgAAAASki9E9asWWMFBQVcO5TqIOrk5LIdYZCUlGQNGjSwatWq7fcxCBAAAAAlbHlQeKhataodcsghrkIGlNZn60BbByJR4N20aZP7/DZp0mS/34sAAQAAUALqtqSKmMJDlSpVuHaImwAh+tz+8MMP7nNMgCgDDKIGAAB7Q8sDgvq5pQUiAgZRAwAAAMWxDgQAAEACOPLII61p06bWqlUra9GihT3//PP7fM0555zjZpvcl3feecfmz59f+PsXX3xhl156qcWC0LJF447+1q1bS/y62rVru65E+/Lggw+W+cxMJUWAAAAASBBjx461hQsX2pQpU2zQoEG2aNGiiPtPnjzZhY6SVtJbt27t3isWRDtAlLWHHnoo5gIEXZgQ876bv84WvLXSdm7fZUGUUrmite7RxI5qVy/aRQEA7MU7f59nuVl5ZXZ9qlavZBc80t73/kcccYQLBitWrLC6devaTTfdZCtXrnSDv2+77Ta78cYbC1stVAFXq0XHjh1dMPjss89s7dq11rlzZ3vxxRddyJg0aZJ98MEHNnr0aLv11lutcePGdscdd7iworvoev3tt99u7733nmVlZdmQIUNc64b85z//sQEDBlhqaqqdffbZNmrUKNeCofcOtWzZMnfcdevWud/79u3ryq6y6fELLrjAPX7JJZfYueeea4ceeugeZevTp0+xY+q1J554on3++eeurL1797aTTz7Z/vGPf7jZiPr162d33nmn21fl0u+aprdy5cr2zDPPWPv2f1z3999/37UG7Nixw7U6/POf/7R27doVvo+urc5T5zBmzBg3S1dRKuc999zj1mHQdSjq7rvvto8++sgNbD7ooINs5MiR7u+nc5dTTz3VDXiePn26G6P73HPPuXJo2tfBgwfbeeedZ+WJAIGYp/CQtTbHgivPvhi/kgABADFM4SH3l7ILECW1ePFiW758ubVs2dIFBlVG3377bdu4caOrTOvxk046aY/XffvttzZ79mxXkVU3qE8//dQFge7du7uQoEq8fPjhh8Vep9Bw3HHHubvlU6dOdWFCr9P7XXfddTZv3jxr1qyZvfLKK7Zly5awZd61a5edf/757hi9evVyj23evDnieYYrWzirV6925/Xbb7+54PLrr7/a3LlzXVDStVEZVeG/6KKLXOX9rLPOso8//tguvvhiW7Vqldvv2muvtTlz5rjz0PXJzc0tPH5eXp4rc61atWzixIl7zG6k66DX6z11XUeMGFHsOihYPPXUU24WpvHjx7vrp+uoAKegotfVqFHD7auy6b0UYhSI9HfU+VWqVMnKCwEiAmZhig2L6/5oc7pk2o7UXZaUHKy5tgt2F1jqjorW4eum1iPahQEARGwhiIXja1yCppZVZfjll192c/2rPrNgwQL3vO7Yq5Ksx8IFCL2+YsWKblOlXIFCd+v3RXfrdVzR/nqdqGuRgoUq3aK7/95d9VAai6GuOl548MYJlAa1WKhSX7NmTWvYsKFrvVAFvH79+oXTmup5LeKmCrqccsopVqdOHdfK8vXXX7tWA+88UlJSrHr16oXH79atmws/9913X9j3966DwoNcf/31Lth51IIydOhQ27Ztm2tV+OWXX/Z6Lt9//71dccUVrvVEfyftq8e8spUHAkQEzMIUG+Yev8J+rRHcFogcy7O5J6yIdjEAABGUpHtRWdK4BFX893caTwUBjyrUahXwQ3e/vePqdbqTvi8KMeq6Iz169CjsnhSOKspFj1nSMQGh5xXuPMOtieB3ytMzzjjDhQC1HKgL0r4UPe6PP/7oul6pi5VaR5YsWWIdOnTY62svu+wye+yxx1wokoMPPrjcx0gQIBDzdqT88eWVtNvskEr/l/aDYFNelhUk/981AABgf26IqlvOI4884lYhVlcmdZMpCVWK1U2ppNTKoYHcal1QV6HXX3/d9d33yqW7+x5V4tVy8uabbxbrwqRWCI250NgMdSnS3XZ1L/Iq0PtbtlAqn+7+Kwho/Mcnn3xi69evd4FM40gefvhh1y2sWZEuTF4rhAasawyGzkkD2NWVqSi1yqgLk/d6tQ5510FlV4tGvXr13DiKYcOGFXttenq628frwqTuV0cddZT7t66nfi9vBAjEjbTfK9sHbR6wIGk/Z5Blp8XWzAsAgPiiAc0333yzHXvssa6C+ve//73Y4F8/rrrqKrvmmmvcgOtbbrnFVej9UJepl156ybUuqJVCFfNq1aoVVoZDWxk04FpdezTAWd2JNIhaA7779+/vulfpHI4++uhi5Q8tW+ggar80yFvhSoOo77rrLtdK8dZbb7ny6nw1fuPKK68sXMFZ4xPatm1b+HqNwUhLS3OtEdOmTXOhw6NuUgoNF154YeFgci9k6JzUqqDz0mOhLTEqi66bwpUGUWsAtcKTrqHe6/DDD7fyllSgTxIiUmrOyMgoHBEfDeoTpwQaRF4lulpOZZvX4R8WJEE+dyDI33uIbeouorvgugtctCsM9v3fsir5AwcOdDMVYU/qphWuK1WsfX5pgQAAAECZ0eBgjc1Q5Vjdjf79739zteMcAQIAAABlRuMDtCFxECAiYBpXAAAAoDgCRARM4woAAAAUR4AA4mRBuTdunWVBlFK5orXu0YSVuAEAiBEECCCGFV2/JveXPAumPPti/EoCBAAAMYIAAcSw1Cp//CealJxkVQ+uZEHz+695pommd25nIT0AAGIFAQKIYRVSK5jtNMtNy7NXbphrQZO7Nc9StlewDl9HZ/0VAIgnRx55pFusrUqVKm6VYy2qpi2Sc845x5555pl9rnOl9Ru0MJpWlpYvvvjCnnzySTc9a7SFli1Ux44d3SJvoQu07a+OpXy8sqQF584991y30F5pIkAAMSytQiUXIHZbgW3cmWWBk/bHNveEFdEuCQDEBVXoW7VqZatXr7bjjjvOTj31VPdzbyZPnuy7kq7jepX01q1bx0R4CFe2eLBr1y638na8it+SAwFwy2Fd7fm1UywnP5jjHzblZVlBstmOFLowAYhtvZY9bZt3biuz49dOSbc3m9/pe/8jjjjCtSqsWLHC3Z2/6aabbOXKlVZQUGC33Xab3XjjjYWtFl4FXHfWFQw+++wzW7t2rXXu3NlefPFFFzImTZpkH3zwgY0ePdpuvfVWa9y4sbsLv3DhQvvhhx/c62+//XZ77733LCsry4YMGeJaN+Q///mPDRgwwFJTU+3ss8+2UaNGuRYMvXcorVCt465bt8793rdvX1f20Lv+3p31Qw89dI+y9enTZ4/jzpw50x555BH79ddf7fzzz7ennnrKkpKS7Omnn7Y333zTdu7caSkpKa7cJ598csSyFDVhwgR7+OGH7e2337ZGjRrZAw884BbKq1mzpp111ln2+uuvu+vjXSNdd5X16quvduXX8TZu3GjJycn24IMP2nnnneeOq7KprDVq1HC/165du/CaadPrdZz169fb9ddfb/fee6/bb/ny5Xbddde5v0GTJk0sNzfXygIBIgLWgUC0da7Z0m1B1X7OIMtO2x7tYgDAPik8xFJL8eLFi11lsmXLli4wKEyokqvK6oknnugeD3fH/ttvv7XZs2e7CnWLFi3s008/dUGge/furgKsCrV8+OGHxV6nCqtaOh566CGbOnWqCxN6nd5PFdp58+ZZs2bN7JVXXrEtW7bs9a68Kvc6Rq9evdxjmzdvjnie4coWztKlS+2TTz5x59WhQwcXGi6//HK76qqr7M47/whm8+fPd119dN38lOXpp5+2iRMn2qxZs6xWrVr2/vvvu0Dx1VdfWbVq1dx5h16jo48+2h5//HH3e7t27dw+ChUKd/p7fP7559awYUPbl61bt7q/jcqk4HLttdda/fr13fkolChU6DOgQKjzLG0EiAhYBwIAAPhtIYiF41966aVuDETVqlXt5ZdfdnehdUN0wYIF7nndsb/ooovcY+EChF6vrjXaVClXoPDuyEdSuXJld1zR/nqdVylXsFB4kN69e+9xF9+TmZlp27dvL6ywu/OuXdtKg+7Yq4VB25VXXunOXxVrVfbVMqFQo3NWGX7//Xf77rvvIpZl8ODBVqdOHdcKoHP3Wjl69Ohh6el//K1UiVcY83jvLdu2bbMvv/zSBSvR3+mUU06xjz/+2FeA8EKByqT9v//+e/e+ahHyxjsce+yx7phlgQABAABwgErSvag8xkBEou4xe+NVhqVChQruTrwfGrztHVevy8/P3+drVIm/++673b9V8Y40KFmV+6LHVOX+QKisGmiu0KNKfps2bey3336z6tWrW17evrsNt2vXzqZPn+6Chlpq9vYeRSnUqatSpDJ5Qq9h6Pn6/TtF+lsfiL2fBQAAABKiR8XIkSPdvzdt2uS6Mml8Q0kcdNBBrgtOSamVY9GiRe7OvmhMgCruXrl0x1zb3//+d9fNSpVsdS/yeN2GNOZCYzNEd9t1p74kZdP7qvuSWhfeeOMN996qlKsshx9+uNtn6NChhftHKovo+qmFR2MW1JIgZ5xxhuvClJ2d7caa6Pm9UWvBCSec4Lp0yapVq9w5adB76Pnq75WTk2P7outw/PHH26uvvup+X7JkSbHrVJoIEAAAAAlMA4M1IFhdWk4//XRXWdcd9JJQ3/px48a5CupLL73k+3XqMqX91bqglhH1y9f4AG9wcGgrgwZcq1Ktsmqchirk0r9/f9dSoMcHDhxYrPx+yta8eXNr3769e70q6ZdddpmrcKsrUtu2bd24EA3y9lMWj44zZswYN6BbXZE0KFrjJnSeatHQOYY7T48GW6vFSMfWMVR2L8xoal2NI1HIUDcrjbHwQ+FhxIgRdswxx7iB1RrvURaSChSREJFSc0ZGhkvv+5onuayor5zXpy6oA2mr5VS2eR3+Ee3ioBzxtw+2IH/vIbbpzrXugh911FHFupJg3/8ta8YnBQAFmkQ+14KCArvrrrtci8fw4cN9v17dltQlKdY/v4yBAAAAQJlR1yDdaVflWHf9dec9UV199dVuylZV0jXjkqbBTUQECAAAAJSZQYMGuS0IJk6caEHAGAgAAID9QC9wBPVzSwsEAABACWg+f02PqRmNDjnkkDKbKhPBk1/GYyAUHvS51WdWn+P9RYCIgJWoAQBAKFXwGjRoYGvWrHH93YHSsnv37ohrRZQGhQd9fg8kqBAgImAlagAAEI6mItXqwVpbACgtWu8hLS3NypJaHg60lYMAASDmFewusDdunWVBlFK5orXu0cSOalcv2kUBEEKVsLKechPBsnPnzriYGpgAASBmFe1WnPtLngVTnn0xfiUBAgAQMwgQceC7+evs83GZlr9jtwVRwSWsdRhUqVX++IrKTcuzV26aa0FseUndUdE6fN3UekS7MAAA/H8EiDiw4K2Vtm397xZ0THIRPAdVqWpbtmdbQbK51ciDKMfybO4JK6JdDAAAChEg4sDiuj/anC6ZtiN1lyUlB2+quNwqecXuRiM4bjmsqz2/dorl5Aez+9KmvCwXnnak7Ip2UQAAKESNLA7MPX6F/Vojx4JOd6MRLJ1rtnRbULWfMyiwLS8AgNhFgIgD3t3HpN1mh1SqbkGUVqGSuxsNAACA6CJAxJG03yvbB20eiHYxAAAAEGBlu9QdAAAAgIRCgAAAAADgG12YIpgxY4bbsrOz/V9RAAAAIIERICLo1KmT2zIzMy0jI6P8/ioAAABAjKILEwAAAADfCBAAAAAAfCNAAAAAAPCNAAEAAADANwIEAAAAAN8IEAAAAAB8I0AAAAAA8I11IAAgxhXsLrA3bp1lQVMhNdnaXtrUjmpXL9pFAQAUQYAAgBiVlPTHz9y0PBt2yQcWNKk7Klrm/HV2f7te0S4KAKAIAgQAxKjqVaraNttuBclmOel5FjQ5lmezj1lm90e7IACAYggQABCj7mh0nj2/dorl5AcvPGzKy3LBaUfKrmgXBQAQggABADGqc82Wbgui9nMGWXba9mgXAwAQBrMwAQAAAPCNAAEAAADANwIEAAAAAN8IEAAAAAB8YxB1BDNmzHBbdna2/ysKAAAAJDACRASdOnVyW2ZmpmVkZJTfXwUAEOhVuCWlckVr3aMJK3EDiDkECABAzAn6KtzeStxnzP+ZlbgBxBwCBAAg5gR9FW5hJW4AsYoAAQCIyVW4h655334v2GlBxErcAGIZAQIAEHO0AvdJFRtaenq6BRErcQOIZUzjCgAAAMA3AgQAAAAA3wgQAAAAAHwjQAAAAADwjQABAAAAwDcCBAAAAADfCBAAAAAAfCNAAAAAAPCNAAEAAACAAAEAAACg9NECAQAAAMA3AgQAAAAA3wgQAAAAAHyr6H9XAABQnnKqbLf2cwYF7qJX2lXRrq92hl3RtmO0iwIgDAIEAAAxWIHONrOCZLPstO0WNDr3UVmz7AojQACxiAABAECM0d13VaDzKu6yILa6KDgF8dyBeEGAAAAgxqjrTlDvvqvLVhBbXYB4EqgA8cYbb9iECRMsOzvbGjRoYEOHDrWqVatGu1gAAABA3AhMgHj77bfts88+sxdeeMEOPfRQ++6776xixcCcPgAAAFAqAlGDzs/Pt9dee82GDRtmderUcY81atQo2sUCAAAA4k5MBojc3FwbM2aMLV261JYtW2bbtm2zgQMHWteuXffYd8eOHTZq1CibPn2620/BoE+fPtamTZvCfTZt2mR5eXn24Ycf2rhx46xatWp22WWX2XnnnVfOZwYAAADEt5hcSC4rK8tGjx5tq1evtsaNG0fc99FHH3WhoHPnztavXz9LTk62/v3726JFi4oFCI17+Omnn9y+Dz/8sI0YMcK+/vrrcjgbAAAAIHHEZICoVauWTZw40caPH28333zzXvdTC8XMmTPthhtusL59+1r37t3t2Weftbp169rw4cML96tUqZL7ec0117h/q5XizDPPtPnz55fL+QAAAACJIiYDRGpqqgsR+/LRRx9ZhQoVXHDwKCB069bNlixZYhs2bHCP/elPf7KUlBRLSkoq3K/ovwEAAADEcYDwa+XKlW461rS0tGKPN2/e3P1ctWqV+1mlShU77bTT7NVXX3VjJn744QebNWuWnXTSSWGPu3nzZsvMzCzc1JUKAAAAQIwOovZry5YtYVsqvMcUBDx//etf7fHHH3cDp6tXr27XX3+9tWzZMuxxJ02a5MZghBvcrYHa0RTt9weA8qLvXAQb/89D0ORG+XsvPT098QOEZlZS16RwXaC854tekMGDB/s6rrpEtW/fvvB3tUDotVp0zu+FLSvRfn8AKE985wUbf38EUXoc1PXiOkBovMPOnTv3eFzdlLzn90ft2rXdBgAAACCBxkCoq5K6MYXyHiMEAAAAAKUrrlsgtEbEV199ZTk5OcUGUmt6V+/5AzFjxgy3aQ0JAAAAAHHeAtGxY0fLz893g56Ldl+aPHmytWjRwurUqXNAx+/UqZM99thjdtttt5VCaQEAAID4F7MtEBMmTHB3/r3uSPPmzbONGze6f1988cVWrVo1FxJOP/10t6r01q1brX79+jZ16lRbv3693XPPPVE+AwAAACDxxGyAGDt2rAsCnjlz5rhNunTp4gKEDBo0yLU0TJs2zQWOhg0buulaW7VqFbWyAwAAAIkqZgPEuHHjfO2nmZb69u3rNgAAAABlK67HQAAAAAAoXzHbAhELmIUJAAAAKI4AsY9ZmLRlZmZaRkZGpF0BAACAQKALEwAAAADfCBAAAAAAfCNAAAAAAPCNMRARMIgaAAAAKI4AEQGDqAEAAIDi6MIEAAAAwDcCBAAAAADfCBAAAAAAfCNAAAAAAPCNAAEAAADAN2ZhioBpXAEAAIDiCBARMI0rAAAAUBxdmAAAAACUfQvEt99+a8uXL7eOHTtaWlqaeywvL8+GDRtm8+bNs0qVKtlll11m559//v6+BQAACKicKtut/ZxB0S5GVFTaVdGur3aGXdG2Y7SLApRugHj11Vdt8eLFds455xQ+NmLECJs0aZJVqVLFsrKy7JlnnrHDDjvM2rRps79vAwAAAlZ5zjazgmSz7LTtFkQ6/1FZs+wKI0AgwQLEsmXL7Pjjj7ekpCT3+65du2zKlCnWvHlze+6552zbtm3Wp08fe+uttwgQAADAF915V+U5r+KuwLa8KDwF9fyR4AFCLQyHHnpo4e/qzpSTk+O6LKn7krb27dvb/PnzS6usAAAgwanbTpDvvKvbVlBbXhCAQdQVKlSwnTt3Fv6+cOFC1xqhVglP9erVXdAAAAAAEPAWiLp169pXX31V+Pvs2bOtXr167nHPpk2bXIiIV6wDAQAAAJRSgOjSpYsNHz7cbrzxRktJSXGzMl111VXF9vnuu++sQYMGFq9YBwIAAAAopS5MF110kZvCNTMz083G1K5dO7vyyisLn//+++9t1apVdsIJJ+zvWwAAAABIlBaI1NRUe+ihh9zAaY19qFq1arHna9asaaNGjSrWpQkAAABAQFsgNGh6w4YNbhG50PAgNWrUsPT0dNcKAQAAACDgAeKOO+5w6z5EMm3aNLcfAAAAgIAHiIKCAl/7eAvNAQAAAAhwgPBjzZo1rosTAAAAgAAOon7ssceK/T537lxbv379Hvvl5+fbxo0bbdGiRW52JgAAAAABDBBFxzyoa5IGSO9tkLSeb9asmd16660HXkoAAAAA8Rcgxo4dWzi24bLLLrMePXrYJZdcssd+ycnJbgamKlWqWDxjJWoAAADgAAJE0TUdBgwYYH/+858Tep0HVqIGAAAASmkhua5du+7vSwEAAAAELUB4li5dasuXL7fs7GzbvXt32LEQvXv3PtC3AQAAABDPAeK3336zQYMG2TfffBNxTQgCBAAAAJA49jtADBs2zBYvXmytWrWys88+2w499FCrUKFC6ZYOAAAAQGIEiE8//dSaN29uzz77LKtNAwAAAAGx3ytR5+XlWcuWLQkPAAAAQIDsd4Bo3Lhx2FWoAQAAACSu/Q4Q11xzjc2bN8+WLFlSuiUCAAAAkHhjIH755Rc76aSTrF+/fta5c2dr0qSJpaWlhd1Xg6wBAAAABDhAPProo278g6ZwnTJlitv0e1F6To8RIAAAAICAB4gBAwaUbkkAAAAAJG6A6Nq1qyW6GTNmuE2rbAMAAAA4gAARBJ06dXJbZmamZWRkRLs4AAAAQPwGiA0bNvjet06dOvv7NgAAAAASIUD07NnT1yJy2mf27Nn7+zYAAAAAEiFAnHXWWWEDhMYLfPvtt7Zu3Tpr1aqV1a1b90DLCAAAACDeA8SgQYP2+pymbx0zZoy9+eabds899+zvWwAAAABIlJWoI1HLRK9eveyoo46yF154oSzeAgAAAECiBAhP06ZN7csvvyzLtwAAAACQKAHi559/tvz8/LJ8CwAAAADxvA7E7t27bdOmTTZ16lSbN2+enXDCCaX9FgAAAADiLUCcdtppEadx1UDq9PR0u+WWW/b3LQAAAAAkSoBo2bJl2AChxxQcmjVrZuecc47VrFnzQMsIAAAAIN4DxJAhQ0q3JAAAAACCPYgaAAAAQGIplUHUixcvtpUrV1pubq5VrVrVmjRpYscee2xpHBoAAABAogQIBYfHHnvMTdfqDZz2xkU0aNDABgwYYMccc0zplBQAAABA/AaI77//3u6++27bvn27tW7d2o4//nirVauW/fLLL/bVV1/Z559/7p5/8cUX7cgjjyzdUgMAAACIrwAxevRo27lzpz3xxBPWrl27Ys9dccUV9tlnn9nAgQPdfg8++GBplBUAAABAvAaIhQsXWseOHfcIDx49rucXLFhg8WrGjBluy87OjnZRAAAAgPgOEDk5OVavXr2I++h57RevOnXq5LbMzEzLyMiIdnEAAACA+J3GVeMdlixZEnGfpUuXuv0AAAAABDxAtG/f3nVjeumllywvL6/Yc/r95ZdfdoOpTznllNIoJwAAAIB47sLUu3dv+/TTT+3111+3SZMmWfPmza1mzZr266+/2vLly23r1q122GGHuf0AAAAABDxAVK9e3YYPH+6maZ05c6bNnz+/8LnU1FTr2rWr3XTTTXbQQQeVVlkBAAAAxPNCcjVq1HCLxWm9h9WrVxeuRH3EEUdYxYqlssg1AAAAgBhS4lr+q6++6haPu+666wpDgn42atSocB+tDzFy5EirUqWKXXnllaVbYgAAAADxMYj6iy++cIOj1S0pUgtDSkqK20cDrL/88svSKCcAAACAeAsQ06ZNs/T0dLvooov2ue+FF17o9p0yZcqBlA8AAABAvAaIb775xk488UQ3SHpftE/r1q1t8eLFB1I+AAAAAPEaIDZv3uymZvVLK1Fv2bJlf8oFAAAAIN4DRHJysu3atcv3/tpXrwEAAACQGEpUu69Vq5Z9//33vvfXvrVr196fcgEAAACI9wBx3HHHuVmV1q1bt899tY/2bdmy5YGUDwAAAEC8BgjNrKRuSffff79t3bp1r/tlZWXZAw88YPn5+Xb++eeXRjkBAAAAxNtCck2bNrUePXrY+PHj7eqrr3bh4Pjjj7dDDjmkcJD1ggUL7N1333UBo2fPnu41AAAAAAK6EvUtt9zipmh988037bXXXnNbUQUFBW7gtFag7tOnT2mWFQAAAEC8BYikpCS74YYbrFu3bjZ58mS3NsQvv/zinjv44IPt2GOPta5du1r9+vXLorwAAAAA4ilAeBQQMjIySrc0AAAAsILdBfbGrbMCeSVSKle01j2a2FHt6kW7KCjtAAEAAIDSlZT0x8/ctDwbdskHgby8qTsq2hnzf7b72/WKdlGwFwQIAACAGFG9SlXbZtutINksJz3PgijH8mz2Mcvs/mgXBHsVmADRr18/W7p0qVWoUKFwTYsnn3wy2sUCAAAodEej8+z5tVMsJz+Y4WFTXpYLTztSdkW7KIggMAFC+vfvb126dIl2MQAAAMLqXLOl24Kq/ZxBlp22PdrFQGkuJAcAAAAg2GKyBSI3N9fGjBnjuhwtW7bMtm3bZgMHDnTTw4basWOHjRo1yqZPn+72a9SokVt/ok2bNnvsO3ToULc1adLErWehfQEAAADEeQtEVlaWjR492lavXm2NGzeOuO+jjz5q48aNs86dO7txDlrETl2VFi1aVGy/m266ycaOHWtvvfWWtW7d2v72t7+5oAIAAAAgzgNErVq1bOLEiTZ+/Hi7+eab97qfWihmzpzpFrbr27evde/e3Z599lmrW7euDR8+vNi+LVq0sKpVq1qlSpXs8ssvd/9esmRJOZwNAAAAkDhiMkCkpqa6ELEvH330kZtVScHBo4CgVbIVDjZs2BBxRe2CgoJSKzMAAAAQBDE5BsKvlStXWoMGDSwtLa3Y482bN3c/V61aZXXq1HFjI5YvX24tW7Z0wUGtG3pMrRLhbN682bZs2VL4u7pSAQAAAIjzAKFKfriWCu8xBQHJz8+3ESNG2I8//mgVK1Z04yoef/xxq1atWtjjTpo0yY3BCKUxEwoe0RTt9weA8sI4NSDYgljnyY3y+Nz09PTEDxB5eXmWkpIStguU97zUqFHDRo4c6fu46hLVvn37Yi0QgwcPduMm/F7YshLt9weA8sR3HhBcQf3vPz0OzjuuA4TGO+zcuTPs1K7e8/ujdu3abgMAAAAQB4Oo/VJXpaJjFTzeY4QAAAAAoHTFdYDQWIY1a9ZYTk7OHtO7es8DAAAAKD1xHSA6duzoBkhr0HPR7kuTJ092MyxpBqYDMWPGDBswYIBbvRoAAABADI+BmDBhgmVnZxd2R5o3b55t3LjR/fviiy92MygpJJx++uluhqWtW7da/fr1berUqbZ+/Xq75557DrgMnTp1cltmZqZlZGQc8PEAAACAeBezAWLs2LEuCHjmzJnjNunSpUvhFKyDBg1yLQ3Tpk1zgaNhw4ZuitZWrVpFrewAAABAoorZADFu3Dhf+2mmpb59+7oNAAAAQNmK6zEQAAAAAMpXzLZAxAINotamrlEAAAAACBARMYgaAAAAKI4uTAAAAAB8I0AAAAAA8I0AAQAAAMA3AgQAAAAA35iFKQJmYQIAAACKI0BEwCxMAAAAQHF0YQIAAADgGwECAAAAgG8ECAAAAAC+ESAAAAAA+EaAAAAAAOAbszBFwDSuAAAAQHEEiAiYxhUAAAAoji5MAAAAAHwjQAAAAADwjQABAAAAwDcCBAAAAADfCBAAAAAAfCNAAAAAAPCNaVwjYB0IAAAAoDgCRASsAwEAAAAURxcmAAAAAL4RIAAAAAD4RoAAAAAA4BsBAgAAAIBvBAgAAAAAvhEgAAAAAPhGgAAAAADgGwECAAAAgG8sJBcBK1EDAAAAxREgImAlagAAAKA4ujABAAAA8I0AAQAAAMA3AgQAAAAA3wgQAAAAAHwjQAAAAADwjQABAAAAwDcCBAAAAADfCBAAAAAAfCNAAAAAAPCNAAEAAADAt4r+dw2eGTNmuC07OzvaRQEAAABiAgEigk6dOrktMzPTMjIyyu+vAgAAAMQoujABAAAA8I0AAQAAAMA3AgQAAAAA3wgQAAAAAHwjQAAAAADwjQABAAAAwDcCBAAAAADfCBAAAAAAfCNAAAAAAPCNAAEAAADANwIEAAAAAN8IEAAAAAB8I0AAAAAA8I0AAQAAAMA3AgQAAAAA3wgQAAAAAHyryLXauxkzZrgtOzubywQAAAAQICLr1KmT2zIzMy0jI4MPDAAAAAKPLkwAAAAAfCNAAAAAAPCNAAEAAADANwIEAAAAAN8IEAAAAAB8I0AAAAAA8I0AAQAAAMA3AgQAAAAA3wgQAAAAAHwjQAAAAADwjQABAAAAwDcCBAAAAADfCBAAAAAAfCNAAAAAAPCNAAEAAADANwIEAAAAAN8IEAAAAAB8I0AAAAAA8I0AAQAAAMC3iv53BQAAAMpewe4Ce+PWWYG71BVSk63tpU3tqHb1LJYRIAAAABATkpL++JmblmfDLvnAgiZ1R0XLnL/O7m/Xy2JZ4ALEN998Y7fccotdd9111rt372gXBwAAAP9f9SpVbZttt4Jks5z0vMBdlxzLs9nHLLP7LbYFKkDs3r3bhg0bZs2aNYt2UQAAABDijkbn2fNrp1hOfvDCw6a8LBecdqTsslgXqADx7rvvWvPmzS0nJyfaRQEAAECIzjVbui2I2s8ZZNlp2y0exOQsTLm5ufbyyy/b3Xffbd26dbMOHTrYlClTwu67Y8cOGz58uF144YXWqVMnu/HGG+3zzz/fY7+srCwbP36867oEAAAAIIEChCr7o0ePttWrV1vjxo0j7vvoo4/auHHjrHPnztavXz9LTk62/v3726JFi4rtN3LkSOvRo4elp6eXcekBAACAxBWTAaJWrVo2ceJE12Jw880373W/pUuX2syZM+2GG26wvn37Wvfu3e3ZZ5+1unXrulYJz4oVK2z58uV27rnnltMZAAAAAIkpJsdApKamuhCxLx999JFVqFDBBQdPpUqVXLenESNG2IYNG6xOnTq2cOFC++mnn+ziiy92+2RnZ7vXrV271gYOHFim5wIAAAAkkpgMEH6tXLnSGjRoYGlpacUe10BpWbVqlQsQChhnnnlm4fNDhgyxevXq2RVXXBH2uJs3b7YtW7YU/q6uVAAAAADiPECokh+upcJ7TEFAKleu7LairRRVqlTZ63iISZMmuTEY4QZ3b9u2zaIp2u8PAOVF37kAEETbolTf8ztWOK4DRF5enqWkpITtAuU9H86gQYMiHlctFu3bty/WAjF48GCrWrVq1AdhR/v9AaA88Z0HIIjSY7y+F9cBQi0JO3fuDDu1q/f8/qhdu7bbAAAAAMTBLEx+qatS0bEKHu8xQgAAAABQuuI6QGiNiDVr1uyxsrSmd/WeBwAAAFB64roLU8eOHW3MmDFu0HOvXr0Kuy9NnjzZWrRo4WZgOhAzZsxw22+//Rbd2ZjWZlty5Tyz7bssMzMzOmUAgCgMotbYMwAIhLWxUd874ogjik0+FE5SQUFBgcWgCRMmuPUa1B3pnXfesQ4dOliTJk3cc1rPoVq1au7fDzzwgM2ZM8d69uxp9evXt6lTp9qyZcvsmWeesVatWpVKWaZPn+4GUQMAAACJbOTIkda0adP4DBAKBOvXrw/73NixY906Dt5MS6NGjXKVfAWOhg0bWp8+faxt27alVpatW7faf//7Xxdkbr/99gM61tChQ+22224r0Wu8WaDuvfdelwpR/vbn7xZPYv38olW+8nrfsnif0jpmaRyH7734FOvfC4l+fnzvRe+aBf177wgfLRAx24Vp3LhxvvbTTEt9+/Z1W1mpUaOGdenSxWbNmrXPRLYvajnZ32PoD3qg74/y/7vFg1g/v2iVr7zetyzep7SOWRrH4XsvPsX690Kinx/fe9G7ZnzvJfgg6vLWqVOnmDgGyl+i/91i/fyiVb7yet+yeJ/SOibfe8EV698LiX5+fO9F75rxvbdvMduFCf9HA2kyMjJ89UkDgETA9x6AoMmMo/oeLRBxst7FNddc434CQBDwvQcgaGrFUX2PFggAAAAAvtECAQAAAMA3AgQAAAAA3wgQCUCrbz/22GN2ySWX2Nlnn2033XSTffPNN9EuFgCUqSeffNIuuOAC973Xu3dvmzdvHlccQCB88803dtppp9m//vWvqLw/YyASwO+//+4W1+vatasdcsghNnv2bHv22WfdY1WrVo128QCgTGjRJS0qmpqaasuWLbM777zTxowZY9WrV+eKA0hYu3fvduufaSLVv/zlL+4GSnmjBSIBVKlSxY3ar1OnjiUnJ9uZZ55pFStWtJ9++inaRQOAMqPFNRUeJCkpyXbu3GmbN2/migNIaO+++641b948qqtVx+xK1IksNzfX3SVbunSpu2u2bds2GzhwoGtBCNc9adSoUTZ9+nS3X6NGjaxPnz7Wpk2bvR5fwUH71q9fv4zPBACi+7339NNP2+TJk91rTjrpJGvYsCF/EgAJ+72XlZVl48ePt+HDh9vQoUMtWmiBiAL98UePHu2a3xs3bhxx30cffdTGjRtnnTt3tn79+rkWhv79+9uiRYvC7p+Xl2eDBw+2K664wi3FDgCJ/L2nbkvTpk2zZ555xv2PVi0RAJCo33sjR460Hj16WHp6ukUTASIKtEDIxIkTXYK8+eab97qfEuvMmTPthhtucH3dunfv7sY21K1b1yXPULt27bL777/ftTyoSxMAJPr3nlSoUMFOPPFEW7BggX366adleBYA4F9pf++tWLHCli9fbueee27U/wwEiChQn10/qwx+9NFH7n+M+iB5KlWqZN26dbMlS5bYhg0big2oUcuD7r4NGjSIu3AAEv57L1R+fr79/PPPpVZmAIil772FCxe6buoXX3yxm4Fu1qxZ9sYbb7jWi/LGGIgYtnLlSmvQoIGlpaUVe1wDZ2TVqlVu4LQ89dRTtmXLFvdTA6gBIJG/97Kzs11rQ/v27d3/pOfOnWtfffWVu4MHAIn4vde9e3c3UY5nyJAhbiY6dVsvb9Q0Y5gCQbjk6j3mzTayfv16e++999z/RIum1yeeeMJatmxZjiUGgPL53lNrq773NPZBUxmq6+Z9991nTZo04U8AICG/9ypXruy2oq0UmokzGuMhCBAxTAOiU1JS9njcm7ZQz4v6yM2ZM6fcywcA0fre05265557jj8AgMB874VSl/VoYQxEDFOy1Lzm4ab68p4HgETC9x6AoKkUh/U9AkQMU9OVmrVCeY/Vrl07CqUCgLLD9x6AoKkVh/U9AkQM05zBa9assZycnD2m+/KeB4BEwvcegKBpHIf1PQJEDOvYsaOblnDSpEnFmrO06mqLFi0KZ2ACgETB9x6AoOkYh/U9BlFHyYQJE9w0hF7z1Lx582zjxo3u35rfV6tI60Nz+umn24gRI2zr1q1ulpGpU6e6WZfuueeeaBUdAPYL33sAgmZCgtb3kgo0/x3KXc+ePd0HI5yxY8e6eX29kfejRo2y6dOnuw9gw4YNrU+fPta2bdtyLjEAHBi+9wAETc8Ere8RIAAAAAD4xhgIAAAAAL4RIAAAAAD4RoAAAAAA4BsBAgAAAIBvBAgAAAAAvhEgAAAAAPhGgAAAAADgGwECAAAAgG8ECAAAAAC+ESAAAAAA+EaAAAAckClTpliHDh3cT5Senj17uuvqbatXry587quvvnKPvfzyyzFxyT/77LNiZe3Xr1+0iwSgDBEgACSsdevWFVZoLrjgAtu1a1fY/X744YfC/VRpA2JFtWrV7JprrnFb9erVy/S9fvzxR/ffwJVXXrnPfUeOHOn2fe2119zv9evXLywngMRHgACQ8CpUqGC//PKLzZ8/P+zz77//viUnJ7sNiLUAcd1117mtRo0aZfpehx9+uB133HEuSCxevHiv++3evdumTp3q/rvq2rWre6xBgwaF5QSQ+Pi/JYCEd8wxx7iK2OTJk/d4Tq0SH3zwgZ144olWsWLFqJQPiBXdunUrDNV789///tc2bdpkbdu2tdq1a5dj6QDECv5vCSDhVapUyc444wxXKfr111+tZs2ahc99+umnrnXitttus6+//jrs6wsKClz40Ou/++47y8/PtyOPPNJ1i/IqXJ7NmzfbpEmTXCVr7dq1lpOTY7Vq1bKTTjrJrr322mLvLdnZ2TZ27Fj78MMPbePGjZaUlOT2Uei5/vrrrW7dum6/f/zjH+6ur/atV69esWOoH/zo0aPtueees+OPP76wj/ztt9/uupS0adPGXnnlFVu+fLl7vzlz5pT4vOS3336zESNG2Ny5cy03N9eOOuooX91dwtHf4fXXX7dPPvnEnXfVqlWtZcuW7g52w4YNi+3rdSvTOarrjK6VyvKnP/3JnV/Hjh33OP7OnTvt7bfftunTp9tPP/3krmuTJk3ssssus1NOOaXYvt61HTNmjLs2uh7625155pk2aNAgt8/ChQvtpZdeshUrVlhqaqoLnH379rVHHnnEPeddU5VP3XoefPBB95kLpWM//vjjdsMNN+z3tYtEf9+BAwfaokWL3DiEiy++2D2uv5fOT9dO56ZzaN68ufXu3du1Onh0LfU5mj17tvv8VKlSZY/38IJ4uM8IgGCgBQJAIJxzzjmugjxt2rQ9KnQHHXSQnXrqqWFfp0r2//zP/7hK39atW61Tp0527rnn2u+//+4ee/7554vtrxCiSr5CgCqgF110kR122GH2zjvv2M033+wqeEWPfffdd9u//vUvV4bzzjvPbarozps3z9asWXPA5/3NN9+4iqAq0Dq2V6kt6Xlt377dVUgVjnQ+l1xyiavAq6KsSmlJ/Pzzz9anTx8bP368O5aukQKWQpeu0dKlS8O2FN111132+eef22mnnWadO3d2FeEHHnjAva6oHTt2uOvqnYMqul26dLH169e7QDBhwoSw5Xr22WddqGnWrJk7Py/I6Ph//etfXQDTe+s6btiwwW699dZif0/Rc+oK995774V9Dz1etOtPaVJ4VRDW9bv//vsLw4PClq6rAlh6erqdf/75bvyCwpA+GwqEHgUGfW71OVCICJWVleU+m/p8/+Uvfyn1cwAQH2iBABAILVq0cHfMNVOQ7kLLli1b3OwxuuOuO7J7q/DNmDHDBRBVSr1uTrrDfd9997mwoMp306ZN3eMnnHCCTZw40d1RL0p3uHWnW3fFr776aveY7vqrsqfwojvZoZXgvQ36LokvvvjCBgwY4Mp/IOf1xhtvuPKqgvy3v/2t8DhnnXWWe31J6FzV6vPUU0+5bjAeXZeMjAx74oknXGU3tHKsiv2QIUMsJSXFPaYQoYr9uHHjih1HgUwtMLq7rhYNhSfvLvwdd9xhL7zwggsCod1vvv32Wxs1apTVqVOn8DGFTpVT/f51Z77o3XqdR2ggVYuRWnwUOjSIv2hr0ffff29Llixxf2+1SpUmtbIoYCksKAC2bt26WDDSe/fv39+FxKKtQLreTz75pLt+aqnzAte7777rWhpCPzfq7qfPiP7udPkDgosWCACBoYqRKlLeHW5V6lVBDK0kFaUKv+7KqqJatMKkSqwqX6KKuEd3ZkPDg6jClZaWZgsWLNjjOa/iVpQCTbjjlNSf//znsOdX0vNSRVnPhQ6SVcVT3Xn80l1vtYroehSt9ItaNFTBVVDRFkp3/L3wIHpfVdjVMuBRRV+tPZoVqGh4EF1PhQpVgD/66KM9jt+rV69i4UE0mFgtF7rbXjQ8iFpR1JoQSnf41cITOo7Aa5UoWokvDcuWLbNbbrnFtRIp5BQND2pdUkuCgm3o++qzqnPWPkU/l17YVjcoBZOivKl6I/03AyDx0QIBIDDUjeXFF190FTtVklQZUnchbeGoQqaKrO5U//vf/97jeYUP0aw1Ralyqq4+qiyri4u3n3cn3XPEEUdYo0aNXEVd4wB0Z7pVq1auPKU1I5Tu2h/oeWkch+6ma3xEuDvnqliHC0bheOFNd7/DrWHgvad+Fh0LoUHw6u4U6pBDDnF39Yu+ftu2be7cNO4jlCrLRd+nKI0JCLVq1arCcwylsHHooYe6a1PUySef7Mqlz5fGvShkKLRoPIb2b9eunZUWVfLVWqQZmtRSohBWlMKV/p56/3DX2+smpzUminZJUtgeNmyYa4W48cYb3WOZmZm2cuVKO/roo91nAUBwESAABIYqWaokzZo1y04//XRXiVQf8L1RRVR3kjXjTGiXmtAKuUcDVdVFRu+lriyqSHotDG+99ZaryHl051/dS1TRVejw+uzrtRoXcNVVV4W9w10SBx988AGflwKEhA4Aj/Qee6MuNt7gdW17oz74RSlAhKPro1aHoucmamnS5udv5gl3ft65720KVZ17aIBQmVQB17VVFzl95jTIWuMHNOi7NKcLVoVe10qftXABy7veakmJNDVr6PVQ2P7nP//pWp68lhYGTwPwECAABIoqdqrMPfroo66bkCpKe6MuR6JxAJpdZ180ZuHVV191d+l1t7dohVQV9jfffHOP12hxMPXLV5DRXeAvv/zSdS/S6xUwvJl6vEpn0daM0EquXyU9L29/tRqEo/EMJX1vna83yLc0ed2+NMZBg8RLomh3p9Dyei0Xfs9d3YU0G5PGEihAeGuNlPbMRQqaatXS8R9++GE3fqVolzSv/Jdeeqnr5uSXApNmq1L3J4UgdYtSS5m6vYWbXQpAsDAGAkCgqN+9WgV0911dhjQrTaTKqLoZqWLv3dmORHeY1WVJXTxC72arK0leXl7Eyqu6hahC+L//+7/uMc12E3oHvmgXqKJ3oUuipOelSqgGA2v2JA08D9eNxi+vm1DRbkelSeel8qq7TWkMQm/cuLH7Ge7uvbqdaQtHXZU0s5QWL9Rr1cVLrQShYywOlD43Ghytwe2q7Cs0FT1vdWHTPvtzvYuuCaGZmvRZUctdaYzNARDfCBAAAkVdMTR7jjbNxb8vms5T3Ts0U01otxrRVKJeFxaFBnVX0tiHol1CVPHS4NZQel1o95eid/qLzgzlVby9QaweTaGqdQhKqiTnJRr0HK4fvWYb8jv+QTT2RNvMmTPdFkrdkfbnfDy6+65BzBr4rC5h4UKExn/srTUl1LHHHusq/VqvQoO/i9KMTeFahDzdu3d3z2uqWbVAqZJfFhQQNBOW3k8hQi0R3nmrNUyVfpVdLWAqR7hxKeG6dKnVQUFIXc0005Ww9gMAoQsTgMDRXdlwg4vDUaVMd281Y5PuJKtSpUqZKqAaQ+HNua879OqioilhNahVg2fbt2/vuhepC4gqoaHThmqA7r333uvCgVof1J9eLSMff/yxO1aPHj0K91V3Es0spAChu94aaO11efLudJdESc5LNFuPxmmoS47GFmjRN5VDFVYNGo40niGUjqtuWw899JAbF6JzUfDS8VTRVUtO0RmgSkqzLynEab0HXReVVV1y1Hqj8KDrPnz48L2O6QgNnJoeVYuzqczqvqPrpJCj46mFQtO/hqPB0polSmFGf9uyXDdBIULl1OdGs1ApKCi4KFDdeeedbjYlnbPGNKiFTC1a+qypZUwDqTX1cOXKlYsdU8fSbEsay6GZng4//HAXqACAAAEA+6iYafExVdI1DafuROuOvSqfDRo0cKsRF53GVDPWaFE4VfRVkfMWlFOg0ADaojQG4fLLL3eVUVXA1f1JFU0dTxV2VfQ8qmA//fTTbmYc3fFXBV938ocOHerKVNIAUdLzUt93vZcG1qo7iyromupTC8kpJJUkQGiwr+7eK2jpWLpWqqyqYq7KfriVpUtCLTdqWVHXG1WYFXzUeqJzU1BTC0XoateR6BqpW5laXxSY9LfQtdG5q/uQN84glM5JY2w0LkYLx5X1ugn6m2paXv1UIFD5tOnzqIH9GlujCQQUztTSo8+aApCmttVYnHBUbq2roUBC6wMAT1JBuPZMAAAQkRam88KIglU499xzjwt3mi5Xwawkevbs6X563YfihVa51nTEWvQPQGJiDAQAABGoZUZhoSiNbdBdfQ2M12D8cH744QcXHtQ9rKThwaPuT6qQa1OXtVilbnpeOQEkProwAQAQgcYIaBVsb60FhQnNPKWAoG5codPRfvDBB27MgcaXSGjXtZIMdFe3Ns/euhnFAo3PKXqe3tgZAImJLkwAAESgNSA0AFljVTTIXK0Pmp1ILQ9a7C90KuB+/fq5gKGB8xpfoIHIAJBICBAAAAAAfGMMBAAAAAACBAAAAIDSRwsEAAAAAN8IEAAAAAB8I0AAAAAA8I0AAQAAAMA3AgQAAAAA3wgQAAAAAMyv/wcx7gKtKxrkHAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 800x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "mock_hist = Histogram.open(binned_output_file)\n",
    "background_hist = Histogram.open(binned_bkg_output_file)\n",
    "\n",
    "mock_energy_counts = mock_hist.project(\"Em\")\n",
    "background_energy_counts = background_hist.project(\"Em\")\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8, 4))\n",
    "\n",
    "mock_energy_counts.draw(\n",
    "    ax=ax,\n",
    "    label=\"Pointing-cut mock data\",\n",
    "    linewidth=2,\n",
    ")\n",
    "background_energy_counts.draw(\n",
    "    ax=ax,\n",
    "    label=\"Pointing-cut background\",\n",
    "    linewidth=2,\n",
    ")\n",
    "\n",
    "ax.set_xscale(\"log\")\n",
    "ax.set_yscale(\"log\")\n",
    "ax.set_xlabel(\"Measured energy [keV]\")\n",
    "ax.set_ylabel(\"Counts\")\n",
    "ax.legend()\n",
    "ax.grid(alpha=0.3)\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2b62264",
   "metadata": {},
   "source": [
    "Save the pointing cut orientation file if needed. The time cut can be performed and used on the orientation file in analysis directy as shown above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bdcb6fc4",
   "metadata": {},
   "outputs": [],
   "source": [
    "orientation_output_file = \"DC4_final_530km_3_month_with_slew_15sbins_GalacticEarth_SAA_NGC4151_cut\"\n",
    "pointing_cut_orientation.write_fits(orientation_output_file, overwrite=True)\n",
    "print(f\"Wrote: {orientation_output_file}\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "cosipy312",
   "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}