Spectrum estimation and modeling (gammapy.spectrum)

Introduction

gammapy.spectrum holds functions and classes related to 1D region based spectral analysis. This includes also simulation tools.

The basic of 1D spectral analysis are explained in this talk. A good reference for the forward-folding on-off likelihood fitting methods is Section 7.5 “Spectra and Light Curves” in [Naurois2012], in publications usually the reference [Piron2001] is used. A standard reference for the unfolding method is [Albert2007].

Getting Started

The following code snippet demonstrates how to load an observation stored in OGIP format and fit a spectral model.

import astropy.units as u
from gammapy.datasets import gammapy_extra
from gammapy.spectrum import SpectrumObservation, SpectrumFit, models

filename = '$GAMMAPY_EXTRA/datasets/hess-crab4_pha/pha_obs23592.fits'
obs = SpectrumObservation.read(filename)

model = models.PowerLaw(
    index=2 * u.Unit(''),
    amplitude=1e-12*u.Unit('cm-2 s-1 TeV-1'),
    reference=1*u.TeV,
)
fit = SpectrumFit(obs_list=obs, model=model)
fit.fit()
fit.est_errors()
print(fit.result[0])

It will print the following output to the console:

Fit result info
---------------
Model: PowerLaw
ParameterList
Parameter(name='index', value=2.1473880540790522, unit=Unit(dimensionless), min=0, max=None, frozen=False)
Parameter(name='amplitude', value=2.7914083679020973e-11, unit=Unit("1 / (cm2 s TeV)"), min=0, max=None, frozen=False)
Parameter(name='reference', value=1.0, unit=Unit("TeV"), min=None, max=None, frozen=True)

Covariance: [[  6.89132245e-03   1.12566759e-13   0.00000000e+00]
 [  1.12566759e-13   7.26865610e-24   0.00000000e+00]
 [  0.00000000e+00   0.00000000e+00   0.00000000e+00]]

Statistic: 46.051 (wstat)
Fit Range: [  5.99484250e+08   1.00000000e+11] keV

For more advanced use cases please go to the tutorial notebooks:

Content

Reference/API

gammapy.spectrum Package

Spectrum estimation and modeling methods (1-dimensional, with an energy axis).

Functions

check_chi2() Execute this function after fitting to see if the
chi2asym_err_func(data) Compute statistical error per bin from the data.
chi2asym_stat_func(data, model[, staterror, ...]) Define asymmetric chi-square errors.
compute_flux_points_dnde(flux_points, model) Compute differential flux points quantities.
cosmic_ray_flux(energy[, particle]) Cosmic ray flux at Earth.
diffuse_gamma_ray_flux(energy[, component]) Diffuse gamma ray flux.
integrate_spectrum(func, xmin, xmax[, ...]) Integrate 1d function using the log-log trapezoidal rule.
load_chi2asym_stat() “Load and set the chi2asym statistic

Classes

CountsPredictor(model[, aeff, edisp, ...]) Calculate npred
CountsSpectrum(energy_lo, energy_hi[, data, ...]) Generic counts spectrum.
CrabSpectrum([reference]) Crab spectral model.
FluxPointEstimator(obs, groups, model) Flux point estimator.
FluxPoints(table) Flux point object.
FluxPointsFitter([stat, optimizer, ...]) Fit a set of flux points with a parametric model.
LogEnergyAxis(energy[, mode]) Log energy axis.
PHACountsSpectrum(energy_lo, energy_hi[, ...]) OGIP PHA equivalent
PHACountsSpectrumList List of PHACountsSpectrum
SEDLikelihoodProfile(table) SED likelihood profile.
SpectrumButterfly([data, masked, names, ...]) Spectral model butterfly class.
SpectrumEnergyGroup(energy_group_idx, ...) Spectrum energy group.
SpectrumEnergyGroupMaker(obs) Energy bin groups for spectral analysis.
SpectrumEnergyGroups([initlist]) List of SpectrumEnergyGroup objects.
SpectrumExtraction(obs_list, bkg_estimate[, ...]) Creating input data to 1D spectrum fitting
SpectrumFit(obs_list, model[, stat, ...]) Spectral Fit
SpectrumFitResult(model[, fit_range, ...]) Class representing the result of a spectral fit
SpectrumObservation(on_vector[, aeff, ...]) 1D spectral analysis storage class
SpectrumObservationList([initlist]) List of SpectrumObservation.
SpectrumObservationStacker(obs_list) Stack SpectrumObervationList
SpectrumResult(model, points) Class holding all results of a spectral analysis
SpectrumSimulation(livetime, source_model, aeff) Simulate SpectrumObservation.
SpectrumStats(**kwargs) Spectrum stats.

gammapy.spectrum.models Module

Spectral models for Gammapy.

Classes

SpectralModel Spectral model base class.
PowerLaw(index, amplitude, reference) Spectral power-law model.
PowerLaw2(amplitude, index, emin, emax) Spectral power-law model with integral as norm parameter
ExponentialCutoffPowerLaw(index, amplitude, ...) Spectral exponential cutoff power-law model.
ExponentialCutoffPowerLaw3FGL(index, ...) Spectral exponential cutoff power-law model used for 3FGL.
PLSuperExpCutoff3FGL(index_1, index_2, ...) Spectral super exponential cutoff power-law model used for 3FGL.
LogParabola(amplitude, reference, alpha, beta) Spectral log parabola model.
TableModel(energy, values[, amplitude, ...]) A model generated from a table of energy and value arrays.
AbsorbedSpectralModel(spectral_model, ...[, ...]) Spectral model with EBL absorption.
Absorption(energy_lo, energy_hi, param_lo, ...) Gamma-ray absorption models.

gammapy.spectrum.powerlaw Module

Power law spectrum helper functions.

This gammapy.spectrum.powerlaw module contains a bunch of helper functions for computations concerning power laws that are standalone and in their separate namespace (i.e. not imported into the gammapy.spectrum namespace). Their main purpose is to serve as building blocks for implementing other code in Gammapy that depends on power laws, e.g. interpolation or extrapolation or integration of spectra or spectral cubes.

End users will rarely need to use the functions here, the gammapy.spectrum.models.PowerLaw class, which is a gammapy.spectrum.models.SpectralModel sub-class provides the common functionality with a convenient API, that can be somewhat slower though than the functions here, because it uses Quantity objects.

TODO: probably it doesn’t make sense to keep power_law_evaluate here ... it’s a direct duplication of the staticmethod PowerLaw.evaluate.

Examples

All the functions contain the powerlaw_ prefix and are pretty long, so we suggest importing them directly and using them like this:

>>> from gammapy.spectrum.powerlaw import power_law_evaluate
>>> power_law_evaluate(energy=42, norm=4.2e-11, gamma=2.3, energy_ref=1)
7.758467093729267e-15

Functions

power_law_evaluate(energy, norm, gamma, ...) Differential flux at a given energy.
power_law_pivot_energy(energy_ref, f0, ...) Compute pivot (a.k.a.
power_law_df_over_f(e, e0, f0, df0, dg, cov) Compute relative flux error at any given energy.
power_law_flux([I, g, e, e1, e2]) Compute differential flux for a given integral flux.
power_law_energy_flux(I[, g, e, e1, e2]) Compute energy flux between e1 and e2 for a given integral flux.
power_law_integral_flux([f, g, e, e1, e2]) Compute integral flux for a given differential flux.
power_law_g_from_f(e, f[, de]) Spectral index at a given energy e for a given function f(e)
power_law_g_from_points(e1, e2, f1, f2) Spectral index for two given differential flux points
power_law_I_from_points(e1, e2, f1, f2) Integral flux in energy bin for power law
power_law_f_from_points(e1, e2, f1, f2, e) Linear interpolation
power_law_f_with_err([I_val, I_err, g_val, ...]) Wrapper for f so the user doesn’t have to know about
power_law_I_with_err([f_val, f_err, g_val, ...]) Wrapper for f so the user doesn’t have to know about
power_law_compatibility(par_low, par_high) Quantify spectral compatibility of power-law measurements in two energy bands.