Synopsis
Convert a XSPEC save file to Sherpa commands *experimental*
Syntax
convert_xspec_script infile outfile Unlike most CIAO contributed scripts, there is no parameter file. If infile or outfile are - then stdin or stdout are used, respectively. The supported command-line flags can be found using -h/--help: --models to include models created by convert_xspec_user_script (can be used multiple times) --clean if the script should start with a call to clean --clobber will overwrite existing files --verbose will change the amount of screen output --version prints the script version --copyright prints the script copyright
Description
The convert_xspec_script tool will attempt to convert a XSPEC save file (often ending in .xcm) into Sherpa commands. The resulting file can then be run from sherpa or modified for further analysis. The results should be considered as the start of a Sherpa analysis, as further work may be needed.
This script is *experimental* and will not work for all XSPEC commands or datasets. Please contact the CXC HelpDesk if you find there are scripts which this script will not convert, or where the results do not appear to match XSPEC.
Examples
Example 1
unix% convert_xspec_script savespec.xcm example.py model phabs*powerlaw
The script displays each model line it processes. With an input file (savexspec.xcm) containing:
method leven 10 0.01
abund angr
xsect vern
cosmo 70 0 0.73
xset delta 0.01
systematic 0
model phabs*powerlaw
0.0115491 0.001 0 0 100000 1e+06
1.98713 0.01 -3 -2 9 10
0.000185474 0.01 0 0 1e+20 1e+24
bayes off
then the script will produce:
from sherpa.astro.ui import *
set_method('levmar')
set_xsabund('angr')
set_xsxsect('vern')
set_xscosmo(70, 0, 0.73)
# model phabs*powerlaw
m1 = create_model_component('xsphabs', 'm1')
m2 = create_model_component('xspowerlaw', 'm2')
set_par(m1.nH, 0.0115491)
set_par(m2.PhoIndex, 1.98713)
set_par(m2.norm, 0.000185474, max=1e+20)
# Set up the model expressions
#
set_source(1, m1 * m2)
Example 2
unix% convert_xspec_script bestfit.xcm bestfit.py model phabs*apec unix% sherpa ... sherpa In [1]: %run bestfit sherpa In [2]: plot_fit(1, ylog=True, alpha=0.5) sherpa In [2]: plot_fit(2, overplot=True, alpha=0.5)
The example above shows how the converted script can be used in Sherpa. In this particular case the XSPEC save file was:
unix% cat bestfit.xcm
statistic chi
data 1:1 obs1.pi
data 1:2 obs2.pi
ignore 1:1-3,76-301 2:1-2,110-376
method leven 10 0.01
abund angr
xsect vern
cosmo 70 0 0.73
xset delta 0.01
systematic 0
model phabs*apec
0.252632 0.001 0 0 100000 1e+06
0.85796 0.01 0.008 0.008 64 64
0.359129 0.001 0 0 5 5
0.0175 -0.01 -0.999 -0.999 10 10
0.000618978 0.01 0 0 1e+20 1e+24
bayes off
and the converted file is
unix% cat bestfit.py
from sherpa.astro.ui import *
from sherpa_contrib.xspec import xcm
set_stat('chi2datavar')
load_pha(1, 'obs1.pi', use_errors=True)
load_pha(2, 'obs2.pi', use_errors=True)
xcm.ignore(1, 1, 3)
xcm.ignore(1, 76, 301)
xcm.ignore(2, 1, 2)
xcm.ignore(2, 110, 376)
set_method('levmar')
set_xsabund('angr')
set_xsxsect('vern')
set_xscosmo(70, 0, 0.73)
xcm.subtract(1)
xcm.subtract(2)
# model phabs*apec
m1 = create_model_component('xsphabs', 'm1')
m2 = create_model_component('xsapec', 'm2')
set_par(m1.nH, 0.252632)
set_par(m2.kT, 0.85796)
set_par(m2.Abundanc, 0.359129)
set_par(m2.redshift, 0.0175, frozen=True)
set_par(m2.norm, 0.000618978, max=1e+20)
# Set up the model expressions
#
set_source(1, m1 * m2)
set_source(2, m1 * m2)
Note that the conversion uses several routines from sherpa_contrib.xspec.xcm, which provide some Sherpa-like functions which simplify the conversion of the XSPEC script.
Example 3
CIAO 4.13 adds support for XSPEC convolution models, so the model expression
model phabs*cflux*powerlaw
will be converted to
# model phabs*cflux*powerlaw
m1 = create_model_component('xsphabs', 'm1')
m2 = create_model_component('xscflux', 'm2')
m3 = create_model_component('xspowerlaw', 'm3')
# Set up the model expressions
#
set_source(1, m1 * m2(m3))
and
model cflux*phabs*powerlaw
will create
# model cflux*phabs*powerlaw
m1 = create_model_component('xscflux', 'm1')
m2 = create_model_component('xsphabs', 'm2')
m3 = create_model_component('xspowerlaw', 'm3')
# Set up the model expressions
#
set_source(1, m1(m2 * m3))
Example 4
unix% echo "model phabs(apec)" | convert_xspec_script - - --verbose 0
from sherpa.astro.ui import *
# model phabs(apec)
m1 = create_model_component('xsphabs', 'm1')
m2 = create_model_component('xsapec', 'm2')
# Set up the model expressions
#
set_source(1, m1 * (m2))Convert the input from stdin - in this case the string "model phabs(apec)" - and write the output to stdin. Using the "--verbose 0" option means that messages from the script are not included in the screen output.
Example 5
Added in 4.15.2 is support for the MDEFINE command. Here we show the output from a simple additive model of two power laws:
mdefine dplaw E**p1 + f*E**p2 model dplaw
Converting this file will create the following Python script:
from sherpa.astro.ui import *
# mdefine dplaw E**p1 + f*E**p2 : ADD
# parameters: p1, f, p2, norm
def model_dplaw(pars, elo, ehi):
p1 = pars[0]
f = pars[1]
p2 = pars[2]
norm = pars[3]
elo = np.asarray(elo)
ehi = np.asarray(ehi)
E = (elo + ehi) / 2
de = ehi - elo
return norm * (E**p1 + f*E**p2) * de
# model dplaw
load_user_model(model_dplaw, 'm1')
add_user_pars('m1', ['p1', 'f', 'p2', 'norm'])
# Parameter settings
# Set up the model expressions
#
set_source(1, m1)
Using XSPEC user models
If the MODEL or MDEFINE definitions include XSPEC user models then:
- the model must be converted to Sherpa with the script convert_xspec_user_model,
- and the --models option is used to tell the script what module, or modules, to be loaded.
For any model that was created by convert_xspec_user_model with the call
convert_xspec_user_model modname lmodel.dat
add a
--models modname
option to the call to xspec_user_script. This argument can be repeated if several sets of models are needed; for example
--models mod1 --models mod2
.
MDEFINE
Support for the MDEFINE command has been added in the 4.15.2 release. The results should be reviewed carefully as there has been limited testing, in particular differences between the various XSPEC language features and Python.
Calling XSPEC functions
MDEFINE expressions can call various unary and binary functions, such as ATAN2, EXP, HEAVISIDE, and SMAX. Although support is provided for these models, they are not guaranteed to behave as they do in XSPEC, so please check the results carefully!
Example
With the following XCM script, based on the XSPEC MDEFINE examples,
mdefine junk a*e+b*log(e)/sin(e) mdefine junk3 0.2+B*e : mul mdefine mymod junk3(p1)*junk(p2,p3) model mymod
running convert_xspec_script will display
model mymod Unable to find parameter value for m1.p1 - skipping other parameters Found MDEFINE command: please consult 'ahelp convert_xspec_script' - junk ADD : check definition of model_junk() - junk3 MUL : check definition of model_junk3() - mymod ADD : check definition of model_mymod()
and create the output file
from sherpa.astro.ui import *
# mdefine junk a*e+b*log(e)/sin(e) : ADD
# parameters: a, b, norm
def model_junk(pars, elo, ehi):
a = pars[0]
b = pars[1]
norm = pars[2]
elo = np.asarray(elo)
ehi = np.asarray(ehi)
E = (elo + ehi) / 2
de = ehi - elo
return norm * (a*E+b*np.log10(E)/np.sin(E)) * de
# mdefine junk3 0.2+B*e : MUL
# parameters: B
def model_junk3(pars, elo, ehi):
B = pars[0]
elo = np.asarray(elo)
ehi = np.asarray(ehi)
E = (elo + ehi) / 2
return 0.2+B*E
# mdefine mymod junk3(p1)*junk(p2,p3) : ADD
# parameters: p1, p2, p3, norm
def model_mymod(pars, elo, ehi):
p1 = pars[0]
p2 = pars[1]
p3 = pars[2]
norm = pars[3]
elo = np.asarray(elo)
ehi = np.asarray(ehi)
def junk3(*args):
pars = list(args)
out = model_junk3(pars, elo, ehi)
return out
def junk(*args):
pars = list(args)
pars.append(1.0) # model is additive
out = model_junk(pars, elo, ehi)
out /= de # model is additive
return out
E = (elo + ehi) / 2
de = ehi - elo
return norm * (junk3(p1)*junk(p2,p3)) * de
# model mymod
load_user_model(model_mymod, 'm1')
add_user_pars('m1', ['p1', 'p2', 'p3', 'norm'])
# Parameter settings
# Set up the model expressions
#
set_source(1, m1)
Known Problems
This script is intended to simplify analysis in Sherpa, but is not guaranteed to produce the same results. In part this is due to missing functionality, but other differences include
- not all parts of the various XSPEC expression languages - such as MODEL and MDEFINE - are suppoted,
- differences in how Sherpa and XSPEC work,
- possible differences in your ~/.xspec/Xspec.init file (as of CIAO 4.12 Sherpa now uses this file),
- and differences between the XSPEC version and that of CIAO (12.12.1c for CIAO 4.15).
The script is also designed to parse the output of the SAVE command from XSPEC. It can handle hand-edited files but it is more likely to fail or work incorrectly. For instance, XSPEC commands must be included in full, so a line like
mo (phabs)apec
will not work, but
model (phabs)apec
will.
Tied parameters
Comlplicated tie expressions, in particular those involving functions, are not guaranteed to be converted correctly.
Unsupported models
Only additive, multiplicative, and convolution models are recognized. Only models included in XSPEC 12.12.1c (used in CIAO 4.15) are available, although the --models argument, added in the 4.15.2 release, can be used to add XSPEC user models that have been converted with convert_xspec_user_model.
Unsupported commands
Some commands are ignored by the script, but others will cause the message
SKIPPING '...'
to be displayed while converting the file.
Multiple files
There is currently no support for files containing multiple data sets (e.g. PHA-II format files). In particular the {} syntax is not support for commands like DATA and RESPONSE.
Further help
Please contact the CXC HelpDesk if you have a problem with the script.
Changes in the scripts 4.16.0 (December 2023) release
Initial support of the ENERGIES command
The ENERGIES command extends the energy range used to evaluate models. There is now preliminary support for this command, although not all modes are supported (in particular the EXTEND LOW/HIGH command fails), and this command has only been lightly tested. Please contact the CXC HelpDesk if there appears to be problems.
Changes in the scripts 4.15.2 (June 2023) release
Allow scripts that reference user models
If the XCM scripts references one or XSPEC user models, then the --models argument must be used to tell the script about these models (this assumes that convert_xspec_user_script was used to convert them). So if
convert_xspec_user_model javier lmodel.dat
was used to convert the extra models, then add
--models javier
to the call to convert_xspec_script.
Initial support of the MDEFINE command
XSPEC allows user models to be created with the MDEFINE command. The Sherpa approach to creating user models is somewhat different, which requires converting from the XSPEC expression into one supported by Sherpa. This is the first attempt at supporting such models, and so the results should be checked against XSPEC. Convolution models are currently unsupported and models that call functions (such as ATAN2, DIM, or HEAVISIDE) may not work as expected because of differences between XSPEC and Python. Please contact the CXC HelpDesk if there appear to be problems or differences to XSPEC when using the output of this script.
Changes to the notice and ignore commands
To match changes made in CIAO 4.15, calls to ignore and notice in the XCM script will now display the change in the filter, such as reporting
dataset 1: 0:15.01 -> 0.32:15.01 Energy (keV) dataset 1: 0.32:15.01 -> 0.32:9.92 Energy (keV)
It turns out that Sherpa and XSPEC handle PHA files with invalid grouping data (the GROUPING and QUALITY columns) differently. When running the output of convert_xspec_script, warning messages will be displayed when the filters used by XSPEC do not match those used by Sherpa, saying
WARNING: Spectrum 1: the grouping scheme does not match OGIP standards. The filtering may not exactly match XSPEC.
Improvements in model parsing
The XSPEC model language is different to Sherpa, for example a model can be expressed as
model constant*phabs( ( zpowerlw )etable{mytorus_Ezero_v00.fits} )which used to get converted to a model expression like
set_source(1, m1 * m2 * ( * (m3) * m4))
but will now create
set_source(1, m1 * m2 * ((m3) * m4))
Changes in the scripts 4.15.1 (January 2023) release
Better handling of un-supported models. Please contact the CXC HelpDesk if you have XCM files with this problem. The script should once again be able to handle multiple groups (creating a model instance for each group).
Changes in the scripts 4.14.0 (December 2021) release
Updated for changes in Sherpa. Support XSPEC table models using atable, mtable, and etable expressions. Input and output can no be taken from stdin and stdout respectively by using the "-" name. Added the --verbose flag.
Changes in the scripts 4.13.0 (December 2020) release
The script is new in this release.
Notes
This script is not an official part of the CIAO release but is made available as "contributed" software via the CIAO scripts page. Please see the installation instructions page for help on installing the package.