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       : 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     :

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)

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

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

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

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.59497706e-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 11.413 seconds)

Gallery generated by Sphinx-Gallery

On this page
Show Source