next up [*] [*]
Next: Simultaneous Fitting: Examples from Up: Walks through XSPEC Previous: Brief Discussion of XSPEC

Fitting Models to Data: An Example from EXOSAT

The 6-s X-ray pulsar 1E1048.1-5937 was observed by EXOSAT in June 1985 for 20 ks. In this example, we'll conduct a general investigation of the spectrum from the Medium Energy (ME) instrument, i.e. follow the same sort of steps as the original investigators (Seward, Charles & Smale,1986). The ME spectrum and corresponding response matrix were obtained from the HEASARC On-line service.

Once installed, XSPEC is invoked by typing xspec, as in this example.

%xspec

               Xspec 11.1.0    21:50:14 08-Jul-2001

 For documentation, notes, and fixes see  http://xspec.gsfc.nasa.gov/

 Plot device not set, use "cpd" to set it

 Type "help" or "?" for further information


XSPEC>data s54405.pha
 Net count rate (cts/s) for file   1   3.783    +/-  0.1367
   using response (RMF) file...       s54405.rsp
   1 data set is in use

The data command tells the program to read the data as well as the response file that is named in the header of the data file. In general, XSPEC commands can be truncated, provided they remain unambiguous. Since the default extension of a data file is .pha, and the abbreviation of data to the first two letters is unambiguous, the above command can be abbreviated to da s54405, if desired. To obtain help on the data command, or on any other command, type help command followed by the name of the command.

One of the first things most users will want to do at this stage--even before fitting models--is to look at their data. The plotting device should be set first, with the command cpd (change plotting device). Here, we use the pgplot X-Window server, /xw.

XSPEC> cpd /xw

To see a list of alternative plot devices, type cpd ? There are more than 20 different things that can be plotted, all related in some way to the data, the model, the fit and the instrument. To see them, type:

XSPEC> plot ? 
 Choose from the following `plot' sub-commands:
 data      counts    ldata     residuals chisq     delchi    ratio       
 summary   model     emodel    eemodel   contour   efficien  ufspec      
 eufspec   eeufspec  dem       insensitv sensitvty genetic   foldmodel   
 icounts
Insert selection: (data)

The most fundamental is the data plotted against instrument channel (data); next most fundamental, and more informative, is the data plotted against channel energy. To do this plot, use the XSPEC command setplot energy. Figure A shows the result of the commands:

XSPEC> setplot energy
XSPEC> plot data


\begin{figure}\begin{center}
\mbox{
{\psfig{figure=figures/figa.eps,width=4.0in}...
...8.1--5937 available from the HEASARC on-line service.}}
\end{center}\end{figure}

People familiar with EXOSAT ME data will recognize the spectrum to be soft, absorbed and without an obvious bright iron emission line.

We are now ready to fit the data with a model. Models in XSPEC are specified using the model command, followed by an algebraic expression of a combination of model components. There are two basic kinds of model components: additive which represent X-Ray sources of different kinds. After being convolved with the instrument response, the components prescribe the number of counts per energy bin (e.g., a bremsstrahlung continuum); and multiplicative models components, which represent phenomena that modify the observed X-Radiation (e.g. reddening or an absorption edge). They apply an energy-dependent multiplicative factor to the source radiation before the result is convolved with the instrumental response. More generally, XSPEC allows three types of modifying components: convolutions and mixing models in addition to the multiplicative type. Since there must be a source, there must be least one additive component in a model, but there is no restriction on the number of modifying components. To see what components are available, type model ?:

XSPEC>model ?

 Possible additive models are :
 apec     bbody    bbodyrad bexrav   bexriv   bknpower bmc      bremss
 c6mekl   c6pmekl  c6pvmkl  c6vmekl  cemekl   cevmkl   cflow    compbb
 compLS   compST   compTT   cutoffpl disk     diskbb   diskline diskm
 disko    diskpn   equil    gaussian gnei     grad     grbm     laor
 lorentz  meka     mekal    mkcflow  nei      npshock  nteea    pegpwrlw
 pexrav   pexriv   plcabs   powerlaw posm     pshock   raymond  redge
 refsch   sedov    srcut    sresc    step     vapec    vbremss  vequil
 vgnei    vmeka    vmekal   vmcflow  vnei     vnpshock vpshock  vraymond
 vsedov   zbbody   zbremss  zgauss   zpowerlw atable

 Possible multiplicative models are :
 absori   constant cabs     cyclabs  dust     edge     expabs   expfac
 highecut hrefl    notch    pcfabs   phabs    plabs    redden   smedge
 spline   SSS_ice  TBabs    TBgrain  TBvarabs uvred    varabs   vphabs
 wabs     wndabs   xion     zedge    zhighect zpcfabs  zphabs   zTBabs
 zvarabs  zvfeabs  zvphabs  zwabs    zwndabs  mtable   etable

 Possible mixing models are :
 ascac    projct

 Possible convolution models are :
 gsmooth  lsmooth  reflect

 Possible pile-up models are :
 pileup

For information about a specific component, type help Models followed by the name of the component):

XSPEC>help Models raymond

     XSPEC_11.1_commands
          Models
               raymond

An emission spectrum from hot diffuse gas based on the model calculations of
Raymond and Smith including line emission from several elements. This model
interpolates on a grid of spectra for different temperatures. The grid is
logarithmically spaced with 80 temperatures ranging from 0.008 to 80 keV.

 par1   = plasma temperature in keV
 par2   = Metal abundances (He fixed at cosmic) The elements 
          included are C, N, O, Ne, Mg, Si, S, Ar, Ca, Fe, Ni and 
          their relative abundances are set by the abund command.
 par3   = redshift, z
 K      = 10**-14 / (4 pi (D_L/(1+z))**2) Int n_e n_H dV, where D_L is 
          the luminosity distance to the source (cm), n_e is the electron 
          density (cm**-3), and n_H is the hydrogen density (cm**-3)


help>

Given the quality of our data, as shown by the plot, we'll choose an absorbed power law, specified as follows :

XSPEC> model phabs(powerlaw)

Or, abbreviating unambiguously:

XSPEC> mo pha(po)

The user is then prompted for the initial values of the parameters. Entering <return> or / in response to a prompt uses the default values. We could also have set all parameters to their default values by entering /* at the first prompt. As well as the parameter values themselves, users also may specify step sizes and ranges (value, delta, min, bot, top, and max values), but here, we'll enter the defaults:

XSPEC>mo pha(po)
  Model:  phabs[1]( powerlaw[2] )
Input parameter value, delta, min, bot, top, and max values for ...
Current:          1     0.001         0         0     1E+05     1E+06
phabs:nH>/*
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
  Model:  phabs[1]( powerlaw[2] )
  Model Fit Model Component  Parameter  Unit     Value
  par   par comp
    1    1    1   phabs      nH       10^22      1.000     +/-      0.
    2    2    2   powerlaw   PhoIndex            1.000     +/-      0.
    3    3    2   powerlaw   norm                1.000     +/-      0.
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
 Chi-Squared =     4.8645994E+08 using   125 PHA bins.
 Reduced chi-squared =      3987376.     for    122 degrees of freedom
 Null hypothesis probability =    0.

Note that most of the other numerical values in this section are have changed from those produced by earlier versions. This is because the default photoelectric absorption cross-sections have changed since XSPEC V.10. To recover identical results to earlier versions, use xsect obcm. bcmc is the default; see help xsect for details).

The Current statistic is $\chi^2$ and is huge for the initial, default values--mostly because the power law normalization is two orders of magnitude too large. This is easily fixed using the renorm command :

XSPEC> renorm
 Chi-Squared =      852.1660     using   125 PHA bins.
 Reduced chi-squared =      6.984967     for    122 degrees of freedom
 Null hypothesis probability =    0.

We are not quite ready to fit the data (and obtain a better $\chi^2$), because not all of the 125 PHA bins should be included in the fitting: some are below the lower discriminator of the instrument and therefore do not contain valid data; some have imperfect background subtraction at the margins of the pass band; and some may not contain enough counts for $\chi^2$ to be strictly meaningful. To find out which channels to discard (ignore in XSPEC terminology), users must consult mission-specific documentation, which will inform them of discriminator settings, background subtraction problems and other issues. For the mature missions in the HEASARC archives, this information already has been encoded in the headers of the spectral files as a list of ``bad'' channels. Simply issue the command:

XSPEC> ignore bad
 Chi-Squared =      799.7109     using    85 PHA bins.
 Reduced chi-squared =      9.752572     for     82 degrees of freedom
 Null hypothesis probability =    0.
XSPEC> setplot chan
XSPEC> plot data


\begin{figure}\begin{center}
\mbox{
{\psfig{figure=figures/figb.eps,width=4.0in}...
... bad} on the {\sl EXOSAT} ME spectrum
1E1048.1--5937.}}
\end{center}\end{figure}

We can see that 40 channels are bad--but do we need to ignore any more? These channels are bad because of certain instrument properties: other channels still may need to be ignored because of the shape and brightness of the spectrum itself. Figure B shows the result of plotting the data in channels (using the commands setplot channel and plot data). We see that above about channel 33 the S/N becomes small. We also see, comparing Figure B with Figure A, which bad channels were ignored. Although visual inspection is not the most rigorous method for deciding which channels to ignore (more on this subject later), it's good enough for now, and will at least prevent us from getting grossly misleading results from the fitting. To ignore channels above 33:

XSPEC> ignore 34-**
 Chi-Squared =      677.6218     using    31 PHA bins.
 Reduced chi-squared =      24.20078     for     28 degrees of freedom
 Null hypothesis probability =    0.

The same result can be achieved with the command ignore 34-125. The inverse of ignore is notice, which has the same syntax.

We are now ready to fit the data. Fitting is initiated by the command fit. As the fit proceeds, the screen displays the status of the fit for each iteration until either the fit converges to the minimum $\chi^2$, or the user is asked whether the fit is to go through another set of iterations to find the minimum. The default number of iterations is ten.

XSPEC>fit
 Chi-Squared    Lvl  Fit param # 1     2           3
   204.136     -3     7.9869E-02   1.564      4.4539E-03
   84.5658     -4     0.3331       2.234      1.0977E-02
   30.2511     -5     0.4422       2.174      1.1965E-02
   30.1202     -6     0.4648       2.196      1.2264E-02
   30.1189     -7     0.4624       2.195      1.2244E-02
 ---------------------------------------------------------------------------
 Variances and Principal axes :
               1     2     3
 4.14E-08 |  0.00 -0.01  1.00
 8.70E-02 | -0.91 -0.41 -0.01
 2.32E-03 | -0.41  0.91  0.01
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
  Model:  phabs[1]( powerlaw[2] )
  Model Fit Model Component  Parameter  Unit     Value
  par   par comp
    1    1    1   phabs      nH       10^22     0.4624     +/-  0.2698
    2    2    2   powerlaw   PhoIndex            2.195     +/-  0.1288
    3    3    2   powerlaw   norm               1.2244E-02 +/-  0.2415E-02
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
 Chi-Squared =      30.11890     using    31 PHA bins.
 Reduced chi-squared =      1.075675     for     28 degrees of freedom
 Null hypothesis probability = 0.358

The fit is good: reduced $\chi^2$ is 1.075 for 31 - 3 = 28 degrees of freedom. The null hypothesis probability is the probability of getting a value of $\chi^2$ as large or larger than observed if the model is correct. If this probability is small then the model is not a good fit. The matrix of principal axes printed out at the end of a fit provides an indication of whether parameters are correlated (at least local to the best fit). In this example the powerlaw norm is not correlated with any other parameter while the column and powerlaw index are slightly correlated.

To see the fit and the residuals, we use the command

XSPEC>plot data resid

The result is shown in Figure C.

\begin{figure}
\begin{center}
\mbox{
{\psfig{figure=figures/figc.eps,width=4.0in...
...g absorbed power-law model; the residuals of the fit.}}
\end{center}\end{figure}

The screen output shows the best-fitting parameter values, as well as approximations to their errors. These errors should be regarded as indications of the uncertainties in the parameters and should not be quoted in publications. The true errors, i.e. the confidence ranges, are obtained using the error command:

XSPEC>error 1 2 3
 Parameter   Confidence Range (     2.706)
     1    3.254765E-02    0.932778
     2     1.99397         2.40950
 *WARNING*:RENORM: No variable model to allow renormalization
 *WARNING*:RENORM: No variable model to allow renormalization
 *WARNING*:RENORM: No variable model to allow renormalization
 *WARNING*:RENORM: No variable model to allow renormalization
 *WARNING*:RENORM: No variable model to allow renormalization
 *WARNING*:RENORM: No variable model to allow renormalization
 *WARNING*:RENORM: No variable model to allow renormalization
 *WARNING*:RENORM: No variable model to allow renormalization
 *WARNING*:RENORM: No variable model to allow renormalization
 *WARNING*:RENORM: No variable model to allow renormalization
     3    8.942987E-03    1.711637E-02

Here, the numbers 1, 2, 3 refer to the parameter numbers in the Model par column of the screen output. For the first parameter, the column of absorbing hydrogen atoms, the 90% confidence range is $3.0\times 10^{20} < N_{H} < 8.6\times 10^{21}$cm-2. This corresponds to an excursion in $\chi^2$ of 2.706. The reason these ``better'' errors are not given automatically as part of the fit output is that they entail further fitting. When the model is simple, this does not require much CPU, but for complicated models the extra time can be considerable. The warning message is generated because there are no free normalizations in the model while the error is being calculated on the normalization itself. In this case, the warning may safely be ignored.

What else can we do with the fit? One thing is to derive the flux of the model. The data by themselves only give the instrument-dependent count rate. The model, on the other hand, is an estimate of the true spectrum emitted. In XSPEC, the model is defined in physical units independent of the instrument. The command flux integrates the current model over the range specified by the user:

XSPEC> flux 2 10
 Model flux  3.5496E-03 photons ( 2.2492E-11 ergs)cm**-2 s**-1 (  2.000- 10.000)

Here, we have chosen the standard X-ray range of 2-10 keV and find that the energy flux is $2.2\times 10^{-11}$ erg cm-2 s-1. Note that flux will integrate only within the energy range of the current response matrix. If the model flux outside this range is desired--in effect, an extrapolation beyond the data--then the command dummyrsp should be used. This command sets up a dummy response that covers the range required. For example, if we want to know the flux of our model in the ROSAT PSPC band of 0.2-2 keV, we enter:

XSPEC>dummy 0.2 2.
 Chi-Squared =      3583.779     using    31 PHA bins.
 Reduced chi-squared =      127.9921     for     28 degrees of freedom
 Null hypothesis probability =    0.
XSPEC>flux 0.2 2.
 Model flux  4.5306E-03 photons ( 9.1030E-12 ergs)cm**-2 s**-1 (  0.200-  2.000)

The energy flux, at $9.1\times 10^{-12}$ erg cm-2 s-1, is lower in this band but the photon flux is higher. To get our original response matrix back we enter :

XSPEC> response
 Chi-Squared =      30.11890     using    31 PHA bins.
 Reduced chi-squared =      1.075675     for     28 degrees of freedom
 Null hypothesis probability = 0.358

The fit, as we've remarked, is good, and the parameters are constrained. But unless the purpose of our investigation is merely to measure a photon index, it's a good idea to check whether alternative models can fit the data just as well. We also should derive upper limits to components such as iron emission lines and additional continua, which, although not evident in the data nor required for a good fit, are nevertheless important to constrain. First, let's try an absorbed black body:

XSPEC>mo pha(bb)
  Model:  phabs[1]( bbody[2] )
Input parameter value, delta, min, bot, top, and max values for ...
Current:          1     0.001         0         0     1E+05     1E+06
phabs:nH>/*
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
  Model:  phabs[1]( bbody[2] )
  Model Fit Model Component  Parameter  Unit     Value
  par   par comp
    1    1    1   phabs      nH       10^22      1.000     +/-      0.
    2    2    2   bbody      kT       keV        3.000     +/-      0.
    3    3    2   bbody      norm                1.000     +/-      0.
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
 Chi-Squared =     3.3142067E+09 using    31 PHA bins.
 Reduced chi-squared =     1.1836453E+08 for     28 degrees of freedom
 Null hypothesis probability =    0.
XSPEC>fit
 Chi-Squared    Lvl  Fit param # 1     2           3
   1420.96      0     0.2116       2.987      7.7666E-04
   1387.72      0     4.4419E-02   2.975      7.7003E-04
   1376.39      0     4.3009E-03   2.963      7.6354E-04
   1371.67      0     1.8192E-03   2.951      7.5730E-04
   1367.04      0     5.8100E-04   2.939      7.5130E-04
   1362.47      0     1.9059E-04   2.926      7.4548E-04
   1357.84      0     6.1455E-05   2.913      7.3982E-04
   1353.14      0     2.5356E-05   2.900      7.3429E-04
   1348.36      0     6.7582E-06   2.887      7.2885E-04
   1343.50      0     2.0479E-06   2.874      7.2350E-04
 Number of trials exceeded - last iteration delta =    4.861
 Continue fitting? (Y)y 
 ...
 
   113.954      0         0.      0.8907      2.7865E-04
   113.954     -1         0.      0.8905      2.7859E-04
   113.954      4         0.      0.8905      2.7859E-04
 ---------------------------------------------------------------------------
 Variances and Principal axes :
               2     3
 2.88E-04 | -1.00  0.00
 8.45E-11 |  0.00 -1.00
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
  Model:  phabs[1]( bbody[2] )
  Model Fit Model Component  Parameter  Unit     Value
  par   par comp
    1    1    1   phabs      nH       10^22         0.     +/-  -1.000
    2    2    2   bbody      kT       keV       0.8905     +/-  0.1696E-01
    3    3    2   bbody      norm               2.7859E-04 +/-  0.9268E-05
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
 Chi-Squared =      113.9542     using    31 PHA bins.
 Reduced chi-squared =      4.069792     for     28 degrees of freedom
 Null hypothesis probability = 2.481E-12

Note that when more than 10 iterations are required for convergence the user is asked whether or not to continue at the end of each set of 10. Saying no at these prompts is a good idea if the fit is not converging quickly. Conversely, to avoid having to keep answering the question, i.e., to increase the number of iterations before the prompting question appears, begin the fit with, say fit 100. This command will put the fit through 100 iterations before pausing.

Plotting the data and residuals again with

XSPEC> plot data resid

we obtain Figure D.

\begin{figure}
\begin{center}
\mbox{
{\psfig{figure=figures/figd.eps,width=4.0in...
...t the continuum is obviously {\sl
not} a black body.}}
\end{center}\end{figure}

The black body fit is obviously not a good one. Not only is $\chi^2$ large, but the best-fitting NH is rather low. Inspection of the residuals confirms this: the pronounced wave-like shape is indicative of a bad choice of overall continuum (see Figure D). Let's try thermal bremsstrahlung next:

XSPEC>mo pha(br)
  Model:  phabs[1]( bremss[2] )
Input parameter value, delta, min, bot, top, and max values for ...
Current:          1     0.001         0         0     1E+05     1E+06
phabs:nH>/*
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
  Model:  phabs[1]( bremss[2] )
  Model Fit Model Component  Parameter  Unit     Value
  par   par comp
    1    1    1   phabs      nH       10^22      1.000     +/-      0.
    2    2    2   bremss     kT       keV        7.000     +/-      0.
    3    3    2   bremss     norm                1.000     +/-      0.
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
 Chi-Squared =     4.5311800E+07 using    31 PHA bins.
 Reduced chi-squared =      1618279.     for     28 degrees of freedom
 Null hypothesis probability =    0.
XSPEC>fit
 Chi-Squared    Lvl  Fit param # 1     2           3
   113.305     -3     0.2441       6.557      6.8962E-03
   40.4519     -4     0.1173       5.816      7.7944E-03
   36.0549     -5     4.4750E-02   5.880      7.7201E-03
   33.4168     -6     1.8882E-02   5.868      7.7476E-03
   32.6766     -7     7.8376E-03   5.864      7.7495E-03
   32.3192     -8     2.7059E-03   5.862      7.7515E-03
   32.1512     -9     2.3222E-04   5.861      7.7525E-03
   32.1471    -10     1.0881E-04   5.861      7.7523E-03
 ---------------------------------------------------------------------------
 Variances and Principal axes :
               1     2     3
 2.29E-08 |  0.00  0.00  1.00
 3.18E-02 |  0.95  0.31  0.00
 8.25E-01 |  0.31 -0.95  0.00
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
  Model:  phabs[1]( bremss[2] )
  Model Fit Model Component  Parameter  Unit     Value
  par   par comp
    1    1    1   phabs      nH       10^22     1.0881E-04 +/-  0.3290
    2    2    2   bremss     kT       keV        5.861     +/-  0.8651
    3    3    2   bremss     norm               7.7523E-03 +/-  0.8122E-03
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
 Chi-Squared =      32.14705     using    31 PHA bins.
 Reduced chi-squared =      1.148109     for     28 degrees of freedom
 Null hypothesis probability = 0.269

Bremsstrahlung is a better fit than the black body--and is as good as the power law--although it shares the low NH. With two good fits, the power law and the bremsstrahlung, it's time to scrutinize their parameters in more detail.

First, we reset our fit to the powerlaw (output omitted):

XSPEC>mo pha(po)

From the EXOSAT database on HEASARC, we know that the target in question, 1E1048.1-5937, has a Galactic latitude of 24 arcmin, i.e., almost on the plane of the Galaxy. In fact, the data base also provides the value of the Galactic NH based on 21-cm radio observations. At $4\times 10^{22}$ cm-2, it is higher than the 90 percent-confidence upper limit from the power-law fit. Perhaps, then, the power-law fit is not so good after all. What we can do is fix (freeze in XSPEC terminology) the value of NH at the Galactic value and refit the power law. Although we won't get a good fit, the shape of the residuals might give us a clue to what is missing. To freeze a parameter in XSPEC, use the command freeze followed by the parameter number, like this:

XSPEC> freeze 1
 Number of variable fit parameters =    2

The inverse of freeze is thaw:

XSPEC> thaw 1
 Number of variable fit parameters =    3

Alternatively, parameters can be frozen using the newpar command, which allows all the quantities associated with a parameter to be changed. The second quantity, delta, is the step size used to calculate the derivative in the fitting, and, if set to a negative number, will cause the parameter to be frozen. In our case, we want NH frozen at $4\times 10^{22}$ cm-2, so we go back to the power law best fit and do the following :

XSPEC>newpar 1
Current:      0.463     0.001         0         0     1E+05     1E+06
phabs:nH>4,0
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
  Model:  phabs[1]( powerlaw[2] )
  Model Fit Model Component  Parameter  Unit     Value
  par   par comp
    1    1    1   phabs      nH       10^22      4.000     frozen
    2    2    2   powerlaw   PhoIndex            2.195     +/-  0.1287
    3    3    2   powerlaw   norm               1.2247E-02 +/-  0.2412E-02
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
    2 variable fit parameters
 Chi-Squared =      829.3545     using    31 PHA bins.
 Reduced chi-squared =      28.59843     for     29 degrees of freedom
 Null hypothesis probability =    0.

Note the useful trick of giving a value of zero for delta in the newpar command. This has the effect of changing delta to the negative of its current value. If the parameter is free, it will be frozen, and if frozen, thawed. The same result can be obtained by putting everything onto the command line, i.e., newpar 1 4, 0, or by issuing the two commands, newpar 1 4 followed by freeze 1. Now, if we fit and plot again, we get the following model (Fig. E).

XSPEC>fit
...
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
  Model:  phabs[1]( powerlaw[2] )
  Model Fit Model Component  Parameter  Unit     Value
  par   par comp
    1    1    1   phabs      nH       10^22      4.000     frozen
    2    2    2   powerlaw   PhoIndex            3.594     +/-  0.6867E-01
    3    3    2   powerlaw   norm               0.1161     +/-  0.9412E-02
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
 Chi-Squared =      125.5134     using    31 PHA bins.
 Reduced chi-squared =      4.328048     for     29 degrees of freedom
 Null hypothesis probability = 5.662E-14
 XSPEC>plot data resid

\begin{figure}
\begin{center}
\mbox{
{\psfig{figure=figures/fige.eps,width=4.0in...
...ents might be added to the model to
improve the fit.}}
\end{center}\end{figure}

The fit is not good. In Figure E we can see why: there appears to be a surplus of softer photons, perhaps indicating a second continuum component. To investigate this possibility we can add a component to our model. The editmod command lets us do this without having to respecify the model from scratch. Here, we'll add a black body component.

XSPEC>editmod pha(po+bb)
  Model:  phabs[1]( powerlaw[2] + bbody[3] )
Input parameter value, delta, min, bot, top, and max values for ...
Current:          3      0.01    0.0001      0.01       100       200
bbody:kT>2,0
Current:          1      0.01         0         0     1E+24     1E+24
bbody:norm>1.e-5
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
  Model:  phabs[1]( powerlaw[2] + bbody[3] )
  Model Fit Model Component  Parameter  Unit     Value
  par   par comp
    1    1    1   phabs      nH       10^22      4.000     frozen
    2    2    2   powerlaw   PhoIndex            3.594     +/-  0.6867E-01
    3    3    2   powerlaw   norm               0.1161     +/-  0.9412E-02
    4    4    3   bbody      kT       keV        2.000     frozen
    5    5    3   bbody      norm               1.0000E-05 +/-      0.
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
 Chi-Squared =      122.1538     using    31 PHA bins.
 Reduced chi-squared =      4.362635     for     28 degrees of freedom
 Null hypothesis probability = 9.963E-14

Notice that in specifying the initial values of the black body, we have frozen kT at 2 keV (the canonical temperature for nuclear burning on the surface of a neutron star in a low-mass X-ray binary) and started the normalization at zero. Without these measures, the fit might ``lose its way''. Now, if we fit, we get (not showing all the iterations this time):

  ---------------------------------------------------------------------------
  Model:  phabs[1]( powerlaw[2] + bbody[3] )
  Model Fit Model Component  Parameter  Unit     Value
  par   par comp
    1    1    1   phabs      nH       10^22      4.000     frozen
    2    2    2   powerlaw   PhoIndex            4.932     +/-  0.1618
    3    3    2   powerlaw   norm               0.3761     +/-  0.5449E-01
    4    4    3   bbody      kT       keV        2.000     frozen
    5    5    3   bbody      norm               2.3212E-04 +/-  0.3966E-04
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
 Chi-Squared =      55.63374     using    31 PHA bins.
 Reduced chi-squared =      1.986919     for     28 degrees of freedom
 Null hypothesis probability = 1.425E-03

The fit is better than the one with just a power law and the fixed Galactic column, but it is still not good. Thawing the black body temperature and fitting gives us:

XSPEC>thaw 4
 Number of variable fit parameters =    4
XSPEC>fit
...
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
  Model:  phabs[1]( powerlaw[2] + bbody[3] )
  Model Fit Model Component  Parameter  Unit     Value
  par   par comp
    1    1    1   phabs      nH       10^22      4.000     frozen
    2    2    2   powerlaw   PhoIndex            6.401     +/-  0.3873
    3    3    2   powerlaw   norm                1.086     +/-  0.3032
    4    4    3   bbody      kT       keV        1.199     +/-  0.8082E-01
    5    5    3   bbody      norm               2.6530E-04 +/-  0.3371E-04
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
 Chi-Squared =      37.21207     using    31 PHA bins.
 Reduced chi-squared =      1.378225     for     27 degrees of freedom
 Null hypothesis probability = 9.118E-02

This, of course, is a better fit, but the photon index of the power law has ended up extremely and implausibly steep. Looking at this odd model with the command

XSPEC> plot model

we see, in Figure F, that the black body and the power law have changed places, in that the power law provides the soft photons required by the high absorption, while the black body provides the harder photons.

\begin{figure}
\begin{center}
\mbox{
{\psfig{figure=figures/figf.eps,width=4.0in...
...ponents (absorbed
to the same degree) and their sum.}}
\end{center}\end{figure}

We could continue to search for a plausible, well-fitting model, but the data, with their limited signal-to-noise and energy resolution, probably don't warrant it (the original investigators published only the power law fit). There is, however, one final, useful thing to do with the data: derive an upper limit to the presence of a fluorescent iron emission line. First we delete the black body component using delcomp:

XSPEC>delcomp 3
  Model:  phabs[1]( powerlaw[2] )
 Chi-Squared =      1285.487     using    31 PHA bins.
 Reduced chi-squared =      44.32712     for     29 degrees of freedom
 Null hypothesis probability =    0.

Then we thaw NH and refit to recover our original, best fit:

XSPEC>thaw 1
 Number of variable fit parameters =    3
XSPEC>fit
 Chi-Squared    Lvl  Fit param # 1     2           3
   924.178     -2      5.087       5.076      0.4056
   305.507     -2      4.525       3.791      0.1249
   140.460     -2      2.930       3.367      6.5553E-02
   64.4275     -3     0.6068       2.244      1.4635E-02
   30.3738     -4     0.4837       2.201      1.2279E-02
   30.1189     -5     0.4641       2.195      1.2258E-02
   30.1189     -6     0.4637       2.195      1.2255E-02
 ---------------------------------------------------------------------------
 Variances and Principal axes :
               1     2     3
 4.13E-08 |  0.00 -0.01  1.00
 8.69E-02 | -0.91 -0.41 -0.01
 2.31E-03 | -0.41  0.91  0.01
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
  Model:  phabs[1]( powerlaw[2] )
  Model Fit Model Component  Parameter  Unit     Value
  par   par comp
    1    1    1   phabs      nH       10^22     0.4637     +/-  0.2696
    2    2    2   powerlaw   PhoIndex            2.195     +/-  0.1287
    3    3    2   powerlaw   norm               1.2255E-02 +/-  0.2412E-02
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
 Chi-Squared =      30.11890     using    31 PHA bins.
 Reduced chi-squared =      1.075675     for     28 degrees of freedom
 Null hypothesis probability = 0.358

Now, we add a gaussian emission line of fixed energy and width:

XSPEC>editmod pha(po+ga)
  Model:  phabs[1]( powerlaw[2] + gaussian[3] )
Input parameter value, delta, min, bot, top, and max values for ...
Current:        6.5      0.05         0         0     1E+06     1E+06
gaussian:LineE>6.4 0
Current:        0.1      0.05         0         0        10        20
gaussian:Sigma>0.1 0
Current:          1      0.01         0         0     1E+24     1E+24
gaussian:norm>1.e-4
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
  Model:  phabs[1]( powerlaw[2] + gaussian[3] )
  Model Fit Model Component  Parameter  Unit     Value
  par   par comp
    1    1    1   phabs      nH       10^22     0.4637     +/-  0.2696
    2    2    2   powerlaw   PhoIndex            2.195     +/-  0.1287
    3    3    2   powerlaw   norm               1.2255E-02 +/-  0.2412E-02
    4    4    3   gaussian   LineE    keV        6.400     frozen
    5    5    3   gaussian   Sigma    keV       0.1000     frozen
    6    6    3   gaussian   norm               1.0000E-04 +/-      0.
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
 Chi-Squared =      32.75002     using    31 PHA bins.
 Reduced chi-squared =      1.212964     for     27 degrees of freedom
 Null hypothesis probability = 0.205
XSPEC>fit
...
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
  Model:  phabs[1]( powerlaw[2] + gaussian[3] )
  Model Fit Model Component  Parameter  Unit     Value
  par   par comp
    1    1    1   phabs      nH       10^22     0.6562     +/-  0.3193
    2    2    2   powerlaw   PhoIndex            2.324     +/-  0.1700
    3    3    2   powerlaw   norm               1.4636E-02 +/-  0.3642E-02
    4    4    3   gaussian   LineE    keV        6.400     frozen
    5    5    3   gaussian   Sigma    keV       0.1000     frozen
    6    6    3   gaussian   norm               9.6462E-05 +/-  0.9542E-04
  ---------------------------------------------------------------------------
  ---------------------------------------------------------------------------
 Chi-Squared =      29.18509     using    31 PHA bins.
 Reduced chi-squared =      1.080929     for     27 degrees of freedom
 Null hypothesis probability = 0.352

The energy and width have to be frozen because,in the absence of an obvious line in the data, the fit would be completely unable to converge on meaningful values. Besides, our aim is to see how bright a line at 6.4 keV can be and still not ruin the fit. To do this, we fit first and then use the error command to derive the maximum allowable iron line normalization. We then set the normalization at this maximum value with newpar and, finally, derive the equivalent width using the eqwidth command. That is:

XSPEC>err 6
 Parameter   Confidence Range (     2.706)
 Parameter pegged at hard limit         0.
 with delta ftstat=    0.9338
     6          0.        1.530722E-04
XSPEC>new 6  1.530722E-04
    4 variable fit parameters
 Chi-Squared =      34.91923     using    31 PHA bins.
 Reduced chi-squared =      1.293305     for     27 degrees of freedom
 Null hypothesis probability = 0.141
XSPEC>eqwidth 3
 Additive group equiv width for model 3 (gaussian): 839. eV

Things to note:


next up [*] [*]
Next: Simultaneous Fitting: Examples from Up: Walks through XSPEC Previous: Brief Discussion of XSPEC
Ben Dorman
2003-11-28