.. 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 : q1OcCrGu Instrument : HESS Models : Dataset 1: Type : SpectrumDatasetOnOff Name : hSXzXYvI Instrument : HESS Models : Dataset 2: Type : SpectrumDatasetOnOff Name : LKDDjYgz Instrument : HESS Models : Dataset 3: Type : SpectrumDatasetOnOff Name : 9hxY3Db3 Instrument : HESS Models : Dataset 4: Type : SpectrumDatasetOnOff Name : UjP1fiPA Instrument : HESS Models : Dataset 5: Type : SpectrumDatasetOnOff Name : lQSRA7Ot Instrument : HESS Models : Dataset 6: Type : SpectrumDatasetOnOff Name : c9PX_C_N Instrument : HESS Models : Dataset 7: Type : SpectrumDatasetOnOff Name : ybAwpQ_a Instrument : HESS Models : Dataset 8: Type : SpectrumDatasetOnOff Name : s5v8Z8vP Instrument : HESS Models : Dataset 9: Type : SpectrumDatasetOnOff Name : CSfLe__Z Instrument : HESS Models : Dataset 10: Type : SpectrumDatasetOnOff Name : IoZf6jQx Instrument : HESS Models : Dataset 11: Type : SpectrumDatasetOnOff Name : MqPjLs34 Instrument : HESS Models : Dataset 12: Type : SpectrumDatasetOnOff Name : RsQz0Gt0 Instrument : HESS Models : Dataset 13: Type : SpectrumDatasetOnOff Name : 6eTlEIyl Instrument : HESS Models : Dataset 14: Type : SpectrumDatasetOnOff Name : HaFfRWZ5 Instrument : HESS Models : Dataset 15: Type : SpectrumDatasetOnOff Name : S8zlaiAF Instrument : HESS Models : Dataset 16: Type : SpectrumDatasetOnOff Name : TyE7NUb1 Instrument : HESS Models : Dataset 17: Type : SpectrumDatasetOnOff Name : 8gXyFVkI Instrument : HESS Models : Dataset 18: Type : SpectrumDatasetOnOff Name : FR-KmXO3 Instrument : HESS Models : Dataset 19: Type : SpectrumDatasetOnOff Name : Mg8kIb1p Instrument : HESS Models : Dataset 20: Type : SpectrumDatasetOnOff Name : VrByGuzD Instrument : HESS Models : Dataset 21: Type : SpectrumDatasetOnOff Name : CNlS-73t Instrument : HESS Models : Dataset 22: Type : SpectrumDatasetOnOff Name : cgV8UJfe Instrument : HESS Models : Dataset 23: Type : SpectrumDatasetOnOff Name : AnlVnt0Y Instrument : HESS Models : Dataset 24: Type : SpectrumDatasetOnOff Name : 51rmRTu2 Instrument : HESS Models : Dataset 25: Type : SpectrumDatasetOnOff Name : NbrazCLa Instrument : HESS Models : Dataset 26: Type : SpectrumDatasetOnOff Name : HU7jTV7z Instrument : HESS Models : Dataset 27: Type : SpectrumDatasetOnOff Name : rVhLImMw Instrument : HESS Models : Dataset 28: Type : SpectrumDatasetOnOff Name : 6PHh9jJB Instrument : HESS Models : Dataset 29: Type : SpectrumDatasetOnOff Name : kZlkUBHF Instrument : HESS Models : Dataset 30: Type : SpectrumDatasetOnOff Name : kIneTwo_ Instrument : HESS Models : Dataset 31: Type : SpectrumDatasetOnOff Name : V3BK9fjH Instrument : HESS Models : Dataset 32: Type : SpectrumDatasetOnOff Name : fmnAzqko Instrument : HESS Models : Dataset 33: Type : SpectrumDatasetOnOff Name : 1PSVIL6E 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 /Users/mregeard/Workspace/dev/code/gammapy/gammapy/.tox/build_docs/lib/python3.11/site-packages/numpy/core/fromnumeric.py:88: RuntimeWarning: overflow encountered in reduce return ufunc.reduce(obj, axis, dtype, out, **passkwargs) /Users/mregeard/Workspace/dev/code/gammapy/gammapy/.tox/build_docs/lib/python3.11/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 visualize 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 convergence index index_err amplitude amplitude_err 1 / (cm2 s TeV) 1 / (cm2 s TeV) ----------------------- ----------------------- ----------- ------------------ ------------------- ---------------------- ---------------------- 2006-07-29T20:30:00.000 2006-07-29T20:45:00.000 True 4.0086734490094305 0.3531753451180378 1.0215442304174445e-10 1.2923272671327138e-11 2006-07-29T20:45:00.000 2006-07-29T21:00:00.000 True 4.124946197229212 0.22570516613106914 1.2506136938230208e-10 1.1393881761299661e-11 2006-07-29T21:00:00.000 2006-07-29T21:15:00.000 True 3.588361705582432 0.13050003825315337 1.6485973133488616e-10 9.820734773683985e-12 2006-07-29T21:15:00.000 2006-07-29T21:30:00.000 True 3.4272511141842323 0.10213517590250513 1.640138036161798e-10 1.0033685315041575e-11 2006-07-29T21:30:00.000 2006-07-29T21:45:00.000 True 3.488160391427099 0.07290480076152729 2.1619396351099933e-10 1.0859121804265987e-11 2006-07-29T21:45:00.000 2006-07-29T22:00:00.000 True 3.6847561893518432 0.0892728529988859 1.8431926467132017e-10 1.1661916125989783e-11 2006-07-29T22:00:00.000 2006-07-29T22:15:00.000 True 3.5496018522628856 0.08729118437678622 1.587628603617094e-10 9.683531578205622e-12 2006-07-29T22:15:00.000 2006-07-29T22:30:00.000 True 3.685149963186663 0.11537399447253352 1.1147742234725352e-10 9.200164221906703e-12 2006-07-29T22:30:00.000 2006-07-29T22:45:00.000 True 3.622909292257767 0.10841490276258127 1.06498208240016e-10 8.243746306437121e-12 2006-07-29T22:45:00.000 2006-07-29T23:00:00.000 True 3.5642236448958657 0.11889238916653445 1.0362682001464577e-10 8.685882328799652e-12 2006-07-29T23:00:00.000 2006-07-29T23:15:00.000 True 3.417860642833073 0.10445367683973007 1.1266135555905379e-10 8.065714177096826e-12 2006-07-29T23:15:00.000 2006-07-29T23:30:00.000 True 3.80347487867966 0.1845605014026926 5.0351196110366946e-11 6.7547188456794526e-12 2006-07-29T23:30:00.000 2006-07-29T23:45:00.000 True 4.119403076791005 0.20578237105496278 3.516178346389486e-11 5.539061634078695e-12 2006-07-29T23:45:00.000 2006-07-30T00:00:00.000 True 3.963169999400986 0.22661395406431778 4.260261956530079e-11 7.202169952790548e-12 2006-07-30T00:00:00.000 2006-07-30T00:15:00.000 True 3.3827191178859177 0.15244913050690279 8.22378535817987e-11 8.538492113839458e-12 2006-07-30T00:15:00.000 2006-07-30T00:30:00.000 True 3.690344344737323 0.18527001217906722 4.260615821949757e-11 5.5971938698284516e-12 2006-07-30T00:30:00.000 2006-07-30T00:45:00.000 True 3.6869093684972123 0.20376940713235722 4.261614584186258e-11 6.166871726925046e-12 2006-07-30T00:45:00.000 2006-07-30T01:00:00.000 True 4.050046162135841 0.19235626954902682 4.042332754171405e-11 5.871732716149608e-12 2006-07-30T01:00:00.000 2006-07-30T01:15:00.000 True 3.620709504356781 0.15963275377206024 6.19965922666747e-11 6.972194169499312e-12 2006-07-30T01:15:00.000 2006-07-30T01:30:00.000 True 3.817796201656764 0.16226332715359063 5.3910755456212985e-11 6.3738978978922026e-12 2006-07-30T01:30:00.000 2006-07-30T01:45:00.000 True 3.7343463377488475 0.17888971333642376 4.971179459757526e-11 6.397035499771326e-12 2006-07-30T01:45:00.000 2006-07-30T02:00:00.000 True 3.907921576839459 0.1970169979311921 3.6203052100234095e-11 5.282953651237949e-12 2006-07-30T02:00:00.000 2006-07-30T02:15:00.000 True 3.564956893264179 0.182281399862031 4.4538359076922994e-11 5.711541740165e-12 2006-07-30T02:15:00.000 2006-07-30T02:30:00.000 True 3.734296015242402 0.2070829554449924 3.064399383964599e-11 4.579177137463697e-12 .. 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.59497706e-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 11.413 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/main?urlpath=lab/tree/notebooks/dev/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 `_