PyXspec  1.0.3
A Tutorial - Quick Version
Previous: Build and Install PyXspec Next: A Tutorial - Extended Version


This assumes the user already has a basic familiarity with both XSPEC and Python. Everything in PyXspec is accessible by importing the package xspec into your Python script.

Mac Snow Leopard and Lion users: Did you remember to make sure you're running the 32-bit version of Python, as shown in the Running on Mac OS X section of the "Build and Install" page?

PyXspec can be utilized in a Python script or from the command line of the plain interactive Python interpreter. PyXspec does not implement its own command handler, so it is NOT intended to be run as the Python equivalent of a traditional interactive XSPEC session (which is really an enhanced interactive Tcl interpreter). In other words you launch an interactive PyXspec session with:

UNIX>python
>>> import xspec
>>>

rather than:

UNIX>xspec
XSPEC12>

Note that in all the tutorial examples the xspec package name qualifier is left off. You must either include the xspec qualifier:

s = xspec.Spectrum("file1.pha")

or use a variation of the Python import or from...import commands:

from xspec import *
s = Spectrum("file1.pha")

Jumping In - The Really Quick Version

A simple Xspec load-fit-plot Python script may look something like this:

#!/usr/bin/python
from xspec import *

Spectrum("file1.pha")
Model("wabs*pow")
Fit.perform()
Plot.device = "/xs"
Plot("data")

Keeping this template in mind, we'll proceed to fill in the details...

Terminology

This description uses the standard Python object-oriented terminology, distinguishing between classes and objects. Class is used when referring to the type or definition of an object. An object refers to a specific instance of a class and is normally assigned to a variable. For example a user may load 3 data files by creating 3 spectral data objects s1, s2, and s3, which are all instances of the class Spectrum.

The functions and stored data members that make up the definition of a class are referred to as methods and attributes respectively.

The term Standard XSPEC refers to the traditional ways of using XSPEC, either with a Tcl script or an interactive XSPEC session.

Getting Help

There are two ways to get help for programming with PyXspec classes. The first is by viewing the Classes section of this manual. The Classes:Class List subsection is particularly useful as an entry point, as it contains hyperlinks to descriptions of every PyXspec class that is part of the public interface. The second way is to call Python's built-in help([class]) function from the interactive Python shell. Both methods will display essentially the same information, which originates in the class docstrings in the code files.

The 6 Global Objects

An XSPEC session fundamentally consists of loading data, fitting that data to a model, and plotting the results. To manage these operations, PyXspec offers the user 6 global objects: AllChains, AllData, AllModels, Fit, Plot, and Xset. Note that these are NOT the names of classes. They are instantiated objects of the class types shown in Table 1.

Table 1. PyXspec global objects
Object Name Class Role
AllChains ChainManager Monte Carlo Markov Chain container
AllData DataManager Container for all loaded data sets (objects of class Spectrum)
AllModels ModelManager Container for all Model objects
Fit FitManager Manager class for setting properties and running a fit
Plot PlotManager Manager class for performing XSPEC plots
Xset XspecSettings Storage class for Xspec settings

PyXspec instantiates these objects immediately upon the importing of the xspec package. You cannot create any other objects of these class types, as they each allow only 1 instance of their type. (They are singletons in the language of design patterns.)

Operations involving these should ALWAYS be performed through the objects and NOT their class names. These class names should never appear in your code.

Loading And Removing Data

Spectral data files can be loaded in several ways. You can create an object of the Spectrum class by passing it the data file name:

s1 = Spectrum("file1.pha")

which also adds the new object s1 to the AllData container. Or you can simply add the new file directly to the container without retrieving a Spectrum object:

AllData += "file1.pha"

Later you can always obtain a Spectrum object reference to any of the loaded spectra by passing AllData an integer:

s2 = AllData(2) # s2 is a reference to the 2nd loaded spectrum

For more complicated data loading, you have access to the same functionality in Standard XSPEC's data command. Simply pass a string to the AllData object's __call__ method:

AllData("file1 file2 2:3 file3")

Note that only the last example allows you to assign multiple data groups, the 3rd spectrum being assigned to data group 2. Also note that in the last example any previously loaded data sets are removed, thus reproducing the behavior of Standard XSPEC's data command.

Other ways of removing Spectrum objects (ie. data sets) from the container:

AllData -= 3 # Removes the 3rd Spectrum object (the spectrum with index number 3) from the container.
AllData -= s1 # Removes the Spectrum object s1.
AllData -= "*" # Removes all Spectrum objects.
AllData.clear() # Removes all Spectrum objects.

You can check the current state of the AllData container at any time by doing:

AllData.show()

Similarly, to view information about a single Spectrum object:

s2.show()

Defining Models

The basic way of defining an XSPEC model is to create an object of the PyXspec class Model. Simply pass in a string containing a combination of 1 or more XSPEC model components. Since this uses the same syntax as Standard XSPEC's model command, component abbreviations are allowed:

m1 = Model("wa*po + ga")

and to see a complete listing of available XSPEC model components, do:

Model.showList()

When you define a model like this, PyXspec also automatically adds the new object to the global AllModels container. If the model is applied to multiple data groups, object copies are added to the container for each data group.

Similar to the case of spectral data, you can also load models directly into the global container:

# Another way to define a new model and create an object for each data group.
AllModels += "wa*po + ga"
# Retrieve the model object assigned to data group 2.
m2 = AllModels(2)
# Various ways to remove all model objects from the container.
AllModels.clear()
AllModels -= "*"

To display models and their parameters:

# This displays all parameters in all model objects:
AllModels.show()
# While this displays just parameters 1,2,3 and 5:
AllModels.show("1-3, 5")
# This displays a single model object:
m2.show()

For defining mulitple (or named) models and assigning multiple sources, please see the Extended Tutorial section.

Component and Parameter Objects

Model objects contain Component objects and Component objects contain Parameter objects. There are several ways to access and set components and parameters individually (and if you want to change many parameter values at once, it may be faster to use the Model or AllModels setPars methods described in the next section). Examples of individual Component and Parameter object access:

# Component objects are accessible-by-name as Model object attributes*:
comp1 = m1.wabs
comp2 = m1.powerlaw
# Parameter objects are accessible-by-name as Component object attributes:
par4 = m1.gaussian.LineE
# ...and we can modify their values:
par4.values = 3.895
m1.wabs.nH = 5.0
comp2.PhoIndex = 1.5

# Can also get a Parameter object directly from a Model, without going through a Component.
# Just pass the Model the Parameter index number:
par5 = m1(5)
# Note change from Beta version: Parameters are now numbered from 1 in EACH model object.
# Example: If m1 and m2 are the data group 1 and data group 2 objects for a 5-parameter model,
# m2's parameters are now also numbered 1-5 and NOT 6-10.

# Examples of numerical operations allowed with Parameter objects:
par4 += 0.75
par4 *= 2.0
y1 = m1.wabs.nH*100.0
y2 = par4 + par5

(*)For models with duplicate copies of components, see the Extended Tutorial for accessing Component objects by name.

Note that in the above examples, only the parameter's value is being accessed or modified. To change all or part of its FULL list of settings including auxiliary values: value, fit delta, min, bot, top, max, you can set its values attribute to a tuple or list of size 1-6:

par4.values = 4.3, .01, 1e-3
par4.values = [4.3, .01, 1e-3, 1e-2, 100, 200]

Or for greater flexibility you can set it to a string using Standard XSPEC's newpar command syntax:

# This allows you to set new values non-consecutively.
par4.values = "1.0, -.01,,,,150"

A quick way to freeze or thaw a parameter is to toggle its frozen attribute:

par4.frozen = False
par5.frozen = True

To link a parameter to one or more others, set its link attribute to a link expression string as you would have with the newpar command. To remove the link, set link to an empty string or call the parameter's untie method.

par5.link = "2.3 * 4" # Link par 5 to par 4 with a multiplicative constant.
par5.link = "" # Removes the link.
par5.untie() # Also removes the link.

Also ALL linked parameters in a model object can be untied with a single call to the Model class untie method.

To display a parameter's full set of values (including auxiliary values), just print its values attribute:

>>> print par4.values
[6.5, 0.05, 0.0, 0.0, 1000000.0, 1000000.0]
>>>

Setting Multiple Parameters At A Time

You can set multiple parameter values with a single call using the Model or AllModels setPars methods. This may be considerably faster than setting parameters one at a time through the individual Parameter objects as shown in the previous section. With setPars, the model will be recalculated just ONCE after all the changes have been made. But when setting through individual Parameter objects, the model will be recalculated after EACH parameter change.

# For Model object m1, supply 1 or more new parameter values in consecutive order:
m1.setPars(2.5, 1.4, 1.0e3) # This changes pars 1, 2, and 3.

# Can also change paramater auxiliary values by passing a string using the same
# syntax as with Standard XSPEC's newpar command:
m1.setPars(.95, "1.8,,-5,-4,10,10")

# Now set parameters NON-CONSECUTIVELY by passing a Python dictionary object.
# This has been changed from the Beta version. No longer using the p[n] syntax.
# This example changes pars 1, 2, 4, and 6:
m1.setPars(.95, 1.2, {4:9.8, 6:2.0})

Parameters can also be initialized by passing values to the Model object constructor. You do this by setting the Model constructor's setPars keyword argument to a tuple, list, or dictionary (or just a single value or string if only setting the first parameter):

# Supply values for parameters 1 and 3, use defaults for the rest.
m = Model("wa*ga", setPars={1:1.5, 3:.2})
# Supply values for 1 and 2, use defaults for the rest.
m = Model("wa*ga", setPars=(1.5, 0.7))
# Supply value only for 1.
m = Model("wa*ga", setPars=1.5)

Finally, if you wish to set multiple parameters that belong to different model objects, you must use the AllModels container's setPars method. This follows the same syntax rules as the single Model setPars, except that you also supply the Model objects as arguments:

# Change pars 1 and 3 in m1, and pars 1 and 2 in m2:
AllModels.setPars(m1, {1:6.4, 3:1.78}, m2, 3.5, 0.99)

Fitting

Once data and models are loaded, fitting is performed by calling the perform method of the Fit global object:

Fit.perform()

Some of the more frequently modified fit settings are the type of statistic to minimize and the maximum number of fit iterations to perform. These settings are attributes of Fit:

Fit.nIterations = 100
Fit.statMethod = "cstat"
Fit.statMethod = "chi"

Please see the class reference guide and the Extended Tutorial for Fit's complete functionality.

To display the fit results at any time:

Fit.show()

Plotting

In Standard XSPEC, plot settings are adjusted using the setplot command while the plot is displayed through the plot command. In PyXspec, all plot settings and functionality is handled through the global Plot object. A device must be set before any plots can be displayed, and this done through the device attribute:

Plot.device = "/xs"

The device can also be set to print to an output file in several formats. The list of possible devices is given by the cpd command in the Standard XSPEC manual.

A typical setting to adjust is the X-axis units. You can choose to plot channel numbers, or select from various energy and wavelength units. The strings can also be abbreviated. Examples:

Plot.xAxis = "channel"
Plot.xAxis = "MeV"
Plot.xAxis = "Hz"
Plot.xAxis = "angstrom"

The displays of individual additive components or background spectra is toggled by setting their attributes to a bool:

Plot.add = True
Plot.background = False

Similarly log/linear settings for data plots (when using energy or wavelength units):

Plot.xLog = True
Plot.yLog = False

The current plot settings are displayed with:

Plot.show()

To actually display a plot, send 1 or more string arguments to the Plot __call__ method:

# Single panel plots
Plot("data")
Plot("model")
Plot("ufspec")
# Multi panel plots
Plot("data chisq")
Plot("data","model","resid")
# Call Plot with no arguments to repeat the previously entered Plot command
Plot()

After displaying a plot, you can get an array of the plotted values by calling one of Plot's retrieval methods. All of these functions take an optional plot group number argument for the case of multiple plot groups, and all return the plot values in a Python list.

Plot("data")
xVals = Plot.x()
yVals = Plot.y()
yVals2 = Plot.y(2) # Gets values for data in the second plot group
modVals = Plot.model()
# To get a background array, Plot.background must be set prior to plot
Plot.background = True
Plot("data")
bkg = Plot.backgroundVals()
# Retrieve error arrays
xErrs = Plot.xErr()
yErrs = Plot.yErr()