{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Spectrum simulation for CTA\n",
    "\n",
    "A quick example how to simulate and fit a spectrum for the [Cherenkov Telescope Array (CTA)](https://www.cta-observatory.org).\n",
    "\n",
    "We will use the following classes:\n",
    "\n",
    "* [gammapy.spectrum.SpectrumObservation](https://docs.gammapy.org/dev/api/gammapy.spectrum.SpectrumObservation.html)\n",
    "* [gammapy.spectrum.SpectrumSimulation](https://docs.gammapy.org/dev/api/gammapy.spectrum.SpectrumSimulation.html)\n",
    "* [gammapy.spectrum.SpectrumFit](https://docs.gammapy.org/dev/api/gammapy.spectrum.SpectrumFit.html)\n",
    "* [gammapy.irf.load_cta_irfs](https://docs.gammapy.org/dev/api/gammapy.irf.load_cta_irfs.html)"
   ]
  },
  {
   "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.irf import EnergyDispersion, EffectiveAreaTable\n",
    "from gammapy.spectrum import SpectrumSimulation, SpectrumFit\n",
    "from gammapy.spectrum.models import PowerLaw\n",
    "from gammapy.irf import load_cta_irfs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define simulation parameters parameters\n",
    "livetime = 1 * u.h\n",
    "offset = 0.5 * u.deg\n",
    "# Energy from 0.1 to 100 TeV with 10 bins/decade\n",
    "energy = np.logspace(-1, 2, 31) * u.TeV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define spectral model\n",
    "model = PowerLaw(\n",
    "    index=2.1,\n",
    "    amplitude=2.5e-12 * u.Unit(\"cm-2 s-1 TeV-1\"),\n",
    "    reference=1 * u.TeV,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load IRFs\n",
    "filename = (\n",
    "    \"$GAMMAPY_DATA/cta-1dc/caldb/data/cta/1dc/bcf/South_z20_50h/irf_file.fits\"\n",
    ")\n",
    "cta_irf = load_cta_irfs(filename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "aeff = cta_irf[\"aeff\"].to_effective_area_table(offset=offset, energy=energy)\n",
    "aeff.plot()\n",
    "plt.loglog()\n",
    "print(cta_irf[\"aeff\"].data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "edisp = cta_irf[\"edisp\"].to_energy_dispersion(\n",
    "    offset=offset, e_true=energy, e_reco=energy\n",
    ")\n",
    "edisp.plot_matrix()\n",
    "print(edisp.data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Simulate data\n",
    "sim = SpectrumSimulation(\n",
    "    aeff=aeff, edisp=edisp, source_model=model, livetime=livetime\n",
    ")\n",
    "sim.simulate_obs(seed=42, obs_id=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.obs.peek()\n",
    "print(sim.obs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Spectral analysis\n",
    "\n",
    "Now that we have some simulated CTA counts spectrum, let's analyse it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fit data\n",
    "fit = SpectrumFit(obs_list=sim.obs, model=model, stat=\"cash\")\n",
    "fit.run()\n",
    "result = fit.result[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "energy_range = [0.1, 100] * u.TeV\n",
    "model.plot(energy_range=energy_range, energy_power=2)\n",
    "result.model.plot(energy_range=energy_range, energy_power=2)\n",
    "result.model.plot_error(energy_range=energy_range, energy_power=2);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "* Change the observation time to something longer or shorter. Does the observation and spectrum results change as you expected?\n",
    "* Change the spectral model, e.g. add a cutoff at 5 TeV, or put a steep-spectrum source with spectral index of 4.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Start the exercises here!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What next?\n",
    "\n",
    "In this tutorial we simulated and analysed the spectrum of source using CTA prod 2 IRFs.\n",
    "\n",
    "If you'd like to go further, please see the other tutorial notebooks."
   ]
  }
 ],
 "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
}