commit docs dev version

docs/dev/_downloads/49945b472a79d59ea792f8e113366e99/Time_resolved_spectroscopy.ipynb

@@ -0,0 +1,284 @@
+{
+  "cells": [
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "\n# A time resolved spectroscopy estimator\n\n## Aim\n\nPerform spectral fits of a blazar in different time bins to investigate\nspectral changes during flares.\n\n## Context\n\nThe `~gammapy.estimators.LightCurveEstimator` in Gammapy (see\n:doc:`light curve notebook </tutorials/analysis-time/light_curve>`,\nand\n:doc:`light curve for flares notebook </tutorials/analysis-time/light_curve_flare>`.)\nfits the amplitude in each time/energy bin, keeping the spectral shape\nfrozen. However, in the analysis of flaring sources, it is often\ninteresting to study not only how the flux changes with time but how the\nspectral shape varies with time.\n\n## Proposed approach\n\nThe main idea behind doing a time resolved spectroscopy is to\n\n-  Select relevant `~gammapy.data.Observations` from the\n   `~gammapy.data.DataStore`\n-  Define time intervals in which to fit the spectral model\n-  Apply the above time selections on the data to obtain new\n   `~gammapy.data.Observations`\n-  Perform standard data reduction on the above data\n-  Define a source model\n-  Fit the reduced data in each time bin with the source model\n-  Extract relevant information in a table\n\nHere, we will use the PKS 2155-304 observations from the\n[H.E.S.S. first public test data release](https://www.mpi-hd.mpg.de/hfm/HESS/pages/dl3-dr1/)_.\n\nWe use time intervals of 15 minutes duration to explore spectral\nvariability.\n\n## Setup\n\nAs usual, we\u2019ll start with some general imports\u2026\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "import logging\nimport numpy as np\nimport astropy.units as u\nfrom astropy.coordinates import Angle, SkyCoord\nfrom astropy.table import QTable\nfrom astropy.time import Time\nfrom regions import CircleSkyRegion\n\n# %matplotlib inline\nimport matplotlib.pyplot as plt\n\nlog = logging.getLogger(__name__)\n\nfrom gammapy.data import DataStore\nfrom gammapy.datasets import Datasets, SpectrumDataset\nfrom gammapy.makers import (\n    ReflectedRegionsBackgroundMaker,\n    SafeMaskMaker,\n    SpectrumDatasetMaker,\n)\nfrom gammapy.maps import MapAxis, RegionGeom, TimeMapAxis\nfrom gammapy.modeling import Fit\nfrom gammapy.modeling.models import (\n    PowerLawSpectralModel,\n    SkyModel,\n)\n\nlog = logging.getLogger(__name__)"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "### Data selection\n\nWe select all runs pointing within 2 degrees of PKS 2155-304.\n\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "data_store = DataStore.from_dir(\"$GAMMAPY_DATA/hess-dl3-dr1/\")\ntarget_position = SkyCoord(329.71693826 * u.deg, -30.2255890 * u.deg, frame=\"icrs\")\nselection = dict(\n    type=\"sky_circle\",\n    frame=\"icrs\",\n    lon=target_position.ra,\n    lat=target_position.dec,\n    radius=2 * u.deg,\n)\nobs_ids = data_store.obs_table.select_observations(selection)[\"OBS_ID\"]\nobservations = data_store.get_observations(obs_ids)\nprint(f\"Number of selected observations : {len(observations)}\")"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "The flaring observations were taken during July 2006. We define\n15-minute time intervals as lists of `~astropy.time.Time` start and stop\nobjects, and apply the intervals to the observations by using\n`~gammapy.data.Observations.select_time`\n\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "t0 = Time(\"2006-07-29T20:30\")\nduration = 15 * u.min\nn_time_bins = 25\ntimes = t0 + np.arange(n_time_bins) * duration\n\ntime_intervals = [Time([tstart, tstop]) for tstart, tstop in zip(times[:-1], times[1:])]\nprint(time_intervals[-1].mjd)\nshort_observations = observations.select_time(time_intervals)\n\n# check that observations have been filtered\nprint(f\"Number of observations after time filtering: {len(short_observations)}\\n\")\nprint(short_observations[1].gti)"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## Data reduction\n\nIn this example, we perform a 1D analysis with a reflected regions\nbackground estimation. For details, see the\n:doc:`/tutorials/analysis-1d/spectral_analysis` tutorial.\n\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "energy_axis = MapAxis.from_energy_bounds(\"0.4 TeV\", \"20 TeV\", nbin=10)\nenergy_axis_true = MapAxis.from_energy_bounds(\n    \"0.1 TeV\", \"40 TeV\", nbin=20, name=\"energy_true\"\n)\n\non_region_radius = Angle(\"0.11 deg\")\non_region = CircleSkyRegion(center=target_position, radius=on_region_radius)\n\ngeom = RegionGeom.create(region=on_region, axes=[energy_axis])\n\ndataset_maker = SpectrumDatasetMaker(\n    containment_correction=True, selection=[\"counts\", \"exposure\", \"edisp\"]\n)\nbkg_maker = ReflectedRegionsBackgroundMaker()\nsafe_mask_masker = SafeMaskMaker(methods=[\"aeff-max\"], aeff_percent=10)\n\ndatasets = Datasets()\n\ndataset_empty = SpectrumDataset.create(geom=geom, energy_axis_true=energy_axis_true)\n\nfor obs in short_observations:\n    dataset = dataset_maker.run(dataset_empty.copy(), obs)\n\n    dataset_on_off = bkg_maker.run(dataset, obs)\n    dataset_on_off = safe_mask_masker.run(dataset_on_off, obs)\n    datasets.append(dataset_on_off)"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "This gives us list of `~gammapy.datasets.SpectrumDatasetOnOff` which can now be\nmodelled.\n\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "print(datasets)"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## Modeling\n\nWe will first fit a simple power law model in each time bin. Note that\nsince we are using an on-off analysis here, no background model is\nrequired. If you are doing a 3D FoV analysis, you will need to model the\nbackground appropriately as well.\n\nThe index and amplitude of the spectral model is kept free. You can\nconfigure the quantities you want to freeze.\n\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "spectral_model = PowerLawSpectralModel(\n    index=3.0, amplitude=2e-11 * u.Unit(\"1 / (cm2 s TeV)\"), reference=1 * u.TeV\n)\nspectral_model.parameters[\"index\"].frozen = False\n\n\nsky_model = SkyModel(spatial_model=None, spectral_model=spectral_model, name=\"pks2155\")\nprint(sky_model)"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## Time resolved spectroscopy algorithm\n\nThe following function is the crux of this tutorial. The ``sky_model``\nis fit in each bin and a list of ``fit_results`` stores the fit\ninformation in each bin.\n\nIf time bins are present without any available observations, those bins\nare discarded and a new list of valid time intervals and fit results are\ncreated.\n\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "def time_resolved_spectroscopy(datasets, model, time_intervals):\n    fit = Fit()\n    valid_intervals = []\n    fit_results = []\n    index = 0\n    for t_min, t_max in time_intervals:\n        datasets_to_fit = datasets.select_time(time_min=t_min, time_max=t_max)\n\n        if len(datasets_to_fit) == 0:\n            log.info(\n                f\"No Dataset for the time interval {t_min} to {t_max}. Skipping interval.\"\n            )\n            continue\n\n        model_in_bin = model.copy(name=\"Model_bin_\" + str(index))\n        datasets_to_fit.models = model_in_bin\n        result = fit.run(datasets_to_fit)\n        fit_results.append(result)\n        valid_intervals.append([t_min, t_max])\n        index += 1\n\n    return valid_intervals, fit_results"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "We now apply it to our data\n\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "valid_times, results = time_resolved_spectroscopy(datasets, sky_model, time_intervals)"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "To view the results of the fit,\n\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "print(results[0])"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Or, to access the fitted models,\n\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "print(results[0].models)"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "To better visualize the data, we can create a table by extracting some\nrelevant information. In the following, we extract the time intervals,\ninformation on the fit convergence and the free parameters. You can\nextract more information if required, eg, the `total_stat` in each\nbin, etc.\n\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "def create_table(time_intervals, fit_result):\n\n    t = QTable()\n\n    t[\"tstart\"] = np.array(time_intervals).T[0]\n    t[\"tstop\"] = np.array(time_intervals).T[1]\n    t[\"convergence\"] = [result.success for result in fit_result]\n    for par in fit_result[0].models.parameters.free_parameters:\n        t[par.name] = [\n            result.models.parameters[par.name].value * par.unit for result in fit_result\n        ]\n        t[par.name + \"_err\"] = [\n            result.models.parameters[par.name].error * par.unit for result in fit_result\n        ]\n\n    return t\n\n\ntable = create_table(valid_times, results)\nprint(table)"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "### Visualising the results\n\nWe can plot the spectral index and the amplitude as a function of time.\nFor convenience, we will convert the times into a `~gammapy.maps.TimeMapAxis`.\n\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "time_axis = TimeMapAxis.from_time_edges(\n    time_min=table[\"tstart\"], time_max=table[\"tstop\"]\n)\n\nfix, axes = plt.subplots(2, 1, figsize=(8, 8))\naxes[0].errorbar(\n    x=time_axis.as_plot_center, y=table[\"index\"], yerr=table[\"index_err\"], fmt=\"o\"\n)\naxes[1].errorbar(\n    x=time_axis.as_plot_center,\n    y=table[\"amplitude\"],\n    yerr=table[\"amplitude_err\"],\n    fmt=\"o\",\n)\n\naxes[0].set_ylabel(\"index\")\naxes[1].set_ylabel(\"amplitude\")\naxes[1].set_xlabel(\"time\")\nplt.show()"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "To get the integrated flux, we can access the model stored in the fit\nresult object, eg\n\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "integral_flux = (\n    results[0]\n    .models[0]\n    .spectral_model.integral_error(energy_min=1 * u.TeV, energy_max=10 * u.TeV)\n)\nprint(\"Integral flux in the first bin:\", integral_flux)"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "To plot hysteresis curves, ie the spectral index as a function of\namplitude\n\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "plt.scatter(table[\"amplitude\"], table[\"index\"], c=time_axis.center.value)\nplt.plot(table[\"amplitude\"], table[\"index\"], linewidth=0.5)\nplt.xlabel(\"amplitude\")\nplt.ylabel(\"index\")\nplt.colorbar()\nplt.show()"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## Exercises\n\n1. Quantify the variability in the spectral index\n2. Rerun the algorithm using a different spectral shape, such as a\n   broken power law.\n3. Compare the significance of the new model with the simple power law.\n   Take note of any fit non-convergence in the bins.\n\n\n"
+      ]
+    }
+  ],
+  "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.9.20"
+    }
+  },
+  "nbformat": 4,
+  "nbformat_minor": 0
+}