hist_FGL.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from tqdm import *\n",
    "from astropy.io import fits\n",
    "from astropy.table import Table\n",
    "import healpy as hp\n",
    "\n",
    "from NPTFit import create_mask as cm \n",
    "from pFGL import plot_FGL\n",
    "\n",
    "%matplotlib inline\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Analysis mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "nside = 128\n",
    "mask = cm.make_mask_total(nside=nside, band_mask = True, band_mask_range = 2, mask_ring = True, inner = 0, outer = 30)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### [3/4]FGL SCDs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4feb64e96a824708add7b3c95471d55b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=228), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "# 3FGL \n",
    "\n",
    "p3FGL = plot_FGL(catalog='3FGL')\n",
    "\n",
    "x_counts_3FGL, y_counts_3FGL, error_L_3FGL, error_H_3FGL, x_errors_L_3FGL, x_errors_H_3FGL = p3FGL.return_counts(\n",
    "                flux_min = 1e-13,\n",
    "                flux_max = 1e-6,\n",
    "                flux_bins = 17,\n",
    "                mask = mask);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "81c4cb9001a34735ba0958a8f0ce9d86",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=532), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "# 4FGL \n",
    "\n",
    "p4FGL = plot_FGL(catalog='4FGL')\n",
    "\n",
    "x_counts_4FGL, y_counts_4FGL, error_L_4FGL, error_H_4FGL, x_errors_L_4FGL, x_errors_H_4FGL = p4FGL.return_counts(\n",
    "                flux_min = 1e-13,\n",
    "                flux_max = 1e-6,\n",
    "                flux_bins = 17,\n",
    "                mask = mask);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Source classes without extragalactic association"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remove sources with extragalactic associations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Source classes with extragalactic association, to remove\n",
    "eg_assoc = ['CSS', 'BLL', 'FSRQ', 'AGN', 'NLSY1', 'RDG', 'SEY', 'BCU', 'GAL', 'SBG', 'SSRQ','css', 'bll', 'fsrq', 'agn', 'nlsy1', 'rdg', 'sey', 'bcu', 'gal', 'sbg', 'ssrq']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "p4FGL_no_eg = plot_FGL(catalog='4FGL')\n",
    "p4FGL_no_eg.cat = p4FGL_no_eg.cat[[src['CLASS'].strip() in eg_assoc for src in p4FGL_no_eg.cat]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "cbbeac8953b042c3aff0457dd30e88de",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=133), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "x_counts_4FGL_noeg, y_counts_4FGL_noeg, error_L_4FGL_noeg, error_H_4FGL_noeg, x_errors_L_4FGL_noeg, x_errors_H_4FGL_noeg = p4FGL_no_eg.return_counts(\n",
    "                flux_min = 1e-13,\n",
    "                flux_max = 1e-6,\n",
    "                flux_bins = 17,\n",
    "                mask = mask);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8VOXZ//HPBbIIKMqiacTIokYR0UJwA7X4SEER60JRf1FRa1O0tvIrVtpHfRYtrYJtfdq60cXWn9FKEX3EHUsKioosIorYqmxCBBFBkQiG5Pr9MWdCgEkyy5nMTOb7fr3yypkzM/d9hUPmyr2c+zZ3R0REJNu0ynQAIiIisShBiYhIVlKCEhGRrKQEJSIiWUkJSkREspISlIiIZCUlKBERyUpKUCIikpX2yXQAqTKzocBXQIG7P5bpeEREJBxZ1YIyswH1jsea2RlmdmMTbzvP3ecB/cysTXojFBGR5pIVCcrM2pjZOcAjweMRQLW7zwa2m9npjbx9VfD9YKBLWgMVEZFmkxVdfO5eDTxlZhOCU0OBZ4PjVcCpwBwzGwacUu+tVcAfg5bXZuDjPcs2szKgDKBjx44DjzrqqLT8DCIissuiRYs+cffuqZSRFQkqhq7Al8FxFUHLyN1nAbPqvzAYg/oCeMZjrHzr7lPN7CNg1Ne+9rWBCxcuTGvgIiICZrY61TKyoosvBgdaB8etaSROd69w9wXBOFRDr5np7mWdO3cOOUwREUmXbE1Qm4COwXEn4NNUCjOzUWY29bPPPks5MBERaR7ZmqAqgN7BcW9gTiqFqQUlIpJ7siZBmdnFQJ9gUsNCoI2ZnQ184e4pJSi1oEREco9ly466wT1MO4mMOdXEmvCQqpKSEtckCRGR9DOzRe5ekkoZWdOCcvdqj9gZdnJSC0pEJPdkTYJKJ41BiYjknrxIUGpBiYjknrxIUGpBiYjknrxIUCIiknvyIkGpi09EJPfkRYJSF5+ISO7JiwQlIiK5RwlKRESyUl4kKI1BiYjknrxIUBqDEhHJPXmRoEREJPcoQYmISFZSghIRkaykBCUiIlkpLxKUZvGJiOSevEhQmsUnIpJ78iJBiYhI7lGCEhGRrKQEJSIiWUkJSkREspISlIiIZKV9Mh1AqszsJMCA3u5enul4REQkHBltQZnZgHrHY83sDDO7sYn3dDezu8zswODUGe7+KtDZzPZLZ7wiItJ8MpKgzKyNmZ0DPBI8HgFUu/tsYLuZnd7Qe919I7CZXbEfHHzfARwao64yM1toZgs3btwY5o8hIpJ31i9ezPrFi5ulrox08bl7NfCUmU0ITg0Fng2OVwGnAnPMbBhwSr23Vrn7lD2K+yr43hGojFHXVGAqQElJiYfyA4iISNplyxhUV+DL4LgK6ALg7rOAWfVfaGZdgJ7AQOAF4AkzGwxsdfctsQo3s1HAqMMPPzwtwYuISPiyJUE50Do4bk3jXY+bgauITIzA3ecF5+c1+A4REck52TLNfBORLjqATsCnDb3Qd6mNt3CtxSciknuyJUFVAL2D497AnDAL12rmIiK5J2MJyswuBvqYWRmwEGhjZmcDX7h7qAlKLSgRyUqVlZEvianJMSgzuwoYQmScqMGXBc+/7O4PxFn3Y8CjRMacatz9d3G+L2GaJCEiknvimSSx1N3/FE9hZlYSb8XBVHOAnfG+J1nuPhOYWVJS8t101yUi0pKt3x65n7SgGepqsovP3RfGW1girxUREWlMkwnKzM40s9+b2fHB47L0hxUuTZIQEck98UyS+D7wY+BSMzsDOD69IYVPkyRERHJPPAlqo7tvcfcbgG8Cg9Ick4hIftiwAS68ENavz3QkWSmeBPV09MDdfwI8mL5w0kNdfCKSle66C+bPh1tvzXQkWcnc41s/1cy6ufsnaY4nrUpKSnzhQs3jEJEEhX2vUu/esGPH3ufbtYMVK8Kpo7AwnHLqKS8vZ8r46/mfTzZxQ48ejL/9dkpLS2O+1swWuXvcM7tjSeRG3bimmouIZEQu3fT66qtw3nnQvn3k8b77wvnnw2uvZTauRpSXl1NWVkbZJ5sYAly5di1lZWWUl6dvn9hEFou1tEUhIpLNwm6NFBZCQUGkFdWuXeR7QQEcH94ctLD3bLrw0kup31a6Fri2qortl17K+qOPDrWuqERaUDm7l5LGoEQk62zYAJdfDjNnwrhxWT9RojdQDmwLHm8DHgJ6pbHOvGhBaSUJEck6M2bs6pIcNiz04gsGDAi1vLaHHcbnq1fTnsjmfe2Bz4F2hx0Wel1RibSgfpqWCEQSUF5ezqAePZhjxqBDD01r/7fkGE3ZTqtJkyZR2KoV9wEnAfcBh7RqxaRJk9JWZ9wJyt3fBjCzsfXPm1lnM/tvMxsYdnDSDD76CE4/PSd+qaODtFeuW9dsg7SSQzRlO61KS0v54sEHmXTwQbxlxpTDDuOLBx9scBZfGOKeZl73BrOfAiXABHdfZWa3Ao8A/dz9b2mIMTSaZh7DtdfC/ffD974H99wTatFhD9IeMHAg7WOc3w5sWbQolDrS1VUhewhztl2OTtkGdv07pKv8NFjyyvMAHH/K8EZfF8Y082S2fP8A+B1wBrCKyBjZ+8DBqQSSTi1mu410/lLfe2/kK8xf6pD1BqYA5xHZfnkb8DiRdbjeyGBcLVIufXC++mqk1fTcc7B9e2TK9ogR8B//kenImrRk1VsAHJ8L/84ZkEyCOhy4DNhpZkcQ+dzoAOwfZmBh0iSJGKK/1M8/D19+mZZf6pYwSJu3NmyItK4ffzwy/TlsYX4gN8OUbcmMZBLUb4D+wFKgL3A2kSnxb4YYl8SSrl/q9u1z4pd60qRJdLr8cu6rrWUqUAb0SPMgbd6qP54TctdvWkSnbJeWwhNPRMZWJeclk6DOBfoRWdX8YeBkd/9FqFFJ89iwIXL/RVkZTJ2a9b/UpaWllANTbrqJNWvWMKWoiEmTJqV1kDYn5GLXb9hdWmmesi2Zkcg086jP3P3fgYXuvoUcvoE3782YAXffDccdF/k+Y0amI2pSaWkpq1atora2llWrVmV9ciovL6dHjx60atWKnj17Zv+MwxxcgieX1XbvTm337pkOI2sl04IqMbOtQGczO45Id9/z4YYlkvui0+L3r6qiArho9WrKyiL7fYaaWNPQ9es7dkDbtlgOdP1Ky5VMgvofYCJwLLAYmBxqRCIZVBlid9nEiROpqqpiCjAE+A/g+1VVTJw4kaFDh4ZWT2HY3WUbNlA1ejRfXnAB3V56Keu7fqPWB/fyFWhGXIuRcIIKuvXqVpUws68BX4QZlEimtArxhuX31q1j33qPrw2+vly3js/CvDE6DeM5W6P3sF18cbhliyQgrgRlZjcC7WI81Qo4ncg9UQkzswHuvjg4Hgt8CJS4e4OtMjPrDtwE/Le7b97zcTJxiEQVhDilelBBAePXr9/rvq3/KShgQTqmboekvLycn9xwA+s2bKAohyai1Gbxv6kkJ94WVDsiO+k6kW3fXwF2Buc3JVqpmbUBhgO/BIrNbARQ7e6zzayfmZ3u7nNivdfdN5rZZoIJHns+FklJiK2R8XfeSdXYsbSvqam7b6uqdWvG33ln1t4AGx03q6qqAmB1usbNROIQV4Jy99uix2b2nru/U+9x10Qrdfdq4CkzmxCcGgo8GxyvAk4F5pjZMOCUem+tcvcpidRlZmVEbpmhqKgo0VBFklZaWsqau+6ifPlyfrVtGz/q1IkRRx1FUcgf9OkYN6uvKhfGzaRFSmaSxHFBUloDFAJHAv9IMY6uRBYHAKgCugC4+yxgVv0XmlkXoCcwEHhhz8d7FuzuU83sI2BU27ZttaCtNKuiBQs4c2slZwKF+2X/h3JDyS7MJCgSr2QS1O+AscHXm8CvQ4jDgdbBcWsa767bDFzFrv2p9ny8d+Fa6khasDBbI0VFRaxevTrmebV6pLklPG7j7rXu/oC7/wCYBiwMxo1Oaeq9jdhEZBwZoBPwaSP1R9XGehyLdtSVTJg8bzIVKyt2O1exsoLJ87L3zoxJkybRoUOH3c516NBBy0lJRqQ0scDdPwfOcfe33f2VFIqqILLoLMH3mBMkkuXuM929rHPnzmEWm/PKy8vp2bNnzqxyEP3Ar9xaSeXWSJdTNn/gDyocxJjpY5i3Zh4QiXXM9DEMKhyU4cgaVlpaytSpUznkkEMwMw477DCmTp2aExMk1m9bz/pt2b+vmcQvmS6+3bj7OqibmXdsdNp4U8zsYqBPMInhMeAIMzsb+KKhGXzJajHbbYQoF2drRT/w7zn7HgYXDa77wJ82elpodUQTXxiKuxVzz9n3UPZUGaP7jmbG8hncN/I+irsVh1pP2GNbpaWldRMi1K0nmZRSgjKz1kQ2L2xLZN+4RHaNewx4lMiYU427/y6VWBrTUsagcm621n7hFBPVXB/4YRpcNJjRfUfzh8V/YPyJ4xlcNDjTIYnkjIQTlJm1IjJjrh3wFbDI3XcmWk4w1Rwi91OllVpQe8vV2Vrp/sAPuzVSsbKCGctnMP7E8Tz01kOcW3wuQ3uFN11bpCWLO0GZ2UBgX6AaWFwvwWS9ltKC0myt3PrAj3ZB3jfyPgYXDebc4nPruiSzNeaobP4/IJl10LHHNltdcU2SCJLT2+7+srvPz6XkBJrFF8vJN5xMu6N2X72q3VHtOPmGkzMUUdOiH/jTvz2dX4/4NdNGT2PM9DF7zZTLFgsqFzBt9LS6Vt7QXkOZNnoaCyoXZDiylqe8vJyR3xhJSXFJTkz4kfiYe/5s51RSUuILFy7MdBhZoWJlBd966Ft0eLoDH7/+MQedcBBVI6v430v/N2v/up88bzKDCgftFl/FygoWVC7gxsE3ZjCylic6ppcLNxfvOeEHIlPjc2H24eL3InPKBhwxIMORxC/e/xtmtsjdS1KpK+EEZWYT3f2OVCrNlFxPUGFPBJi3Zt5eEw7CHNPJhQ83iS3dCSrMsc4TTjiBdevW7XX+kEMO4fXXXw+ljnR1eSpBNa7JMSgzqz+H14hs9Z5TCUqTJGLTDDNpCXJ1wo80LZ5JEp+7+9XRB2Z2bxrjSYsWM0lCM8ykhdCEH4lHPJMk9lzj5KZ0BNIifPQRnH46hLkZXZrk2oQDkYZoeaaWq8kE5e4rAcysW/C4wXXy8t5tt8HLL8Ott2Y6kiZFZ5hFW0yaYSa5Kro8U0FhQc4tzySNi3uShJk96e7npjmetKg3BvXd9957L3IyzP7p3r1hx469z7drBytWhFePuiukmeTSLL6oXJxwkIsxN+ckiUQWi21wO4tsl/bFYl99Fc47D9q3jzzed184/3x47bX01CcikgcSWeqoZd0wFWZrpLAQCgoirah27SLfCwrg+OPDq0NEJM/kRQuqWWzYANdcA/Pnw7hxOTFRQkQkmyXSgvpp2qJoCWbM2HV8992Zi0NEpIWIO0G5+9sAZnYhsAw4APgm8Iy75+7yDCIikpWS2VF3C/Av4F7gTqBPqBGlgRaLFRHJPckkqPbAcOB5d68isv1GVtOW7yLxmTxv8l43a1esrGDyvMkZikjyWTIJajMwAJgSbNt+dLghiUimDCocxJjpY5i3Zh6wa8WRQYWDMhyZ5KOEd9R191eAV4KHfw03HBFJVJir3Bd3K+aes+/Za5X74m7FodaTSzcAS+Yk04KqY2b7m9kbZtbPzE4JKygRyZz6q9xf3v9yrXIvGZNwC6o+d//czM5x9703YxGRZqFV7qWlSqkFBRBNTmbWxsyafUEpMzvazErM7OqmXy0ijYmOOd038j5+PPjHWuVeMirVLr7WZnaimZ1KZCPDJQm+f0C947FmdoaZNbp3t5l1N7O7zOzA4NSAoN4jzCylFqFIvouuch/t1tMq95JJCX+gm1krYCDQDvgKWOTuOxMsow2Rqeq/BIrNbARQ7e6zg/Gs0919Tqz3uvtGM9tMkFzdvTxITNsSjUNEdnfj4Mjfh/UnRAztNVRdfJIRcScoMxsI7EvkvqfF7p70/U/Be58yswnBqaHAs8HxKuBUYI6ZDQPqT76ocvcpMYo8B/iNmbV295o94i4DyiCyw2Y6TJ43mUGFgyjuVgxExgQqVlawoHJB3S+8iIgkJq4uviA5ve3uL7v7/FSSUwO6Al8Gx1VAFwB3n+Xu/13va4qZdQF6EmnFYWYXAWcCvwBq9yzY3ae6e4m7l3Tv3j3ksCN074iISPjiakG5+6I0x+FA6+C4NY0nzs3AVQSrq7v7o8CjjRVeb8PCunO6d0REJLs12YIys2Fm9nszOz54XJaGODYBHYPjTkCD28r7Lnu1ljJJ946IiIQrnhbUtcCVwM1B91o6duGrAHoHx72BmBMkkuXuM4GZJSUl342e070jIiLZLZ4xqI3uvsXdbyCyvUYoAyvBOn59ghbZQqCNmZ0NfNHQDL4U6krraubRMafp357Or0f8WveOiGRAQUEBBQUFmQ5DQhRPgnoawMy6uftPgAdDqvsx4DDgT8Cn7v47d3/G3e8Nqfw66V7NPHrvSLTFpHtHRERS12QXn7v/b3D4J+Bcd/9tGBXXmwmY9nuXYk2SCFOsqeS6d0REJDWJrCRhaYsizbQflIhI7kkkQXnaohAREdlDXrSgtOW7iEjuSSRB/TRtUaSZuvhERHJP3AnK3d+Odd7MOsY6n48qKyuprAxv5QgRkXyW8HYbZvY1M7vIzC4zs8uAUGb1pVNzdPGVl5dzwgkn0KNHD3r27El5eXna6hIRyQfJ7Af1X0SWI7LgK+v3YEp3F195eTllZWWsW7cOd2f16tWUlZUpSYmIpCCZ5PKQu78UfWBmOXk3aphdcRMnTqSqqmq3c1VVVUycOJGhQ8O7F6qwUIvFikj+SCZB9TGzO4B3g8f9gBPCCyn3NJTsNB4lIpK8ZBJUd+AyIhsXAnwrvHDSI9ZKEmG2RoqKili9enXM82r1iIgkJ5kxqCVEFpBd4+5rgNdDjil06R6DmjRpEh06dNjtXIcOHZg0aVJa6hMRyQfJtKD+E+hiZuuDx0cBed1MKC0tBSJjUZWVlRQVFTFp0qS68yK5RptjSjZIJkH9DHjB3WvNrBNwXMgx5aTS0tK6CRHq1hMRSV0yCaoY6A9MBjqQ562n+pSYRETCk0yCesfdZwG4+8dmlvWLyKZ7uw0REQlfMpMkBpnZ/zGzk82slEhrKqtpLT4RkdyTTIK6A9gXKAU6AreGGpGIiAhJdPG5ew3wxzTEIiIiUieZFlQdM9vfzN4ws35mdkpYQYmIiKS00Ku7f25m57j7urACEhERgRRbUADR5GRmbcxsQOohiYiIpNiCMrPWQAnQFtgOLAojqARjKCay/ccgd7+vuesXEZH0SGbDwlZmNsjMhgADgUXu/pK7L3D32gTLGlDveKyZnWFmNzbxnu5mdpeZHRicOhL4F9A7wR9FRESyWNwtKDMbSGR6eTWw2N2rm3hLY2W1AYYDvwSKzWwEUO3us4MJF6e7+5xY73X3jWa2mSC5uvtMM+sHxNySXkREclNcCSpITm+7+44wKg2S21NmNiE4NRR4NjheBZwKzDGzYUD92YFV7j5lj9i+AbwMfM/M5rr7qj2eLwPKILL9hYiI5Ia4EpS7Nzi2ZGa93H1linF0Bb4MjquALkG9s4BZe9TXBehJpHvxBaA1kfUBPwE+jBH7VDP7CBjVtm3bgSnGKSIizSThSRJBd9x5QBvAgL7ASSnG4UQSDcH3xsbGNgNXBXXj7n8Pzi9rsHD3mcDMkpKS76YYp4iINJNkZvEdB0wh3B11NxFZNgkiM/I+beiF7h5dnDbuRWq1WKyISO5J5j6oNwh/R90Kds3C6w3EnCCRrLrFYjt0gMrKyJeIiGS1eCdJvARsJdJq2QfoaWbR1SOS2lHXzC4G+gSTGB4DjjCzs4EvGprBl6xoC+qInj3DLFZERNIo3i6+29z9hfonzKyLu39qZoOTrPsx4FEiY0417v67JMtp0m5jUNpUUEQkJ8TVxVc/OZnZmWa2DHjCzN4imKyQKHev9oid9caV0sLMRpnZ1M8++yyd1YiISIiSGYM6Gujn7qcBXycyaSKracNCEZHck0yCejfa4nH3ncBSgGA1h6ykFpSISO5JZpr5L8zsJqCWyCKxnc1sI9Ar+Mo6ug9KRCT3JJOgfuPuD+550swuCCEeERERII4uPjM7pP7jWMkpOD9jz9dmC3XxiYjknnhaUGPN7HB2X7kh1sw9B94HfhFGYGFSF5+ISO5pMkG5+8+bIxAREZH6Ut7yXUREJB3yIkFFx6A2bdpEZWUllVqLT0Qk6+VFgoreqLv//vtnOhQREYlTXAnKzP6S7kCaQ5s2bSgsLKRQ6/GJiGS9eFtQ/aMHZvZCYy8UEREJQ7wJqv4U8+7pCERERKS+eFeSKDCzK4A3SXL18mxQXVtN5dbIBInC/dTNJyKSzeJNUP8FlABXAj2CbTaWBV/vuPtj6QkvHNENC3v27pnpUEREJE7xdvF9CPzc3U93927AWcCDwFfA+ekKLizRWXxdD+xK4X6Faj2JiOSAeFtQFwK3mdnBwLvAEiLdfU8Dd6YpNhERyWPx7qh7tbuXAL8E/gWsBIYCrwOr0xeeiIjkq0S327jS3et20DWze4AfhxuSiIhI4itJfG5mA6MP3H0RcGS4IYmIiCTegroKeMjM3gEWAccC1aFHlSAz2xe41t1/melYREQkHAm1oNz9PeAU4BngYGA5cHaylZvZgHrHY83sDDO7sYn3dDezu8zswHqni4FOycYhIiLZJ961+E42MwNw9xp3/5u73+Lud7n7pkQrNbM2ZnYO8EjweARQ7e6zge1mdnpD73X3jcDmaOxmdiiwJtEYREQku8XbxTcWuNvM/gU8Bzzn7uuTrdTdq4GnzGxCcGoo8GxwvAo4FZhjZsOItNiiqtx9yh7FFQLbgUPNrLO777avu5mVAWUARUVFyYYsIiLNLK4E5e7jAMzsKCI36f7ZzDoDFUQS1jx3r0khjq7Al8FxFdAlqHcWMKv+C82sC9ATGAi84O7zzaxnUEb9NQOjsU8FpgKUlJTs9byIiGSneLv4igDc/V13/7W7jwDOAF4Gvg3MTzEOB1oHx62biGszkckaL9Y7txq4ANga6w3RDQs/++yzWE+LiEgWireL7zkzO4jIpIi3gKXB95fd/ZkQ4tgEdAyOOwGfNvRCd4+2gjzGORERaSHiXUmiL5Gxnh8CrwGHA7cA75rZyhDiqAB6B8e9gTkhlFknuhZf586dwyxWRETSKO5p5u7+lbu/ATxOpEtvPZFxozeTqdjMLgb6BJMYFgJtzOxs4At3DzVBqYtPRCT3xNXFZ2bFwEjgHCIbFs4CyoEyd/8qybofAx4lMuZU4+6/S7KcJrn7TGBmSUnJd9NVh4iIhCveMajlwBvA7cCT7r4j1YqDqeYAO1MtqynR/aAOP/zwdFclIiIhibeL7xrgVeA64EMzW25m08zsFjM7L33hhUNjUCIiuSfe+6Dur//YzHoA/YmsxXch8ET4oYmISD6LdwzqLOBNd68EcPe1wFoia/Jlvca6+Kqrq1m7di3bt29v/sAkJe3bt6dHjx60adMm06GISBqkuqPuEmB5iqtIpF1jkyTWrl3LfvvtR8+ePQmWG5Qc4O5s2rSJtWvX0qtXr0yHIyJpkPc76m7fvp2uXbsqOeUYM6Nr165q+Yq0YHmxo25Ts/iUnHKTrptIy5YXO+qGPYuvsrKSysrKUMoKU01NDbW1tdTW1lJTk9W9riIiTWoRO+rmsq1bt1JRUUG7du2YPXs2P//5z3nuuecYN24cEyZM4OCDD+bdd9/lnXfe4W9/+xsAzz77LB988AGFhYVs376d5cuXc9ttt7F8+XImTJhAQUEB119/PQMGDIhZ5yOPPMK4ceP47W9/S21tLUuXLuXqq6+mb9++vPjii3zwwQd069aNpUuXcs0111BQUNCc/yQiIkCCCcrd3zOzU4isHN6fyA28/56OwLJVeXk5EydOpLKykqKiIiZNmkRpaWnS5c2ePZulS5dyyy23cPfdd/POO+8wcuRIDj/8cC677DK6du0KwC9/GdnN/tlnn2X27NlMmRLZFqu2tpbhw4cD0K9fP0455RT69evXYHICuOSSS/j5z3/O5ZdfDsD777/PpZdeymuvvca0adOYOnUqAB988IFaYiKSMYm2oAhm7P0t+MoJYa0kUV5eTllZGVVVVQCsXr2asrIygKST1DnnnMOwYcNwd6qqqujbt+9uz7s7W7Zs4YgjjgBg4sSJPP3003XPt2rVih/+8IdJ1R31ta99jeXLlwPwz3/+ky1btnDAAQdw/vnn07FjxybeLSKSHomOQeWkVMegomNOEydOrEtOUVVVVXUtqmS0bt2aL7/8kttvv53rr7+e1q1b7/Z8ZWUlS5cu5dxzz2Xbtm2sWbOGQw89dLfXjBo1Kqm6oxYsWMDo0aMBuOiiizj66KO56qqraNOmDQcccEBKZYuIJKvJBGVmw8zs92Z2fPC4LP1hZaeGklCqEya6du3KT3/6U+6++25WrVpVd37GjBlcf/31u722traW6PZX27Zt4y9/+QuXXHIJc+fOjVl2bW1tzPM7duzg+eef5/HHH2f16tXcd999AFx77bX8/e9/56ijjmLYsGFs27YtpZ9NRCRZ8bSgriUylfxSMzsDOD69IWWfwsJCCgsLKSoqivl8UVERhYWFSZU9a9Ys3njjDQCOOeYYnnlm1+IcF1xwAVOnTqV9+/YAdOzYkZ49e9YlsY4dOzJ27FjeeOMNTjvttJjlP/nkkzHPuzvDhw/n/PPPZ+zYsXWrMSxZsoS+ffty4403cv311zNnTqg7n4iIxC2eBLXR3be4+w3AN4FBaY4pa02aNIkOHTrsdq5Dhw5MmjQp6TLnzJnDpk2bANiyZQs9evTY7fkuXbpw4okn1j2+4447uP322+se19TU0Ng+V++//35C8fz1r3+tO27Xrh3du3dP6P0iImGJZ5JE3Yi8u//EzH6QxnjSIqxJEtGJEGHO4rvuuut48skneeCBB+jevTujRo3ixRdfZMWKFTz44IN8+9vf3i1pnXXWWbg7t956K0VFRbslyCVLljB37lw++ugjtm/fzrxYficuAAARPUlEQVR589i4ceNedc6YMYMNGzbw5z//mcsuu2y3ca/a2loefvhhampq+PLLLxk0KG//HhGRDLPoeEajLzL7i7uPbYZ40qqkpMQXLly427nly5dz9NFHJ1ROdMwp2W69puzcuZPWrVtTU1NDq1ataNUq3Lks1dXV7LPPPuzcuTPnF1pN5vpJy1S5Nfi93C89v5fpsPi9xQAMOKLh20KyTbz/zma2KFgiL2nxTjPvX6/SF9z9m6lUmuvSlZii9tlnn92+hy2alHI9OYlIyxbvn+b1m1kalBARkbSL90/0AjO7gsgWG1qhU0RE0i7eBPVfQAlwJdDDzN4ClgVf77j7Y+kJT0RE8lW8W75Prf94jy3fzwOUoEREJFRJjcJn05bvZjaYyI6/s939qWapNLpyRJonSyQqumqEmVFTU5O2SRYiIs0ho59gZjbA3RcHx2OBD4ESd5/cyHu6AzcB/+3umwHc/UfNEW+63XLLLdx222289NJL/OAHP2DkyJEce+yxfPjhh8ycObNuOaPXX3+duXPn0rNnT3bu3Ml7773Hddddx5YtW7jlllvYvHkzEydObHB1iaeffpprrrmGH/3oR3Tp0oVly5YxatQohgwZwqJFi5g3bx49evRg7dq1nHTSSZxwwgnN+c8gIgJkKEGZWRtgOJEt5IvNbARQ7e6zzayfmZ3u7jHX2HH3jWa2mV0zENeZWW/gOHd/vFl+gA0b4Npr4fHHIaS9kubPn8/LL78MwKmnnsrXv/51LrjgAgYOjOwPGb0X6q233uKuu+7i4YcfrnvvmDFjqK2tpVevXgwbNoxPPvmkweQEMHLkSPr06VO3ncf27ds54ogjWLNmDffff3/ddhvTpk2julrbfYlIZmRkNXN3rw6646KrrA6td7wKOBXqFqr9z3pfsbaXHwysAfqaWYc9nzSzMjNbaGYLY62qkJS77oL58+HWW8MpDxpdlHXTpk307duX6upqbrrpJq655prdnr/22mtp27Zt0nW3b9+e7du3s2XLFtauXcu6desAGD58uJY6EpGMyZZBiq7Al8FxFdAFwN1nAbPqv9DMugA9gYHAC8BCIpM13nf33ffCiJQx1cw+Aka1bdt24J7PxyU65tS7N+zYsev8vfdGvtq1gxUrkh6Teu211xpcUqi2tpZZs2Zx8cUXA5GW1p57Rn3jG99Iqt6oVatWcdxxx3HggQfyve99jwEDBnDGGWfwk5/8hOOOOy6lskVEkpUt+0E5EF0QrjWNx7WZyNbzLwK4+z/d/Q13f7TBwlPcD6rOq6/CeedBsLo4++4L558Pr72WUrFbt25lv/322+v8888/z/e//3121E+K7JoM4e6Ul5dzxRVX8NhjsSdSNraU1QsvvMCTTz7J3//+d5544gkAvvWtbzF//nwGDx7MRRddxMqVK5P9sUREUpItCWoTEN26tRPwaUMv9F1ib3QUg5mNMrOpja363ajCwsjX178eGXPasSPSatqxI/L4+OOTbj3Nnz+/wUkIw4cP59577+Xggw+uO3fiiSfy9ttvA5HZeqWlpaxYsaLR7TYaSlKnnXYa5557Lt/5znfo1KkTEFlwtmfPnlx33XXceeedPPvss0n9XCIiqcqWBFUB9A6OewOhbkIUWgsKIhMkLr8cZs6EceNg/fqUivvwww957rnnePTRR/n444+pqKjY6zUjRoyoO/7Zz37Gr371q92SzsaNGxtcUPa9997DLP7FP7Tdhohki4yNQZnZxUCfYIfex4AjzOxs4IuGZvClUFco220AMGPGrjGpYcNSLi661frDDz/Mjh072GeffXjttddYunQp06dPp1OnThQXF9e9vn///tx8883cfPPNHH744XTo0IGbb76Z1q1b8/777/PMM8+wc+dOCgsLefPNN5k7dy433HDDbnVGt/N4+OGHGTdu3G7di/vuuy8PP/ww7s66deuYMGFCyj+jiEgy4tpuIy0VR6aa7yQy5lTjzRBIWNttpONG3egWGDU1NUBkWnn0n6T+fk1hqL+dR+vWrRNqYWUbbbchUdpuo3lk43YboXP36A02O9NdV6gtKEjLChLRrS+aY/WHdG/nISIShmwZg0qrUMegRESkWeRFgkp5Fp+IiDS7vEhQakGJiOSevEhQIiKSe/IiQYXdxVe5tbJuJks2qa2tpba2Fndn5860zz0REUmrvJjG5e4zgZklJSXfzXQsjcmX7TaWLFnChAkT6NOnD//2b//Gli1buP/++/nrX//KkUceyVdffcVvfvMbCgoK2H///amurmbdunX88Ic/DDUOEclueZGgwjJ53mQGFQ6iuNuuG2crVlawoHIBNw6+MaWy82m7jeOPP57TTjuNfv36ceGFFwKRBW/Xr1/PkUceydVXX8348eMZMCByb8jbb79dl5xFJH/kRRdfWAYVDmLM9DHMWzMPiCSnMdPHMKgw9krkicjn7TbWr19PcXEx27dvZ/HixVRWVtYlJ4B+/fpx5plnpj0OEckueZGgUh2Dio45FXcr5p6z76HsqTL+8x//yei/jeaes++huFtxSmNS8Wy3cdZZZ9GmTZsGt9uItRp6vGJtt3HJJZewatUqjjzyyAbf9+qrr9KvXz/+8Y9/8Ic//IHf//73QGT9vz/96U/MmjWLP//5z03W/8ILL9T9HPPnz4+5MsSoUaOS++FEJGflRYIKc5r54KLBjO47mj8s/gOX97+cwUWDUy4zV7fbOPnkkznooIM46KCDuOKKK5g+fToA1113HZdffjnDhg1j7dq1LFq0KOb7X375Ze644w5eeeUVgLpWYPTnA1i2bBm33nor48eP5+OPP24wFhFpefIiQaWqcL/Cuq9/fvJPZiyfwfgTx/PQWw/xz0/+WfdcMnJ9uw13p2/fvuyzzz589dVXAFRWVtYto3TQQQex5/qHUUOGDGHixIl1C+YuW7aME088kWXLltW95phjjmHIkCFUV1dz0EEHNRqLiLQsSlAJiI453TfyPn48+MdMGz2NMdPHULFy7y0y4pUr220k0j3apUuXuskVK1as4Nhjj2309WeeeSY7d+5kzpw5DBgwgEMOOYSXXnqp7vnGfj4Rabk0iy8BCyoXMG30tLpZfEN7DWXa6GksqFzA0F5DkyozF7bbWLt2LSNHjuTNN9/crZwlS5awcuVKZs6cCUTGst58803uuece7r33XoqKijj00EM55ZRT9nrf3Llz2bBhAzt37mTbtm08/fTTXHLJJQA88MADTJ48mTfffJMuXbrQo0ePvcoQkZYvY9ttZEJY222kY1n/bN9uY+nSpfTv3z/UOMKg7TYkStttNI+82G6jOYW93UY6fgGyfbuNzz//PF3hiIjElBcd+1osNnVDhgzJdAgikmfyIkE1JZ+6OVsSXTeRli3vE1T79u3ZtGmTPuxyjLuzadMm2rdvn+lQRCRN8mIMqjHRRVE3btyY6VAkQe3bt6dHjx6ZDkNE0iTvE1SbNm3o1atXpsMQEZE95HyCMrMCYAiw3d2fynQ8IiISjoyOQZnZgHrHY83sDDNrdN8KM+tuZneZ2YHBqTJgBhpPExFpUTLyoW5mbczsHOCR4PEIoNrdZwPbzez0ht7r7huBzeyKvRdwfJpDFhGRZpaRLj53rwaeMrMJwamhQHRV0lXAqcAcMxsG1F/jpsrdp+xR3IfuvtjMhpvZi+5eVf9JMysj0sqCSPJbRtM6A40tPtfY8w09F+t8rHPdgE/iiDEdmvq501lOvO/J1LXRdUntdfqdCbecXLguxaTK3TP2BVQE3/8AnBgcnwn8qpH3dAEeAL4ZPC4BTgLOj6O+qXHG1ejrGnu+oedinW/g3MIMXo+4/n3SUU62Xxtdl+y8Lvl6bfLlumTLJAkHogvOtabxrsfNwFWAAbh77L0cYpsZ0usae76h52Kdjzee5hJWPMmUo2vTMF2X+Otpbpm6NnlxXTK6WKyZVbj7UDO7Hfi7u88yswuBY9z91owFlkFmttBTXGBRwqfrkr10bbJTGNclW2a+VQC9g+PewJwMxpJpUzMdgMSk65K9dG2yU8rXJWMtKDO7GJgM/Ax4DLgEWAEc5u73ZiQoERHJGplMUG2AnUTGnGo8k32NIiKSdfJqw0IREckd2TIGJSIishslKBERyUrZch+UxMnMTiJyD1hvdy/PdDz5zsyGAl8BBe7+WKbjkQgz6wW0BQ5297mZjkd2MbNLgb7AvsBv3X1FQ69VC6oZhbQ47hnu/irQ2cz2S2e8+SaZ6wOc5+7zgH7BxB8JWZLX5WQiq85oR8s0SvLavODu/w480VhyAiWoZhHy4rgHB993AIemL+r8kcr1IbJ2JESuS5e0BppnUrwurwHXE0lUErIUP9M+NrP+wPtN1aMuvmbg4S6O+1XwvSNQmZ6I80sq1wf4Y/BX5Gbg4+aJOD+keF1WuvvFZnatme3r7l82V9z5IITPtDPc/a6m6lGCyoyuQPQXporgL293nwXMqv9CM+sC9AQGAi8AT5jZYGCru29proDzTCLXZyjwBfCM7uVLu0SuyynBHw5vKzk1i7ivTSCu4QklqMxIZXHcecH5eQ2+Q1IV9/Vx94pmiUggsevySrNEJFGJfKbh7rfFU6gSVGZsItJFB9AJ+LShF9b7q1x/nTefuK+PNCtdl+yVlmujSRKZocVxs5uuT3bSdcleabk2SlDNJFgct0+ww+9CoI2ZnQ184e76RcswXZ/spOuSvZrj2mgtvmaixXGzm65PdtJ1yV7NcW2UoEREJCupi09ERLKSEpSIiGQlJSgREclKSlAiIpKVlKBERCQrKUGJiEhWUoISSYCZ3WRmWbnNiZmdZ2a/N7P/NbNvhlhubzP7o5lNr3dukJl9J6w6RGJRghIJmNn3zOwjM1sSfP2/GC873N0/bOD9w83s5T3OlZhZs/yeufsT7v5d4ArgorDqd/cV7v6dPc4tAIakUq5IU7RYrMgu/YGb3f2PybzZ3Z83s/+7x+m3gQFEloIJhZkdC/xij9NXuXt0P6qbgbuTqT+OskWajRKUyC7HAn9KsYzdlmZx9+1mtm+KZe5egftbwDl7njczA24HnnX3xcnU31DZIpmgLj6RXY4BHgi6915s6EVmNtbMFgVdej8Mdg2Nam1mI8zsJ2Z2VnBup5nF/GPQzPYJxrWGm9n44NxZZjYvOHejmV1pZqeZ2VQzO6SR+H8AnAmMNrNx9c7HrN/MBpvZOWY2zMy+3cjP29XM7gO+bmY/baR+kVBpLT4RIJj48Jy7H9PE6x5w9yvNbK67nxaM77wNHOPubmb/Ak4AugE/dvfvmVlb4Hh3fz1GeVcD2939ITP7CTDF3WvM7F3gpKCch9z9pOC129z9kQR/tpj1m9nNwBvAc0Cxu7+TYLkPuPuVibxHJBFqQYlE9AeWJfD6GgB3ryWyOVvX4Hylu28BqoG2wWu+Ato1UM7XgZXB625395rg/Pp65bwbnKsrMxGN1H8fcBqwGDgg0XJF0k0JSiTiWBJLUG0h0kUHHOjunwTnG+qSqA62JyB4X+fg8C2Cjd7M7GAz2y9GOWF0c8Sqf5i7TwROAUKbli4SFiUokYhjgUS6uDqb2enAOOAaADMbAhxtZicD3wKOMbM+wesXAQOD1/UA5gbn/wgUm9kY4BvuvjVGOcea2TFExpeG1ktiiYhVfx8zGw2MAp5JokyRtNIYlEgC6o1BVbj70ATfe5q7zw2O+7v70vRE2Tz1awxK0k0tKJEEmdmFQG8z65/gW3eYWXQsaP+Qw8qF+kUSogQlkpj3gdfd/bAkWiALiXQl4u4vN/HadAitfjMbBGTiZ5A8oi4+ERHJSmpBiYhIVlKCEhGRrKQEJSIiWUkJSkREspISlIiIZCUlKBERyUpKUCIikpWUoEREJCv9f7CO3Kh9QsqqAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot\n",
    "\n",
    "plt.errorbar(x_counts_3FGL,x_counts_3FGL**2*y_counts_3FGL,xerr=[x_errors_L_3FGL,x_errors_H_3FGL],\n",
    "             yerr=x_counts_3FGL**2*np.array([error_L_3FGL,error_H_3FGL]), fmt='o', \n",
    "             color='black', label='3FGL PS', elinewidth=0.2)\n",
    "\n",
    "plt.errorbar(x_counts_4FGL,x_counts_4FGL**2*y_counts_4FGL,xerr=[x_errors_L_4FGL,x_errors_H_4FGL],\n",
    "             yerr=x_counts_4FGL**2*np.array([error_L_4FGL,error_H_4FGL]), fmt='*', \n",
    "             color='red', label='4FGL PS', elinewidth=0.2)\n",
    "\n",
    "plt.errorbar(x_counts_4FGL_noeg,x_counts_4FGL_noeg**2*y_counts_4FGL_noeg,xerr=[x_errors_L_4FGL_noeg,x_errors_H_4FGL_noeg],\n",
    "             yerr=x_counts_4FGL_noeg**2*np.array([error_L_4FGL_noeg,error_H_4FGL_noeg]), fmt='x', \n",
    "             color='green', label='4FGL PS, no EG', elinewidth=0.2)\n",
    "\n",
    "plt.xscale(\"log\")\n",
    "plt.yscale(\"log\")\n",
    "\n",
    "plt.xlim(4e-11,1e-7)\n",
    "plt.ylim(1e-16,1e-9)\n",
    "\n",
    "plt.xlabel(\"$F$ [ph\\,cm$^{-2}$\\,s$^{-1}$]\")\n",
    "plt.ylabel(\"$F^2\\,dN/dF$ [ph\\,cm$^{-2}$\\,s$^{-1}$\\,deg$^{-2}$]\")\n",
    "\n",
    "plt.legend(loc='lower left')\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}