Extended Tutorial

This assumes the user is familiar with the basics of PyXspec as explained in the Quick Tutorial.

Data

Background, Response, and Arf

When a Spectrum object is created from a spectral data file, PyXspec also reads the file's BACKFILE, RESPFILE, and ANCRFILE keywords and will load the corresponding background, response, and arf files. The spectrum's Background and Response objects are then available as attributes of the Spectrum class, while the arf file name becomes an attribute of the Response class:

s1 = Spectrum("file1")
b1 = s1.background
r1 = s1.response
arfFileName = r1.arf

Note

You never create Background and Response objects directly. They are accessible only through the Spectrum class attributes.

These attributes may also be used to add, change, or remove auxiliary files to an existing Spectrum object:

# Add or replace files:
s1.background = "newBackground.pha"
s1.response.arf = "newArf.pha"
# Removal examples:
s1.response = None
s1.background = ""

Background and Spectrum store their original file names in their fileName attribute. This means that while you SET the Spectrum.background object by assigning it a file name (as shown above), to GET the file name you must access its fileName attribute:

bkgFileName = s1.background  # Wrong!!! This returns the entire Background
                                # object, not a string.
bkgFileName = s1.background.fileName  # Correct

Response stores its RMF and optional ARF file names in its rmf and arf attributes respectively:

rmfFileName = r1.rmf
arfFileName = r1.arf

Background objects have some of the same attributes as Spectrum objects, such as areaScale, exposure, energies, and values. The Spectrum object's values array (actually a tuple) does NOT include contributions from the background. Those are stored separately in the associated Background object.

The Spectrum class also provides a multiresponse array attribute for assigning multiple detectors (or sources) to a spectrum. The standard 0-based Python array indices corresponding to the 1-based XSPEC source numbers:

# Set a response for source 2
s1.multiresponse[1] = "resp2.rsp"
# Get the response object for source 2
r2 = s1.multiresponse[1]
# Remove the response from source 2
s1.multiresponse[1] = None
# This is the same as doing s1.response = "resp1.rsp"
s1.multiresponse[0] = "resp1.rsp"

The rule is: when doing single-source analysis (typical of most XSPEC sessions) use the response attribute, otherwise use the multiresponse array.

Ignore/Notice

To ignore channels for a SINGLE spectrum, call the Spectrum object's ignore method passing a string following the same syntax as for Standard XSPEC's ignore command:

s1.ignore("20-30 50-**")
s1.ignore("**-5")

Similarly, to notice channels in a single spectrum:

s1.notice("10-30,80-**")
s1.notice("all")

As with Standard XSPEC, if the x-axis plot units are set to energies or wavelengths, ignore and notice will accept floating-point input assumed to be in those same units:

Plot.xAxis = "nm"
# Ignore channel bins corresponding to 15.0 to 20.0 nm wavelengths:
s1.ignore("15.-20.")

The currently noticed channel ranges are displayed for each spectrum in the AllData.show() output. You can also get a list of the individual noticed channel numbers from Spectrum's noticed attribute:

>>> s1.noticed
[3,4,5,7,8,10]

To apply ignore and notice commands to ALL loaded spectra, call the methods from the global AllData object. To apply to a subset of loaded spectra, add a range specifier to the left of the colon:

# These apply to all loaded spectra
AllData.ignore("100-120, 150-200")
AllData.notice("all")
AllData.ignore("bad")
# These apply to a subset of loaded spectra
AllData.ignore("1-3: 60-65")
AllData.notice("2-**:50-60")

Models

Model With Multiple Data Groups

When a model is defined and spectra are assigned to multiple data groups, PyXspec will generate a Model object copy for each data group (assuming the spectra also have responses attached). So if:

m1 = Model("phabs*ga")
AllData("file1 2:2 file2")

then there are 2 Model objects for the model definition phabs*gaussian. The variable m1 is set to the object belonging to data group 1, and to get the object for data group 2 do:

m2 = AllModels(2)

m1 and m2 will each have the same set of Component and Parameter objects.

Parameters can be accessed directly by index from the Model objects, and these indices are numbered from 1 to nParameters for ALL data group copies. So for the phabs*ga example above:

p = m1(2) # Returns the 2nd ('LineE') parameter from the model for data group 1.
p = m2(2) # Returns the 2nd ('LineE') parameter from the model for data group 2.
p = m2(6) # Wrong!!

Defining Multiple Models

Beginning with XSPEC12, it became possible to assign multiple sources to spectra, and each source may have its own model function definition. To keep track of multiple model definitions, XSPEC requires that you assign them names. In PyXspec, the model name and source number are supplied as additional arguments to the Model __init__ function:

# Define a model named "alpha" assigned to source 1
m_1_1 = Model("phabs*po","alpha")
# Define a model named "beta" assigned to source 2
m_2_1 = Model("const*bbody","beta",2)
# (In both of these cases, the returned object belongs to data group 1)

Note

As with Standard XSPEC, to define a model for source numbers > 1 you first must load a detector response for the source. See "Background, Response, and Arf" in the previous section.

Note that in all previous examples in this tutorial, we have been using unnamed models which were assigned to source 1. Named models and source numbers may also be defined directly into the AllModels container by passing in a tuple:

# Define a model named "defn1" assigned to source 1
AllModels += ("phabs*po", "defn1")
# Define a model named "defn2" assigned to source 2
AllModels += ("const*bbody", "defn2", 2)
# This replaces "defn1" with an unnamed model for source 1
AllModels += "phabs*gaussian"

and from which Model objects can be retrieved:

# Get the "defn2" Model object for data group 1
m_2_1 = AllModels(1,"defn2")
# ...and for data group 2
m_2_2 = AllModels(2,"defn2")

To view all current source number and model assignments, see the AllModels.sources attribute, which displays a dictionary of the [source number]:[model name] pairs.

To remove model definitions:

# Remove all data group copies of "defn2"
AllModels -= "defn2"
# Remove all data group copies of the unnamed model (defined above
#    as "phabs*gaussian")
AllModels -= ""
# Remove all copies of ALL model definitions
AllModels.clear()

Component And Parameter Access Part 2

When PyXspec constructs a Model object, it immediately adds to it an attribute of type Component for every component in the model expression. The attribute has the same (full) name as the component in the original expression, allowing you to access it like:

m = Model("phabs*pow")
c2 = m.powerlaw

However when a model contains multiple copies of the same component, this type of access becomes ambiguous. So to distinguish among copies, for any component making its 2nd or more appearance (from left to right), PyXspec will append _n to the attribute name where n refers to the component's position in the expression (again from left to right). Or to put it more simply:

m = Model("phabs*po + po")
# This gets the leftmost powerlaw component
pow1 = m.powerlaw
# This gets the rightmost, which is the 3rd component in the expression.
pow2 = m.powerlaw_3

The Model object also stores an attribute which is a just a list of the names of its constituent Component attributes:

>>> m.componentNames
['phabs', 'powerlaw', 'powerlaw_3']

This may be useful for example if writing a loop to access each of a model's components. Similarly Component objects have a parameterNames attribute, listing the names of their constituent Parameter attributes:

>>> m.powerlaw.parameterNames
['PhoIndex', 'norm']

Gain Parameters (Response Models)

Response Models differ from the regular kind in that they act on a Response rather than directly calculate a flux. At present there is only one kind of Response Model in Xspec, and this is gain. gain is a built-in attribute of all Response objects, and is of the class type RModel. It has 2 parameters for adjusting the energies of a Response: slope and offset. Gain parameters are initially off by default, but may be turned on simply by setting either one. For example:

s = Spectrum("file1")
# The spectrum's response has a gain attribute that is not in use,
# which is the equivalent of having a slope fixed at 1.0 and offset = 0.0.
r = s.response
# Setting either the slope or offset turns the gain on for this response.
# Both slope and offset will now be fit parameters.
r.gain.slope = 1.05
# The previous setting leaves the offset at 0.0.  Now we'll change it.
r.gain.offset = .05
# You can set slope and offset at the same time using Response's setPars method.
r.setPars(.99, .03)

slope and offset are Parameter objects and therefore have the same interface as regular model parameters:

# Modify the parameter's auxilliary values
r.gain.offset = ".08,,.01,.01,.5,.5"
# Set a parameter link
r.gain.offset.link = ".005*1"

To remove the response fit parameters and return the Response back to its original state, call the gain.off() method:

# This deletes the slope and offset parameters.
# Any references to them become invalid.
r.gain.off()

Flux Calculations

To perform a Standard XSPEC flux or lumin calculation, call the AllModels methods calcFlux or calcLumin respectively:

AllModels.calcFlux(".3 1.0")
AllModels.calcFlux(".1 10.0 err")
AllModels.calcLumin(".1 10. .05 err")

As in Standard XSPEC the results will be stored with the currently loaded spectra:

>>> s1 = AllData(1)
>>> s1.flux
(5.7141821510911499e-14, 0.0, 0.0, 4.0744161672429196e-05, 0.0, 0.0)
>>> s1.lumin
(30.972086553634504, 0.0, 0.0, 0.056670019567301853, 0.0, 0.0)

unless there are no spectra, in which case the results are stored with the model object:

>>> AllModels(1).flux
(5.6336924399373855e-10, 0.0, 0.0, 0.05929616315253175, 0.0, 0.0)

Local Models in C/C++/Fortran

In Standard XSPEC, local model libraries are built with the initpackage comamnd and then loaded with lmod. The AllModels container supplies both of these functions for doing the same thing in PyXspec:

AllModels.initpackage("myLocalMods","lmodel.dat")
AllModels.lmod("myLocalMods")

By default this looks in the directory set by the LOCAL_MODEL_DIRECTORY variable in your ~/.xspec/Xspec.init start-up file. You can override this by giving these functions an absolute or relative path as a dirPath keyword argument (see ModelManager for details).

Local Models in Python

You can also write model functions in Python and insert them into the XSPEC library with the AllModels addPyMod method. You simply define a function with 3 arguments for energies, parameters, and flux. For example a powerlaw model function might look like:

def lpow(engs, params, flux):
   for i in range(len(engs)-1):
      pconst = 1.0 - params[0]
      val = math.pow(engs[i+1],pconst)/pconst - math.pow(engs[i],pconst)/pconst
      flux[i] = val

XSPEC will pass tuples containing the energy and parameter values to your function. For the flux array, it will pass a list pre-sized to nE-1, where nE is the size of the energies array. Your model function should fill in this list with the proper flux values. (For additional optional arguments to your model function, please see the documentation for the addPyMod method in ModelManager.)

The second thing you must define is a tuple containing the parameters' information strings, one string for each parameter in your model. This is equivalent to the parameter strings you would define in a model.dat file when adding local models in standard XSPEC, and it requires the same format. (See Appendix C of the XSPEC manual for more details.) So with the powerlaw function above which takes just 1 parameter, you might define a tuple as:

powInfo = ("phoIndex   \"\"  1.1  -3.  -2.  9.  10.  0.01",)

Note the need for the trailing comma when there's just one parameter string. This is to let Python know that powInfo is a tuple type and not a string.

Once you've defined your function and parameter information, simply call:

AllModels.addPyMod(lpow, powInfo, 'add')

The 3rd argument tells XSPEC the type of your model function ('add', 'mul', or 'con'). After this call your function will be added to the list of available model components, which you can see by doing Model.showList(). Your model will show up with the same name as your original Python function ('lpow'), and is ready for use in future model definitions.

Fitting

Error

The error command is implemented through Fit (see FitManager), and the results are stored with the chosen Parameter object(s). The error attribute stores a tuple containing the low and high range values for the parameter, and the 9-letter status string to report problems incurred during the error calculation.

Example: Estimate the 90% confidence range for the 4th parameter

>>> Fit.error("2.706 4")
>>> par4 = AllModels(1)(4)
>>> par4.error
(0.11350354517707145, 0.14372981075906774, 'FFFFFFFFF')

Query

During a Fit.perform() operation, the default is to query the user whenever the fit has run the maximum number of iterations, as set by the Fit.nIterations attribute. You can change this behavior with the query attribute:

# When nIterations is reached, continue the fit without stopping to query.
Fit.query = "yes"
# Stop fit at nIterations and do not query.
Fit.query = "no"
# Query the user when nIterations is reached.
Fit.query = "on"

Steppar

The Standard XSPEC steppar command is also implemented through the global Fit object. You supply it with a string following the same steppar command syntax rules. For example:

# Step parameters 1 and 2 through the given range values
# over a 10x10 2-D grid.
Fit.steppar("1 20. 30. 10 2 .05 .08 10")

Fakeit

PyXspec provides access to standard XSPEC's fakeit command, which is for creating spectra with simulated data. It is called through the AllData.fakeit method:

AllData.fakeit(nSpectra=1, settings=None, applyStats=True, filePrefix="")

Note

If AllData.fakeit is run when spectra are currently loaded, it will follow the same rules as the standard XSPEC fakeit function: It will REMOVE ALL pre-existing spectra and replace each one with a simulated spectrum (even if nSpectra is less than the number originally loaded).

As those familiar with standard fakeit know, the user is normally prompted for quite a bit of additional information needed to generate the fakeit files. However the goal here is to have NO additional prompting, and that requires that all information must be entered as arguments to the AllData fakeit method call. This is done by passing objects of the FakeitSettings class to AllData.fakeit, as we'll show further below.

Note

Unless stated otherwise, assume all spectra are OGIP type-1 (1 spectrum per file).

For the simplest of cases, you don't need to create any FakeitSettings objects. Just pass in the number of fake spectra you'd like to create:

# Create 3 fake spectra using only default settings.
AllData.fakeit(3)

The fakeit function will then create a default FakeitSettings object for each of the 3 spectra. By default, a FakeitSettings object will have empty strings for all of its attributes, and these are handled differently depending on whether the fake spectrum is replacing a currently loaded spectrum, or creating one from scratch.

From Existing Spectra

When replacing an existing spectrum, FakeitSettings attributes with empty strings will simply take their value from the original spectrum. Also note that the response and arf settings for the original spectrum CANNOT be modified for the fakeit spectrum. If a name is filled in for either of these attributes, it will be ignored. If you wish to modify these, you can make the change to the original spectrum prior to calling fakeit. [The one exception is when the original spectrum has no response, in which case the response attribute MUST be filled in.] If the fileName attribute is empty, XSPEC will generate a default output name derived from the original file name.

From Scratch

When creating from scratch, an empty string implies "none" for the arf and background, 1.0 for exposure and correction, and XSPEC's default dummy response for the response attribute. If the fileName attribute is empty, XSPEC will generate a default output file name based on the response name, and it will include an auto-incremented index to prevent multiple output files from overwriting each other.

FakeitSettings Objects

To create a fake spectrum with anything other than default settings, you must supply a FakeitSettings object for that spectrum. The FakeitSettings attributes are: response, arf, background, exposure, correction, backExposure, and fileName. All are string types, though exposure, backExposure, and correction can also be entered as floats. Attributes can be set upon object construction, or anytime afterwards:

fs1 = FakeitSettings("response1.rsp", exposure = 1500.0)
fs1.background = "back1.pha"

A new FakeitSettings object can also be made by copying an existing one:

fs2 = FakeitSettings(fs1)

And now pass the objects to the fakeit method, either in a list, dictionary, or as a single object:

# Apply settings to fakeit spectra 1 and 2:
AllData.fakeit(2,[fs1,fs2])
# Apply setting to fakeit spectrum 1, use defaults for spectrum 2:
AllData.fakeit(2, fs1)
# Apply settings to fakeit spectra 2 and 4, use defaults for 1 and 3:
settingsDict = {2:fs1, 4:fs2}
AllData.fakeit(4, settingsDict)
# Create 4 fakeit spectra from the same settings object:
settingsList = 4*[fs1]
AllData.fakeit(4, settingsList)

The remaining 2 arguments to the AllData.fakeit function are for choosing whether to apply statistical fluctuations (default = True), and whether to add an optional prefix string to the names of all output files.

OGIP Type-2 Files

With OGIP type-2 files, multiple spectra may be placed in a single file. The important thing to recognize when generating type-2 fakeit files is that the exposure, correction, backExposure, and fileName attributes apply to the output files and not the individual spectra. Therefore these settings will be ignored for all but the first spectrum in a file. For example:

# Start with 4 spectra loaded, in 2 type-2 files:
AllData("myDataFile1.pha{1-2} myDataFile2.pha{7-8}")
# Create settings for the 4 fake spectra that will be generated from these:
fs1 = FakeitSettings(background="back1.pha", exposure=250.)
# The exposure setting in fs2 will be ignored!!!
fs2 = FakeitSettings(background="back2.pha", exposure=100.)
fs3 = FakeitSettings(fileName="myFakeitFile_2.pha")
fs4 = FakeitSettings(fs3)
# The following change will be ignored!!!
fs4.fileName = "myFakeitFile_3.pha"
# Now generate the fakeit files:
AllData.fakeit(4, [fs1,fs2,fs3,fs4])

The above will generate 4 fakeit spectra, placed in 2 type-2 files. The exposure setting for spectrum 2 and the fileName setting for spectrum 4 will be ignored. Those values are only set by spectra 1 and 3.

For more fakeit details and examples, please check FakeitSettings and the fakeit method description in DataManager.

Monte Carlo Markov Chains (MCMC)

All MCMC operations are handled either by objects of class Chain, or the global AllChains container object (see the ChainManager class description). To create a new chain based on the current fit parameters, simply create a Chain object by passing it an output file name:

c1 = Chain("chain1.fits")

The above call creates the file chain1.fits, performs an MCMC run using the default burn, fileType, length, proposal, rand, and temperature values, and automatically places the new object in the AllChains container. These default settings are stored as attributes of AllChains:

# Ensure that new chains will burn the first 100 iterations, will
# have length 1000, and will use the proposal "gaussian fit"
AllChains.defBurn = 100
AllChains.defLength = 1000
AllChains.defProposal = "gaussian fit"
c2 = Chain("chain2.fits")

You can also override the AllChains default settings by passing additional arguments to Chain upon construction:

# Length will be 2000 for this chain, use defaults for all other settings.
c3 = Chain("chain3.fits", runLength = 2000)

The new chain objects will then store their own settings as attributes:

>>> c2.burn
100
>>> c2.runLength
1000
>>> c3.runLength
2000

All of a chain object's attributes will be displayed when calling its show() method.

To append a new run to an existing chain object, call the object's run() method. The appending run will use the object's current attribute settings, and not the AllChains default settings:

# This will append a run of length 3000 to the c3 chain object, and with a
# Metropolis-Hastings temperature of 50.0:
c3.runLength = 3000
c3.temperature = 50.0
c3.run()

Display the new chain length:

>>> c3.totalLength
5000

To overwrite rather than append to an existing chain object, call run with its append argument set to False:

# This erases the results of any previous runs for object c3.
c3.run(False)

>>> c3.totalLength
3000

New chains are loaded into AllChains by default, but you can unload or reload them using the AllChains arithmetic operators:

# Chain c2 may be unloaded by passing its chain index number
AllChains -= 2
# OR by passing the object itself
AllChains -= c2
# 2 ways to remove ALL chains
AllChains -= '*'
AllChains.clear()

# Reload an existing chain object
AllChains += c2
# Load a chain from an existing chain file
AllChains += "earlierChain.fits"
# Create a new chain, to be stored in file "chain4.fits"
AllChains += "chain4.fits"

As with Standard XSPEC, unloading a chain will leave the chain's file intact. It merely removes the chain from XSPEC's calculations. To display information about the currently loaded chains, call AllChains.show().

You may also get a chain object from the container at any time by passing it an index number:

# Retrieve a chain object for the 4th chain in the container
c4 = AllChains(4)

Plotting

All of the plotting options available in Standard XSPEC's setplot command are now implemented as attributes of the Plot object. Some of these are mentioned in the Quick Tutorial, and please see the PlotManager class reference for the complete guide.

One setting of particular interest is the commands attribute. This is a tuple of user-entered PLT command strings which are added to XSPEC's auto-generated commands when performing a plot, and is modified through Plot's addCommand and delCommand methods. For example, to enter a PLT command to place an additional label at the specified coordinates on the plot:

Plot.addCommand("label 1 pos 10 .05 \"Another Label\"")

To view the currently loaded commands:

print(Plot.commands)

and to remove the 3rd command from the tuple:

Plot.delCommand(3)

XSPEC Settings

Most of the internal switches set through Standard XSPEC's xset command are now set through attributes of the global Xset object. (See XspecSettings for the full class description.)

Examples:

Xset.abund = "angr"
Xset.cosmo = "50 .5 0."
Xset.xsect = "bcmc"

Xset also provides the methods addModelString and delModelString to set the <string name>,<string value> pairs which are used by certain models. The <string name> argument is case-insensitive:

Xset.addModelString("APECROOT","1.3.1")
Xset.addModelString("APECTHERMAL","yes")
Xset.delModelString("neivers")

The entire collection of <name>,<value> pairs may be set or retrieved with the Xset.modelStrings attribute:

# Replace all previous entries with a new dictionary
Xset.modelStrings = {"neivers":"1.1", "apecroot":"1.3.1"}
# Clear out all entries:
Xset.modelStrings = {}

Xset.show() will display all of the current settings including the current <string name>,<string value> pairs.

Logging And XSPEC Output

The Xset object provides attributes and methods for controlling output chatter level and for creating log files:

# Get/Set the console chatter level
ch = Xset.chatter
Xset.chatter = 10
# Get/Set the log chatter level
lch = Xset.logChatter
Xset.logChatter = 20

# Create and open a log file for XSPEC output
# This returns a Python file object
logFile = Xset.openLog("newLogFile.txt")
# Get the Python file object for the currently opened log
logFile = Xset.log
# Close XSPEC's currently opened log file.
Xset.closeLog()

Exceptions And Error Handling

PyXspec utilizes the standard Python try/except/raise mechanism for handling and reporting errors. Currently only exception objects of the class Exception are ever raised. In the future other (more specific) error classes may be used, but they should always be derived from Exception. So you can catch all PyXspec exceptions with code such as:

try:
    # Only 4 spectra are currently loaded
     s = xspec.AllData(5)
except Exception as msg:
     print(str(msg))

which will print the error message

Error: Spectrum index number is out of range: 5

PyXspec raises errors in a variety of situations, such as for invalid input argument syntax, or for input which is invalid within the context of the call (as in the example above). It can also raise exceptions if you try to rebind a class attribute when such modification is not permitted.

Adding Attributes To PyXspec Objects

A particularly novel feature of Python (in comparison with say C++) is that it allows you to create new attributes "on the fly". The attributes don't have to have been part of the original class definition:

class C:
    pass

x = C()
x.pi = 3.1416

The downside of course is that spelling or case sensitive errors become much harder to detect. For example, with PyXspec's Plot object:

Plot.yLog = True  # Correct
Plot.ylog = True  # Wrong!

In the second case, standard Python will simply add a new attribute named ylog to Plot, and this will have no effect on the actual plot since PyXspec is only looking for an attribute named yLog.

So operating under the assumption that this downside outweighs the benefits, we've decided to disable the ability to add new attributes to PyXspec class objects. A misspelling or case error will instead raise an Exception object. And since some users may genuinely wish to add their own attributes to PyXspec classes, this default behavior may be overridden by toggling the Xset.allowNewAttributes flag:

s = Spectrum("dataFile.pha")
s.myNewIndex = 10  # Error: Will raise an exception
Xset.allowNewAttributes = True

s.myNewIndex = 10 # OK
#  .
#  .
#  .  In this section you can now add new attributes to any PyXspec object,
#  .  but attribute spelling errors will go undetected.
#  .
#  .

Xset.allowNewAttributes = False

Using With Other Packages

One of the primary benefits of PyXspec is that it makes it much easier to use XSPEC data and results in 3rd party packages. For example you can bypass XSPEC's built-in plotting functions in favor of a Python plotting library such as Matplotlib:

#!/usr/bin/python

from xspec import *
import matplotlib.pyplot as plt

# PyXspec operations:
s = Spectrum("file1.pha")
m = Model("phabs*po")
Fit.perform()
Plot.device = '/null'
Plot('data')
# Get coordinates from plot:
chans = Plot.x()
rates = Plot.y()
folded = Plot.model()

# Plot using Matplotlib:
plt.plot(chans, rates, 'ro', chans, folded)
plt.xlabel('channels')
plt.ylabel('counts/cm^2/sec/chan')
plt.savefig('myplot')

The above code produces a Matplotlib plot of the spectral data and folded model vs. channels (similar to what you get with Standard XSPEC's "plot data" command). It makes use of the Plot object's x and y methods to return the coordinates of the spectrum plot, and the model method to get the folded model values. These are then forwarded to Matplotlib's plotting command.