{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Spectral analysis with Gammapy" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Introduction\n", "\n", "This notebook explains in detail how to use the classes in [gammapy.spectrum](https://docs.gammapy.org/dev/spectrum/index.html) and related ones. Note, that there is also [spectrum_pipe.ipynb](spectrum_pipe.ipynb) which explains how to do the analysis using a high-level interface. This notebook is aimed at advanced users who want to script taylor-made analysis pipelines and implement new methods.\n", "\n", "Based on a datasets of 4 Crab observations with H.E.S.S. (simulated events for now) we will perform a full region based spectral analysis, i.e. extracting source and background counts from certain \n", "regions, and fitting them using the forward-folding approach. We will use the following classes\n", "\n", "Data handling:\n", "\n", "* [gammapy.data.DataStore](https://docs.gammapy.org/dev/api/gammapy.data.DataStore.html)\n", "* [gammapy.data.DataStoreObservation](https://docs.gammapy.org/dev/api/gammapy.data.DataStoreObservation.html)\n", "* [gammapy.data.ObservationStats](https://docs.gammapy.org/dev/api/gammapy.data.ObservationStats.html)\n", "* [gammapy.data.ObservationSummary](https://docs.gammapy.org/dev/api/gammapy.data.ObservationSummary.html)\n", "\n", "To extract the 1-dim spectral information:\n", "\n", "* [gammapy.spectrum.SpectrumObservation](https://docs.gammapy.org/dev/api/gammapy.spectrum.SpectrumObservation.html)\n", "* [gammapy.spectrum.SpectrumExtraction](https://docs.gammapy.org/dev/api/gammapy.spectrum.SpectrumExtraction.html)\n", "* [gammapy.background.ReflectedRegionsBackgroundEstimator](https://docs.gammapy.org/dev/api/gammapy.background.ReflectedRegionsBackgroundEstimator.html)\n", "\n", "\n", "For the global fit (using Sherpa and WSTAT in the background):\n", "\n", "* [gammapy.spectrum.SpectrumFit](https://docs.gammapy.org/dev/api/gammapy.spectrum.SpectrumFit.html)\n", "* [gammapy.spectrum.models.PowerLaw](https://docs.gammapy.org/dev/api/gammapy.spectrum.models.PowerLaw.html)\n", "* [gammapy.spectrum.models.ExponentialCutoffPowerLaw](https://docs.gammapy.org/dev/api/gammapy.spectrum.models.ExponentialCutoffPowerLaw.html)\n", "* [gammapy.spectrum.models.LogParabola](https://docs.gammapy.org/dev/api/gammapy.spectrum.models.LogParabola.html)\n", "\n", "To compute flux points (a.k.a. \"SED\" = \"spectral energy distribution\")\n", "\n", "* [gammapy.spectrum.SpectrumResult](https://docs.gammapy.org/dev/api/gammapy.spectrum.SpectrumResult.html)\n", "* [gammapy.spectrum.FluxPoints](https://docs.gammapy.org/dev/api/gammapy.spectrum.FluxPoints.html)\n", "* [gammapy.spectrum.SpectrumEnergyGroupMaker](https://docs.gammapy.org/dev/api/gammapy.spectrum.SpectrumEnergyGroupMaker.html)\n", "* [gammapy.spectrum.FluxPointEstimator](https://docs.gammapy.org/dev/api/gammapy.spectrum.FluxPointEstimator.html)\n", "\n", "Feedback welcome!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Setup\n", "\n", "As usual, we'll start with some 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": [ "# Check package versions\n", "import gammapy\n", "import numpy as np\n", "import astropy\n", "import regions\n", "import sherpa\n", "\n", "print(\"gammapy:\", gammapy.__version__)\n", "print(\"numpy:\", np.__version__)\n", "print(\"astropy\", astropy.__version__)\n", "print(\"regions\", regions.__version__)\n", "print(\"sherpa\", sherpa.__version__)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import astropy.units as u\n", "from astropy.coordinates import SkyCoord, Angle\n", "from astropy.table import vstack as vstack_table\n", "from regions import CircleSkyRegion\n", "from gammapy.data import DataStore, Observations\n", "from gammapy.data import ObservationStats, ObservationSummary\n", "from gammapy.background.reflected import ReflectedRegionsBackgroundEstimator\n", "from gammapy.utils.energy import EnergyBounds\n", "from gammapy.spectrum import (\n", " SpectrumExtraction,\n", " SpectrumObservation,\n", " SpectrumFit,\n", " SpectrumResult,\n", ")\n", "from gammapy.spectrum.models import (\n", " PowerLaw,\n", " ExponentialCutoffPowerLaw,\n", " LogParabola,\n", ")\n", "from gammapy.spectrum import (\n", " FluxPoints,\n", " SpectrumEnergyGroupMaker,\n", " FluxPointEstimator,\n", ")\n", "from gammapy.maps import Map" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Configure logger\n", "\n", "Most high level classes in gammapy have the possibility to turn on logging or debug output. We well configure the logger in the following. For more info see https://docs.python.org/2/howto/logging.html#logging-basic-tutorial" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Setup the logger\n", "import logging\n", "\n", "logging.basicConfig()\n", "logging.getLogger(\"gammapy.spectrum\").setLevel(\"WARNING\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Load Data\n", "\n", "First, we select and load some H.E.S.S. observations of the Crab nebula (simulated events for now).\n", "\n", "We will access the events, effective area, energy dispersion, livetime and PSF for containement correction." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "datastore = DataStore.from_dir(\"$GAMMAPY_DATA/hess-dl3-dr1/\")\n", "obs_ids = [23523, 23526, 23559, 23592]\n", "observations = datastore.get_observations(obs_ids)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Define Target Region\n", "\n", "The next step is to define a signal extraction region, also known as on region. In the simplest case this is just a [CircleSkyRegion](http://astropy-regions.readthedocs.io/en/latest/api/regions.CircleSkyRegion.html#regions.CircleSkyRegion), but here we will use the ``Target`` class in gammapy that is useful for book-keeping if you run several analysis in a script." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "target_position = SkyCoord(ra=83.63, dec=22.01, unit=\"deg\", frame=\"icrs\")\n", "on_region_radius = Angle(\"0.11 deg\")\n", "on_region = CircleSkyRegion(center=target_position, radius=on_region_radius)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Create exclusion mask\n", "\n", "We will use the reflected regions method to place off regions to estimate the background level in the on region.\n", "To make sure the off regions don't contain gamma-ray emission, we create an exclusion mask.\n", "\n", "Using http://gamma-sky.net/ we find that there's only one known gamma-ray source near the Crab nebula: the AGN called [RGB J0521+212](http://gamma-sky.net/#/cat/tev/23) at GLON = 183.604 deg and GLAT = -8.708 deg." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "exclusion_region = CircleSkyRegion(\n", " center=SkyCoord(183.604, -8.708, unit=\"deg\", frame=\"galactic\"),\n", " radius=0.5 * u.deg,\n", ")\n", "\n", "skydir = target_position.galactic\n", "exclusion_mask = Map.create(\n", " npix=(150, 150), binsz=0.05, skydir=skydir, proj=\"TAN\", coordsys=\"GAL\"\n", ")\n", "\n", "mask = exclusion_mask.geom.region_mask([exclusion_region], inside=False)\n", "exclusion_mask.data = mask\n", "exclusion_mask.plot()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Estimate background\n", "\n", "Next we will manually perform a background estimate by placing [reflected regions](https://docs.gammapy.org/dev/background/reflected.html) around the pointing position and looking at the source statistics. This will result in a [gammapy.background.BackgroundEstimate](https://docs.gammapy.org/dev/api/gammapy.background.BackgroundEstimate.html) that serves as input for other classes in gammapy." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "background_estimator = ReflectedRegionsBackgroundEstimator(\n", " observations=observations,\n", " on_region=on_region,\n", " exclusion_mask=exclusion_mask,\n", ")\n", "\n", "background_estimator.run()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# print(background_estimator.result[0])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(8, 8))\n", "background_estimator.plot(add_legend=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Source statistic\n", "\n", "Next we're going to look at the overall source statistics in our signal region. For more info about what debug plots you can create check out the [ObservationSummary](https://docs.gammapy.org/dev/api/gammapy.data.ObservationSummary.html#gammapy.data.ObservationSummary) class." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "stats = []\n", "for obs, bkg in zip(observations, background_estimator.result):\n", " stats.append(ObservationStats.from_observation(obs, bkg))\n", "\n", "print(stats[1])\n", "\n", "obs_summary = ObservationSummary(stats)\n", "fig = plt.figure(figsize=(10, 6))\n", "ax1 = fig.add_subplot(121)\n", "\n", "obs_summary.plot_excess_vs_livetime(ax=ax1)\n", "ax2 = fig.add_subplot(122)\n", "obs_summary.plot_significance_vs_livetime(ax=ax2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Extract spectrum\n", "\n", "Now, we're going to extract a spectrum using the [SpectrumExtraction](https://docs.gammapy.org/dev/api/gammapy.spectrum.SpectrumExtraction.html) class. We provide the reconstructed energy binning we want to use. It is expected to be a Quantity with unit energy, i.e. an array with an energy unit. We use a utility function to create it. We also provide the true energy binning to use." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "e_reco = EnergyBounds.equal_log_spacing(0.1, 40, 40, unit=\"TeV\")\n", "e_true = EnergyBounds.equal_log_spacing(0.05, 100.0, 200, unit=\"TeV\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Instantiate a [SpectrumExtraction](https://docs.gammapy.org/dev/api/gammapy.spectrum.SpectrumExtraction.html) object that will do the extraction. The containment_correction parameter is there to allow for PSF leakage correction if one is working with full enclosure IRFs. We also compute a threshold energy and store the result in OGIP compliant files (pha, rmf, arf). This last step might be omitted though." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ANALYSIS_DIR = \"crab_analysis\"\n", "\n", "extraction = SpectrumExtraction(\n", " observations=observations,\n", " bkg_estimate=background_estimator.result,\n", " containment_correction=False,\n", ")\n", "extraction.run()\n", "\n", "# Add a condition on correct energy range in case it is not set by default\n", "extraction.compute_energy_threshold(method_lo=\"area_max\", area_percent_lo=10.0)\n", "\n", "print(extraction.spectrum_observations[0])\n", "# Write output in the form of OGIP files: PHA, ARF, RMF, BKG\n", "# extraction.run(observations=observations, bkg_estimate=background_estimator.result, outdir=ANALYSIS_DIR)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Look at observations\n", "\n", "Now we will look at the files we just created. We will use the [SpectrumObservation](https://docs.gammapy.org/dev/api/gammapy.spectrum.SpectrumObservation.html) object that are still in memory from the extraction step. Note, however, that you could also read them from disk if you have written them in the step above. The ``ANALYSIS_DIR`` folder contains 4 ``FITS`` files for each observation. These files are described in detail [here](https://gamma-astro-data-formats.readthedocs.io/en/latest/spectra/ogip/index.html). In short, they correspond to the on vector, the off vector, the effectie area, and the energy dispersion." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# filename = ANALYSIS_DIR + '/ogip_data/pha_obs23523.fits'\n", "# obs = SpectrumObservation.read(filename)\n", "\n", "# Requires IPython widgets\n", "# _ = extraction.spectrum_observations.peek()\n", "\n", "extraction.spectrum_observations[0].peek()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Fit spectrum\n", "\n", "Now we'll fit a global model to the spectrum. First we do a joint likelihood fit to all observations. If you want to stack the observations see below. We will also produce a debug plot in order to show how the global fit matches one of the individual observations." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model = PowerLaw(\n", " index=2, amplitude=2e-11 * u.Unit(\"cm-2 s-1 TeV-1\"), reference=1 * u.TeV\n", ")\n", "\n", "joint_fit = SpectrumFit(obs_list=extraction.spectrum_observations, model=model)\n", "joint_fit.run()\n", "joint_result = joint_fit.result" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ax0, ax1 = joint_result[0].plot(figsize=(8, 8))\n", "ax0.set_ylim(0, 20)\n", "print(joint_result[0])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Compute Flux Points\n", "\n", "To round up our analysis we can compute flux points by fitting the norm of the global model in energy bands. We'll use a fixed energy binning for now." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Define energy binning\n", "stacked_obs = extraction.spectrum_observations.stack()\n", "\n", "e_min, e_max = stacked_obs.lo_threshold.to_value(\"TeV\"), 30\n", "ebounds = np.logspace(np.log10(e_min), np.log10(e_max), 15) * u.TeV\n", "\n", "\n", "seg = SpectrumEnergyGroupMaker(obs=stacked_obs)\n", "seg.compute_groups_fixed(ebounds=ebounds)\n", "\n", "print(seg.groups)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fpe = FluxPointEstimator(\n", " obs=stacked_obs, groups=seg.groups, model=joint_result[0].model\n", ")\n", "flux_points = fpe.run()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "flux_points.table_formatted" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we plot the flux points and their likelihood profiles. For the plotting of upper limits we choose a threshold of TS < 4." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "flux_points.table[\"is_ul\"] = flux_points.table[\"ts\"] < 4\n", "ax = flux_points.plot(\n", " energy_power=2, flux_unit=\"erg-1 cm-2 s-1\", color=\"darkorange\"\n", ")\n", "flux_points.to_sed_type(\"e2dnde\").plot_likelihood(ax=ax)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The final plot with the best fit model, flux points and residuals can be quickly made like this" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "spectrum_result = SpectrumResult(\n", " points=flux_points, model=joint_result[0].model\n", ")\n", "ax0, ax1 = spectrum_result.plot(\n", " energy_range=joint_fit.result[0].fit_range,\n", " energy_power=2,\n", " flux_unit=\"erg-1 cm-2 s-1\",\n", " fig_kwargs=dict(figsize=(8, 8)),\n", ")\n", "\n", "ax0.set_xlim(0.4, 50)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Stack observations\n", "\n", "And alternative approach to fitting the spectrum is stacking all observations first and the fitting a model to the stacked observation. This works as follows. A comparison to the joint likelihood fit is also printed." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "stacked_obs = extraction.spectrum_observations.stack()\n", "\n", "stacked_fit = SpectrumFit(obs_list=stacked_obs, model=model)\n", "stacked_fit.run()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "stacked_result = stacked_fit.result\n", "print(stacked_result[0])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "stacked_table = stacked_result[0].to_table(format=\".3g\")\n", "stacked_table[\"method\"] = \"stacked\"\n", "joint_table = joint_result[0].to_table(format=\".3g\")\n", "joint_table[\"method\"] = \"joint\"\n", "total_table = vstack_table([stacked_table, joint_table])\n", "print(\n", " total_table[\"method\", \"index\", \"index_err\", \"amplitude\", \"amplitude_err\"]\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercises\n", "\n", "Some things we might do:\n", "\n", "- Fit a different spectral model (ECPL or CPL or ...)\n", "- Use different method or parameters to compute the flux points\n", "- Do a chi^2 fit to the flux points and compare\n", "\n", "TODO: give pointers how to do this (and maybe write a notebook with solutions)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Start exercises here" ] } ], "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.6.0" } }, "nbformat": 4, "nbformat_minor": 2 }