SPEX is a spectral fitting program used to fit high-resolution X-ray spectra and works with a command-line interface. A first-time user should get familiar with the syntax of the commands to be able to work with it. This page provides some basic commands and threads to fit X-ray spectra.
Since SPEX uses its own spectrum and response file format, you may need to convert your OGIP spectra and responses to SPEX format first. This can be done with the trafo program in the SPEX package.
Before you start, make sure that the SPEX environment is properly set, e.g. the SPEX90 environment variable is set and the spexdist.(c)sh script has been sourced:
user@linux:~> echo $SPEX90 /opt/spex
In case you already have spectra in SPEX format, you can skip the trafo step and move on to the SPEX section immediately.
Before trafo is started, you need the OGIP spectra and responses first. Please read the documentation regarding your dataset for more information on how to create OGIP files. The minimum that trafo needs is a spectrum (.pi or .pha) and a response matrix (.rmf or .rsp). If you want to subtract background spectra or if you have an additional arf file, then please also collect these files in your working directory.
NOTE: There is also an (experimental) alternative for the trafo program called ogip2spex. This program is part of the SPEX Python tools (PYSPEX). This program takes the necessary input from the command line and is easy for scripting. More information about ogip2spex can be found on the Pyspex Github documentation page.
It is best to run trafo in the directory where your OGIP files are. It is an interactive program, so it will ask the user for information when the program is run.
user@linux:~> trafo Program trafo: transform data to SPEX 2.0 format This is version 1.04, of trafo Are your data in OGIP format : type=1 Old (Version 1.10 and below) SPEX format: type=2 New (Version 2.00 and above) SPEX format: type=3 Enter the type:
The first questions are quite straightforward. In the case of OGIP spectra, the type is always 1. In principle, it is possible to put more than one spectrum in a spo and res file, but for most simple cases transforming a single spectrum is sufficient.
Enter the type: 1 Enter the number of spectra you want to transform: 1 Enter the maximum number of response groups per energy per spectrum: 100000
The maximum number of response groups should just be a large number in nearly all cases. In special cases, it can put a limit on the number of response groups that will be saved in the .res file.
Optimizing the response matrix
The following feature is present in trafo since version 1.02 (SPEX version 2.02). It allows the user to re-arrange the response matrix to make more optimal use of parallel processing. There are three options: 1. Keep the matrix as provided. 2. Try to re-arrange the matrix into contiguous groups. The program tries to identify physically distinctive components and avoids overlapping data. 3. Split the matrix into N equal-sized components. This is particularly useful for grating spectra (RGS) and allows for efficient matrix multiplication on multi-core processors. Any power of 2 between 8 and 32 should provide a fast response matrix. In the terminal, this option is provided in the following way:
How should the matrix be partioned? Option 1: keep as provided (1 component, no re-arrangements) Option 2: rearrange into contiguous groups Option 3: split into N roughly equal-sized components Enter your preferred option (1,2,3): 1
Option number 1 is the safest option to choose, but also the slowest. Option 2 and 3 can provide a significant increase in performance, but results should be carefully checked. More information about re-arranging response matrices can be found in the SPEX Manual.
Reading the spectra
Then, trafo asks for the filenames of the source and background spectra. First, provide the file name of the source spectrum. trafo will return some of the basic properties of the spectral file, like exposure time and values of the most important FITS keywords.
Enter filename spectrum to be read: PN-source.pi Exposure time (s): 2.10571992E+04 Assuming Poissonian Errors Areascal: 1.00000000E+00 Backscal: 1.00000000E+00 No BACKFILE keyword found Corrscal: 1.00000000E+00 No CORRFILE keyword found No RESPFILE keyword found No ANCRFILE keyword found No background specified in pha-file.
A background spectrum can be provided (optional), which will be subtracted from the source spectrum. If a background file is already specified in the FITS header of the source spectrum, this question will not be asked.
Read nevertheless a background file? (y/n) [no]: y Enter filename background spectrum to be read: PN-background.pi Exposure time (s): 2.10572832E+04 Assuming Poissonian Errors Areascal: 1.00000000E+00 Backscal: 1.00000000E+00 No BACKFILE keyword found Corrscal: 1.00000000E+00 No CORRFILE keyword found No RESPFILE keyword found No ANCRFILE keyword found
Bad channels and grouping
Depending on the instrument used, there is a chance that the spectrum contains bad channels. This is especially true for grating spectra. Sometimes the background spectrum can have a different number of bad channels than the source spectrum. It is therefore important that a particular bad channel in either of the two spectra is ignored. In this example, there are no bad channels, so either yes or no will do.
Checking data quality and grouping ... Ogip files have quality flags. Quality 0 means okay Your spectrum file has 0 bins with bad quality Your background file has 0 bins with bad quality Your combined file has 0 bins with bad quality Shall we ignore bad channels? (y/n) [no]:y
If grppha has been used on the spectrum, trafo will also ask whether the spectra should be binned according to the groups defined in the PHA file.
Important note: We do not recommend the use of grppha for binning spectra. For spectra with Poisson statistics (most X-ray spectra), it is much better to use C-statistics and use an optimal binning algorithm in SPEX based on the spectral resolution of the instrument.
Read response and effective area files
In the next step, the response matrix is read. Sometimes, the response matrices start at channel 0, which can be somewhat confusing. Especially when some arrays start at channel 0 and others at channel 1. If both data sets start at zero, it is best to shift the channel numbers with 1 unit. For most instruments this is fine, however, there are situations when this does not apply. In that case, please check your energy grid by loading a delta line component in SPEX and check the energy of the line manually. Then, compare the output with a delta line defined in XSPEC.
Determining background subtracted spectra ... No response matrix file specified in pha-file. Enter filename response matrix to be read: PN.rmf Reading response matrix ... Warning, ebounds data started at channel 0 Warning, response data started at channel 0 Possible response conflict; check xspec/spex with delta line! Enter shift to response array (1 recommended, but some cases may be 0):1 No effective area file specified in pha-file.
Sometimes, also an effective area file needs to be provided separately:
Read nevertheless an effective area file? (y/n) [no]: y Enter filename arf-file to be read: PN.arf Reading effective area ... Determining zero response data ... Total number of channels with zero response: 373 Original number of data channels : 4096 Channels after passing mask and omitting zero response channels: 3723 Rebinning data where necessary ... Rebinning response where necessary ... old number of response elements: 435950 new number of response elements: 435950 old number of response groups : 1481 new number of response groups : 1481 Correcting for effective area ... Determine number of components ... Found 1 components Enter any shift in bins (0 recommended): 0 order will not be swapped ...
If there are bins with zero response, then they are excluded from the resulting file. Also here a shift in bins can be set, but the recommended value is 0.
Writing res and spo files
The final step is writing the spectra in SPEX format. The file names should be provided without an extension. The .spo and .res extension will be added automatically.
Enter filename spectrum to be saved (without .spo): PN Enter filename response to be saved (without .res): PN Final number of response elements: 435950
The PN.spo and PN.res file have been saved in the current directory.
When you have a .res and .spo file, it is time to try SPEX. SPEX can read in the spectra and responses from the files and fit the spectrum with a model. In the example below, we show how a simple fit is done. If you want to reproduce it, you can download PN.spo and PN.res.
The SPEX program is started by entering spex in a linux terminal window. In the following sections, we describe one run of the program. To start SPEX do this:
user@linux:~> spex Welcome user to SPEX version 3.05.00 SPEX>
First, we have to load the SPEX data files. This is done using the command data. It is a general thing in SPEX that filename extensions are not typed explicitly when issuing a command. If you have a file called PN.spo and PN.res then you type:
SPEX> data PN PN
The responsefile (.res) is entered first and then the file containing the spectrum (.spo). You can avoid confusion by giving the same filename to both .res and .spo files. Remember that the order of the words in the commands is very important!
Plotting the data
If the data command was successful, we can now have a look at the spectra. SPEX offers a lot of different plot commands through the PGPLOT package (known from, for example, Fortran programs). Below is an example of a log-log plot for XMM-Newton PN data:
SPEX> plot dev xs SPEX> plot type data SPEX> plot x log SPEX> plot y log SPEX> plot rx 0.2:10. SPEX> plot ry 1E-4:10. SPEX> plot
The sequence above opens a PGPLOT window (dev xs) and tells SPEX that the plot will contain data with logarithmic axes (plot x log and plot y log). Then we set the ranges of the X and Y axis with plot rx 0.2:10 and plot ry 1E-4:10., respectively. The final plot command puts the plot on screen:
Ignoring and rebinning
High-resolution X-ray spectra from Chandra and XMM-Newton are usually oversampled (e.g. the energy bins are much smaller than the spectral resolution) and contain a lot more channels then is useful. It is necessary to remove energy intervals that contain bad data and to rebin your spectrum. In the next example, we exclude energies below 0.2 keV and above 10 keV. Then we optimally bin the spectrum with the obin command, which determines the binning based on the spectral resolution detected in the response.
SPEX> ignore 0.0:0.2 unit kev SPEX> ignore 10.:100. unit kev SPEX> obin 0.2:10. un kev SPEX> plot
The final plot command shows the result of the obin command. The spectrum appears to be much less noisy then before:
Defining a model
Now we have a clean and rebinned spectrum that is ready to fit. Before we can start fitting, we first need to define a model. In this case, we take a simple absorbed power law. The absorptionmodel is absm and the power-law model is called pow. They can be loaded with the com command that stands for component:
SPEX> com abs You have defined 1 component. SPEX> com pow You have defined 2 components. SPEX> com rel 2 1
The last command (com rel 2 1) tells SPEX that component 2, the powerlaw, is absorbed by component 1.
Fitting the data
Now we have to estimate the initial parameters. This takes some experience, or trial and error, to figure out. For this spectrum, we assume a hydrogen column density (NH) value of 1E+20 cm-2. We set the power-law normalisation to 1000 for now. If you plot the spectrum again, is the red model line visible? If yes, then the initial guess is good enough.
SPEX> par 1 1 nh value 1E-4 SPEX> par 1 2 norm value 1E+3 SPEX> calc SPEX> plot
We are ready to fit the data! SPEX has a nice feature to look at the progress of the fit. To activate this feature you have to give the command fit print 1. If your initial parameters were acceptable, you can see the model converge to the data in the plot window after you entered the fit command. When the fit is done, then the parameters and C-stat are printed on screen. If the C-stat value is close to the number of degrees of freedom, then your fit is acceptable. Sometimes more runs of the command fit are necessary after changing some initial parameters. This is especially true when using complex models. Again this is a game of trial and error.
You also might want to fix or free certain parameters to see if they can be constrained. In SPEX fixing is f (frozen) and freeing is t (thawn). For example, the NH parameter can be frozen with the command: par 1 1 nh stat f.
The commands below start the fitting process:
SPEX> fit print 1 SPEX> fit -------------------------------------------------------------------------------------------------- sect comp mod acro parameter with unit value status minimum maximum lsec lcom lpar 1 1 absm nh Column (1E28/m**2) 1.9518682E-04 thawn 0.0 1.00E+20 1 1 absm fcov Covering fraction 1.000000 frozen 0.0 1.0 1 2 pow norm Norm (1E44 ph/s/keV) 99.62738 thawn 0.0 1.00E+20 1 2 pow gamm Photon index 2.000923 thawn -10. 10. 1 2 pow dgam Photon index break 0.000000 frozen -10. 10. 1 2 pow e0 Break energy (keV) 1.0000000E+10 frozen 0.0 1.00E+20 1 2 pow b Break strength 0.000000 frozen 0.0 10. 1 2 pow type Type of norm 0.000000 frozen 0.0 1.0 1 2 pow elow Low flux limit (keV) 2.000000 frozen 1.00E-20 1.00E+10 1 2 pow eupp Upp flux limit (keV) 10.00000 frozen 1.00E-20 1.00E+10 1 2 pow lum Luminosity (1E30 W) 2.565488 frozen 0.0 1.00E+20 Instrument 1 region 1 has norm 1.00000E+00 and is frozen --------------------------------------------------------------------------------------------- sect comp mod acro parameter with unit sect comp mod acro parameter with unit correl. 1 1 absm nh Column (1E28/m**2) <-> 1 2 pow norm Norm (1E44 ph/s/keV) 0.786 1 1 absm nh Column (1E28/m**2) <-> 1 2 pow gamm Photon index 0.789 1 2 pow norm Norm (1E44 ph/s/keV) <-> 1 2 pow gamm Photon index 0.557 -------------------------------------------------------------------------------- Fluxes and restframe luminosities between 2.0000 and 10.000 keV sect comp mod photon flux energy flux nr of photons luminosity (phot/m**2/s) (W/m**2) (photons/s) (W) 1 2 pow 3.15846 2.037560E-15 3.980407E+45 2.565460E+30 Fit method : Classical Levenberg-Marquardt Fit statistic : C-statistic C-statistic : 245.96 Expected C-stat : 226.73 +/- 21.34 Chi-squared value : 255.59 Degrees of freedom: 223 W-statistic : 243.47
With the fit print 1 command, the plot is automatically update and now shows a model that fits the data points:
When the fit is acceptable, you might want to know the uncertainties on your fitted parameters. Errors are determined one-by-one by fixing the parameter to some value and calculate the values where delta C-stat = 1 with respect to the best fit. This corresponds to the 1 sigma error (68%). The error command is simply:
SPEX> error 1 1 nh parameter C-stat Delta Delta value value parameter C-stat ---------------------------------------------------- 1.882095E-04 247.02 -6.901551E-06 1.06 1.886280E-04 246.90 -6.483009E-06 0.94 1.884233E-04 246.96 -6.687726E-06 1.00 2.020126E-04 247.03 6.901551E-06 1.07 2.015895E-04 246.90 6.478454E-06 0.94 2.017922E-04 246.96 6.681148E-06 1.00 Parameter 1 1 nh : 1.95111E-04 Errors: -6.68773E-06 , 6.68115E-06
In this case, the best fit value and error of the NH is 1.95 +/- 0.07 1020 cm-2.
End the program
SPEX can be stopped by typing quit on the SPEX prompt:
Note: The example above is meant to provide a simple quick fit to a simple dataset. In reality, spectral analysis is much more complicated. See the SPEX Cookbook for more examples.