{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "
\n", "**This is a fixed-text formatted version of a Jupyter notebook.**\n", "\n", " You can contribute with your own notebooks in this\n", " [GitHub repository](https://github.com/gammapy/gammapy-extra/tree/master/notebooks).\n", "\n", "**Source files:**\n", "[source_diffuse_estimation.ipynb](../_static/notebooks/source_diffuse_estimation.ipynb) |\n", "[source_diffuse_estimation.py](../_static/notebooks/source_diffuse_estimation.py)\n", "
\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Simulating and analysing sources and diffuse emission\n", "\n", "## Introduction\n", "\n", "This notebook shows how to do the analysis presented in this paper using Gammapy:\n", "\n", "* Ellis Owen et al. (2015). *\"The $\\gamma$-ray Milky Way above 10 GeV: Distinguishing Sources from Diffuse Emission\"* ([ADS](http://adsabs.harvard.edu/abs/2015arXiv150602319O))\n", "\n", "The following parts of Gammapy are used:\n", "\n", "* [gammapy.image.SkyImage](http://docs.gammapy.org/dev/api/gammapy.image.SkyImage.html)\n", "* [gammapy.detect.KernelBackgroundEstimator](http://docs.gammapy.org/dev/api/gammapy.detect.KernelBackgroundEstimator.html)\n", "* TODO: code to simulate images from a source catalog is where?!?\n", "\n", "\n", "## Setup" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "plt.style.use('ggplot')" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import logging\n", "logging.basicConfig(level=logging.INFO)\n", "from astropy.io import fits\n", "from scipy.ndimage import convolve\n", "from gammapy.stats import significance\n", "# from gammapy.image import SkyImage\n", "from gammapy.detect import KernelBackgroundEstimator" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Simulation\n", "\n", "TODO" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Parameters\n", "TOTAL_COUNTS = 1e6\n", "SOURCE_FRACTION = 0.2\n", "\n", "CORRELATION_RADIUS = 0.1 # deg\n", "SIGNIFICANCE_THRESHOLD = 4.\n", "MASK_DILATION_RADIUS = 0.2 # deg\n", "NUMBER_OF_ITERATIONS = 3\n", "\n", "# Derived parameters\n", "DIFFUSE_FRACTION = 1. - SOURCE_FRACTION" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Load example model images\n", "source_image_true = fits.getdata('sources.fits.gz')\n", "diffuse_image_true = fits.getdata('diffuse.fits.gz')" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Generate example data\n", "source_image_true *= SOURCE_FRACTION * TOTAL_COUNTS / source_image_true.sum()\n", "diffuse_image_true *= DIFFUSE_FRACTION * TOTAL_COUNTS / diffuse_image_true.sum()\n", "total_image_true = source_image_true + diffuse_image_true\n", "\n", "counts = np.random.poisson(total_image_true)\n", "\n", "print('source counts: {0}'.format(source_image_true.sum()))\n", "print('diffuse counts: {0}'.format(diffuse_image_true.sum()))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Analysis\n", "\n", "TODO" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Start with flat background estimate\n", "background=np.ones_like(counts, dtype=float)\n", "images = KernelBackgroundEstimatorData(counts=counts, background=background)\n", "\n", "# CORRELATION_RADIUS\n", "source_kernel = np.ones((5, 5))\n", "background_kernel = np.ones((100, 10))\n", "\n", "kbe = KernelBackgroundEstimator(\n", " images=images,\n", " source_kernel=source_kernel,\n", " background_kernel=background_kernel,\n", " significance_threshold=SIGNIFICANCE_THRESHOLD,\n", " mask_dilation_radius=MASK_DILATION_RADIUS\n", ")\n", "\n", "kbe.run(n_iterations=3)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "kbe.run_iteration()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "kbe._data[1].print_info()" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "## Plot results\n", "\n", "TODO" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "n_iterations = 5\n", "\n", "#plt.clf()\n", "plt.figure(figsize=(10, 10))\n", "\n", "for iteration in range(n_iterations):\n", " filename = 'test{0:02d}_mask.fits'.format(iteration)\n", " mask = fits.getdata(filename)[100:300,:]\n", "\n", " plt.subplot(n_iterations, 2, 2 * iteration + 1)\n", " filename = 'test{0:02d}_background.fits'.format(iteration)\n", " data = fits.getdata(filename)[100:300,:]\n", " plt.imshow(data)\n", " plt.contour(mask, levels=[0], linewidths=2, colors='white')\n", " plt.axis('off')\n", " plt.title(filename)\n", " \n", " plt.subplot(n_iterations, 2, 2 * iteration + 2)\n", " filename = 'test{0:02d}_significance.fits'.format(iteration)\n", " data = fits.getdata(filename)[100:300,:]\n", " plt.imshow(data, vmin=-3, vmax=5)\n", " plt.contour(mask, levels=[0], linewidths=2, colors='white')\n", " plt.axis('off')\n", " plt.title(filename)\n", " #plt.colorbar()\n", "\n", "#plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercises\n", "\n", "- Vary some parameter and check how results change\n", "- Apply same method to 2FHL or some other real dataset?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## What next?\n", "\n", "TODO: summarise\n", "\n", "TODO: pointers to other docs" ] } ], "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.2" }, "nbsphinx": { "orphan": true } }, "nbformat": 4, "nbformat_minor": 0 }