Note

Go to the end to download the full example code. or to run this example in your browser via Binder

Time resolved spectroscopy estimator#

Perform spectral fits of a blazar in different time bins to investigate spectral changes during flares.

Context#

The LightCurveEstimator in Gammapy (see light curve notebook, and 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 Observations from the DataStore

  • Define time intervals in which to fit the spectral model

  • Apply the above time selections on the data to obtain new 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…

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__)

Data selection#

We select all runs pointing within 2 degrees of PKS 2155-304.

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)}")

Number of selected observations : 21

The flaring observations were taken during July 2006. We define 15-minute time intervals as lists of Time start and stop objects, and apply the intervals to the observations by using select_time

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)

[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)

Data reduction#

In this example, we perform a 1D analysis with a reflected regions background estimation. For details, see the Spectral analysis tutorial.

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)

This gives us list of SpectrumDatasetOnOff which can now be modelled.

print(datasets)

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     :

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.

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)

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

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.

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

We now apply it to our data

valid_times, results = time_resolved_spectroscopy(datasets, sky_model, time_intervals)

/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)

To view the results of the fit,

print(results[0])

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.

Or, to access the fitted models,

print(results[0].models)

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

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.

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)

         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

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 TimeMapAxis.

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()
time resolved spectroscopy

To get the integrated flux, we can access the model stored in the fit result object, eg

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)

Integral flux in the first bin: [3.39200283e-11 4.59497703e-12] 1 / (cm2 s)

To plot hysteresis curves, ie the spectral index as a function of amplitude

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()
time resolved spectroscopy

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.

Total running time of the script: (0 minutes 22.502 seconds)

Gallery generated by Sphinx-Gallery