tutorials/analysis_3d.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3D analysis\n",
    "\n",
    "This tutorial shows how to run a 3D map-based analysis using three example observations of the Galactic center region with CTA."
   ]
  },
  {
   "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 os\n",
    "from pathlib import Path\n",
    "import numpy as np\n",
    "import astropy.units as u\n",
    "from astropy.coordinates import SkyCoord\n",
    "from gammapy.data import DataStore\n",
    "from gammapy.irf import EnergyDispersion, make_mean_psf, make_mean_edisp\n",
    "from gammapy.maps import WcsGeom, MapAxis, Map, WcsNDMap\n",
    "from gammapy.cube import MapMaker, MapEvaluator, PSFKernel, MapDataset\n",
    "from gammapy.cube.models import SkyModel, SkyDiffuseCube, BackgroundModel\n",
    "from gammapy.spectrum.models import PowerLaw, ExponentialCutoffPowerLaw\n",
    "from gammapy.image.models import SkyGaussian, SkyPointSource\n",
    "from gammapy.utils.fitting import Fit\n",
    "from regions import CircleSkyRegion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!gammapy info --no-system"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prepare modeling input data\n",
    "\n",
    "### Prepare input maps\n",
    "\n",
    "We first use the `DataStore` object to access the CTA observations and retrieve a list of observations by passing the observations IDs to the `.get_observations()` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define which data to use and print some information\n",
    "data_store = DataStore.from_dir(\"$GAMMAPY_DATA/cta-1dc/index/gps/\")\n",
    "data_store.info()\n",
    "print(\n",
    "    \"Total observation time (hours): \",\n",
    "    data_store.obs_table[\"ONTIME\"].sum() / 3600,\n",
    ")\n",
    "print(\"Observation table: \", data_store.obs_table.colnames)\n",
    "print(\"HDU table: \", data_store.hdu_table.colnames)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select some observations from these dataset by hand\n",
    "obs_ids = [110380, 111140, 111159]\n",
    "observations = data_store.get_observations(obs_ids)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we define a reference geometry for our analysis, We choose a WCS based gemoetry with a binsize of 0.02 deg and also define an energy axis: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "energy_axis = MapAxis.from_edges(\n",
    "    np.logspace(-1.0, 1.0, 10), unit=\"TeV\", name=\"energy\", interp=\"log\"\n",
    ")\n",
    "geom = WcsGeom.create(\n",
    "    skydir=(0, 0),\n",
    "    binsz=0.02,\n",
    "    width=(10, 8),\n",
    "    coordsys=\"GAL\",\n",
    "    proj=\"CAR\",\n",
    "    axes=[energy_axis],\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `MapMaker` object is initialized with this reference geometry and a field of view cut of 4 deg:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "maker = MapMaker(geom, offset_max=4.0 * u.deg)\n",
    "maps = maker.run(observations)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The maps are prepared by calling the `.run()` method and passing the `observations`. The `.run()` method returns a Python `dict` containing a `counts`, `background` and `exposure` map:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(maps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is what the summed counts image looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "counts = maps[\"counts\"].sum_over_axes()\n",
    "counts.smooth(width=0.1 * u.deg).plot(stretch=\"sqrt\", add_cbar=True, vmax=6);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the background image:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "background = maps[\"background\"].sum_over_axes()\n",
    "background.smooth(width=0.1 * u.deg).plot(\n",
    "    stretch=\"sqrt\", add_cbar=True, vmax=6\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And this one the exposure image:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "exposure = maps[\"exposure\"].sum_over_axes()\n",
    "exposure.smooth(width=0.1 * u.deg).plot(stretch=\"sqrt\", add_cbar=True);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also compute an excess image just with  a few lines of code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "excess = counts - background\n",
    "excess.smooth(5).plot(stretch=\"sqrt\", add_cbar=True);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a more realistic excess plot we can also take into account the diffuse galactic emission. For this tutorial we will load a Fermi diffuse model map that represents a small cutout for the Galactic center region:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diffuse_gal = Map.read(\"$GAMMAPY_DATA/fermi-3fhl-gc/gll_iem_v06_gc.fits.gz\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Diffuse image: \", diffuse_gal.geom)\n",
    "print(\"counts: \", maps[\"counts\"].geom)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the geometry of the images is completely different, so we need to apply our geometric configuration to the diffuse emission file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "coord = maps[\"counts\"].geom.get_coord()\n",
    "\n",
    "data = diffuse_gal.interp_by_coord(\n",
    "    {\n",
    "        \"skycoord\": coord.skycoord,\n",
    "        \"energy\": coord[\"energy\"]\n",
    "        * maps[\"counts\"].geom.get_axis_by_name(\"energy\").unit,\n",
    "    },\n",
    "    interp=3,\n",
    ")\n",
    "diffuse_galactic = WcsNDMap(maps[\"counts\"].geom, data)\n",
    "print(\"Before: \\n\", diffuse_gal.geom)\n",
    "print(\"Now (same as maps): \\n\", diffuse_galactic.geom)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# diffuse_galactic.slice_by_idx({\"energy\": 0}).plot(add_cbar=True); # this can be used to check image at different energy bins\n",
    "diffuse = diffuse_galactic.sum_over_axes()\n",
    "diffuse.smooth(5).plot(stretch=\"sqrt\", add_cbar=True)\n",
    "print(diffuse)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now multiply the exposure for this diffuse emission to subtract the result from the counts along with the background."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "combination = diffuse * exposure\n",
    "combination.unit = \"\"\n",
    "combination.smooth(5).plot(stretch=\"sqrt\", add_cbar=True);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot then the excess image subtracting now the effect of the diffuse galactic emission."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "excess2 = counts - background - combination\n",
    "\n",
    "fig, axs = plt.subplots(1, 2, figsize=(15, 5))\n",
    "\n",
    "axs[0].set_title(\"With diffuse emission subtraction\")\n",
    "axs[1].set_title(\"Without diffuse emission subtraction\")\n",
    "excess2.smooth(5).plot(\n",
    "    cmap=\"coolwarm\", vmin=-1, vmax=1, add_cbar=True, ax=axs[0]\n",
    ")\n",
    "excess.smooth(5).plot(\n",
    "    cmap=\"coolwarm\", vmin=-1, vmax=1, add_cbar=True, ax=axs[1]\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Prepare IRFs\n",
    "\n",
    "To estimate the mean PSF across all observations at a given source position `src_pos`, we use `make_mean_psf()`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# mean PSF\n",
    "src_pos = SkyCoord(0, 0, unit=\"deg\", frame=\"galactic\")\n",
    "table_psf = make_mean_psf(observations, src_pos)\n",
    "\n",
    "# PSF kernel used for the model convolution\n",
    "psf_kernel = PSFKernel.from_table_psf(table_psf, geom, max_radius=\"0.3 deg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To estimate the mean energy dispersion across all observations at a given source position `src_pos`, we use `make_mean_edisp()`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define energy grid\n",
    "energy = energy_axis.edges * energy_axis.unit\n",
    "\n",
    "# mean edisp\n",
    "edisp = make_mean_edisp(\n",
    "    observations, position=src_pos, e_true=energy, e_reco=energy\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Save maps and IRFs to disk\n",
    "\n",
    "It is common to run the preparation step independent of the likelihood fit, because often the preparation of maps, PSF and energy dispersion is slow if you have a lot of data. We first create a folder:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "path = Path(\"analysis_3d\")\n",
    "path.mkdir(exist_ok=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And then write the maps and IRFs to disk by calling the dedicated `.write()` methods:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# write maps\n",
    "maps[\"counts\"].write(str(path / \"counts.fits\"), overwrite=True)\n",
    "maps[\"background\"].write(str(path / \"background.fits\"), overwrite=True)\n",
    "maps[\"exposure\"].write(str(path / \"exposure.fits\"), overwrite=True)\n",
    "\n",
    "# write IRFs\n",
    "psf_kernel.write(str(path / \"psf.fits\"), overwrite=True)\n",
    "edisp.write(str(path / \"edisp.fits\"), overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Likelihood fit\n",
    "\n",
    "### Reading maps and IRFs\n",
    "As first step we read in the maps and IRFs that we have saved to disk again:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read maps\n",
    "maps = {\n",
    "    \"counts\": Map.read(str(path / \"counts.fits\")),\n",
    "    \"background\": Map.read(str(path / \"background.fits\")),\n",
    "    \"exposure\": Map.read(str(path / \"exposure.fits\")),\n",
    "}\n",
    "\n",
    "# read IRFs\n",
    "psf_kernel = PSFKernel.read(str(path / \"psf.fits\"))\n",
    "edisp = EnergyDispersion.read(str(path / \"edisp.fits\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's cut out only part of the maps, so that we the fitting step does not take so long (we go from left to right one):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmaps = {\n",
    "    name: m.cutout(SkyCoord(0, 0, unit=\"deg\", frame=\"galactic\"), 2 * u.deg)\n",
    "    for name, m in maps.items()\n",
    "}\n",
    "cmaps[\"counts\"].sum_over_axes().plot(stretch=\"sqrt\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Insted of the complete one, which was:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "counts.plot(stretch=\"sqrt\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fit mask\n",
    "\n",
    "To select a certain spatial region and/or energy range for the fit we can create a fit mask:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask = Map.from_geom(cmaps[\"counts\"].geom)\n",
    "\n",
    "region = CircleSkyRegion(center=src_pos, radius=0.6 * u.deg)\n",
    "mask.data = mask.geom.region_mask([region])\n",
    "\n",
    "mask.get_image_by_idx((0,)).plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition we also exclude the range below 0.3 TeV for the fit:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "coords = mask.geom.get_coord()\n",
    "mask.data &= coords[\"energy\"] > 0.3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model fit\n",
    "\n",
    "No we are ready for the actual likelihood fit. We first define the model as a combination of a point source with a powerlaw:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spatial_model = SkyPointSource(lon_0=\"0.01 deg\", lat_0=\"0.01 deg\")\n",
    "spectral_model = PowerLaw(\n",
    "    index=2.2, amplitude=\"3e-12 cm-2 s-1 TeV-1\", reference=\"1 TeV\"\n",
    ")\n",
    "model = SkyModel(spatial_model=spatial_model, spectral_model=spectral_model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Defining a background model\n",
    "\n",
    "Often, it is useful to fit the normalisation (and also the tilt) of the background. To do so, we have to define the background as a model. In this example, we will keep the tilt fixed and the norm free."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "background_model = BackgroundModel(cmaps[\"background\"], norm=1.1, tilt=0.0)\n",
    "background_model.parameters['norm'].frozen = False\n",
    "background_model.parameters['tilt'].frozen = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we set up the `MapDataset` object by passing the prepared maps, IRFs as well as the model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset = MapDataset(\n",
    "    model=model,\n",
    "    counts=cmaps[\"counts\"],\n",
    "    exposure=cmaps[\"exposure\"],\n",
    "    background_model=background_model,\n",
    "    mask=mask,\n",
    "    psf=psf_kernel,\n",
    "    edisp=edisp,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "No we run the model fit:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "fit = Fit(dataset)\n",
    "result = fit.run(optimize_opts={\"print_level\": 1})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result.parameters.to_table()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the background normalisation is very large (2.15). This is because of the large number of sources, and the presence of diffuse emission in the galactic center."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Check model fit\n",
    "\n",
    "We check the model fit by computing a residual image. For this we first get the number of predicted counts:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "npred = dataset.npred()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And compute a residual image:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "residual = cmaps[\"counts\"] - npred"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "residual.sum_over_axes().smooth(width=0.05 * u.deg).plot(\n",
    "    cmap=\"coolwarm\", vmin=-3, vmax=3, add_cbar=True\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also plot the best fit spectrum. For that need to extract the covariance of the spectral parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spec = model.spectral_model\n",
    "\n",
    "# set covariance on the spectral model\n",
    "covariance = result.parameters.covariance\n",
    "spec.parameters.covariance = covariance[2:5, 2:5]\n",
    "\n",
    "energy_range = [0.3, 10] * u.TeV\n",
    "spec.plot(energy_range=energy_range, energy_power=2)\n",
    "ax = spec.plot_error(energy_range=energy_range, energy_power=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Apparently our model should be improved by adding a component for diffuse Galactic emission and at least one second point\n",
    "source. But before we do that in the next section, we will fit the background as a model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Add Galactic diffuse emission to model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use both models at the same time, our diffuse model (the same from the Fermi file used before) and our model for the central source. This time, in order to make it more realistic, we will consider an exponential cut off power law spectral model for the source  (note that we are not constraining the fit with any mask this time). We will again fot the normalisation of the background."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diffuse_model = SkyDiffuseCube.read(\n",
    "    \"$GAMMAPY_DATA/fermi-3fhl-gc/gll_iem_v06_gc.fits.gz\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spatial_model = SkyPointSource(lon_0=\"0.01 deg\", lat_0=\"0.01 deg\")\n",
    "spectral_model = ExponentialCutoffPowerLaw(\n",
    "    index=2 * u.Unit(\"\"),\n",
    "    amplitude=1e-12 * u.Unit(\"cm-2 s-1 TeV-1\"),\n",
    "    reference=1.0 * u.TeV,\n",
    "    lambda_=1 / u.TeV,\n",
    ")\n",
    "model_ecpl = SkyModel(\n",
    "    spatial_model=spatial_model, spectral_model=spectral_model\n",
    ")\n",
    "\n",
    "background_model = BackgroundModel(cmaps[\"background\"], norm=1.0, tilt=0.0)\n",
    "background_model.parameters['norm'].frozen = False\n",
    "background_model.parameters['tilt'].frozen = True\n",
    "\n",
    "model_combined = diffuse_model + model_ecpl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset_combined = MapDataset(\n",
    "    model=model_combined,\n",
    "    counts=cmaps[\"counts\"],\n",
    "    exposure=cmaps[\"exposure\"],\n",
    "    background_model=background_model,\n",
    "    psf=psf_kernel,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "fit_combined = Fit(dataset_combined)\n",
    "result_combined = fit_combined.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(model_ecpl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see we have now two components in our model, and we can access them separately."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Checking normalization value (the closer to 1 the better)\n",
    "print(model_combined)\n",
    "print(\"Background model: {}\\n\".format(background_model))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see that the normalisation of the background has vastly improved"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now plot the residual image considering this improved model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "residual2 = cmaps[\"counts\"] - dataset_combined.npred()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just as a comparison, we can plot our previous residual map (left) and the new one (right) with the same scale:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(15, 5))\n",
    "ax_1 = plt.subplot(121, projection=residual.geom.wcs)\n",
    "ax_2 = plt.subplot(122, projection=residual.geom.wcs)\n",
    "\n",
    "ax_1.set_title(\"Without diffuse emission subtraction\")\n",
    "ax_2.set_title(\"With diffuse emission subtraction\")\n",
    "\n",
    "residual.sum_over_axes().smooth(width=0.05 * u.deg).plot(\n",
    "    cmap=\"coolwarm\", vmin=-2, vmax=2, add_cbar=True, ax=ax_1\n",
    ")\n",
    "residual2.sum_over_axes().smooth(width=0.05 * u.deg).plot(\n",
    "    cmap=\"coolwarm\", vmin=-2, vmax=2, add_cbar=True, ax=ax_2\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally we can check again our model (including now the diffuse emission):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spec_ecpl = model_ecpl.spectral_model\n",
    "ax = spec_ecpl.plot(energy_range=energy_range, energy_power=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Results seems to be better (but not perfect yet). Next step to improve our model even more would be getting rid of the other bright source (G0.9+0.1)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that this notebook aims to show you the procedure of a 3D analysis using just a few observations and a cutted Fermi model. Results get much better for a more complete analysis considering the GPS dataset from the CTA First Data Challenge (DC-1) and also the CTA model for the Galactic diffuse emission, as shown in the next image:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](images/DC1_3d.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The complete tutorial notebook of this analysis is available to be downloaded in [GAMMAPY-EXTRA](https://github.com/gammapy/gammapy-extra) repository at https://github.com/gammapy/gammapy-extra/blob/master/analyses/cta_1dc_gc_3d.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "* Analyse the second source in the field of view: G0.9+0.1 and add it to the combined model."
   ]
  }
 ],
 "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.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}