This is a fixed-text formatted version of a Jupyter notebook.
Fitting gammapy spectra with sherpa¶
Once we have exported the spectral files (PHA, ARF, RMF and BKG) in the OGIP format, it becomes possible to fit them later with gammapy or with any existing OGIP compliant tool such as XSpec or sherpa.
We show here how to do so with sherpa using the high-level user interface. For a general view on how to use stand-alone sherpa, see https://sherpa.readthedocs.io.
Load data stack¶
- We first need to import the user interface and load the data with load_data.
- One can load files one by one, or more simply load them all at once through a DataStack.
%matplotlib inline import matplotlib.pyplot as plt import glob # to list files from os.path import expandvars from sherpa.astro.datastack import DataStack import sherpa.astro.datastack as sh
WARNING: imaging routines will not be available, failed to import sherpa.image.ds9_backend due to 'RuntimeErr: DS9Win unusable: Could not find ds9 on your PATH' WARNING: failed to import sherpa.astro.xspec; XSPEC models will not be available
import sherpa sherpa.__version__
ds = DataStack() ANALYSIS_DIR = expandvars("$GAMMAPY_DATA/joint-crab/spectra/hess/") filenames = glob.glob(ANALYSIS_DIR + "pha_obs*.fits") for filename in filenames: sh.load_data(ds, filename) # see what is stored ds.show_stack()
read ARF file /Users/deil/work/gammapy-tutorials/datasets/joint-crab/spectra/hess/arf_obs23559.fits read RMF file /Users/deil/work/gammapy-tutorials/datasets/joint-crab/spectra/hess/rmf_obs23559.fits read background file /Users/deil/work/gammapy-tutorials/datasets/joint-crab/spectra/hess/bkg_obs23559.fits read ARF file /Users/deil/work/gammapy-tutorials/datasets/joint-crab/spectra/hess/arf_obs23523.fits read RMF file /Users/deil/work/gammapy-tutorials/datasets/joint-crab/spectra/hess/rmf_obs23523.fits read background file /Users/deil/work/gammapy-tutorials/datasets/joint-crab/spectra/hess/bkg_obs23523.fits read ARF file /Users/deil/work/gammapy-tutorials/datasets/joint-crab/spectra/hess/arf_obs23592.fits read RMF file /Users/deil/work/gammapy-tutorials/datasets/joint-crab/spectra/hess/rmf_obs23592.fits read background file /Users/deil/work/gammapy-tutorials/datasets/joint-crab/spectra/hess/bkg_obs23592.fits read ARF file /Users/deil/work/gammapy-tutorials/datasets/joint-crab/spectra/hess/arf_obs23526.fits read RMF file /Users/deil/work/gammapy-tutorials/datasets/joint-crab/spectra/hess/rmf_obs23526.fits read background file /Users/deil/work/gammapy-tutorials/datasets/joint-crab/spectra/hess/bkg_obs23526.fits 1: /Users/deil/work/gammapy-tutorials/datasets/joint-crab/spectra/hess/pha_obs23559.fits OBS_ID: 23559 MJD_OBS: N/A 2: /Users/deil/work/gammapy-tutorials/datasets/joint-crab/spectra/hess/pha_obs23523.fits OBS_ID: 23523 MJD_OBS: N/A 3: /Users/deil/work/gammapy-tutorials/datasets/joint-crab/spectra/hess/pha_obs23592.fits OBS_ID: 23592 MJD_OBS: N/A 4: /Users/deil/work/gammapy-tutorials/datasets/joint-crab/spectra/hess/pha_obs23526.fits OBS_ID: 23526 MJD_OBS: N/A
Define source model¶
We can now use sherpa models. We need to remember that they were designed for X-ray astronomy and energy is written in keV.
Here we start with a simple PL.
# Define the source model ds.set_source("powlaw1d.p1") # Change reference energy of the model p1.ref = 1e9 # 1 TeV = 1e9 keV p1.gamma = 2.0 p1.ampl = 1e-20 # in cm**-2 s**-1 keV**-1 # View parameters print(p1)
powlaw1d.p1 Param Type Value Min Max Units ----- ---- ----- --- --- ----- p1.gamma thawed 2 -10 10 p1.ref frozen 1e+09 -3.40282e+38 3.40282e+38 p1.ampl thawed 1e-20 0 3.40282e+38
Fit and error estimation¶
### Define the statistic sh.set_stat("wstat") ### Define the fit range ds.notice(0.6e9, 20e9) ### Do the fit ds.fit()
Datasets = 1, 2, 3, 4 Method = levmar Statistic = wstat Initial fit statistic = 641.651 Final fit statistic = 135.43 at function evaluation 193 Data points = 112 Degrees of freedom = 110 Probability [Q-value] = 0.0503069 Reduced statistic = 1.23118 Change in statistic = 506.222 p1.gamma 2.6166 p1.ampl 4.34008e-20 WARNING: parameter value p1.ampl is at its minimum boundary 0.0
Note that sherpa does not provide flux points. It also only provides plot for each individual spectrum.
sh.get_data_plot_prefs()["xlog"] = True sh.get_data_plot_prefs()["ylog"] = True ds.plot_fit()
WARNING: The displayed errorbars have been supplied with the data or calculated using chi2xspecvar; the errors are not used in fits with wstat WARNING: The displayed errorbars have been supplied with the data or calculated using chi2xspecvar; the errors are not used in fits with wstat WARNING: The displayed errorbars have been supplied with the data or calculated using chi2xspecvar; the errors are not used in fits with wstat WARNING: The displayed errorbars have been supplied with the data or calculated using chi2xspecvar; the errors are not used in fits with wstat
Errors and confidence contours¶
# Compute confidence intervals ds.conf()
p1.ampl lower bound: -2.08218e-21 p1.gamma lower bound: -0.0643092 p1.ampl upper bound: 2.14828e-21 p1.gamma upper bound: 0.0655594 Datasets = 1, 2, 3, 4 Confidence Method = confidence Iterative Fit Method = None Fitting Method = levmar Statistic = wstat confidence 1-sigma (68.2689%) bounds: Param Best-Fit Lower Bound Upper Bound ----- -------- ----------- ----------- p1.gamma 2.6166 -0.0643092 0.0655594 p1.ampl 4.34008e-20 -2.08218e-21 2.14828e-21
# Compute confidence contours for amplitude and index sh.reg_unc(p1.gamma, p1.ampl)
- Change the energy range of the fit to be only 2 to 10 TeV
- Fit the built-in Sherpa exponential cutoff powerlaw model
- Define your own spectral model class (e.g. powerlaw again to practice) and fit it