.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "tutorials/analysis-time/time_resolved_spectroscopy.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. or to run this example in your browser via Binder .. rst-class:: sphx-glr-example-title .. _sphx_glr_tutorials_analysis-time_time_resolved_spectroscopy.py: Time resolved spectroscopy estimator ==================================== Perform spectral fits of a blazar in different time bins to investigate spectral changes during flares. Context ------- The `~gammapy.estimators.LightCurveEstimator` in Gammapy (see :doc:`light curve notebook `, and :doc:`light curve for flares notebook `.) fits the amplitude in each time/energy bin, keeping the spectral shape frozen. However, in the analysis of flaring sources, it is often interesting to study not only how the flux changes with time but how the spectral shape varies with time. Proposed approach ----------------- The main idea behind doing a time resolved spectroscopy is to - Select relevant `~gammapy.data.Observations` from the `~gammapy.data.DataStore` - Define time intervals in which to fit the spectral model - Apply the above time selections on the data to obtain new `~gammapy.data.Observations` - Perform standard data reduction on the above data - Define a source model - Fit the reduced data in each time bin with the source model - Extract relevant information in a table Here, we will use the PKS 2155-304 observations from the `H.E.S.S. first public test data release `__. We use time intervals of 15 minutes duration to explore spectral variability. Setup ----- As usual, we’ll start with some general imports… .. GENERATED FROM PYTHON SOURCE LINES 47-78 .. code-block:: Python import logging import numpy as np import astropy.units as u from astropy.coordinates import Angle, SkyCoord from astropy.table import QTable from astropy.time import Time from regions import CircleSkyRegion # %matplotlib inline import matplotlib.pyplot as plt log = logging.getLogger(__name__) from gammapy.data import DataStore from gammapy.datasets import Datasets, SpectrumDataset from gammapy.makers import ( ReflectedRegionsBackgroundMaker, SafeMaskMaker, SpectrumDatasetMaker, ) from gammapy.maps import MapAxis, RegionGeom, TimeMapAxis from gammapy.modeling import Fit from gammapy.modeling.models import ( PowerLawSpectralModel, SkyModel, ) log = logging.getLogger(__name__) .. GENERATED FROM PYTHON SOURCE LINES 79-84 Data selection ~~~~~~~~~~~~~~ We select all runs pointing within 2 degrees of PKS 2155-304. .. GENERATED FROM PYTHON SOURCE LINES 84-99 .. code-block:: Python data_store = DataStore.from_dir("$GAMMAPY_DATA/hess-dl3-dr1/") target_position = SkyCoord(329.71693826 * u.deg, -30.2255890 * u.deg, frame="icrs") selection = dict( type="sky_circle", frame="icrs", lon=target_position.ra, lat=target_position.dec, radius=2 * u.deg, ) obs_ids = data_store.obs_table.select_observations(selection)["OBS_ID"] observations = data_store.get_observations(obs_ids) print(f"Number of selected observations : {len(observations)}") .. rst-class:: sphx-glr-script-out .. code-block:: none Number of selected observations : 21 .. GENERATED FROM PYTHON SOURCE LINES 100-105 The flaring observations were taken during July 2006. We define 15-minute time intervals as lists of `~astropy.time.Time` start and stop objects, and apply the intervals to the observations by using `~gammapy.data.Observations.select_time` .. GENERATED FROM PYTHON SOURCE LINES 105-120 .. code-block:: Python t0 = Time("2006-07-29T20:30") duration = 15 * u.min n_time_bins = 25 times = t0 + np.arange(n_time_bins) * duration time_intervals = [Time([tstart, tstop]) for tstart, tstop in zip(times[:-1], times[1:])] print(time_intervals[-1].mjd) short_observations = observations.select_time(time_intervals) # check that observations have been filtered print(f"Number of observations after time filtering: {len(short_observations)}\n") print(short_observations[1].gti) .. rst-class:: sphx-glr-script-out .. code-block:: none [53946.09375 53946.10416667] Number of observations after time filtering: 34 GTI info: - Number of GTIs: 1 - Duration: 461.99999999999545 s - Start: 207521165.184 s MET - Start: 2006-07-29T20:45:00.000 (time standard: UTC) - Stop: 207521627.184 s MET - Stop: 2006-07-29T20:53:47.184 (time standard: TT) .. GENERATED FROM PYTHON SOURCE LINES 121-128 Data reduction -------------- In this example, we perform a 1D analysis with a reflected regions background estimation. For details, see the :doc:`/tutorials/analysis-1d/spectral_analysis` tutorial. .. GENERATED FROM PYTHON SOURCE LINES 128-157 .. code-block:: Python energy_axis = MapAxis.from_energy_bounds("0.4 TeV", "20 TeV", nbin=10) energy_axis_true = MapAxis.from_energy_bounds( "0.1 TeV", "40 TeV", nbin=20, name="energy_true" ) on_region_radius = Angle("0.11 deg") on_region = CircleSkyRegion(center=target_position, radius=on_region_radius) geom = RegionGeom.create(region=on_region, axes=[energy_axis]) dataset_maker = SpectrumDatasetMaker( containment_correction=True, selection=["counts", "exposure", "edisp"] ) bkg_maker = ReflectedRegionsBackgroundMaker() safe_mask_masker = SafeMaskMaker(methods=["aeff-max"], aeff_percent=10) datasets = Datasets() dataset_empty = SpectrumDataset.create(geom=geom, energy_axis_true=energy_axis_true) for obs in short_observations: dataset = dataset_maker.run(dataset_empty.copy(), obs) dataset_on_off = bkg_maker.run(dataset, obs) dataset_on_off = safe_mask_masker.run(dataset_on_off, obs) datasets.append(dataset_on_off) .. GENERATED FROM PYTHON SOURCE LINES 158-161 This gives us list of `~gammapy.datasets.SpectrumDatasetOnOff` which can now be modelled. .. GENERATED FROM PYTHON SOURCE LINES 161-165 .. code-block:: Python print(datasets) .. rst-class:: sphx-glr-script-out .. code-block:: none Datasets -------- Dataset 0: Type : SpectrumDatasetOnOff Name : um4eC0bQ Instrument : HESS Models : Dataset 1: Type : SpectrumDatasetOnOff Name : y9tCHvwm Instrument : HESS Models : Dataset 2: Type : SpectrumDatasetOnOff Name : CsNKgP8z Instrument : HESS Models : Dataset 3: Type : SpectrumDatasetOnOff Name : -ujYvuBF Instrument : HESS Models : Dataset 4: Type : SpectrumDatasetOnOff Name : XYdsfU27 Instrument : HESS Models : Dataset 5: Type : SpectrumDatasetOnOff Name : qc2GbhL- Instrument : HESS Models : Dataset 6: Type : SpectrumDatasetOnOff Name : fM-S1BfV Instrument : HESS Models : Dataset 7: Type : SpectrumDatasetOnOff Name : SuiZXK-n Instrument : HESS Models : Dataset 8: Type : SpectrumDatasetOnOff Name : qAMc1BMv Instrument : HESS Models : Dataset 9: Type : SpectrumDatasetOnOff Name : 5R7k8RuH Instrument : HESS Models : Dataset 10: Type : SpectrumDatasetOnOff Name : WAVmP4u1 Instrument : HESS Models : Dataset 11: Type : SpectrumDatasetOnOff Name : DLpU_BKk Instrument : HESS Models : Dataset 12: Type : SpectrumDatasetOnOff Name : szdLSlTv Instrument : HESS Models : Dataset 13: Type : SpectrumDatasetOnOff Name : ZQA6-yD1 Instrument : HESS Models : Dataset 14: Type : SpectrumDatasetOnOff Name : RUezUiJe Instrument : HESS Models : Dataset 15: Type : SpectrumDatasetOnOff Name : fkesaDeB Instrument : HESS Models : Dataset 16: Type : SpectrumDatasetOnOff Name : v1A0JU9F Instrument : HESS Models : Dataset 17: Type : SpectrumDatasetOnOff Name : 5i9p2jQz Instrument : HESS Models : Dataset 18: Type : SpectrumDatasetOnOff Name : RnLxEhdH Instrument : HESS Models : Dataset 19: Type : SpectrumDatasetOnOff Name : FjYAsVgX Instrument : HESS Models : Dataset 20: Type : SpectrumDatasetOnOff Name : rT_I2uP0 Instrument : HESS Models : Dataset 21: Type : SpectrumDatasetOnOff Name : l2ngbhL4 Instrument : HESS Models : Dataset 22: Type : SpectrumDatasetOnOff Name : WMxhVdAu Instrument : HESS Models : Dataset 23: Type : SpectrumDatasetOnOff Name : n2hS2cZ2 Instrument : HESS Models : Dataset 24: Type : SpectrumDatasetOnOff Name : pY-uwzNj Instrument : HESS Models : Dataset 25: Type : SpectrumDatasetOnOff Name : TMoqyWM6 Instrument : HESS Models : Dataset 26: Type : SpectrumDatasetOnOff Name : rMsbjU0v Instrument : HESS Models : Dataset 27: Type : SpectrumDatasetOnOff Name : w9oWeK0H Instrument : HESS Models : Dataset 28: Type : SpectrumDatasetOnOff Name : 8sLj2dDU Instrument : HESS Models : Dataset 29: Type : SpectrumDatasetOnOff Name : 7mYf8T2K Instrument : HESS Models : Dataset 30: Type : SpectrumDatasetOnOff Name : 4rAc4qYJ Instrument : HESS Models : Dataset 31: Type : SpectrumDatasetOnOff Name : aYqHYBx6 Instrument : HESS Models : Dataset 32: Type : SpectrumDatasetOnOff Name : 04HMFgFz Instrument : HESS Models : Dataset 33: Type : SpectrumDatasetOnOff Name : 19D0x7GV Instrument : HESS Models : .. GENERATED FROM PYTHON SOURCE LINES 166-177 Modeling -------- We will first fit a simple power law model in each time bin. Note that since we are using an on-off analysis here, no background model is required. If you are doing a 3D FoV analysis, you will need to model the background appropriately as well. The index and amplitude of the spectral model is kept free. You can configure the quantities you want to freeze. .. GENERATED FROM PYTHON SOURCE LINES 177-188 .. code-block:: Python spectral_model = PowerLawSpectralModel( index=3.0, amplitude=2e-11 * u.Unit("1 / (cm2 s TeV)"), reference=1 * u.TeV ) spectral_model.parameters["index"].frozen = False sky_model = SkyModel(spatial_model=None, spectral_model=spectral_model, name="pks2155") print(sky_model) .. rst-class:: sphx-glr-script-out .. code-block:: none SkyModel Name : pks2155 Datasets names : None Spectral model type : PowerLawSpectralModel Spatial model type : Temporal model type : Parameters: index : 3.000 +/- 0.00 amplitude : 2.00e-11 +/- 0.0e+00 1 / (cm2 s TeV) reference (frozen): 1.000 TeV .. GENERATED FROM PYTHON SOURCE LINES 189-200 Time resolved spectroscopy algorithm ------------------------------------ The following function is the crux of this tutorial. The ``sky_model`` is fit in each bin and a list of ``fit_results`` stores the fit information in each bin. If time bins are present without any available observations, those bins are discarded and a new list of valid time intervals and fit results are created. .. GENERATED FROM PYTHON SOURCE LINES 200-226 .. code-block:: Python def time_resolved_spectroscopy(datasets, model, time_intervals): fit = Fit() valid_intervals = [] fit_results = [] index = 0 for t_min, t_max in time_intervals: datasets_to_fit = datasets.select_time(time_min=t_min, time_max=t_max) if len(datasets_to_fit) == 0: log.info( f"No Dataset for the time interval {t_min} to {t_max}. Skipping interval." ) continue model_in_bin = model.copy(name="Model_bin_" + str(index)) datasets_to_fit.models = model_in_bin result = fit.run(datasets_to_fit) fit_results.append(result) valid_intervals.append([t_min, t_max]) index += 1 return valid_intervals, fit_results .. GENERATED FROM PYTHON SOURCE LINES 227-229 We now apply it to our data .. GENERATED FROM PYTHON SOURCE LINES 229-233 .. code-block:: Python valid_times, results = time_resolved_spectroscopy(datasets, sky_model, time_intervals) .. rst-class:: sphx-glr-script-out .. code-block:: none /home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/lib/python3.9/site-packages/numpy/core/fromnumeric.py:88: RuntimeWarning: overflow encountered in reduce return ufunc.reduce(obj, axis, dtype, out, **passkwargs) /home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/lib/python3.9/site-packages/numpy/core/fromnumeric.py:88: RuntimeWarning: overflow encountered in reduce return ufunc.reduce(obj, axis, dtype, out, **passkwargs) .. GENERATED FROM PYTHON SOURCE LINES 234-236 To view the results of the fit, .. GENERATED FROM PYTHON SOURCE LINES 236-240 .. code-block:: Python print(results[0]) .. rst-class:: sphx-glr-script-out .. code-block:: none OptimizeResult backend : minuit method : migrad success : True message : Optimization terminated successfully. nfev : 76 total stat : 6.00 CovarianceResult backend : minuit method : hesse success : True message : Hesse terminated successfully. .. GENERATED FROM PYTHON SOURCE LINES 241-243 Or, to access the fitted models, .. GENERATED FROM PYTHON SOURCE LINES 243-247 .. code-block:: Python print(results[0].models) .. rst-class:: sphx-glr-script-out .. code-block:: none DatasetModels Component 0: SkyModel Name : Model_bin_0 Datasets names : None Spectral model type : PowerLawSpectralModel Spatial model type : Temporal model type : Parameters: index : 4.009 +/- 0.35 amplitude : 1.02e-10 +/- 1.3e-11 1 / (cm2 s TeV) reference (frozen): 1.000 TeV .. GENERATED FROM PYTHON SOURCE LINES 248-254 To better visualise the data, we can create a table by extracting some relevant information. In the following, we extract the time intervals, information on the fit convergence and the free parameters. You can extract more information if required, eg, the `total_stat` in each bin, etc. .. GENERATED FROM PYTHON SOURCE LINES 254-277 .. code-block:: Python def create_table(time_intervals, fit_result): t = QTable() t["tstart"] = np.array(time_intervals).T[0] t["tstop"] = np.array(time_intervals).T[1] t["convergence"] = [result.success for result in fit_result] for par in fit_result[0].models.parameters.free_parameters: t[par.name] = [ result.models.parameters[par.name].value * par.unit for result in fit_result ] t[par.name + "_err"] = [ result.models.parameters[par.name].error * par.unit for result in fit_result ] return t table = create_table(valid_times, results) print(table) .. rst-class:: sphx-glr-script-out .. code-block:: none tstart tstop ... amplitude_err ... 1 / (cm2 s TeV) ----------------------- ----------------------- ... ---------------------- 2006-07-29T20:30:00.000 2006-07-29T20:45:00.000 ... 1.2923272203680937e-11 2006-07-29T20:45:00.000 2006-07-29T21:00:00.000 ... 1.1393881668086242e-11 2006-07-29T21:00:00.000 2006-07-29T21:15:00.000 ... 9.820734654694367e-12 2006-07-29T21:15:00.000 2006-07-29T21:30:00.000 ... 1.00336861061002e-11 2006-07-29T21:30:00.000 2006-07-29T21:45:00.000 ... 1.0859117235629056e-11 2006-07-29T21:45:00.000 2006-07-29T22:00:00.000 ... 1.1661916367760008e-11 2006-07-29T22:00:00.000 2006-07-29T22:15:00.000 ... 9.683535101178161e-12 2006-07-29T22:15:00.000 2006-07-29T22:30:00.000 ... 9.20016447648031e-12 2006-07-29T22:30:00.000 2006-07-29T22:45:00.000 ... 8.243742672227972e-12 2006-07-29T22:45:00.000 2006-07-29T23:00:00.000 ... 8.68588096384626e-12 ... ... ... ... 2006-07-30T00:00:00.000 2006-07-30T00:15:00.000 ... 8.53849342472254e-12 2006-07-30T00:15:00.000 2006-07-30T00:30:00.000 ... 5.597194451442478e-12 2006-07-30T00:30:00.000 2006-07-30T00:45:00.000 ... 6.166871893629602e-12 2006-07-30T00:45:00.000 2006-07-30T01:00:00.000 ... 5.871730083414263e-12 2006-07-30T01:00:00.000 2006-07-30T01:15:00.000 ... 6.9721952295357764e-12 2006-07-30T01:15:00.000 2006-07-30T01:30:00.000 ... 6.373892759762174e-12 2006-07-30T01:30:00.000 2006-07-30T01:45:00.000 ... 6.397035933646679e-12 2006-07-30T01:45:00.000 2006-07-30T02:00:00.000 ... 5.282951078740491e-12 2006-07-30T02:00:00.000 2006-07-30T02:15:00.000 ... 5.7115414505410296e-12 2006-07-30T02:15:00.000 2006-07-30T02:30:00.000 ... 4.579178801616083e-12 Length = 24 rows .. GENERATED FROM PYTHON SOURCE LINES 278-284 Visualising the results ~~~~~~~~~~~~~~~~~~~~~~~~ We can plot the spectral index and the amplitude as a function of time. For convenience, we will convert the times into a `~gammapy.maps.TimeMapAxis`. .. GENERATED FROM PYTHON SOURCE LINES 284-306 .. code-block:: Python time_axis = TimeMapAxis.from_time_edges( time_min=table["tstart"], time_max=table["tstop"] ) fix, axes = plt.subplots(2, 1, figsize=(8, 8)) axes[0].errorbar( x=time_axis.as_plot_center, y=table["index"], yerr=table["index_err"], fmt="o" ) axes[1].errorbar( x=time_axis.as_plot_center, y=table["amplitude"], yerr=table["amplitude_err"], fmt="o", ) axes[0].set_ylabel("index") axes[1].set_ylabel("amplitude") axes[1].set_xlabel("time") plt.show() .. image-sg:: /tutorials/analysis-time/images/sphx_glr_time_resolved_spectroscopy_001.png :alt: time resolved spectroscopy :srcset: /tutorials/analysis-time/images/sphx_glr_time_resolved_spectroscopy_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 307-310 To get the integrated flux, we can access the model stored in the fit result object, eg .. GENERATED FROM PYTHON SOURCE LINES 310-319 .. code-block:: Python integral_flux = ( results[0] .models[0] .spectral_model.integral_error(energy_min=1 * u.TeV, energy_max=10 * u.TeV) ) print("Integral flux in the first bin:", integral_flux) .. rst-class:: sphx-glr-script-out .. code-block:: none Integral flux in the first bin: [3.39200283e-11 4.59497703e-12] 1 / (cm2 s) .. GENERATED FROM PYTHON SOURCE LINES 320-323 To plot hysteresis curves, ie the spectral index as a function of amplitude .. GENERATED FROM PYTHON SOURCE LINES 323-339 .. code-block:: Python plt.errorbar( table["amplitude"], table["index"], xerr=table["amplitude_err"], yerr=table["index_err"], linestyle=":", linewidth=0.5, ) plt.scatter(table["amplitude"], table["index"], c=time_axis.center.value) plt.xlabel("amplitude") plt.ylabel("index") plt.colorbar().set_label("time") plt.show() .. image-sg:: /tutorials/analysis-time/images/sphx_glr_time_resolved_spectroscopy_002.png :alt: time resolved spectroscopy :srcset: /tutorials/analysis-time/images/sphx_glr_time_resolved_spectroscopy_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 340-349 Exercises --------- 1. Quantify the variability in the spectral index 2. Rerun the algorithm using a different spectral shape, such as a broken power law. 3. Compare the significance of the new model with the simple power law. Take note of any fit non-convergence in the bins. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 22.502 seconds) .. _sphx_glr_download_tutorials_analysis-time_time_resolved_spectroscopy.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: binder-badge .. image:: images/binder_badge_logo.svg :target: https://mybinder.org/v2/gh/gammapy/gammapy-webpage/v1.3?urlpath=lab/tree/notebooks/1.3/tutorials/analysis-time/time_resolved_spectroscopy.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: time_resolved_spectroscopy.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: time_resolved_spectroscopy.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: time_resolved_spectroscopy.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_