X-ray and Gamma-ray Spectral Fitting

Redback provides a complete pipeline for X-ray and gamma-ray spectral analysis using photon counting statistics. This includes simulation of high-energy observations, loading OGIP-standard data products, spectral fitting with physically motivated models, and flux measurement with credible intervals.

Overview

The high-energy spectral pipeline consists of four components:

  1. SimulateHighEnergyTransient — generate synthetic observations with Poisson statistics

  2. SpectralDataset / CountsSpectrumTransient — load and manage count spectra from OGIP files, simulations, or other formats

  3. Spectral models — physically motivated emission and absorption models in the keV range

  4. fit_model — Bayesian parameter estimation using Cash (C-stat), W-stat, or classic Gaussian (chi-square) likelihoods.

The interface mirrors the rest of redback: create a transient object, choose a model and prior, call fit_model, and inspect the result.

Simulating High-Energy Transients

SimulateHighEnergyTransient generates realistic X-ray/gamma-ray observations from any spectral model. It supports:

  • Time-tagged events (TTE): individual photon arrival times and energies

  • Binned photon counts: count rates in time–energy bins

  • Energy-dependent detector response: effective area as a scalar or callable

  • Background modeling: flat or energy-dependent background rate

  • Proper Poisson statistics: counts drawn from Poisson distributions

The model function must follow redback’s standard (time, frequency, **kwargs) interface, where frequency is energy in Hz (energy_keV × 2.417989e17).

Basic usage

import numpy as np
from redback.simulate_transients import SimulateHighEnergyTransient
from redback.transient_models import spectral_models

keV_to_Hz = 2.417989e17

def band_function_transient(times, **kwargs):
    """Band-function spectrum with a Gaussian time profile."""
    energies_keV = kwargs['frequency'] / keV_to_Hz
    spectrum_mjy = spectral_models.band_function_high_energy(
        energies_keV=energies_keV,
        redshift=0.2,
        log10_norm=-2.5,
        alpha=-1.0,
        beta=-2.3,
        e_peak=200.0,
    )
    time_profile = np.exp(-0.5 * ((times - 6.0) / 2.0) ** 2)
    return spectrum_mjy * time_profile

energy_edges = np.logspace(np.log10(10.0), np.log10(1500.0), 60)

sim = SimulateHighEnergyTransient(
    model=band_function_transient,
    parameters={},
    energy_edges=energy_edges,       # keV bin edges
    time_range=(0.1, 30.0),          # seconds
    effective_area=effective_area,   # callable or scalar (cm^2)
    background_rate=background_rate, # callable or scalar (counts/s/keV)
    time_resolution=0.01,            # seconds
    seed=42,
)

Generating counts

# Energy-integrated light curve
time_bins = np.logspace(np.log10(0.1), np.log10(30.0), 40)
lc = sim.generate_binned_counts(time_bins=time_bins, energy_integrated=True)
# Returns DataFrame with: time_center, count_rate, count_rate_error

# Count spectrum in a 1-second window around the peak
spectrum_bins = np.array([5.5, 6.5])
spectrum = sim.generate_binned_counts(time_bins=spectrum_bins, energy_integrated=False)
# Returns DataFrame with: counts, energy_center, energy_low, energy_high

# Time-tagged events
tte = sim.generate_time_tagged_events(max_events=150000)
# Returns DataFrame with: time, energy_keV, is_background

See examples/simulate_gamma_xray_spectra_ray.py for a full comparison of Band and blackbody models including plots.

Loading OGIP Data

Redback reads standard OGIP spectral files (PHA, RMF, ARF) produced by X-ray and gamma-ray telescope pipelines (Swift/XRT, XMM-Newton, Chandra, Fermi/GBM, etc.).

Both split (separate RMF + ARF) and combined response files (RSP, containing both redistribution matrix and effective area in a single file) are supported. Combined responses are common for gamma-ray instruments such as Fermi/GBM.

Runnable example

Redback includes a Fermi/GBM BGO example spectrum (simulated power-law GRB) in examples/example_data/:

import os
from redback.transient.spectral import CountsSpectrumTransient

data_dir = "examples/example_data"

spec = CountsSpectrumTransient.from_ogip(
    pha=os.path.join(data_dir, "ogip_powerlaw.pha"),
    bkg=os.path.join(data_dir, "ogip_powerlaw.bak"),
    rmf=os.path.join(data_dir, "ogip_powerlaw.rsp"),  # combined RSP
    name="gbm_grb",
)

dataset = spec.dataset
dataset.set_active_interval(200.0, 40000.0)  # 200 keV – 40 MeV (GBM BGO range)

The CountsSpectrumTransient wraps a SpectralDataset which holds the counts, response matrices, background, and metadata. Access it via spec.dataset.

General usage

from redback.transient.spectral import CountsSpectrumTransient

# Split RMF + ARF (typical for X-ray observatories)
spec = CountsSpectrumTransient.from_ogip(
    pha="spectrum.pha",
    bkg="background.pha",  # background spectrum (optional)
    rmf="response.rmf",    # redistribution matrix file
    arf="effective.arf",   # ancillary response file (optional)
    name="my_source",
)

# Combined RSP (typical for gamma-ray instruments)
spec = CountsSpectrumTransient.from_ogip(
    pha="spectrum.pha",
    bkg="background.pha",
    rmf="response.rsp",    # combined RSP passed as rmf
    name="my_source",
)

From a directory

If all OGIP files for a source live in one directory, from_ogip_directory will locate them automatically:

spec = CountsSpectrumTransient.from_ogip_directory(
    directory="data/my_source/",
    name="my_source",
)

Background files are found by trying common suffixes (bk.pha, _bkg.pha, _back.pha, -back.pha, _background.pha).

From a simulation

After running SimulateHighEnergyTransient, build a CountsSpectrumTransient directly, passing the time bins you want to integrate over:

import numpy as np
from redback.transient.spectral import CountsSpectrumTransient

time_bins = np.linspace(0.0, 100.0, 101)  # 1-second bins
spec = CountsSpectrumTransient.from_simulator(sim=sim, time_bins=time_bins, name="sim_source")

Setting the energy interval

Before fitting, restrict the analysis to a well-calibrated energy band:

dataset = spec.dataset
dataset.set_active_interval(0.3, 5.0)  # keV

Spectral Models

All high-energy spectral models take energies_keV as their first argument and return flux density in mJy. They are registered in redback.model_library and can be used by name with fit_model.

Available models

Calling models directly

All models are callable functions for exploration and simulation:

import numpy as np
from redback.transient_models.spectral_models import (
    powerlaw_high_energy,
    band_function_high_energy,
    blackbody_high_energy,
    tbabs_powerlaw_high_energy,
)

energies = np.logspace(-1, 2, 200)  # 0.1–100 keV

# Power law
flux = powerlaw_high_energy(energies, log10_norm=-3.0, alpha=-1.5, redshift=0.0)

# Band function
flux = band_function_high_energy(
    energies, log10_norm=-2.5, alpha=-1.0, beta=-2.3, e_peak=200.0, redshift=0.2
)

# Blackbody
flux = blackbody_high_energy(energies, kT=1.5, r_photosphere_rs=1.0, redshift=0.1)

# TBabs-absorbed power law (requires cross-section table)
flux = tbabs_powerlaw_high_energy(
    energies, log10_norm=-3.0, alpha=-1.5, nh=0.05, redshift=0.0
)

Absorption: TBabs

The tbabs_powerlaw_high_energy model includes interstellar absorption using the Tuebingen-Boulder ISM absorption model (Wilms et al. 2000). The column density parameter nh is in units of 10²² cm⁻². The cross-section table is loaded from redback/tables/xsect/ and cached after the first call.

# Fix nh at the Galactic value (from HI maps) by assigning a DeltaFunction prior
from bilby.core.prior import DeltaFunction
prior['nh'] = DeltaFunction(peak=0.04, name='nh')

Fitting Count Spectra

The fitting interface is the same as for photometric transients. Pass a CountsSpectrumTransient (or SpectralDataset) to fit_model along with a model name, prior, and sampler.

Choosing the fit statistic

The appropriate likelihood depends on whether background data are available:

Statistic

Use when

Notes

cstat

Background negligible or subtracted

Cash (1979); Poisson

wstat

Background spectrum available

Wachter et al. (1979); Poisson

chi2

Counts ≳ 20 per bin (Gaussian approximation)

Standard chi-square

auto

Automatically selects wstat (bkg) or cstat (no bkg)

Recommended default

Fitting from OGIP data

Using the bundled Fermi/GBM example data:

import os
import redback.priors
from redback.transient.spectral import CountsSpectrumTransient
from redback.sampler import fit_model

data_dir = "examples/example_data"

spec = CountsSpectrumTransient.from_ogip(
    pha=os.path.join(data_dir, "ogip_powerlaw.pha"),
    bkg=os.path.join(data_dir, "ogip_powerlaw.bak"),
    rmf=os.path.join(data_dir, "ogip_powerlaw.rsp"),
    name="gbm_grb",
)
spec.dataset.set_active_interval(200.0, 40000.0)

prior = redback.priors.get_priors("powerlaw_high_energy")
prior["redshift"] = 0.0  # fix to zero for unknown redshift

result = fit_model(
    transient=spec,
    model="powerlaw_high_energy",
    prior=prior,
    sampler="nestle",
    statistic="auto",   # selects wstat because background is present
    nlive=500,
)

result.plot_corner()

For X-ray data with ISM absorption (e.g. Swift/XRT, XMM-Newton), use tbabs_powerlaw_high_energy with split RMF/ARF files:

spec = CountsSpectrumTransient.from_ogip(
    pha="source.pha", bkg="background.pha",
    rmf="response.rmf", arf="effective.arf", name="xray_source",
)
spec.dataset.set_active_interval(0.3, 10.0)  # keV

prior = redback.priors.get_priors("tbabs_powerlaw_high_energy")
prior["redshift"] = 0.0

result = fit_model(transient=spec, model="tbabs_powerlaw_high_energy",
                   prior=prior, sampler="nestle", statistic="auto", nlive=500)

See examples/fit_ogip_spectrum.py for the complete runnable workflow including flux measurement.

Fitting simulated data

from redback.simulate_transients import SimulateHighEnergyTransient
from redback.transient.spectral import CountsSpectrumTransient
from redback.sampler import fit_model
import redback.priors
import numpy as np

# ... (set up sim as above) ...
time_bins = np.linspace(0.0, 100.0, 101)
spec = CountsSpectrumTransient.from_simulator(sim=sim, time_bins=time_bins, name="sim")
spec.dataset.set_active_interval(0.3, 5.0)

prior = redback.priors.get_priors("blackbody_high_energy")
prior["redshift"] = 0.1  # fix known redshift

result = fit_model(
    transient=spec,
    model="blackbody_high_energy",
    prior=prior,
    sampler="dynesty",
    statistic="auto",
    nlive=300,
)

See examples/fit_simulated_high_energy_blackbody.py for a complete end-to-end example.

Fitting GRB/prompt spectra

For time-resolved GRB spectroscopy, use SpectralDataset.from_simulator with short time bins to extract a time-integrated count spectrum, then fit with the Band function:

from redback.spectral.dataset import SpectralDataset
import redback.priors

time_bins = np.array([5.5, 6.5])  # 1-second window around peak
dataset = SpectralDataset.from_simulator(sim, time_bins=time_bins)

prior = redback.priors.get_priors("band_function_high_energy")

result = fit_model(
    transient=dataset,
    model="band_function_high_energy",
    prior=prior,
    sampler="dynesty",
    statistic="wstat",
    nlive=500,
)

See examples/fit_gamma_xray_spectra.py for the complete GRB prompt emission example.

Flux Measurement

After fitting, compute absorbed or unabsorbed fluxes with credible intervals from the posterior samples using dataset.compute_band_flux:

from redback.utils import calc_credible_intervals
import numpy as np

posterior = result.posterior
samples = posterior.sample(n=min(len(posterior), 500), random_state=0)

fluxes = []
for _, row in samples.iterrows():
    params = row.to_dict()
    fluxes.append(
        dataset.compute_band_flux(
            model="tbabs_powerlaw_high_energy",
            parameters=params,
            energy_min_keV=0.5,
            energy_max_keV=10.0,
            unabsorbed=False,   # set True to remove TBabs absorption
        )
    )

fluxes = np.asarray(fluxes)
lo, hi, med = calc_credible_intervals(fluxes, interval=0.68)
print(f"Flux (0.5–10 keV): {med:.3e} (+{hi-med:.3e}/-{med-lo:.3e}) erg/cm²/s")

Plotting

Spectral plots are produced via the CountsSpectrumTransient or SpectralDataset objects.

# Plot the raw count spectrum
spec.plot_data(
    filename="spectrum_data.png",
    min_counts=5,
    xscale="linear",
    xlim=(0.3, 5.0),
    plot_background=True,
)

# Plot the best-fit spectrum with residuals
spec.plot_fit(
    model="tbabs_powerlaw_high_energy",
    posterior=result.posterior,
    filename="spectrum_fit.png",
    min_counts=5,
    xscale="linear",
    xlim=(0.3, 5.0),
)

# Corner plot of posterior samples
result.plot_corner(filename="corner.png")

Priors

Pre-defined priors for all high-energy models are included in redback. Load them with redback.priors.get_priors:

import redback.priors

prior = redback.priors.get_priors("tbabs_powerlaw_high_energy")
print(prior)

# Fix a parameter to a known value
prior['redshift'] = 0.0

# Use a custom range
from bilby.core.prior import Uniform
prior['alpha'] = Uniform(minimum=-3.0, maximum=0.0, name='alpha', latex_label=r'$\alpha$')

Available prior files: powerlaw_high_energy, cutoff_powerlaw_high_energy, band_function_high_energy, blackbody_high_energy, comptonized_high_energy, tbabs_powerlaw_high_energy.

Examples

Four ready-to-run example scripts are provided in examples/:

Script

Description

simulate_gamma_xray_spectra_ray.py

Simulate Band and blackbody spectra, generate TTE and binned counts, compare models

fit_simulated_high_energy_blackbody.py

End-to-end: simulate a blackbody source, fit with Bayesian inference, recover parameters

fit_gamma_xray_spectra.py

Simulate a GRB Band spectrum, integrate over 1 s window, fit with wstat

fit_ogip_spectrum.py

Load bundled Fermi/GBM OGIP data (power-law GRB), fit with wstat, compute band flux