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.
%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
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()
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)
### Define the statistic
sh.set_stat("wstat")
### Define the fit range
ds.notice(0.6e9, 20e9)
### Do the fit
ds.fit()
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()
# Compute confidence intervals
ds.conf()
# Compute confidence contours for amplitude and index
sh.reg_unc(p1.gamma, p1.ampl)