3D_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# **A full 3D analysis using the low-level Gammapy API**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Objective: Performing a full 3D anaysis of the extended source [MSH 15-52](http://tevcat.uchicago.edu/?mode=1;id=95)**\n",
    "\n",
    "In practice, we have to:\n",
    "- Prepare the **data access and selection**\n",
    "  - Create a `~gammapy.data.DataStore` poiting to the relevant data \n",
    "  - Apply an observation selection to produce a list of observations, a `~gammapy.data.Observations` object.\n",
    "- Set up the **analyis parameters**\n",
    "  - Define a geometry of the Map we want to produce, with a sky projection and an energy range.\n",
    "    - Create a `~gammapy.maps.MapAxis` for the energy\n",
    "    - Create a `~gammapy.maps.WcsGeom` for the geometry\n",
    "    - Define the exclusion mask\n",
    "    - Choose the correct ~gammapy.datasets.Dataset type and define it\n",
    "- Do the **data reduction**\n",
    "  - Create the necessary makers : \n",
    "    - the map dataset maker : `~gammapy.makers.MapDatasetMaker`\n",
    "    - the [background normalization](https://docs.gammapy.org/dev/user-guide/makers/fov.html) maker, here a `~gammapy.makers.FoVBackgroundMaker`\n",
    "    - and usually the safe range maker : `~gammapy.makers.SafeMaskMaker`\n",
    "  - Perform the data reduction loop. And for every observation:\n",
    "    - Apply the makers sequentially to produce the current `~gammapy.datasets.MapDataset`\n",
    "    - Stack it on the target one.\n",
    "- Make the **modeling and fitting**\n",
    "  - Define the `~gammapy.modeling.models.SkyModel` to apply to the dataset.\n",
    "  - Create a `~gammapy.modeling.Fit` object and run it to fit the model parameters\n",
    "  - Apply a `~gammapy.estimators.FluxPointsEstimator` to compute flux points for the spectral part of the fit.\n",
    "\n",
    "## Setup\n",
    "First, we setup the analysis by performing required imports."
   ]
  },
  {
   "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": [
    "from pathlib import Path\n",
    "import numpy as np\n",
    "import logging\n",
    "from astropy import units as u\n",
    "from astropy.coordinates import SkyCoord\n",
    "from regions import CircleSkyRegion\n",
    "from scipy.stats import norm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gammapy.data import DataStore\n",
    "from gammapy.datasets import MapDataset\n",
    "from gammapy.maps import WcsGeom, MapAxis\n",
    "from gammapy.makers import MapDatasetMaker, SafeMaskMaker, FoVBackgroundMaker\n",
    "from gammapy.modeling.models import (\n",
    "    SkyModel,\n",
    "    PowerLawSpectralModel,\n",
    "    PointSpatialModel,\n",
    "    FoVBackgroundModel,\n",
    "    GaussianSpatialModel,\n",
    "    Models\n",
    ")\n",
    "from gammapy.modeling import Fit\n",
    "from gammapy.estimators import FluxPointsEstimator, ExcessMapEstimator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Optional set-up"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "logging.basicConfig()    \n",
    "log = logging.getLogger(\"1Danalysis\")\n",
    "log.setLevel(logging.WARNING) #INFO, WARNING, DEBUG\n",
    "\n",
    "from astropy.io.fits.verify import VerifyWarning\n",
    "import warnings\n",
    "with warnings.catch_warnings():\n",
    "    warnings.simplefilter('ignore', VerifyWarning)\n",
    "\n",
    "from gammapy.utils import pbar\n",
    "pbar.SHOW_PROGRESS_BAR = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Defining the datastore and selecting observations\n",
    "\n",
    "We first use the `~gammapy.data.DataStore` object to access the observations we want to analyse, here the H.E.S.S. DL3 DR1. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_store = DataStore.from_dir(\"$GAMMAPY_DATA/hess-dl3-dr1\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pos = SkyCoord(228.32083333, -59.08166667, unit=u.deg, frame=\"icrs\")\n",
    "pos"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now define an observation filter to select only the relevant observations. \n",
    "Here we use a cone search which we define with a python dict.\n",
    "\n",
    "We then filter the `ObservationTable` with `~gammapy.data.ObservationTable.select_observations()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "obs_table_filtered = data_store.obs_table.select_sky_circle(center=pos, radius=2 * u.deg)\n",
    "obs_ids = obs_table_filtered[\"OBS_ID\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now retrieve the relevant observations by passing their `obs_id` to the`~gammapy.data.DataStore.get_observations()` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "observations = data_store.get_observations(obs_ids)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(observations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "obs = observations[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "obs.events.peek()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preparing reduced datasets geometry\n",
    "\n",
    "Now we define a reference geometry for our analysis, We choose a WCS based geometry 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_energy_bounds(0.3, 10.0, 15, unit=\"TeV\")\n",
    "\n",
    "geom = WcsGeom.create(\n",
    "    skydir=(pos.ra.value, pos.dec.value),\n",
    "    binsz=0.02,\n",
    "    width=(4, 4),\n",
    "    frame=\"icrs\",\n",
    "    proj=\"CAR\",\n",
    "    axes=[energy_axis],\n",
    ")\n",
    "\n",
    "# Reduced IRFs are defined in true energy (i.e. not measured energy).\n",
    "energy_axis_true = MapAxis.from_energy_bounds(\n",
    "    0.1, 20, 20, unit=\"TeV\", name=\"energy_true\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can define the target dataset with this geometry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stacked = MapDataset.create(\n",
    "    geom=geom, energy_axis_true=energy_axis_true, name=\"msh-stacked\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data reduction\n",
    "\n",
    "### Create the maker classes to be used\n",
    "\n",
    "We first initialize the `Maker` objects that will take care of the data reduction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "offset_max = 2.0 * u.deg\n",
    "maker = MapDatasetMaker()\n",
    "maker_safe_mask = SafeMaskMaker(\n",
    "    methods=[\"offset-max\", \"aeff-max\", \"bkg-peak\"], offset_max=offset_max\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "circle = CircleSkyRegion(\n",
    "    center=SkyCoord(pos.ra, pos.dec), radius=0.4 * u.deg\n",
    ")\n",
    "exclusion_mask = geom.region_mask(regions=[circle], inside=False)\n",
    "maker_fov = FoVBackgroundMaker(method=\"scale\", exclusion_mask=exclusion_mask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "exclusion_mask.plot_interactive()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Perform the data reduction loop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "for obs in observations:\n",
    "    # First a cutout of the target map is produced\n",
    "    cutout = stacked.cutout(\n",
    "        obs.pointing_radec, width=2 * offset_max, name=f\"obs-{obs.obs_id}\"\n",
    "    )\n",
    "    # A MapDataset is filled in this cutout geometry\n",
    "    dataset = maker.run(cutout, obs)\n",
    "    # The data quality cut is applied\n",
    "    dataset = maker_safe_mask.run(dataset, obs)\n",
    "    # fit background model\n",
    "    dataset = maker_fov.run(dataset)\n",
    "    print(\n",
    "        f\"Background norm obs {obs.obs_id}: {dataset.background_model.spectral_model.norm.value:.2f} +/- {dataset.background_model.spectral_model.norm.error:.3f}\"\n",
    "    )\n",
    "    # The resulting dataset cutout is stacked onto the final one\n",
    "    stacked.stack(dataset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(stacked)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inspect the reduced dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stacked.counts.smooth(0.04 * u.deg).plot_interactive(\n",
    "    stretch=\"linear\", add_cbar=True\n",
    ");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset.mask_safe.plot_grid(vmin=0, vmax=1);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save the dataset to disc using `~gammapy.datasets.MapDataset.write()` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = \"msh-stacked-dataset.fits.gz\"\n",
    "stacked.write(filename, overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "region = CircleSkyRegion(\n",
    "    center=SkyCoord(228.5, -59.15, unit=\"deg\"), radius=0.4 * u.deg\n",
    ")\n",
    "stacked.plot_residuals(kwargs_spectral={\"region\":region});"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this stage, we have created a DL4 `MapDataset`. The rest of the notebook is agnostic to how this dataset was created (simulated/H.E.S.S./Fermi-LAT, etc)...\n",
    "\n",
    "If you are interested only in the spectrum, not the morphology, you can do a 1D spectral extraction using the `SpectrumDatasetMaker`. See https://docs.gammapy.org/1.0.1/tutorials/index.html#d-spectral for details."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute an excess and a significance map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "estimator = ExcessMapEstimator(0.04 * u.deg, selection_optional=[])\n",
    "lima_maps = estimator.run(stacked)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "significance_map = lima_maps[\"sqrt_ts\"]\n",
    "excess_map = lima_maps[\"npred_excess\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can plot the excess and significance maps\n",
    "plt.figure(figsize=(10, 10))\n",
    "ax1 = plt.subplot(221, projection=significance_map.geom.wcs)\n",
    "ax2 = plt.subplot(222, projection=excess_map.geom.wcs)\n",
    "\n",
    "ax1.set_title(\"Significance map\")\n",
    "significance_map.plot(ax=ax1, add_cbar=True);\n",
    "\n",
    "ax2.set_title(\"Excess map\")\n",
    "excess_map.plot(ax=ax2, add_cbar=True);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## You can zoom into your region\n",
    "significance_map.cutout(position=SkyCoord(228.5, -59.15, unit=\"deg\"), width=0.8*u.deg).plot(add_cbar=True);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define the model\n",
    "We first define the model, a `SkyModel`, as the combination of a point source `SpatialModel` with a powerlaw `SpectralModel`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spatial_model = GaussianSpatialModel(\n",
    "    lon_0=pos.ra, lat_0=pos.dec, sigma=0.3*u.deg, frame=\"icrs\"\n",
    ")\n",
    "# We limit the model position inside a 1 deg box centered on the reference position of MSH 15-52\n",
    "spatial_model.lon_0.min = spatial_model.lon_0.value - 0.5\n",
    "spatial_model.lon_0.max = spatial_model.lon_0.value + 0.5\n",
    "spatial_model.lat_0.min = spatial_model.lat_0.value - 0.5\n",
    "spatial_model.lat_0.max = spatial_model.lat_0.value + 0.5\n",
    "\n",
    "spectral_model = PowerLawSpectralModel(\n",
    "    index=2.702,\n",
    "    amplitude=4.712e-11 * u.Unit(\"1 / (cm2 s TeV)\"),\n",
    "    reference=1 * u.TeV,\n",
    ")\n",
    "\n",
    "sky_model = SkyModel(\n",
    "    spatial_model=spatial_model, spectral_model=spectral_model, name=\"msh\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we define a global `~gammapy.modeling.model.FoVBackgroundModel` in order to finely adjust the level of residual CR backgroud. This should **not** be forgotten."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bkg_model = FoVBackgroundModel(dataset_name=\"msh-stacked\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we assign these models to our reduced dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stacked.models = Models([sky_model, bkg_model])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(stacked)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fit the model\n",
    "\n",
    "The `~gammapy.modeling.Fit` class is orchestrating the fit, connecting the `stats` method of the dataset to the minimizer. By default, it uses `iminuit`.\n",
    "\n",
    "Its constructor takes a list of dataset as argument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "fit = Fit(optimize_opts={\"print_level\": 1})\n",
    "result = fit.run([stacked])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `FitResult` contains information about the optimization and parameter error calculation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result.success"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The fitted parameters are visible from the `~astropy.modeling.models.Models` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stacked.models.to_parameters_table()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inspecting residuals\n",
    "\n",
    "For any fit it is useful to inspect the residual images. We have a few options on the dataset object to handle this. First we can use `.plot_residuals_spatial()` to plot a residual image, summed over all energies:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbsphinx-thumbnail": {
     "tooltip": "Introduction to Gammapy analysis using the low level API."
    }
   },
   "outputs": [],
   "source": [
    "region = CircleSkyRegion(\n",
    "    center=pos, radius=0.4 * u.deg\n",
    ")\n",
    "stacked.plot_residuals(kwargs_spectral={\"region\":region});"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot the fitted spectrum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Making a butterfly plot \n",
    "\n",
    "The `SpectralModel` component can be used to produce a, so-called, butterfly plot showing the envelope of the model taking into account parameter uncertainties:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spec = sky_model.spectral_model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can actually do the plot using the `plot_error` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "energy_bounds = [1, 10] * u.TeV\n",
    "spec.plot(energy_bounds=energy_bounds, sed_type=\"e2dnde\")\n",
    "ax = spec.plot_error(energy_bounds=energy_bounds, sed_type=\"e2dnde\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Computing flux points\n",
    "\n",
    "We can now compute some flux points using the `~gammapy.estimators.FluxPointsEstimator`. \n",
    "\n",
    "Besides the list of datasets to use, we must provide it the energy intervals on which to compute flux points as well as the model component name. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "energy_edges = [0.5, 1, 2, 4, 10] * u.TeV\n",
    "fpe = FluxPointsEstimator(energy_edges=energy_edges, source=\"msh\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "flux_points = fpe.run(datasets=[stacked])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = spec.plot_error(energy_bounds=[0.3, 20] * u.TeV, sed_type=\"e2dnde\")\n",
    "flux_points.plot(ax=ax, sed_type=\"e2dnde\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dataset simulations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The standard way of running simulations consists of sampling data event-by-event and then bin them into a `Dataset` object. This is demonstrated in the tutorial https://docs.gammapy.org/1.0.1/tutorials/analysis-3d/event_sampling.html\n",
    "\n",
    "But is it possible to directly simulate binned 3D datasets? The answer is yes. Based on the 3D dataset produced previous sections, we can generate fake datasets based on a Poisson fluctuation of the model-predicted counts and the dataset IRFs.  \n",
    "\n",
    "The model-predicted data can be accessed through the `npred` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "npred = stacked.npred()\n",
    "\n",
    "npred.plot_interactive(\n",
    "    stretch=\"sqrt\", add_cbar=True\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now create a copy of the reduced dataset, in order not to erase the real data from memory, and then call the `Dataset.fake` method to obtain a mock data realization based on the dataset `npred`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a dataset copy and assign the models\n",
    "stacked_fake = stacked.copy(name=\"msh-stacked-fake\")\n",
    "bkg_model_fake = FoVBackgroundModel(dataset_name=\"msh-stacked\")\n",
    "stacked_fake.models = Models([sky_model, bkg_model_fake])\n",
    "\n",
    "# Perform the binned simulation\n",
    "stacked_fake.fake(random_state='random-seed')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The fake data will now be stored in the `stacked_fake.counts`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stacked_fake.counts.smooth(0.04 * u.deg).plot_interactive(\n",
    "    stretch=\"linear\", add_cbar=True\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercises\n",
    "\n",
    "## Beginner \n",
    "- Since MSH 1552 is a Galactic source, it can make sense to analyze it using Galactic coordinates instead of ICRS. Try to repeat the data reduction in the Galactic coordinate frame\n",
    "- Modify the data reduction loop to fit also the background model `tilt` parameter, in addition to the default `norm` fit\n",
    "- By default the model fit is performed in the full energy range defined by the dataset `mask_safe`. Try to restrict it by using the `Dataset.mask_fit` property\n",
    "- Repeat the flux points estimation by reoptimizing, in each energy bin, all the model free parameters. This will take longer, but it will provide a more reliable estimation of the source flux in each energy bin\n",
    "- What is the TS-based significance of MSH 1552? You can estimate it by comparing the model likelihood without (null hypothesis) and with the source model \n",
    "\n",
    "## Advanced \n",
    "- Verify that the fake data obtained in the simulation are obtained from a Poisson distribution with the `Dataset.npred` as expected value. (Hint: If this is the case, the significance of the difference between the fake data and the `Dataset.npred` has to be distributed as a standard normal Gaussian)\n",
    "- Is the source extended? You can try to understand this (without any model fitting) by producing a radial flux profiles using annulus regions concentric around the source. (Tutorial reference: https://docs.gammapy.org/1.0.1/tutorials/analysis-3d/flux_profiles.html)\n",
    "- Try to estimate the level of analysis systematics related to our imperfect knowledge of the IRFs. For example, what is the impact of a spectral bias in the background model with respect to the data? (Hint: A simulation-based approach can help you tackle this problem!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.16"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {
     "095494c994ba43aa96803803704a7cdb": {
      "model_module": "@jupyter-widgets/base",
      "model_module_version": "1.2.0",
      "model_name": "LayoutModel",
      "state": {
       "_model_module": "@jupyter-widgets/base",
       "_model_module_version": "1.2.0",
       "_model_name": "LayoutModel",
       "_view_count": null,
       "_view_module": "@jupyter-widgets/base",
       "_view_module_version": "1.2.0",
       "_view_name": "LayoutView",
       "align_content": null,
       "align_items": null,
       "align_self": null,
       "border": null,
       "bottom": null,
       "display": null,
       "flex": null,
       "flex_flow": null,
       "grid_area": null,
       "grid_auto_columns": null,
       "grid_auto_flow": null,
       "grid_auto_rows": null,
       "grid_column": null,
       "grid_gap": null,
       "grid_row": null,
       "grid_template_areas": null,
       "grid_template_columns": null,
       "grid_template_rows": null,
       "height": null,
       "justify_content": null,
       "justify_items": null,
       "left": null,
       "margin": null,
       "max_height": null,
       "max_width": null,
       "min_height": null,
       "min_width": null,
       "object_fit": null,
       "object_position": null,
       "order": null,
       "overflow": null,
       "overflow_x": null,
       "overflow_y": null,
       "padding": null,
       "right": null,
       "top": null,
       "visibility": null,
       "width": "50%"
      }
     },
     "203fa8f78fd848f596a7e354c693d2f0": {
      "model_module": "@jupyter-widgets/controls",
      "model_module_version": "1.5.0",
      "model_name": "DescriptionStyleModel",
      "state": {
       "_model_module": "@jupyter-widgets/controls",
       "_model_module_version": "1.5.0",
       "_model_name": "DescriptionStyleModel",
       "_view_count": null,
       "_view_module": "@jupyter-widgets/base",
       "_view_module_version": "1.2.0",
       "_view_name": "StyleView",
       "description_width": "initial"
      }
     },
     "37a5edfadf444b67a88a631e9cd5ed3a": {
      "model_module": "@jupyter-widgets/controls",
      "model_module_version": "1.5.0",
      "model_name": "SelectionSliderModel",
      "state": {
       "_dom_classes": [],
       "_model_module": "@jupyter-widgets/controls",
       "_model_module_version": "1.5.0",
       "_model_name": "SelectionSliderModel",
       "_options_labels": [
        "1.00e+00 TeV - 1.78e+00 TeV TeV",
        "1.78e+00 TeV - 3.16e+00 TeV TeV",
        "3.16e+00 TeV - 5.62e+00 TeV TeV",
        "5.62e+00 TeV - 1.00e+01 TeV TeV"
       ],
       "_view_count": null,
       "_view_module": "@jupyter-widgets/controls",
       "_view_module_version": "1.5.0",
       "_view_name": "SelectionSliderView",
       "continuous_update": false,
       "description": "Select energy:",
       "description_tooltip": null,
       "disabled": false,
       "index": 0,
       "layout": "IPY_MODEL_095494c994ba43aa96803803704a7cdb",
       "orientation": "horizontal",
       "readout": true,
       "style": "IPY_MODEL_5bbb5dc220274665842a1912afd311b1"
      }
     },
     "5b00e5e638ae422bbb878743039a9bd2": {
      "model_module": "@jupyter-widgets/output",
      "model_module_version": "1.0.0",
      "model_name": "OutputModel",
      "state": {
       "_dom_classes": [],
       "_model_module": "@jupyter-widgets/output",
       "_model_module_version": "1.0.0",
       "_model_name": "OutputModel",
       "_view_count": null,
       "_view_module": "@jupyter-widgets/output",
       "_view_module_version": "1.0.0",
       "_view_name": "OutputView",
       "layout": "IPY_MODEL_dcd41c15345c4232a1e35c43630c62e5",
       "msg_id": "",
       "outputs": []
      }
     },
     "5bbb5dc220274665842a1912afd311b1": {
      "model_module": "@jupyter-widgets/controls",
      "model_module_version": "1.5.0",
      "model_name": "DescriptionStyleModel",
      "state": {
       "_model_module": "@jupyter-widgets/controls",
       "_model_module_version": "1.5.0",
       "_model_name": "DescriptionStyleModel",
       "_view_count": null,
       "_view_module": "@jupyter-widgets/base",
       "_view_module_version": "1.2.0",
       "_view_name": "StyleView",
       "description_width": "initial"
      }
     },
     "699e1b20c76e4add9724d8f32f9925f9": {
      "model_module": "@jupyter-widgets/base",
      "model_module_version": "1.2.0",
      "model_name": "LayoutModel",
      "state": {
       "_model_module": "@jupyter-widgets/base",
       "_model_module_version": "1.2.0",
       "_model_name": "LayoutModel",
       "_view_count": null,
       "_view_module": "@jupyter-widgets/base",
       "_view_module_version": "1.2.0",
       "_view_name": "LayoutView",
       "align_content": null,
       "align_items": null,
       "align_self": null,
       "border": null,
       "bottom": null,
       "display": null,
       "flex": null,
       "flex_flow": null,
       "grid_area": null,
       "grid_auto_columns": null,
       "grid_auto_flow": null,
       "grid_auto_rows": null,
       "grid_column": null,
       "grid_gap": null,
       "grid_row": null,
       "grid_template_areas": null,
       "grid_template_columns": null,
       "grid_template_rows": null,
       "height": null,
       "justify_content": null,
       "justify_items": null,
       "left": null,
       "margin": null,
       "max_height": null,
       "max_width": null,
       "min_height": null,
       "min_width": null,
       "object_fit": null,
       "object_position": null,
       "order": null,
       "overflow": null,
       "overflow_x": null,
       "overflow_y": null,
       "padding": null,
       "right": null,
       "top": null,
       "visibility": null,
       "width": null
      }
     },
     "74c471760f1947f8b2e33e1c58637e33": {
      "model_module": "@jupyter-widgets/controls",
      "model_module_version": "1.5.0",
      "model_name": "VBoxModel",
      "state": {
       "_dom_classes": [
        "widget-interact"
       ],
       "_model_module": "@jupyter-widgets/controls",
       "_model_module_version": "1.5.0",
       "_model_name": "VBoxModel",
       "_view_count": null,
       "_view_module": "@jupyter-widgets/controls",
       "_view_module_version": "1.5.0",
       "_view_name": "VBoxView",
       "box_style": "",
       "children": [
        "IPY_MODEL_37a5edfadf444b67a88a631e9cd5ed3a",
        "IPY_MODEL_a4ea50bb615e41b3a1db656c1edf3a34",
        "IPY_MODEL_5b00e5e638ae422bbb878743039a9bd2"
       ],
       "layout": "IPY_MODEL_699e1b20c76e4add9724d8f32f9925f9"
      }
     },
     "a4ea50bb615e41b3a1db656c1edf3a34": {
      "model_module": "@jupyter-widgets/controls",
      "model_module_version": "1.5.0",
      "model_name": "RadioButtonsModel",
      "state": {
       "_dom_classes": [],
       "_model_module": "@jupyter-widgets/controls",
       "_model_module_version": "1.5.0",
       "_model_name": "RadioButtonsModel",
       "_options_labels": [
        "linear",
        "sqrt",
        "log"
       ],
       "_view_count": null,
       "_view_module": "@jupyter-widgets/controls",
       "_view_module_version": "1.5.0",
       "_view_name": "RadioButtonsView",
       "description": "Select stretch:",
       "description_tooltip": null,
       "disabled": false,
       "index": 0,
       "layout": "IPY_MODEL_cceb2b5d8e8241978dca97767054211d",
       "style": "IPY_MODEL_203fa8f78fd848f596a7e354c693d2f0"
      }
     },
     "cceb2b5d8e8241978dca97767054211d": {
      "model_module": "@jupyter-widgets/base",
      "model_module_version": "1.2.0",
      "model_name": "LayoutModel",
      "state": {
       "_model_module": "@jupyter-widgets/base",
       "_model_module_version": "1.2.0",
       "_model_name": "LayoutModel",
       "_view_count": null,
       "_view_module": "@jupyter-widgets/base",
       "_view_module_version": "1.2.0",
       "_view_name": "LayoutView",
       "align_content": null,
       "align_items": null,
       "align_self": null,
       "border": null,
       "bottom": null,
       "display": null,
       "flex": null,
       "flex_flow": null,
       "grid_area": null,
       "grid_auto_columns": null,
       "grid_auto_flow": null,
       "grid_auto_rows": null,
       "grid_column": null,
       "grid_gap": null,
       "grid_row": null,
       "grid_template_areas": null,
       "grid_template_columns": null,
       "grid_template_rows": null,
       "height": null,
       "justify_content": null,
       "justify_items": null,
       "left": null,
       "margin": null,
       "max_height": null,
       "max_width": null,
       "min_height": null,
       "min_width": null,
       "object_fit": null,
       "object_position": null,
       "order": null,
       "overflow": null,
       "overflow_x": null,
       "overflow_y": null,
       "padding": null,
       "right": null,
       "top": null,
       "visibility": null,
       "width": null
      }
     },
     "dcd41c15345c4232a1e35c43630c62e5": {
      "model_module": "@jupyter-widgets/base",
      "model_module_version": "1.2.0",
      "model_name": "LayoutModel",
      "state": {
       "_model_module": "@jupyter-widgets/base",
       "_model_module_version": "1.2.0",
       "_model_name": "LayoutModel",
       "_view_count": null,
       "_view_module": "@jupyter-widgets/base",
       "_view_module_version": "1.2.0",
       "_view_name": "LayoutView",
       "align_content": null,
       "align_items": null,
       "align_self": null,
       "border": null,
       "bottom": null,
       "display": null,
       "flex": null,
       "flex_flow": null,
       "grid_area": null,
       "grid_auto_columns": null,
       "grid_auto_flow": null,
       "grid_auto_rows": null,
       "grid_column": null,
       "grid_gap": null,
       "grid_row": null,
       "grid_template_areas": null,
       "grid_template_columns": null,
       "grid_template_rows": null,
       "height": null,
       "justify_content": null,
       "justify_items": null,
       "left": null,
       "margin": null,
       "max_height": null,
       "max_width": null,
       "min_height": null,
       "min_width": null,
       "object_fit": null,
       "object_position": null,
       "order": null,
       "overflow": null,
       "overflow_x": null,
       "overflow_y": null,
       "padding": null,
       "right": null,
       "top": null,
       "visibility": null,
       "width": null
      }
     }
    },
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}