{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Flux point fitting in Gammapy\n",
    "\n",
    "\n",
    "## Introduction\n",
    "\n",
    "In this tutorial we're going to learn how to fit spectral models to combined Fermi-LAT and IACT flux points.\n",
    "\n",
    "The central class we're going to use for this example analysis is:  \n",
    "\n",
    "- [gammapy.spectrum.FluxPointsDataset](https://docs.gammapy.org/dev/spectrum/index.html#reference-api)\n",
    "\n",
    "In addition we will work with the following data classes:\n",
    "\n",
    "- [gammapy.spectrum.FluxPoints](https://docs.gammapy.org/dev/api/gammapy.spectrum.FluxPoints.html)\n",
    "- [gammapy.catalog.SourceCatalogGammaCat](https://docs.gammapy.org/dev/api/gammapy.catalog.SourceCatalogGammaCat.html)\n",
    "- [gammapy.catalog.SourceCatalog3FHL](https://docs.gammapy.org/dev/api/gammapy.catalog.SourceCatalog3FHL.html)\n",
    "- [gammapy.catalog.SourceCatalog3FGL](https://docs.gammapy.org/dev/api/gammapy.catalog.SourceCatalog3FGL.html)\n",
    "\n",
    "And the following spectral model classes:\n",
    "\n",
    "- [PowerLaw](https://docs.gammapy.org/dev/api/gammapy.spectrum.models.PowerLaw.html)\n",
    "- [ExponentialCutoffPowerLaw](https://docs.gammapy.org/dev/api/gammapy.spectrum.models.ExponentialCutoffPowerLaw.html)\n",
    "- [LogParabola](https://docs.gammapy.org/dev/api/gammapy.spectrum.models.LogParabola.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup\n",
    "\n",
    "Let us start with the usual IPython notebook and Python imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy import units as u\n",
    "from gammapy.spectrum.models import (\n",
    "    PowerLaw,\n",
    "    ExponentialCutoffPowerLaw,\n",
    "    LogParabola,\n",
    ")\n",
    "from gammapy.spectrum import FluxPointsDataset, FluxPoints\n",
    "from gammapy.catalog import (\n",
    "    SourceCatalog3FGL,\n",
    "    SourceCatalogGammaCat,\n",
    "    SourceCatalog3FHL,\n",
    ")\n",
    "from gammapy.utils.fitting import Fit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load spectral points\n",
    "\n",
    "For this analysis we choose to work with the source 'HESS J1507-622' and the associated Fermi-LAT sources '3FGL J1506.6-6219' and '3FHL J1507.9-6228e'. We load the source catalogs, and then access source of interest by name:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fermi_3fgl = SourceCatalog3FGL()\n",
    "fermi_3fhl = SourceCatalog3FHL()\n",
    "gammacat = SourceCatalogGammaCat(\"$GAMMAPY_DATA/gamma-cat/gammacat.fits.gz\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "source_gammacat = gammacat[\"HESS J1507-622\"]\n",
    "source_fermi_3fgl = fermi_3fgl[\"3FGL J1506.6-6219\"]\n",
    "source_fermi_3fhl = fermi_3fhl[\"3FHL J1507.9-6228e\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The corresponding flux points data can be accessed with `.flux_points` attribute:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "flux_points_gammacat = source_gammacat.flux_points\n",
    "flux_points_gammacat.table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the Fermi-LAT catalogs, integral flux points are given. Currently the flux point fitter only works with differential flux points, so we apply the conversion here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "flux_points_3fgl = source_fermi_3fgl.flux_points.to_sed_type(\n",
    "    sed_type=\"dnde\", model=source_fermi_3fgl.spectral_model\n",
    ")\n",
    "flux_points_3fhl = source_fermi_3fhl.flux_points.to_sed_type(\n",
    "    sed_type=\"dnde\", model=source_fermi_3fhl.spectral_model\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally we stack the flux points into a single `FluxPoints` object and drop the upper limit values, because currently we can't handle them in the fit:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# stack flux point tables\n",
    "flux_points = FluxPoints.stack(\n",
    "    [flux_points_gammacat, flux_points_3fhl, flux_points_3fgl]\n",
    ")\n",
    "\n",
    "# drop the flux upper limit values\n",
    "flux_points = flux_points.drop_ul()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Power Law Fit\n",
    "\n",
    "First we start with fitting a simple [power law](https://docs.gammapy.org/dev/api/gammapy.spectrum.models.PowerLaw.html#gammapy.spectrum.models.PowerLaw)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pwl = PowerLaw(index=2, amplitude=\"1e-12 cm-2 s-1 TeV-1\", reference=\"1 TeV\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After creating the model we run the fit by passing the `'flux_points'` and `'pwl'` objects:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset_pwl = FluxPointsDataset(pwl, flux_points, likelihood=\"chi2assym\")\n",
    "fitter = Fit(dataset_pwl)\n",
    "result_pwl = fitter.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And print the result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(result_pwl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(pwl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally we plot the data points and the best fit model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = flux_points.plot(energy_power=2)\n",
    "pwl.plot(energy_range=[1e-4, 1e2] * u.TeV, ax=ax, energy_power=2)\n",
    "\n",
    "# assign covariance for plotting\n",
    "pwl.parameters.covariance = result_pwl.parameters.covariance\n",
    "\n",
    "pwl.plot_error(\n",
    "    energy_range=[1e-4, 1e2] * u.TeV, ax=ax, energy_power=2\n",
    ")\n",
    "ax.set_ylim(1e-13, 1e-11);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exponential Cut-Off Powerlaw Fit\n",
    "\n",
    "Next we fit an [exponential cut-off power](https://docs.gammapy.org/dev/api/gammapy.spectrum.models.ExponentialCutoffPowerLaw.html#gammapy.spectrum.models.ExponentialCutoffPowerLaw) law to the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ecpl = ExponentialCutoffPowerLaw(\n",
    "    index=1.8,\n",
    "    amplitude=\"2e-12 cm-2 s-1 TeV-1\",\n",
    "    reference=\"1 TeV\",\n",
    "    lambda_=\"0.1 TeV-1\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We run the fitter again by passing the flux points and the `ecpl` model instance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset_ecpl = FluxPointsDataset(ecpl, flux_points, likelihood=\"chi2assym\")\n",
    "fitter = Fit(dataset_ecpl)\n",
    "result_ecpl = fitter.run()\n",
    "print(ecpl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We plot the data and best fit model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = flux_points.plot(energy_power=2)\n",
    "ecpl.plot(energy_range=[1e-4, 1e2] * u.TeV, ax=ax, energy_power=2)\n",
    "\n",
    "# assign covariance for plotting\n",
    "ecpl.parameters.covariance = result_ecpl.parameters.covariance\n",
    "\n",
    "ecpl.plot_error(\n",
    "    energy_range=[1e-4, 1e2] * u.TeV, ax=ax, energy_power=2\n",
    ")\n",
    "ax.set_ylim(1e-13, 1e-11)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Log-Parabola Fit\n",
    "\n",
    "Finally we try to fit a [log-parabola](https://docs.gammapy.org/dev/api/gammapy.spectrum.models.LogParabola.html#gammapy.spectrum.models.LogParabola) model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "log_parabola = LogParabola(\n",
    "    alpha=2, amplitude=\"1e-12 cm-2 s-1 TeV-1\", reference=\"1 TeV\", beta=0.1\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset_log_parabola = FluxPointsDataset(log_parabola, flux_points, likelihood=\"chi2assym\")\n",
    "fitter = Fit(dataset_log_parabola)\n",
    "result_log_parabola = fitter.run()\n",
    "print(log_parabola)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = flux_points.plot(energy_power=2)\n",
    "log_parabola.plot(\n",
    "    energy_range=[1e-4, 1e2] * u.TeV, ax=ax, energy_power=2\n",
    ")\n",
    "\n",
    "# assign covariance for plotting\n",
    "log_parabola.parameters.covariance = result_log_parabola.parameters.covariance\n",
    "\n",
    "log_parabola.plot_error(\n",
    "    energy_range=[1e-4, 1e2] * u.TeV, ax=ax, energy_power=2\n",
    ")\n",
    "ax.set_ylim(1e-13, 1e-11);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "- Fit a `PowerLaw2` and `ExponentialCutoffPowerLaw3FGL` to the same data.\n",
    "- Fit a `ExponentialCutoffPowerLaw` model to Vela X ('HESS J0835-455') only and check if the best fit values correspond to the values given in the Gammacat catalog"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What next?\n",
    "\n",
    "This was an introduction to SED fitting in Gammapy.\n",
    "\n",
    "* If you would like to learn how to perform a full Poisson maximum likelihood spectral fit, please check out the [spectrum pipe](spectrum_pipe.ipynb) tutorial.\n",
    "* If you are interested in simulation of spectral data in the context of CTA, please check out the [spectrum simulation cta](spectrum_simulation_cta.ipynb) notebook.\n",
    "* To learn more about other parts of Gammapy (e.g. Fermi-LAT and TeV data analysis), check out the other tutorial notebooks.\n",
    "* To see what's available in Gammapy, browse the [Gammapy docs](https://docs.gammapy.org/) or use the full-text search.\n",
    "* If you have any questions, ask on the mailing list ."
   ]
  }
 ],
 "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
}