{
"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",
"[cube_analysis_part2.ipynb](../_static/notebooks/cube_analysis_part2.ipynb) |\n",
"[cube_analysis_part2.py](../_static/notebooks/cube_analysis_part2.py)\n",
"
\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Cube analysis with Gammapy (part 2)\n",
"\n",
"## Introduction \n",
"\n",
"In this tutorial we will learn how to compute a morphological and spectral fit simultanously.\n",
"\n",
"This is part 2 of 2 for IACT cube analysis. If you haven't prepared your cubes, PSF and EDISP yet, do part 1 first.\n",
"\n",
"The fitting is done using [Sherpa](http://cxc.harvard.edu/contrib/sherpa47/).\n",
"\n",
"We will use the following classes:\n",
"\n",
"- [gammapy.cube](http://docs.gammapy.org/dev/cube/index.html) where are strored the counts, the background, the exposure and the mean psf of this Crab dataset.\n",
"- [gammapy.irf.EnergyDispersion](http://docs.gammapy.org/dev/api/gammapy.irf.EnergyDispersion.html) where is stored the mean rmf of this Crab dataset.\n",
"- the method [gammapy.cube.SkyCube.to_sherpa_data3d](http://docs.gammapy.org/dev/api/gammapy.cube.SkyCube.html#gammapy.cube.SkyCube.to_sherpa_data3d) to transform the counts cube in Sherpa object.\n",
"- [gammapy.cube.CombinedModel3DInt](http://docs.gammapy.org/dev/api/gammapy.cube.CombinedModel3DInt.html) to combine the spectral and spatial model for the fit if you consider a perfect energy resolution\n",
"- [gammapy.cube.CombinedModel3DIntConvolveEdisp](http://docs.gammapy.org/dev/api/gammapy.cube.CombinedModel3DIntConvolveEdisp.html) to combine the spectral and spatial model for the fit taking into account that the true energy is different than the reconstructed one.\n",
"\n",
"We will use the cubes built on the 4 Crab observations of gammapy-extra.\n",
"\n",
"You could use the cubes we just created with the notebook `cube_analysis_part1.ipynb`"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setup\n",
"\n",
"As always, we start with some notebook setup and imports"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"\n",
"import logging\n",
"logging.getLogger(\"sherpa.fit\").setLevel(logging.ERROR)"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"WARNING: imaging routines will not be available, \n",
"failed to import sherpa.image.ds9_backend due to \n",
"'RuntimeErr: DS9Win unusable: Could not find xpaget on your PATH'\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"WARNING:sherpa.image:imaging routines will not be available, \n",
"failed to import sherpa.image.ds9_backend due to \n",
"'RuntimeErr: DS9Win unusable: Could not find xpaget on your PATH'\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"WARNING: failed to import sherpa.astro.xspec; XSPEC models will not be available\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"WARNING:sherpa.astro.all:failed to import sherpa.astro.xspec; XSPEC models will not be available\n"
]
}
],
"source": [
"import numpy as np\n",
"from astropy import units as u\n",
"from astropy.io import fits\n",
"from astropy.coordinates import SkyCoord, Angle\n",
"from regions import CircleSkyRegion\n",
"\n",
"from gammapy.extern.pathlib import Path\n",
"from gammapy.irf import EnergyDispersion\n",
"from gammapy.cube import SkyCube\n",
"from gammapy.cube.sherpa_ import (\n",
" CombinedModel3DInt,\n",
" CombinedModel3DIntConvolveEdisp,\n",
" NormGauss2DInt,\n",
")\n",
"\n",
"from sherpa.models import PowLaw1D, TableModel\n",
"from sherpa.estmethods import Covariance\n",
"from sherpa.optmethods import NelderMead\n",
"from sherpa.stats import Cash\n",
"from sherpa.fit import Fit"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 3D analysis assuming that there is no energy dispersion (perfect energy resolution)\n",
"\n",
"### Load the different cubes needed for the analysis\n",
"\n",
"We will use the Cubes build on the 4 Crab observations of gammapy-extra. You could use the Cubes we just created with the notebook cube_analysis.ipynb by changing the cube_directory by your local path.\n",
"\n",
"- Counts cube"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"cube_dir = Path('$GAMMAPY_EXTRA/test_datasets/cube')\n",
"counts_3d = SkyCube.read(cube_dir / 'counts_cube.fits')\n",
"# Transformation to a sherpa object\n",
"cube = counts_3d.to_sherpa_data3d(dstype='Data3DInt')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- Background Cube"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"bkg_3d = SkyCube.read(cube_dir / 'bkg_cube.fits')\n",
"bkg = TableModel('bkg')\n",
"bkg.load(None, bkg_3d.data.value.ravel())\n",
"bkg.ampl = 1\n",
"bkg.ampl.freeze()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- Exposure Cube"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"exposure_3d = SkyCube.read(cube_dir / 'exposure_cube.fits')\n",
"i_nan = np.where(np.isnan(exposure_3d.data))\n",
"exposure_3d.data[i_nan] = 0\n",
"# In order to have the exposure in cm2 s\n",
"exposure_3d.data = exposure_3d.data * u.Unit('m2 / cm2').to('')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- PSF Cube"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"psf_3d = SkyCube.read(cube_dir / 'psf_cube.fits')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Setup Combined spatial and spectral model"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Define a 2D gaussian for the spatial model\n",
"spatial_model = NormGauss2DInt('spatial-model')\n",
"\n",
"# Define a power law for the spectral model\n",
"spectral_model = PowLaw1D('spectral-model')\n",
"\n",
"# Combine spectral and spatial model\n",
"coord = counts_3d.sky_image_ref.coordinates(mode=\"edges\")\n",
"energies = counts_3d.energies(mode='edges').to(\"TeV\")\n",
"# Here the source model will be convolve by the psf:\n",
"# PSF * (source_model * exposure)\n",
"source_model = CombinedModel3DInt(\n",
" coord=coord,\n",
" energies=energies,\n",
" use_psf=True,\n",
" exposure=exposure_3d,\n",
" psf=psf_3d,\n",
" spatial_model=spatial_model,\n",
" spectral_model=spectral_model,\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Set starting value"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"center = SkyCoord(83.633083, 22.0145, unit=\"deg\").galactic\n",
"source_model.gamma = 2.2\n",
"source_model.xpos = center.l.value\n",
"source_model.ypos = center.b.value\n",
"source_model.fwhm = 0.12\n",
"source_model.ampl = 1.0"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Fit"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Define the model\n",
"flux_factor = 1e-11\n",
"model = bkg + flux_factor * source_model"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Method = neldermead\n",
"Statistic = cash\n",
"Initial fit statistic = 8614.61\n",
"Final fit statistic = 8028.8 at function evaluation 619\n",
"Data points = 12500\n",
"Degrees of freedom = 12495\n",
"Change in statistic = 585.806\n",
" spatial-model.xpos 184.195 \n",
" spatial-model.ypos -6.16917 \n",
" spatial-model.ampl 6.16628 \n",
" spatial-model.fwhm 0.0765675 \n",
" spectral-model.gamma 2.30648 \n"
]
}
],
"source": [
"# Fit to the counts Cube sherpa object\n",
"fit = Fit(\n",
" data=cube,\n",
" model=model,\n",
" stat=Cash(),\n",
" method=NelderMead(),\n",
" estmethod=Covariance(),\n",
")\n",
"fit_results = fit.fit()\n",
"print(fit_results.format())"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Confidence Method = covariance\n",
"Iterative Fit Method = None\n",
"Fitting Method = neldermead\n",
"Statistic = cash\n",
"covariance 1-sigma (68.2689%) bounds:\n",
" Param Best-Fit Lower Bound Upper Bound\n",
" ----- -------- ----------- -----------\n",
" spatial-model.xpos 184.195 -0.0128654 0.0128654\n",
" spatial-model.ypos -6.16917 -0.0124289 0.0124289\n",
" spatial-model.ampl 6.16628 -0.308859 0.308859\n",
" spatial-model.fwhm 0.0765675 -0.01881 0.01881\n",
" spectral-model.gamma 2.30648 -0.0599405 0.0599405\n"
]
}
],
"source": [
"err_est_results = fit.est_errors()\n",
"print(err_est_results.format())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Add an exlusion mask for the Fit\n",
"For example if you want to exclude a region in the FoV of your cube, you just have to provide a SkyCube with the same dimension than the Counts cube and with 0 in the region you want to exlude and 1 outside. With this SkyCube mask you can select only some energy band for the fit or just some spatial region whatever the energy band or both."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Define the mask\n",
"Here this is a test case, there is no source to exlude in our FOV but we will create a mask that remove some events from the source. You just see at the end that the amplitude fitted in the 3D analysis is lower than the one when you use all the events from the Crab. The principle works even if this is not usefull here, just a testcase."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"exclusion_region = CircleSkyRegion(\n",
" center=SkyCoord(83.60, 21.88, unit='deg'),\n",
" radius=Angle(0.1, \"deg\"),\n",
")\n",
"\n",
"sky_mask_cube = counts_3d.region_mask(exclusion_region)\n",
"sky_mask_cube.data = sky_mask_cube.data.astype(int)\n",
"index_region_selected_3d = np.where(sky_mask_cube.data.value == 1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Add the mask to the data"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Set the counts\n",
"cube = counts_3d.to_sherpa_data3d(dstype='Data3DInt')\n",
"cube.mask = sky_mask_cube.data.value.ravel()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Select only the background pixels of the cube of the selected region to create the background model"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Set the bkg and select only the data points of the selected region\n",
"bkg = TableModel('bkg')\n",
"bkg.load(None, bkg_3d.data.value[index_region_selected_3d].ravel())\n",
"bkg.ampl = 1\n",
"bkg.ampl.freeze()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Add the indices of the selected region of the Cube to combine the model"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# The model is evaluated on all the points then it is compared with the data only on the selected_region\n",
"source_model = CombinedModel3DInt(\n",
" coord=coord,\n",
" energies=energies,\n",
" use_psf=True,\n",
" exposure=exposure_3d,\n",
" psf=psf_3d,\n",
" spatial_model=spatial_model,\n",
" spectral_model=spectral_model,\n",
" select_region=True,\n",
" index_selected_region=index_region_selected_3d,\n",
")\n",
"\n",
"# Set starting values\n",
"source_model.gamma = 2.2\n",
"source_model.xpos = center.l.value\n",
"source_model.ypos = center.b.value\n",
"source_model.fwhm = 0.12\n",
"source_model.ampl = 1.0\n",
"\n",
"# Define the model\n",
"model = bkg + flux_factor * (source_model)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Method = neldermead\n",
"Statistic = cash\n",
"Initial fit statistic = 670.276\n",
"Final fit statistic = 448.967 at function evaluation 553\n",
"Data points = 390\n",
"Degrees of freedom = 385\n",
"Change in statistic = 221.309\n",
" spatial-model.xpos 184.607 \n",
" spatial-model.ypos -5.81938 \n",
" spatial-model.ampl 10.4384 \n",
" spatial-model.fwhm 0.00276986 \n",
" spectral-model.gamma 2.41052 \n"
]
}
],
"source": [
"fit = Fit(\n",
" data=cube,\n",
" model=model,\n",
" stat=Cash(),\n",
" method=NelderMead(),\n",
" estmethod=Covariance(),\n",
")\n",
"fit_results = fit.fit()\n",
"print(fit_results.format())"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Confidence Method = covariance\n",
"Iterative Fit Method = None\n",
"Fitting Method = neldermead\n",
"Statistic = cash\n",
"covariance 1-sigma (68.2689%) bounds:\n",
" Param Best-Fit Lower Bound Upper Bound\n",
" ----- -------- ----------- -----------\n",
" spatial-model.xpos 184.607 ----- -----\n",
" spatial-model.ypos -5.81938 ----- -----\n",
" spatial-model.ampl 10.4384 ----- -----\n",
" spatial-model.fwhm 0.00276986 ----- -----\n",
" spectral-model.gamma 2.41052 ----- -----\n"
]
}
],
"source": [
"err_est_results = fit.est_errors()\n",
"\n",
"print(err_est_results.format())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The fitted flux is less than in the previous example since here the mask remove some events from the Crab"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 3D analysis taking into account the energy dispersion"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Set the counts\n",
"counts_3d = SkyCube.read(cube_dir / \"counts_cube.fits\")\n",
"cube = counts_3d.to_sherpa_data3d(dstype='Data3DInt')\n",
"\n",
"# Set the bkg\n",
"bkg_3d = SkyCube.read(cube_dir / 'bkg_cube.fits')\n",
"bkg = TableModel('bkg')\n",
"bkg.load(None, bkg_3d.data.value.ravel())\n",
"bkg.ampl = 1\n",
"bkg.ampl.freeze()\n",
"\n",
"# Set the exposure\n",
"exposure_3d = SkyCube.read(cube_dir / 'exposure_cube_etrue.fits')\n",
"i_nan = np.where(np.isnan(exposure_3d.data))\n",
"exposure_3d.data[i_nan] = 0\n",
"exposure_3d.data = exposure_3d.data * 1e4\n",
"\n",
"# Set the mean psf model\n",
"psf_3d = SkyCube.read(cube_dir / 'psf_cube_etrue.fits')\n",
"\n",
"# Load the mean rmf calculated for the 4 Crab runs\n",
"rmf = EnergyDispersion.read(cube_dir / 'rmf.fits')"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Setup combined spatial and spectral model\n",
"spatial_model = NormGauss2DInt('spatial-model')\n",
"spectral_model = PowLaw1D('spectral-model')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Add the mean RMF to the Combine3DInt object \n",
"\n",
"The model is evaluated on the true energy bin then it is convolved by the energy dispersion to compare to the counts data in reconstructed energy"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"coord = counts_3d.sky_image_ref.coordinates(mode=\"edges\")\n",
"energies = counts_3d.energies(mode='edges').to(\"TeV\")\n",
"source_model = CombinedModel3DIntConvolveEdisp(\n",
" coord=coord,\n",
" energies=energies,\n",
" use_psf=True,\n",
" exposure=exposure_3d,\n",
" psf=psf_3d,\n",
" spatial_model=spatial_model,\n",
" spectral_model=spectral_model,\n",
" edisp=rmf.data.data,\n",
")\n",
"\n",
"# Set starting values\n",
"center = SkyCoord(83.633083, 22.0145, unit=\"deg\").galactic\n",
"source_model.gamma = 2.2\n",
"source_model.xpos = center.l.value\n",
"source_model.ypos = center.b.value\n",
"source_model.fwhm = 0.12\n",
"source_model.ampl = 1.0\n",
"\n",
"# Define the model\n",
"model = bkg + flux_factor * source_model"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Method = neldermead\n",
"Statistic = cash\n",
"Initial fit statistic = 6415.62\n",
"Final fit statistic = 5848.02 at function evaluation 704\n",
"Data points = 12500\n",
"Degrees of freedom = 12495\n",
"Change in statistic = 567.598\n",
" spatial-model.xpos 184.192 \n",
" spatial-model.ypos -6.17565 \n",
" spatial-model.ampl 6.22776 \n",
" spatial-model.fwhm 0.0712976 \n",
" spectral-model.gamma 2.26903 \n"
]
}
],
"source": [
"fit = Fit(\n",
" data=cube,\n",
" model=model,\n",
" stat=Cash(),\n",
" method=NelderMead(),\n",
" estmethod=Covariance(),\n",
")\n",
"fit_results = fit.fit()\n",
"print(fit_results.format())"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Confidence Method = covariance\n",
"Iterative Fit Method = None\n",
"Fitting Method = neldermead\n",
"Statistic = cash\n",
"covariance 1-sigma (68.2689%) bounds:\n",
" Param Best-Fit Lower Bound Upper Bound\n",
" ----- -------- ----------- -----------\n",
" spatial-model.xpos 184.192 -0.0117762 0.0117762\n",
" spatial-model.ypos -6.17565 -0.0136609 0.0136609\n",
" spatial-model.ampl 6.22776 -0.312187 0.312187\n",
" spatial-model.fwhm 0.0712976 -0.0205746 0.0205746\n",
" spectral-model.gamma 2.26903 -0.0515044 0.0515044\n"
]
}
],
"source": [
"err_est_results = fit.est_errors()\n",
"print(err_est_results.format())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Discussion\n",
"\n",
"Here the Cubes are constructed from dummy data. We don't expect to find the good spectral shape or amplitude since there is a lot of problem with the PFS, rmf etc.. in these dummy data.\n",
"On real data, we find the good Crab position and spectral shape for a power law or a power law with an exponential cutoff...."
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# TODO: there should be some visualisation of results.\n",
"# E.g. a plotted spectral model + butterfly\n",
"# Or a counts, model and residual image"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"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.5.2"
},
"nbsphinx": {
"orphan": true
}
},
"nbformat": 4,
"nbformat_minor": 2
}