{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Astrophysical source population modeling with Gammapy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "The [gammapy.astro.population](https://docs.gammapy.org/dev/astro/population/index.html) package contains some simple Galactic source population models.\n",
    "\n",
    "Here we provide some Python code to compute observable parameter distributions for Galactic gamma-ray source populations.\n",
    "\n",
    "* Observables: Flux, GLON, GLAT\n",
    "* Source classes: Pulsar (PSR), Supernova remnant (SNR), pulsar wind nebula (PWN)\n",
    "\n",
    "References:\n",
    "\n",
    "* Section 6.2 in the Fermi-LAT collaboration paper [\"The First Fermi-LAT Catalog of Sources Above 10 GeV\"](http://adsabs.harvard.edu/abs/2013arXiv1306.6772T)\n",
    "* Axel Donath's bachelor thesis [\"Modelling Galactic gamma-ray source populations\"](http://pubman.mpdl.mpg.de/pubman/item/escidoc:912132:1/component/escidoc:912131/BScThesis_ddonath.pdf), specifically Chapter 4.\n",
    "* Casanova & Dingus (2008), [\"Constraints on the TeV source population and its contribution to the galactic diffuse TeV emission\"](http://adsabs.harvard.edu/abs/2008APh....29...63C)\n",
    "* Strong (2007), [\"Source population synthesis and the Galactic diffuse gamma-ray emission\"](http://adsabs.harvard.edu/abs/2007Ap%26SS.309...35S)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import astropy.units as u\n",
    "from gammapy.utils.random import sample_powerlaw\n",
    "from gammapy.astro import population"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulate positions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Spatial distribution using Lorimer (2006) model\n",
    "n_sources = int(1e5)\n",
    "\n",
    "table = population.make_base_catalog_galactic(\n",
    "    n_sources=n_sources,\n",
    "    rad_dis=\"L06\",\n",
    "    vel_dis=\"F06B\",\n",
    "    max_age=1e6 * u.yr,\n",
    "    spiralarms=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulate luminosities"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Several source population models, e.g. the 1FHL paper or Strong (2007), use power-law luminosity functions.\n",
    "\n",
    "Here we implement the \"reference model\" from the 1FHL catalog paper section 6.2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Source luminosity (ph s^-1)\n",
    "\n",
    "luminosity = sample_powerlaw(x_min=1e34, x_max=1e37, gamma=1.5, size=n_sources)\n",
    "table[\"luminosity\"] = luminosity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute observable parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "table = population.add_observed_parameters(table)\n",
    "table.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Check output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The simulation is done, you could save the simulated catalog to a file.\n",
    "\n",
    "Here we just plot a few distributions to check if the results look OK."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(table[\"x\"][:1000], table[\"y\"][:1000]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(table[\"GLON\"], bins=100);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(table[\"GLAT\"], bins=100, log=True);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(table[\"GLON\"][:1000], table[\"GLAT\"][:1000]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(table[\"distance\"], bins=100, log=True);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(np.log10(table[\"luminosity\"]), bins=100, log=True);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(np.log10(table[\"flux\"]), bins=100, log=True);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: plot GLON, GLAT, FLUX distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "TODO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Start exercises here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What next?\n",
    "\n",
    "TODO: summarise what was done here briefly.\n",
    "\n",
    "TODO: add some pointers to other documentation."
   ]
  }
 ],
 "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.5.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}