examples/Example8_Galactic_Center_nonPoissonian.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example 8: Application of NPTFit to the Galactic Center Excess"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It was found in Example 3 that a non-zero value of the GCE template is preferred by a fit in the galactic center. In this example we will test the point source interpretation of this excess by including, in addition to the Poissonian templates considered, non-Poissonian point source templates of various morphologies.\n",
    "\n",
    "Here we simply perform the run, a detailed analysis of the results can be found in the next example. A Python script version of this notebook can be found as `Example8_Galactic_Center_nonPoissonian.py`, which can be run faster on multiple processors with MPI (see example in `Example8_Galactic_Center_Batch.batch`.\n",
    "\n",
    "**NB:** Even with `nlive=100`, this notebook can take up to an hour to complete. This highlights that for realistic non-Poissonian runs, running on multiple cores becomes necessary. We show an explicit application of this in Example 10.\n",
    "\n",
    "**NB:** This example makes use of the Fermi Data, which needs to already be installed. See Example 1 for details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import relevant modules\n",
    "\n",
    "%matplotlib inline\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "from NPTFit import nptfit # module for performing scan\n",
    "from NPTFit import create_mask as cm # module for creating the mask\n",
    "from NPTFit import dnds_analysis # module for analysing the output\n",
    "from NPTFit import psf_correction as pc # module for determining the PSF correction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 1: Set up the Scan"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first need to \n",
    "1. Set up an instance of `NPTF` from `npfit.py`\n",
    "2. Load in the data and exposure maps\n",
    "3. Set up and load the mask used for the scan\n",
    "4. Load in the spatial templates\n",
    "\n",
    "These are done identically to Example 3, and we refer to that notebook for details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = nptfit.NPTF(tag='GCE_Example')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "fermi_data = np.load('fermi_data/fermidata_counts.npy').astype(np.int32)\n",
    "fermi_exposure = np.load('fermi_data/fermidata_exposure.npy')\n",
    "n.load_data(fermi_data, fermi_exposure)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "pscmask=np.array(np.load('fermi_data/fermidata_pscmask.npy'), dtype=bool)\n",
    "analysis_mask = cm.make_mask_total(band_mask = True, band_mask_range = 2,\n",
    "                                   mask_ring = True, inner = 0, outer = 30,\n",
    "                                   custom_mask = pscmask)\n",
    "n.load_mask(analysis_mask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "dif = np.load('fermi_data/template_dif.npy')\n",
    "iso = np.load('fermi_data/template_iso.npy')\n",
    "bub = np.load('fermi_data/template_bub.npy')\n",
    "gce = np.load('fermi_data/template_gce.npy')\n",
    "dsk = np.load('fermi_data/template_dsk.npy')\n",
    "\n",
    "n.add_template(dif, 'dif')\n",
    "n.add_template(iso, 'iso')\n",
    "n.add_template(bub, 'bub')\n",
    "n.add_template(gce, 'gce')\n",
    "n.add_template(dsk, 'dsk')\n",
    "\n",
    "# Remove the exposure correction for PS templates\n",
    "rescale = fermi_exposure/np.mean(fermi_exposure)\n",
    "n.add_template(gce/rescale, 'gce_np', units='PS')\n",
    "n.add_template(dsk/rescale, 'dsk_np', units='PS')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 2: Add Models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "n.add_poiss_model('dif', '$A_\\mathrm{dif}$', fixed=True, fixed_norm=12.85)\n",
    "n.add_poiss_model('iso', '$A_\\mathrm{iso}$', [0,2], False)\n",
    "n.add_poiss_model('gce', '$A_\\mathrm{gce}$', [0,2], False)\n",
    "n.add_poiss_model('bub', '$A_\\mathrm{bub}$', [0,2], False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This time we add a non-Poissonian template correlated with the Galactic Center Excess and also one spatially distributed as a thin disk. The latter is designed to account for the unresolved point sources attributed to the disk of the Milky Way (known sources in the 3FGL are masked)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "n.add_non_poiss_model('gce_np',\n",
    "                      ['$A_\\mathrm{gce}^\\mathrm{ps}$','$n_1^\\mathrm{gce}$','$n_2^\\mathrm{gce}$','$S_b^{(1), \\mathrm{gce}}$'],\n",
    "                      [[-6,1],[2.05,30],[-2,1.95],[0.05,40]],\n",
    "                      [True,False,False,False])\n",
    "n.add_non_poiss_model('dsk_np',\n",
    "                      ['$A_\\mathrm{dsk}^\\mathrm{ps}$','$n_1^\\mathrm{dsk}$','$n_2^\\mathrm{dsk}$','$S_b^{(1), \\mathrm{dsk}}$'],\n",
    "                      [[-6,1],[2.05,30],[-2,1.95],[0.05,40]],\n",
    "                      [True,False,False,False])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 3: Configure Scan with PSF correction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loading the psf correction from: /zfs/nrodd/NPTFRemakeExamples/psf_dir/gauss_128_0.181_10_50000_1000_0.01.npy\n"
     ]
    }
   ],
   "source": [
    "pc_inst = pc.PSFCorrection(psf_sigma_deg=0.1812)\n",
    "f_ary, df_rho_div_f_ary = pc_inst.f_ary, pc_inst.df_rho_div_f_ary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The number of parameters to be fit is 11\n"
     ]
    }
   ],
   "source": [
    "n.configure_for_scan(f_ary, df_rho_div_f_ary, nexp=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Step 4: Perform the Scan"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As noted above, we take a small value of `nlive` simply to ensure the run finishes in a reasonable time on a single core."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "n.perform_scan(nlive=100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This can take **up to an hour to run**. The output of this run will be analyzed in detail in the next example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"https://imgs.xkcd.com/comics/compiling.png\"/>"
      ],
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import Image\n",
    "Image(url = \"https://imgs.xkcd.com/comics/compiling.png\")"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}