Last modified: December 2023

URL: https://cxc.cfa.harvard.edu/ciao/ahelp/convert_xspec_script.html
AHELP for CIAO 4.16

convert_xspec_script

Context: Tools::Utilities

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:

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

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.