Last modified: 11 Dec 2018

URL: http://cxc.harvard.edu/sherpa/threads/lines_acishetg/

Measuring Line Parameters with an HETG/ACIS-S Spectrum

Sherpa Threads (CIAO 4.11 Sherpa v1)


Overview

Synopsis:

After having created or downloaded a set of PHA2 and response files (RMFs and ARFs), for a HETG/ACIS-S observation, perform a simple fit to one of the line features present in the spectrum. This thread takes the user through a calculation of the error bars on the line normalizations, positions, and widths, and shows how to calculate the line equivalent widths.

For those wishing to perform their analysis with XSPEC, the PHA2 file must first be split into individual type 1 PHA files, and then grouped to increase the signal-to-noise in each spectral bin. The Grouping a Grating Spectrum CIAO thread shows how to do this, and the procedure is also included within this thread.

Run this thread if:

You are working with a HETG/ACIS-S data set, and want to begin a simple analysis of the data.

Related Links:

Last Update: 11 Dec 2018 - reviewed for CIAO 4.11, revised screen output


Contents


Data Preparation

This thread makes use of archival data products for Vela X-1, downloaded from TGcat. They were obtained by performing a "Search by Name" on Vela X-1 (using the SIMBAD name resolver), then choosing "viewfile table" for the target, and clicking the download arrow next to each filename. These files were unzipped and placed in the directory in which the data analysis is to be performed.

Note that we are not considering any background files. For HETG/ACIS-S observations, the background is suppressed by the process of order sorting. That is, an event at a given location along the grating arms must fall into a narrow range of energies (as determined by the CCD) in order for it to be considered a valid dispersed photon. Most background events do not meet these "order sorting" criteria. For this, and many other HETG/ACIS-S observations, the first simple analysis can be performed without reference to the background. (Detailed analysis of low S/N spectra likely would require proper consideration of the background. See the Fitting Grating Data thread for an example of analysis including background data).

Those using either Sherpa or ISIS to perform their analysis can use the data files directly in the analysis; skip to the "Load the Spectrum and Responses" section.

XSPEC users, however, must split and group the spectra before the analysis as described here.


Load the Spectrum and Responses

Start a Sherpa session:

unix% sherpa
-----------------------------------------------------
Welcome to Sherpa: CXC's Modeling and Fitting Package
-----------------------------------------------------
CIAO 4.11 Sherpa version 1 Wednesday, December 5, 2018

Python 3.5.4 (default, Oct 15 2018, 13:47:46) 
Type 'copyright', 'credits' or 'license' for more information
IPython 6.5.0 -- An enhanced Interactive Python. Type '?' for help.

IPython profile: sherpa

sherpa>

We begin by reading the PHA2 file that contains the grating spectra. This file contains 12 spectra, 3 spectral orders each for both positive and negative order HEG and MEG spectra. They are stored in the order: HEG -3, -2, -1, 1, 2, 3; MEG -3, -2, -1, 1, 2, 3. By default, Sherpa will read in all 12 spectral files in this order.

We are only interested, however, in the first-order spectra (positive and negative); therefore, we use a Data Model-type filter to select the third and fourth, ninth, and tenth rows of the PHA2 file (HEG -/+1 and MEG -/+1, respectively).

sherpa> load_pha("pha2[#row=3,4,9,10]")    
WARNING: systematic errors were not found in file 'pha2[#row=3,4,9,10]'
statistical errors were found in file 'pha2[#row=3,4,9,10]' 
but not used; to use them, re-read with use_errors=True
read background_up into a dataset from file pha2[#row=3,4,9,10]
read background_down into a dataset from file pha2[#row=3,4,9,10]
Multiple data sets have been input: 1-4

Alternatively, one could use a filter to select all negative and postitive first-order spectra via the TG_M keyword:

sherpa> load_pha("pha2[TG_M=-1,1]")

Note that associated background spectra (two each for each spectrum, one from either side of the gratings arm) will be read in with the source spectra. However, unless the subtract command is called, or unless one creates a background model, they will not be part of the fitting process. For the remainder of this thread, we ignore these background spectra.

We now read in the associated ARF files and assign them to their corresponding data sets, then do the same for the RMF files.

sherpa> load_arf(1,"heg_-1.arf")
sherpa> load_arf(2,"heg_1.arf")
sherpa> load_arf(3,"meg_-1.arf")
sherpa> load_arf(4,"meg_1.arf")

sherpa> load_rmf(1,"heg_-1.rmf") 
sherpa> load_rmf(2,"heg_1.rmf")
sherpa> load_rmf(3,"meg_-1.rmf")
sherpa> load_rmf(4,"meg_1.rmf")

Group and Filter the Data

Sherpa allows one to bin the spectra during the analysis session; see the Sherpa thread Changing the grouping scheme of a data set within Sherpa for details. We bin each of the spectra to have a minimum of 20 counts per energy channel.

sherpa> group_counts(1,20)
sherpa> group_counts(2,20) 
sherpa> group_counts(3,20)
sherpa> group_counts(4,20)

For the purposes of this analysis, we are interested in fitting the line near 6.4 keV. We restrict the energy range to 5-7 keV for all four spectra using the notice command, since the same energy filter is to be applied to all data set IDs; the notice_id command is available for filtering data sets individually.

sherpa> notice(5.,7.)

sherpa> show_filter()

Data Set Filter: 1
5.0096-6.9904 Energy (keV)

Data Set Filter: 2
4.9721-6.9910 Energy (keV)

Data Set Filter: 3
4.9999-9.5411 Energy (keV)

Data Set Filter: 4
5.0154-9.5411 Energy (keV)

The spectra can all be plotted in the same figure using the plot command. We select all four plots to be considered "current" so that changing the axis will be applied to all of them, and then we choose the y-axes to be logarithmic.

sherpa> plot("data",1,"data",2,"data",3,"data",4)
sherpa> current_plot("all")
sherpa> log_scale(Y_AXIS)
sherpa> print_window("all",["format","eps"])

The plot, shown in Figure 1, is also saved in EPS format to the file "all.eps".

Figure 1: Plot of the +1 and -1 HEG and MEG spectra

Figure 1: Plot of the +1 and -1 HEG and MEG spectra

The HEG (upper) and MEG (lower) plots with -1 order on the left and and the +1 order on the right. The plots are displayed with a logarithmic scale and have been filtered to the range 5.0-7.0 keV. Note that the MEG data includes a data point above 9 keV, as the lower bound of this energy bin falls below the 7 keV cutoff of the notice/notice_id command.

The MEG spectra do not have very good statistics in the region of interest; therefore, we delete these data sets from the current session before proceeding with the fits.

sherpa> delete_data(3)
sherpa> delete_data(4)

The data sets currently being considered in the analysis can be viewed with the show_data, show_bkg, and show_all commands:

sherpa> show_data()
Data Set: 1
Filter: 5.0096-6.9904 Energy (keV)
Bkg Scale: 0.248829
Noticed Channels: 7597-7888
name           = pha2[#row=3,4,9,10]
channel        = Float64[8192]
counts         = Float64[8192]
staterror      = None
syserror       = None
bin_lo         = Float64[8192]
bin_hi         = Float64[8192]
grouping       = Float64[8192]
quality        = Float64[8192]
exposure       = 83150.1366728
backscal       = 1.0
areascal       = 1.0
grouped        = True
subtracted     = False
units          = energy
rate           = True
plot_fac       = 0
response_ids   = [1]
background_ids = [1, 2]

RMF Data Set: 1:1
name     = heg_-1.rmf
detchans = 8192
energ_lo = Float64[8192]
energ_hi = Float64[8192]
n_grp    = UInt64[8192]
f_chan   = UInt32[8192]
n_chan   = UInt32[8192]
matrix   = Float64[841124]
offset   = 1
e_min    = Float64[8192]
e_max    = Float64[8192]
ethresh  = 1e-10

ARF Data Set: 1:1
name     = heg_-1.arf
energ_lo = Float64[8192]
energ_hi = Float64[8192]
specresp = Float64[8192]
bin_lo   = Float64[8192]
bin_hi   = Float64[8192]
exposure = 83149.4231809
ethresh  = 1e-10

Background Data Set: 1:1
Filter: 5.3913,8.8601 Energy (keV)
Noticed Channels: 7435-8150
name           = pha2[#row=3,4,9,10]
channel        = Float64[8192]
counts         = Float64[8192]
staterror      = None
syserror       = None
bin_lo         = Float64[8192]
bin_hi         = Float64[8192]
grouping       = Float64[8192]
quality        = Float64[8192]
exposure       = 83150.1366728
backscal       = 4.0188284
areascal       = None
grouped        = True
subtracted     = False
units          = energy
rate           = True
plot_fac       = 0
response_ids   = []
background_ids = []

Background Data Set: 1:2
Filter: 4.8200,9.1407 Energy (keV)
Noticed Channels: 7273-8192
name           = pha2[#row=3,4,9,10]
channel        = Float64[8192]
counts         = Float64[8192]
staterror      = None
syserror       = None
bin_lo         = Float64[8192]
bin_hi         = Float64[8192]
grouping       = Float64[8192]
quality        = Float64[8192]
exposure       = 83150.1366728
backscal       = 4.0188284
areascal       = None
grouped        = True
subtracted     = False
units          = energy
rate           = True
plot_fac       = 0
response_ids   = []
background_ids = []

Data Set: 2
Filter: 4.9721-6.9910 Energy (keV)
Bkg Scale: 0.248829
Noticed Channels: 7588-7891
name           = pha2[#row=3,4,9,10]
channel        = Float64[8192]
counts         = Float64[8192]
staterror      = None
syserror       = None
bin_lo         = Float64[8192]
bin_hi         = Float64[8192]
grouping       = Float64[8192]
quality        = Float64[8192]
exposure       = 83150.1366728
backscal       = 1.0
areascal       = 1.0
grouped        = True
subtracted     = False
units          = energy
rate           = True
plot_fac       = 0
response_ids   = [1]
background_ids = [1, 2]

RMF Data Set: 2:1
name     = heg_1.rmf
detchans = 8192
energ_lo = Float64[8192]
energ_hi = Float64[8192]
n_grp    = UInt64[8192]
f_chan   = UInt32[8192]
n_chan   = UInt32[8192]
matrix   = Float64[841124]
offset   = 1
e_min    = Float64[8192]
e_max    = Float64[8192]
ethresh  = 1e-10

ARF Data Set: 2:1
name     = heg_1.arf
energ_lo = Float64[8192]
energ_hi = Float64[8192]
specresp = Float64[8192]
bin_lo   = Float64[8192]
bin_hi   = Float64[8192]
exposure = 83150.6981615
ethresh  = 1e-10

Background Data Set: 2:1
Filter: 7.7908 Energy (keV)
Noticed Channels: 7035-8192
name           = pha2[#row=3,4,9,10]
channel        = Float64[8192]
counts         = Float64[8192]
staterror      = None
syserror       = None
bin_lo         = Float64[8192]
bin_hi         = Float64[8192]
grouping       = Float64[8192]
quality        = Float64[8192]
exposure       = 83150.1366728
backscal       = 4.0188284
areascal       = None
grouped        = True
subtracted     = False
units          = energy
rate           = True
plot_fac       = 0
response_ids   = []
background_ids = []

Background Data Set: 2:2
Filter: 3.7843,6.8976 Energy (keV)
Noticed Channels: 6538-8018
name           = pha2[#row=3,4,9,10]
channel        = Float64[8192]
counts         = Float64[8192]
staterror      = None
syserror       = None
bin_lo         = Float64[8192]
bin_hi         = Float64[8192]
grouping       = Float64[8192]
quality        = Float64[8192]
exposure       = 83150.1366728
backscal       = 4.0188284
areascal       = None
grouped        = True
subtracted     = False
units          = energy
rate           = True
plot_fac       = 0
response_ids   = []
background_ids = []

Defining the Source Models

We will first fit a broad continuum to the spectra, without any lines. The XSPEC power-law model (xspowerlaw) is chosen for the continuum. It is assigned as the source for data set 1 and given the model name powr; the model name is then used to assign the same source model to data set 2.

sherpa> set_source(1,xspowerlaw.powr)
sherpa> set_source(2,powr)

sherpa> show_source()
Model: 1
xspowerlaw.powr
   Param         Type          Value          Min          Max      Units
   -----         ----          -----          ---          ---      -----
   powr.phoindex thawed           1           -2            9           
   powr.norm     thawed           1            0        1e+24           

Model: 2
xspowerlaw.powr
   Param         Type          Value          Min          Max      Units
   -----         ----          -----          ---          ---      -----
   powr.phoindex thawed            1           -2            9           
   powr.norm     thawed            1            0        1e+24           

Fitting

The fit statistic is set to be χ2 with the data counts acting as the variance (i.e., a Gaussian approximation to simple Poisson statistics). Both data sets are fit simultaneously.

sherpa> set_stat("chi2datavar")

sherpa> fit()
WARNING: data set 1 has associated backgrounds, but they have not been subtracted, nor have background models been set
WARNING: data set 2 has associated backgrounds, but they have not been subtracted, nor have background models been set
Data Sets              = 1, 2
Method                = levmar
Statistic             = chi2datavar
Initial fit statistic = 5.98398e+08
Final fit statistic   = 321.803 at function evaluation 77
Data points           = 56
Degrees of freedom    = 54
Probability [Q-value] = 9.16795e-40
Reduced statistic     = 5.95931
Change in statistic   = 5.98397e+08
   powr.PhoIndex   -1.36276     +/- 0.336331    
   powr.norm      1.98057e-05  +/- 1.18241e-05 

At any time, the results of the last fit can be displayed using the show_fit command. The fit yields the following results:

sherpa> show_fit()
Optimization Method: LevMar
name    = levmar
ftol    = 1.19209289551e-07
xtol    = 1.19209289551e-07
gtol    = 1.19209289551e-07
maxfev  = None
epsfcn  = 1.19209289551e-07
factor  = 100.0
verbose = 0

Statistic: Chi2DataVar
Chi Squared with data variance.

    The variance in each bin is estimated from the data value in that
    bin. See also `Chi2Gehrels`, `Chi2XSpecVar` and `Chi2ModVar`.

    If the number of counts in each bin is large, then the shape of
    the Poisson distribution from which the counts are sampled tends
    asymptotically towards that of a Gaussian distribution, with
    variance

    sigma(i)^2 = N(i,S) + [A(S)/A(B)]^2 N(i,B)

    where N is the number of on-source (and off-source) bins included
    in the fit. The background term appears only if an estimate of the
    background has been subtracted from the data.

    See Also
    --------
    Chi2Gehrels, Chi2ModVar, Chi2XspecVar

    

Fit:Datasets              = 1, 2
Method                = levmar
Statistic             = chi2datavar
Initial fit statistic = 5.98398e+08
Final fit statistic   = 321.803 at function evaluation 77
Data points           = 56
Degrees of freedom    = 54
Probability [Q-value] = 9.16795e-40
Reduced statistic     = 5.95931
Change in statistic   = 5.98397e+08
   powr.PhoIndex   -1.36276     +/- 0.336331    
   powr.norm      1.98057e-05  +/- 1.18241e-05 

The results of these fits can then be plotted. Here we choose to display just the data and fits, not the residuals, and use logarithmic y-axes.

sherpa> plot("fit",1,"fit",2)
sherpa> current_plot("all")
sherpa> log_scale(Y_AXIS)
sherpa> print_window("pl",["format","eps"])

The plot, shown in Figure 2, is also saved in EPS format to the file "pl.eps".

Figure 2: Continuum Fit to the Data

Figure 2: Continuum Fit to the Data

The -1 HEG order is in the upper plot and the +1 HEG order is in the lower plot. The plots are displayed with a logarithmic scale and have been filtered to the range 5.0-7.0 keV.


Adding Lines to the Model and Fitting Again

We start by creating a Gaussian model component for the feature at 6.4 keV, named g1 and set the Gaussian normalization to 10-4.

The create_model_component command is used to establish the model component before setting the complex source expressions.

sherpa> create_model_component("xsgaussian","g1")
           <XSgaussian model instance 'xsgaussian.g1'>
sherpa> g1.norm=1.e-4

We create a model with a power-law continuum plus the Gaussian line components, examine the source definition for one of the data sets (they are identical), then fit the data.

sherpa> set_source(1,powr+g1)
sherpa> set_source(2,powr+g1)

sherpa> show_source(1)
Model: 1
(xspowerlaw.powr + xsgaussian.g1)
   Param         Type          Value          Min          Max      Units
   -----         ----          -----          ---          ---      -----
   powr.phoindex thawed     -1.36276           -2            9           
   powr.norm     thawed  1.98057e-05            0        1e+24           
   g1.linee      thawed          6.5            0        1e+06        keV
   g1.sigma      thawed          0.1            0           10        keV
   g1.norm       thawed       0.0001            0        1e+24           

sherpa> fit()
WARNING: data set 1 has associated backgrounds, but they have not been subtracted, nor have background models been set
WARNING: data set 2 has associated backgrounds, but they have not been subtracted, nor have background models been set
Datasets              = 1, 2
Method                = levmar
Statistic             = chi2datavar
Initial fit statistic = 366.776
Final fit statistic   = 53.3692 at function evaluation 70
Data points           = 56
Degrees of freedom    = 51
Probability [Q-value] = 0.383285
Reduced statistic     = 1.04645
Change in statistic   = 313.407
   powr.PhoIndex   -0.912574    +/- 0.35764     
   powr.norm      4.18777e-05  +/- 2.64444e-05 
   g1.LineE       6.3945       +/- 0.00115733  
   g1.Sigma       0.00956428   +/- 0.00172542  
   g1.norm        0.000184658  +/- 1.1385e-05  

Figure 3: Continuum Fit to the Data

Figure 3: Continuum Fit to the Data

The -1 HEG order is in the upper plot and the +1 HEG order is in the lower plot. The fit is overlaid in red. The plots are displayed with a logarithmic scale and have been filtered to the range 5.0-7.0 keV.


Examine the Fit Results

Plotting out our results (Figure 3), we see that we have a good fit to the spectra, including a line component. We can now explore the errors in the line parameters using the conf() command to run the Sherpa confidence routines. As is common in X-ray astronomy, we search for the 90% confidence level values for one interesting parameter (i.e., the variation of the parameter of interest that when refitting the remaining unfrozen parameters yields a change in the χ2 value of 2.706). This corresponds to an error bar of 1.64σ (1σ being the 68% confidence interval), and set this as the default error value in the confidence routines via the set_conf_opt command. We then run conf for the three parameters of each of the Gaussian model components: line energy, line width, and line normalization.

sherpa> set_conf_opt("sigma",1.64)

sherpa> conf(g1.linee,g1.sigma,g1.norm)
WARNING: data set 1 has associated backgrounds, but they have not been subtracted, nor have background models been set
WARNING: data set 2 has associated backgrounds, but they have not been subtracted, nor have background models been set
g1.LineE lower bound:	-0.00189519
g1.norm lower bound:	-1.86834e-05
g1.Sigma lower bound:	-0.00318074
g1.norm upper bound:	1.86903e-05
g1.LineE upper bound:	0.00189851
g1.Sigma upper bound:	0.00283821
Datasets              = 1, 2
Confidence Method     = confidence
Iterative Fit Method  = None
Fitting Method        = levmar
Statistic             = chi2datavar
confidence 1.64-sigma (89.8995%) bounds:
   Param            Best-Fit  Lower Bound  Upper Bound
   -----            --------  -----------  -----------
   g1.LineE           6.3945  -0.00189519   0.00189851
   g1.Sigma       0.00956484  -0.00318074   0.00283821
   g1.norm       0.000184658 -1.86834e-05  1.86903e-05

We see that the line position is determined to approximately 0.002 keV, i.e., an equivalent velocity of 80 km/s. This is essentially the limit of the HEG spectral resolution. The line position is also 0.005 keV lower in energy than expected, equivalent to a 230 km/sec redshift.

The line width is <0.012 keV (i.e., <560 km/sec), but is likely >0.006 keV (i.e., >280 km/sec). Again, these are near the limits of the capabilities of the HEG detector, but far better than the limits one could place with CCD detectors.

Finally, we see that the line normalization is determined to within +/-10%.

We can now plot the data, fits, and residuals (i.e., four separate plots), all on one figure. We select the first and second plots (i.e., the data and fits, plots "plot1" and "plot2") to have logarithmic y-axes. We then print the resulting figure (Figure 4) to an EPS file, sherpa_one_line.eps:

sherpa> plot("fit",1,"fit",2,"delchi",1,"delchi",2) 
sherpa> current_plot("plot1")                           
sherpa> log_scale(Y_AXIS) 
sherpa> current_plot("plot2")
sherpa> log_scale(Y_AXIS)
sherpa> print_window("sherpa_one_line",["format","eps"])

Figure 4: Fit and Residuals to the Spectral Line

Figure 4: Fit and Residuals to the Spectral Line

The HEG -1 order (left) and +1-order (right) displayed with a logarithmic scale on the data plots and linear scale on the residual plots.

To explore whether adding another, broader line component changes the results any, we create a second Gaussian component, and call it g2. We then define a model consisting of a power-law plus two Gaussian components, and apply it to both data sets. We change the normalization of the second Gaussian to be somewhat smaller, then fit the data.

sherpa> create_model_component("xsgaussian","g2")
           <XSgaussian model instance 'xsgaussian.g2'>
sherpa> g2.norm=1.e-4
sherpa> set_model(1,powr+g1+g2)
sherpa> set_model(2,powr+g1+g2)

sherpa> fit()
WARNING: data set 1 has associated backgrounds, but they have not been subtracted, nor have background models been set
WARNING: data set 2 has associated backgrounds, but they have not been subtracted, nor have background models been set
Datasets              = 1, 2
Method                = levmar
Statistic             = chi2datavar
Initial fit statistic = 157.78
Final fit statistic   = 47.2091 at function evaluation 434
Data points           = 56
Degrees of freedom    = 48
Probability [Q-value] = 0.505171
Reduced statistic     = 0.983524
Change in statistic   = 110.569
   powr.PhoIndex   1.70777      +/- 2.38088     
   powr.norm      0.00327119   +/- 0.0127288   
   g1.LineE       6.39445      +/- 0.00116623  
   g1.Sigma       0.00889812   +/- 0.0018528   
   g1.norm        0.000180085  +/- 1.1658e-05  
   g2.LineE       6.59739      +/- 0.121619    
   g2.Sigma       0.43752      +/- 0.225138    
   g2.norm        0.000152171  +/- 0.000148234 

sherpa> conf()
WARNING: data set 1 has associated backgrounds, but they have not been subtracted, nor have background models been set
WARNING: data set 2 has associated backgrounds, but they have not been subtracted, nor have background models been set
g1.LineE lower bound:	-0.00190849
g1.norm lower bound:	-1.92304e-05
g1.LineE upper bound:	0.00191562
g1.Sigma lower bound:	-0.00346968
powr.PhoIndex lower bound:	-3.06562
g1.norm upper bound:	1.92329e-05
g2.Sigma -: WARNING: The confidence level lies within (1.092544e-01, 1.123116e-01)
g2.Sigma lower bound:	-0.326737
g1.Sigma upper bound:	0.00295318
powr.PhoIndex upper bound:	13.0699
g2.LineE -: WARNING: The confidence level lies within (6.311673e+00, 6.316105e+00)
g2.LineE lower bound:	-0.283498
powr.norm -: WARNING: The confidence level lies within (5.090753e-05, 5.232004e-05)
powr.norm lower bound:	-0.00321957
g2.Sigma upper bound:	1.27247
g2.LineE upper bound:	0.873918
powr.norm +: WARNING: The confidence level lies within (7.242173e-02, 7.241075e-02)
powr.norm upper bound:	0.0691451
g2.norm -: WARNING: The confidence level lies within (1.880199e-05, 1.873642e-05)
g2.norm lower bound:	-0.000133401
g2.norm upper bound:	0.000982883
Datasets              = 1, 2
Confidence Method     = confidence
Iterative Fit Method  = None
Fitting Method        = levmar
Statistic             = chi2datavar
confidence 1.64-sigma (89.8995%) bounds:
   Param            Best-Fit  Lower Bound  Upper Bound
   -----            --------  -----------  -----------
   powr.PhoIndex      1.70777     -3.06562      13.0699
   powr.norm      0.00327119  -0.00321957    0.0691451
   g1.LineE          6.39445  -0.00190849   0.00191562
   g1.Sigma       0.00889812  -0.00346968   0.00295318
   g1.norm       0.000180085 -1.92304e-05  1.92329e-05
   g2.LineE          6.59739    -0.283498     0.873918
   g2.Sigma          0.43752    -0.326737      1.27247
   g2.norm       0.000152171 -0.000133401  0.000982883

The inclusion of the second line component does not greatly affect the conclusions about the parameters of the first line component. The line centroid is the same with comparable error bars. The uncertainty on the line width has increased slightly, as has the uncertainty on the line normalization.

We plot the results for this fit in a manner identical to the one Gaussian line fit (Figure 5).

sherpa> plot("fit",1,"fit",2,"delchi",1,"delchi",2) 
sherpa> current_plot("plot1")                           
sherpa> log_scale(Y_AXIS) 
sherpa> current_plot("plot2")
sherpa> log_scale(Y_AXIS)
sherpa> print_window("sherpa_two_line",["format","eps"])

Figure 5: HEG +/-1 Order Two Gaussian Component Fit

Figure 5: HEG +/-1 Order Two Gaussian Component Fit

The HEG -1-order (left) and +1-order (right) displayed with a logarithmic-scale on the data plots and linear-scale on the residual plots.

Given that the extra component adds little to the description of the data, we drop it from the model and refit the data with a single line, resulting in the same parameters from before.

sherpa> set_model(1,powr+g1)
sherpa> set_model(2,powr+g1)
sherpa> fit()
WARNING: data set 1 has associated backgrounds, but they have not been subtracted, nor have background models been set
WARNING: data set 2 has associated backgrounds, but they have not been subtracted, nor have background models been set
Datasets              = 1, 2
Method                = levmar
Statistic             = chi2datavar
Initial fit statistic = 121.647
Final fit statistic   = 53.3692 at function evaluation 99
Data points           = 56
Degrees of freedom    = 51
Probability [Q-value] = 0.383285
Reduced statistic     = 1.04645
Change in statistic   = 68.3405
   powr.PhoIndex   -0.912574    +/- 0.35764     
   powr.norm      4.18778e-05  +/- 2.64445e-05 
   g1.LineE       6.3945       +/- 0.00115733  
   g1.Sigma       0.00956428   +/- 0.00172542  
   g1.norm        0.000184658  +/- 1.1385e-05  

Calculate the Equivalent Width of a Line and Continuum Flux

Rather than merely presenting the line normalizations, we can instead calculate the line equivalent width. The line equivalent width is an often used measure for the line strength relative to the underlying continuum. Sherpa has a function for evaluating the equivalent width, eqwidth. The required inputs are the source model absent the line feature of interest and the source model with the line feature of interest.

Optionally, the data set number for which the equivalent width will be calculated can also be given as an input, but here as we apply the same model and parameters to both data sets, this is not necessary. There are also optional lo and hi arguments for restricting the calculation of the equivalent width to a subset of the full energy or wavelength range of a data set.

sherpa> eqwidth(powr,powr+g1)
           0.811011851650702

The units of the equivalent width are the same as the units of the x-axis, in this case, keV. Thus the line equivalent width is ~811 eV, which indicates a very strong line. That is, there is as nearly as much flux in the line as there is in the continuum over the 6.0-7.0 keV region.

We can verify this by calculating fluxes directly, using the calc_photon_flux and calc_energy_flux commands. We first calculate these fluxes from the full model (continuum plus line) over just the 6.0-7.0 keV range. The returned fluxes are in units of photons/cm2/sec and ergs/cm2/sec, respectively.

sherpa> calc_photon_flux(6.,7.)
           0.00041605434631919572

sherpa> calc_energy_flux(6.,7.)
           4.3076192427398961e-12

Deleting the power-law continuum component, we calculate the flux values for just the line component. We do this by changing the model with set_model to only the Gaussian component without refitting.

sherpa> set_model(1,g1)
sherpa> set_model(2,g1)
sherpa> calc_photon_flux(6.,7.)
           0.00018465793791707486
sherpa> calc_energy_flux(6.,7.)
           1.8918448000041687e-12

Over this one keV bandwidth, the line contributes a fraction, ~1.84/(4.16-1.84)~0.79, to the photon flux. This is approximately consistent with the equivalent width calculation of 811 eV.


Scripting It

The file fit.py is a Python script which performs the primary commands used above; it can be executed by typing exec(open("fit.py").read()) on the Sherpa command line.

The Sherpa script command may be used to save everything typed on the command line in a Sherpa session:

sherpa> script(filename="sherpa.log", clobber=False)

(Note that restoring a Sherpa session from such a file could be problematic since it may include syntax errors, unwanted fitting trials, et cetera.)

The CXC is committed to helping Sherpa users transition to new syntax as smoothly as possible. If you have existing Sherpa scripts or save files, submit them to us via the CXC Helpdesk and we will provide the CIAO/Sherpa 4.11 syntax to you.


Summary

As shown above, working with gratings spectra is fundamentally no different than working with CCD spectra. Aside from the fact that gratings data will often begin with a PHA2 file containing multiple spectra, rather than a type 1 PHA file containing a single spectrum, the analysis paths can be very much the same for both types of spectra. (If one is performing analysis in either Sherpa or ISIS, the PHA2 file can be used directly, whereas for XSPEC one has to use the intermediate step of creating type 1 PHA files in order to be able to bin the data.) The data are read in, response matrices are assigned, models are defined and fit, and error bars on the parameters are determined. Line fluxes and equivalent widths can be calculated easily. Gratings spectra can even be slightly less complicated than CCD spectra in that ACIS-S/HETG spectra often do not require a background. (Again, order sorting of the spectra often ensures that the background is negligible.)


History

01 Apr 2009 new for Sherpa 4.1
11 Jan 2010 updated for CIAO 4.2
13 Jul 2010 updated for CIAO 4.2 Sherpa v2: removal of S-Lang version of thread.
22 Jan 2012 reviewed for CIAO 4.4 (no changes)
13 Dec 2012 updated for CIAO 4.5: group commands no longer clear the existing data filter
10 Dec 2013 updated for CIAO 4.6: updated fit results and screen output
27 Feb 2015 updated for CIAO 4.7, no content change
14 Dec 2015 updated for CIAO 4.8, no content change
10 Nov 2016 updated for CIAO 4.9, no content change, updated fit results.
30 May 2018 updated for CIAO 4.10, no content change
11 Dec 2018 reviewed for CIAO 4.11, revised screen output