Simulating IXO/CAT Spectra with Sherpa

Sherpa Threads (CIAO 4.5 Sherpa v1)

Overview

Synopsis:

This thread illustrates the use of the Sherpa fake_pha command to simulate the spectrum of a point source that would be produced with the CAT spectrometer aboard the International X-Ray Observatory (IXO), based on current IXO/CAT response information. IXO, scheduled to launch in 2021, represents the next-generation X-ray mission, as its instrumental capabilities will surpass those of the current X-ray observatories, namely Chandra, XMM, RXTE, and Suzaku. Here, we use Sherpa to simulate several IXO/CAT grating spectra, using the appropriate IXO response files and source model expression.

Last Update: 15 Dec 2011 - reviewed for CIAO 4.4: a work-around for a save_pha bug was added

Getting started - downloading calibration response files for simulations

In order to simulate an IXO X-ray spectrum with Sherpa, we must define an instrument response with the appropriate IXO "rsp" response file(s) ("rsp" implies that both the RMF and ARF are included). IXO grating and no-grating response files are available for download on the IXO/CAT: Response Files page of the MIT IXO/CAT Development website, maintained by David P. Huenemoerder; see also the Response Matrices page of the Goddard Space Flight Center website.

In this thread, we simulate the zeroth through fifth order IXO/CAT grating responses using the files CAT_D54_m0.rsp, CAT_D54_R3000_Lm1.rsp, CAT_D54_R3000_Lm2.rsp, CAT_D54_R3000_Lm3.rsp, CAT_D54_R3000_Lm4.rsp, and CAT_D54_R3000_Lm5.rsp (54 degrees of sectors tiled per side, resolving power 3000, Lorentzian instrumental profile, effective area ~3000 cm²). The zeroth order response includes the direct beam from the open inner annulus of the grating array; this portion falls on the microcalorimeter detector. Note that a full IXO/CAT instrumental configuration involves all five (1-5) diffraction responses for grating spectroscopy, and the single no-grating response (not used in this thread) for imaging spectroscopy.

Figure 1: Effective Area Comparison

[Print media version: Effective Area Comparison]

Figure 1: Effective Area Comparison

Plot of the effective areas of the IXO/CAT high-resolution spectrometer compared to those of current X-ray observatories* (Chandra HEG and MEG and XMM RGS).

*http://space.mit.edu/home/dph/ixo/comparison.html

We will use the Sherpa fake_pha command to simulate IXO/CAT 1-D spectral data sets based on the defined instrument responses and source model expression. This command has four required arguments: data set ID, ARF, RMF, and exposure time. For details on the functionality of fake_pha, with other examples of its usage, see the other Sherpa simulation threads.

Establishing the Source Model Expression for the Simulation

In this thread, we simulate the IXO/CAT zeroth, first, second, third, fourth, and fifth order line spectra of the coronally active binary star system (RS CVn) Sigma Geminorum ("Sigma Gem"), using a source model expression consisting of a two-temperature APEC thermal plasma model (xsvapec).

The source model chosen for the simulation is defined with the Sherpa set_source command, as follows, where it is assigned to six data set IDs, one for each of the IXO/CAT grating orders to be simulated:

set_source("sim_0", xsvapec.apec1 + xsvapec.apec2)
set_source("sim_1", apec1 + apec2)
set_source("sim_2", apec1 + apec2)
set_source("sim_3", apec1 + apec2)
set_source("sim_4", apec1 + apec2)
set_source("sim_5", apec1 + apec2)

set_par(apec1.norm, 0.025, frozen=True)
set_par(apec2.norm, 0.095, frozen=True)

set_par(apec1.kT, 0.78, frozen=True)
set_par(apec2.kT, 3.01, frozen=True)

set_par(apec1.O, 0.5, frozen=True)
set_par(apec2.O, 0.5, frozen=True)

set_par(apec1.Ne, 1.5, frozen=True)
set_par(apec2.Ne, 1.5, frozen=True)

set_par(apec1.Mg, 0.4, frozen=True)
set_par(apec2.Mg, 0.4, frozen=True)

set_par(apec1.Si, 0.4, frozen=True)
set_par(apec2.Si, 0.4, frozen=True)

set_par(apec1.S, 0.2, frozen=True)
set_par(apec2.S, 0.2, frozen=True)

set_par(apec1.Fe, 0.3, frozen=True)
set_par(apec2.Fe, 0.3, frozen=True)

This model gives a good approximation to the Chandra ACIS-S/HETG spectrum of Sigma Gem in reproducing the continuum and general distribution of lines. In the section "Fitting the Simulated Data", the simulated spectrum is fit with a polynomial continuum model plus several Gaussian line profiles for calculating equivalent widths and line fluxes.

Defining the Instrument Response for the Simulation

Now that we have defined a source model for the simulations, the only thing left to do to begin the simulation process is load the IXO/CAT instrument responses. For this step, we use the zeroth through fifth order response files downloaded from the MIT IXO/CAT development website:

sherpa> ixo_rsp0=unpack_rmf("CAT_D54_m0.rsp")
sherpa> ixo_rsp1=unpack_rmf("CAT_D54_R3000_Lm1.rsp") 
sherpa> ixo_rsp2=unpack_rmf("CAT_D54_R3000_Lm2.rsp")
sherpa> ixo_rsp3=unpack_rmf("CAT_D54_R3000_Lm3.rsp")
sherpa> ixo_rsp4=unpack_rmf("CAT_D54_R3000_Lm4.rsp")
sherpa> ixo_rsp5=unpack_rmf("CAT_D54_R3000_Lm5.rsp")

The response data is loaded into variables 'ixo_rsp*' with the unpack_rmf command. These variables will be used to assign the instrument response for each grating order to the corresponding faked data set in the next section of this thread, "Running the Simulation with fake_pha." (Note that the effective area is included in each response; since the fake_pha command requires both an 'arf' and 'rmf' input, we will specify later that 'arf=None' and 'rmf=ixo_rsp'.)

Running the Simulation with fake_pha

Simulating the zeroth through fifth IXO/CAT grating spectra with Sherpa involves convolving the chosen source model expression with the corresponding IXO/CAT responses, and applying Poisson noise to the counts predicted by the model. These steps are run by the fake_pha command, which has four required arguments: data set ID, ARF, RMF, and exposure time. We will assign the IXO response stored in each of the 'ixo_rsp*' variables to the fake_pha 'rmf' argument and set 'arf' to 'None' (since the ARF is included in the response), and the 'exposure' to 100000 s. We will repeat this step six times, once for each spectral order, each time assigning a different grating order response to a different data set ID (but using the same source model for each).

We are now ready to run the simulation for each of the six spectral orders with fake_pha, using the two-temperature APEC thermal plasma source model defined in the previous section:

sherpa> fake_pha("sim_0", arf=None, rmf=ixo_rsp0, exposure=100000)
sherpa> fake_pha("sim_1", arf=None, rmf=ixo_rsp1, exposure=100000)
sherpa> fake_pha("sim_2", arf=None, rmf=ixo_rsp2, exposure=100000)
sherpa> fake_pha("sim_3", arf=None, rmf=ixo_rsp3, exposure=100000)
sherpa> fake_pha("sim_4", arf=None, rmf=ixo_rsp4, exposure=100000)
sherpa> fake_pha("sim_5", arf=None, rmf=ixo_rsp5, exposure=100000)

These commands associate the data set IDs "sim_0", "sim_1", "sim_2", "sim_3", "sim_4", and "sim_5" with simulated IXO/CAT zeroth, first, second, third, fourth, and fifth order line spectra, based on the assumed exposure time, IXO/CAT instrument responses, and source model expression. Poisson noise is added to the modeled data.

We may use the show_data command to inspect the basic properties of the new data sets; observe that the "name" field reads "faked", and the "exposure" field has been set to the chosen exposure time:

sherpa> show_data()
Data Set: sim_0
Filter: 1-22000 Channel
Noticed Channels: 1-22000
name           = faked
channel        = Int32[22000]
counts         = Float64[22000]
staterror      = None
syserror       = None
bin_lo         = None
bin_hi         = None
grouping       = None
quality        = None
exposure       = 100000
backscal       = None
areascal       = None
grouped        = False
subtracted     = False
units          = channel
rate           = True
plot_fac       = 0
response_ids   = [1]
background_ids = []

RMF Data Set: sim_0:1
header   = 
 CHANTYPE      = PHA
 CREATOR       = marfrmf 3.2.6
 EFFAREA       = 1.0
 FILTER        = 
 HDUCLAS1      = RESPONSE
 HDUCLAS2      = EBOUNDS
 HDUCLAS3      = FULL
 HDUCLASS      = OGIP
 HDUDOC        = OGIP memos CAL/GEN/92-002 & 92-002a
 HDUVERS       = 1.3.0
 HDUVERS1      = 1.0.0
 HDUVERS2      = 1.3.0
 INSTRUME      = XMS
 LO_THRES      = 0.0
 NUMELT        = 321425
 NUMGRP        = 22000
 RMFVERSN      = 1992a
 TELESCOP      = ConX
name     = CAT_D54_m0.rsp
detchans = 22000
energ_lo = Float64[22000]
energ_hi = Float64[22000]
n_grp    = Int16[22000]
f_chan   = UInt32[22000]
n_chan   = UInt32[22000]
matrix   = Float64[321425]
offset   = 1
e_min    = Float64[22000]
e_max    = Float64[22000]

[output for the rest of the data sets was omitted for brevity]

The simulated data sets may be visualized with the plot command, as shown below. To plot in wavelength space, we use the set_analysis command. To save the plot, we use print_window.

sherpa> set_analysis("wave")

sherpa> prefs = get_data_plot_prefs()
sherpa> prefs["linestyle"] = chips_solid
sherpa> prefs["symbolstyle"] = chips_none

sherpa> plot("data", "sim_0", "data", "sim_1", "data", "sim_2", "data","sim_3", "data", "sim_4", "data", "sim_5")
sherpa> split(6)
sherpa> current_plot("all")
sherpa> linear_scale()

sherpa> current_plot("plot1")
sherpa> set_plot_xlabel("")
sherpa> set_plot_ylabel("")
sherpa> set_plot_title("Simulated IXO/CAT grating orders 0-5")
sherpa> add_label(35., 700, "CAT 0th order + direct beam")
sherpa> set_label(["color","green"])

sherpa> current_plot("plot2")
sherpa> set_plot_xlabel("")
sherpa> set_plot_ylabel("")
sherpa> set_plot_title("")
sherpa> add_label(66., 16, "CAT 1st order")
sherpa> set_label(["color","green"])

sherpa> current_plot("plot3")
sherpa> set_plot_xlabel("")
sherpa> set_plot_title("")
sherpa> add_label(33., 40, "CAT 2nd order")
sherpa> set_label(["color","green"])

sherpa> current_plot("plot4")
sherpa> set_plot_xlabel("")
sherpa> set_plot_ylabel("")
sherpa> set_plot_title("")
sherpa> add_label(22., 190, "CAT 3rd order")
sherpa> set_label(["color","green"])

sherpa> current_plot("plot5")
sherpa> set_plot_xlabel("")
sherpa> set_plot_ylabel("")
sherpa> set_plot_title("")
sherpa> add_label(16.5, 145, "CAT 4th order")
sherpa> set_label(["color","green"])

sherpa> current_plot("plot6")
sherpa> set_plot_ylabel("")
sherpa> set_plot_title("")
sherpa> add_label(13.3, 75, "CAT 5th order")
sherpa> set_label(["color","green"])

sherpa> print_window("ixo_sim_data.ps")

The resulting plot is shown Figure 2.

Figure 2: Plot of simulated IXO/CAT source spectra

[Plot of simulated IXO/CAT 0th through 5th order grating spectra]

[Print media version: Plot of simulated IXO/CAT 0th through 5th order grating spectra]

Figure 2: Plot of simulated IXO/CAT source spectra

Simulated IXO/CAT grating spectra of coronally active binary Sigma Gem. The zeroth order response includes the direct beam from the open inner annulus of the grating array; this portion falls on the microcalorimeter detector. The effective area of each of the CAT orders (1-5) is ~3000 cm² and the resolving power is ~3000 in the 0.3-1 keV (12.4-41.3 Angstroms) range, and the combined effective area is ~3000 cm² in the 1-2 keV (6.2-12.4 Angstroms) range. The zeroth+direct beam effective area lies well above 3000 cm² for most of the 0.3-10 keV (1.24-41.3 Angstroms) range (see the "Project Requirements" plots on the MIT IXO/CAT: Performance page for details).

Note that had we not applied a customized, properly normalized source model expression to the faked data, we would have had to re-normalize the simulated data to match the known flux (or counts) of Sigma Gem. This process is illustrated in the threads "Simulating Chandra ACIS-S Spectra with Sherpa" and "Simulating 1-D Data: the Sherpa fake_pha Command."

Writing the Simulated Data to Output Files

We may use the save_pha command to write each simulated data set to a PHA file, with a header containing the exposure time value and path to the corresponding response file:

# NOTE: the 'get_data' lines below are necessary to work around a CIAO
# 4.4 bug with the save_pha() function; normally, one would
# just need to use save_pha() to write a simulated PHA data set to file.

sherpa> get_data("sim_0").header = {}            # Make header an empty dictionary
sherpa> get_data("sim_0").name = "ixo_sim0.pha"  # Update name for save_pha

sherpa> get_data("sim_1").header = {}               
sherpa> get_data("sim_1").name = "ixo_sim1.pha"  

sherpa> get_data("sim_2").header = {}              
sherpa> get_data("sim_2").name = "ixo_sim2.pha"  

sherpa> get_data("sim_3").header = {}              
sherpa> get_data("sim_3").name = "ixo_sim3.pha"  

sherpa> get_data("sim_4").header = {}              
sherpa> get_data("sim_4").name = "ixo_sim4.pha"  

sherpa> get_data("sim_5").header = {}              
sherpa> get_data("sim_5").name = "ixo_sim5.pha"  

sherpa> save_pha("sim_0", "ixo_sim0.pha")
sherpa> save_pha("sim_1", "ixo_sim1.pha")
sherpa> save_pha("sim_2", "ixo_sim2.pha")
sherpa> save_pha("sim_3", "ixo_sim3.pha")
sherpa> save_pha("sim_4", "ixo_sim4.pha")
sherpa> save_pha("sim_5", "ixo_sim5.pha")

sherpa> !dmlist ixo_sim0.pha header | grep EXPOSURE
0002 EXPOSURE             100000                         Int4 

sherpa> !dmlist ixo_sim0.pha header | grep RESPFILE
0001 RESPFILE             CAT_D54_m0.rsp                 String

Fitting the Simulated Data

The simulated data sets may be filtered and fit as any other data set in Sherpa. For example, we can filter the simulated X-ray spectra of Sigma Gem to include only the counts in a region which includes emission lines, e.g., to fit line profiles and determine line fluxes:

sherpa> print get_analysis()
wavelength

sherpa> ignore()
sherpa> notice(17.5, 18.83)

sherpa> show_filter()
Data Set Filter: sim_5
15.2848 Wavelength (Angstrom)

Data Set Filter: sim_4
17.4997-18.8300 Wavelength (Angstrom)

Data Set Filter: sim_1
24.0016 Wavelength (Angstrom)

Data Set Filter: sim_0
17.5009-18.8331 Wavelength (Angstrom)

Data Set Filter: sim_3
17.4997-18.8297 Wavelength (Angstrom)

Data Set Filter: sim_2
17.4997-18.8308 Wavelength (Angstrom)

Then, we can fit the simulated data with a source model expression. In this example, we choose to fit orders 0, 2, 3, and 4 with a polynomial plus several Gaussians to model the continuum and emission lines, respectively.

sherpa> create_model_component("polynom1d","poly")
sherpa> set_par(poly.a0, 0.002, 0.0, 0.04, frozen=False)

sherpa> create_model_component("gauss1d", "g1")
sherpa> create_model_component("gauss1d", "g2")
sherpa> create_model_component("gauss1d", "g3")
sherpa> create_model_component("gauss1d", "g4")
sherpa> create_model_component("gauss1d", "g5")

sherpa> set_par(g1.pos, 17.6, 17.55, 17.65, frozen=False)
sherpa> set_par(g2.pos, 18.63, 18.58, 18.68, frozen=False)
sherpa> set_par(g3.pos, 18.69, 18.64, 18.74, frozen=False)
sherpa> set_par(g4.pos, 18.73, 18.68, 18.78,  frozen=False)
sherpa> set_par(g5.pos, 17.75, 17.70, 17.80,  frozen=False)

sherpa> for num in range(1, 5):
            set_par(g%i.fwhm %num, 1.e-4, frozen=True) 
            set_par(g%i.ampl %num, 4., frozen=False)

The individual polynomial and Gaussian model components are now combined into a single, multi-component source model using the set_source command. (Note: we chose to define the model components separately first with the create_model_component command, and then combine them into a single, multi-component model with set_source; the initial step is not necessary but is shown for clarity, i.e., you can just use set_source to define a multi-component source model expression):

sherpa> set_source("sim_0", poly + g1 + g2 + g3 + g4 + g5)
sherpa> set_source("sim_2", poly + g1 + g2 + g3 + g4 + g5)
sherpa> set_source("sim_3", poly + g1 + g2 + g3 + g4 + g5)
sherpa> set_source("sim_4", poly + g1 + g2 + g3 + g4 + g5)

Here we have used the set_source command to redefine the source model associated with the faked data sets. Initially, we used a two-temperature APEC thermal plasma model to produce the simulated spectra, now we choose to fit the resulting spectra with a polynomial plus five Gaussians to model the continuum and emission lines between 17.5 and 18.83 Angstroms in the spectrum of Sigma Gem.

We may now fit the model to the simulated spectra using the Neldermead-Simplex optimization method and Sherpa's default fit statistic, Chi2Gehrels, and then plot the fit results with the plot command.

sherpa> show_method()

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

sherpa> show_stat()

Statistic: Chi2Gehrels


sherpa> set_method("simplex")

sherpa> fit("sim_0", "sim_2", "sim_3", "sim_4")
Datasets              = 'sim_0', 'sim_2', 'sim_3', 'sim_4'
Method                = neldermead
Statistic             = chi2gehrels
Initial fit statistic = 22568.1
Final fit statistic   = 4396.81 at function evaluation 7291
Data points           = 3864
Degrees of freedom    = 3853
Probability [Q-value] = 1.48859e-09
Reduced statistic     = 1.14114
Change in statistic   = 18171.3
   poly.c0        0.00207708
   g1.pos         17.6228
   g1.ampl        0.762622
   g2.pos         18.6257
   g2.ampl        0.0920954
   g3.pos         18.6909
   g3.ampl        0.203393
   g4.pos         18.732
   g4.ampl        0.0939529
   g5.pos         17.7317
   g5.ampl        0.0816013


sherpa> plot("fit", "sim_0", "fit", "sim_2", "fit","sim_3", "fit", "sim_4")
sherpa> split(4)
sherpa> current_plot("all")
sherpa> linear_scale()

sherpa> current_plot("plot1")
sherpa> set_plot_xlabel("")
sherpa> set_plot_ylabel("")
sherpa> set_plot_title("")
sherpa> add_label(18., 40., "CAT 0th order + direct beam")
sherpa> set_label(["color","green"])

sherpa> current_plot("plot2")
sherpa> set_plot_xlabel("")
sherpa> set_plot_ylabel("")
sherpa> set_plot_title("")
sherpa> add_label(18., 3., "CAT 2nd order")
sherpa> set_label(["color","green"])

sherpa> current_plot("plot3")
sherpa> set_plot_xlabel("")
sherpa> set_plot_title("")
sherpa> add_label(18., 15., "CAT 3rd order")
sherpa> set_label(["color","green"])

sherpa> current_plot("plot4")
sherpa> set_plot_ylabel("")
sherpa> set_plot_title("")
sherpa> add_label(18., 4., "CAT 4th order")
sherpa> set_label(["color","green"])

sherpa> print_window("ixo_sim_fits.ps")

The resulting plot is shown in Figure 3.

Figure 3: Plot of fit to simulated IXO/CAT source spectra

[Print media version: Plot of fit to simulated IXO spectrum]

Figure 3: Plot of fit to simulated IXO/CAT source spectra

Plot of fit to simulated IXO/CAT grating spectra of coronally active binary Sigma Gem, using a polynomial plus Gaussian source model.

Examining the Fit Results

Plotting our results (Figure 3), we see that we have a good fit to the spectra, including the line components. We can now explore the errors in the model parameters using the confidence command; the 68% confidence level values are returned by default (1-sigma), and may be changed with the set_conf_opt command. In this example, we run the confidence routines on the four data sets simultaneously, for all thawed model parameters (note that this type of confidence fit can take time).

sherpa> conf("sim_0","sim_2","sim_3","sim_4")
WARNING: New minimum statistic found while computing confidence limits
WARNING: New best-fit parameters:
Method                = neldermead
Statistic             = chi2gehrels
Initial fit statistic = 4396.81
Final fit statistic   = 4126.95 at function evaluation 5578
Data points           = 3864
Degrees of freedom    = 3853
Probability [Q-value] = 0.0011167
Reduced statistic     = 1.0711
Change in statistic   = 269.859
   poly.c0        0.00206652
   g1.pos         17.6226
   g1.ampl        0.763434
   g2.pos         18.6265
   g2.ampl        0.0989297
   g3.pos         18.6911
   g3.ampl        0.23027
   g4.pos         18.732
   g4.ampl        0.117179
   g5.pos         17.7317
   g5.ampl        0.0764028

g4.ampl lower bound:    -0.00461964
g5.ampl lower bound:    -0.0023044
g3.ampl lower bound:    -0.00464142
g2.pos lower bound:     -7.14814e-05
g5.pos -: WARNING: The confidence level lies within (1.773166e+01, 1.773168e+01)
g5.pos lower bound:     -3.4881e-05
g1.pos -: WARNING: The confidence level lies within (1.762250e+01, 1.762253e+01)
g1.pos lower bound:     -3.87605e-05
g2.ampl -: WARNING: The confidence level lies within (9.388294e-02, 9.640631e-02)
g2.ampl lower bound:    -0.00378506
g1.ampl -: WARNING: The confidence level lies within (7.570599e-01, 7.602470e-01)
g1.ampl lower bound:    -0.00478058
poly.c0 -: WARNING: The confidence level lies within (2.064413e-03, 2.065469e-03)
poly.c0 lower bound:    -1.58289e-06
g3.pos -: WARNING: The confidence level lies within (1.868864e+01, 1.869115e+01)
g3.pos lower bound:     -0.00125111
g2.pos upper bound:     3.31521e-05
g4.pos -: WARNING: The confidence level lies within (1.873124e+01, 1.873157e+01)
g4.pos lower bound:     -0.00061419

g1.pos upper bound:     2.90703e-05

g5.ampl upper bound:    0.00418329
g5.pos upper bound:     9.61042e-05
g3.ampl upper bound:    0.00334795
g2.ampl upper bound:    0.00433705
g4.pos upper bound:     0.00010125
g4.ampl upper bound:    0.00276796

poly.c0 +: WARNING: The confidence level lies within (2.067580e-03, 2.067579e-03)
poly.c0 upper bound:    1.05578e-06
g3.pos +: WARNING: The confidence level lies within (1.869121e+01, 1.869121e+01)
g3.pos upper bound:     6.59767e-05
g1.ampl +: WARNING: The confidence level lies within (7.684760e-01, 7.684729e-01)
g1.ampl upper bound:    0.00504046
Datasets              = 'sim_0', 'sim_2', 'sim_3', 'sim_4'
Confidence Method     = confidence
Fitting Method        = neldermead
Statistic             = chi2gehrels
confidence 1-sigma (68.2689%) bounds:
   Param            Best-Fit  Lower Bound  Upper Bound
   -----            --------  -----------  -----------
   poly.c0        0.00206652 -1.58289e-06  1.05578e-06
   g1.pos            17.6226 -3.87605e-05  2.90703e-05
   g1.ampl          0.763434  -0.00478058   0.00504046
   g2.pos            18.6265 -7.14814e-05  3.31521e-05
   g2.ampl         0.0989297  -0.00378506   0.00433705
   g3.pos            18.6911  -0.00125111  6.59767e-05
   g3.ampl           0.23027  -0.00464142   0.00334795
   g4.pos             18.732  -0.00061419   0.00010125
   g4.ampl          0.117179  -0.00461964   0.00276796
   g5.pos            17.7317  -3.4881e-05  9.61042e-05
   g5.ampl         0.0764028   -0.0023044   0.00418329

We can now visualize the data, fit, and residuals with the plot_fit_delchi command. Here, we examine the fit and delta chi residuals for the 3rd spectral order:

sherpa> plot_fit_delchi("sim_3")
sherpa> set_plot_title("Model fit to IXO/CAT 3rd grating order")

The resulting plot of the fit to the simulated source spectrum, along with the residuals divided by the uncertainties, is shown in Figure 4.

Figure 4: Plot of fit to simulated IXO/CAT grating spectrum, with delta chi residuals.

Figure 4: Plot of fit to simulated IXO/CAT grating spectrum, with delta chi residuals.

Delta chi residuals from the fit to the IXO/CAT third order grating spectrum.

The plot may be saved to a PostScript (or other format) file with the print_window command:

sherpa> print_window("ixo_sim_fit_delchi.ps")

Comparison with the Chandra mission

A comparison of the predicted IXO/CAT grating line spectra and the corresponding Chandra ACIS-S/HETG spectrum in Figure 5 reveals some of the scientific gains promised by IXO.

Figure 5: Comparison of predicted IXO/CAT RS CVn binary source spectra and Chandra/HETG source spectrum

The IXO/CAT second, third, and fourth order grating responses (and first and fifth, CAT_D54_R3000_Lm*.rsp) each have an effective area of ~3000 cm² and resolving power of ~3000 in the 0.3-1.0 keV (12.4-41.3 Angstroms) range, and a combined effective area of ~3000 in the 1-2 keV (6.2-12.4 Angstroms) range, while the Chandra ACIS-S/HETG first order response has an effective area of ~100 cm² and resolving power of ~1000 in this range. The IXO/CAT zeroth+direct beam effective area lies well above 3000 cm² for most of the 0.3-10 keV (1.24-41.3 Angstroms) range.

After loading a saved Sherpa session in which the positive first order Chandra ACIS/HETG spectrum of Sigma Gem was fit, we can compare the total counts, line flux, and equivalent width of a selected emission line in the IXO/CAT grating spectra to that in the Chandra/HETG first order spectrum.

The Sherpa eqwidth command evaluates the equivalent width of a line given the source model absent the line feature of interest in the first argument, and the source model with the line feature of interest in the second argument. The 'lo' and 'hi' arguments may be used to restrict the calculation of the equivalent width to a subset of the full energy or wavelength range of a data set.

We can calculate fluxes directly using the calc_energy_flux command (or calc_photon_flux), which returns fluxes in units of ergs/cm2/sec from the full model (continuum plus line) over just the chosen wavelength range (re-defining the source model by removing just the continuum component, without re-fitting, will return the flux from just the line region). The calc_data_sum function returns the total number of counts in the spectral range specified.

sherpa> restore("chandra_session")

sherpa> show_data(1)                          #Chandra ACIS/HETG +1
Filter: 17.4912-17.9560 Wavelength (Angstrom)
Noticed Channels: 1382-1612
{header omitted for brevity}
name           = pha2[TG_M=-1,1]
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       = 62827.7722802
backscal       = 1.0
areascal       = 1.0
grouped        = True
subtracted     = False
units          = wavelength
rate           = True
plot_fac       = 0
response_ids   = [1]
background_ids = [1, 2]


sherpa> calc_data_sum(17.55, 17.65, id=1)        #Chandra ACIS/HETG +1
           31.0

sherpa> calc_data_sum(17.55, 17.65, id="sim_0")  #IXO/CAT
           429140.0

sherpa> calc_data_sum(17.55, 17.65, id="sim_2")
           15200.0

sherpa> calc_data_sum(17.55, 17.65, id="sim_3")
           69053.0

sherpa> calc_data_sum(17.55, 17.65, id="sim_4")
           12959.0



sherpa> calc_energy_flux(17.55, 17.65, id=1)    #Chandra ACIS/HETG +1
        6.0715112709209197e-13

sherpa> eqwidth(poly, poly+g1, lo=17.55, hi=17.65, id=1)
        0.11786604246590549
            

sherpa> calc_energy_flux(17.55, 17.65, id="sim_0")
        3.1542577202444669e-13

sherpa> eqwidth(poly, poly+g1, lo=17.55, hi=17.65, id="sim_0")
        0.039035857783371569


sherpa> calc_energy_flux(17.55, 17.65, id="sim_2")
        3.2779064180803955e-13

sherpa> eqwidth(poly, poly+g1, lo=17.55, hi=17.65, id="sim_2")
        0.039039465746696034


sherpa> calc_energy_flux(17.55, 17.65, id="sim_3")
        3.2654314616391017e-13

sherpa> eqwidth(poly, poly+g1, lo=17.55, hi=17.65, id="sim_3")
        0.039120462115528443

sherpa> calc_energy_flux(17.55, 17.65, id="sim_4")
        3.2591620508079683e-13

sherpa> eqwidth(poly, poly+g1, lo=17.55, hi=17.65, id="sim_4")
        0.039079462666832203

Figure 6: Comparison of predicted IXO/CAT RS CVn binary source spectra and actual Chandra source spectrum

Simulated IXO/CAT grating line spectra versus Chandra ACIS-S/HETG positive first order line spectrum. Note the larger effective area and resolving power of IXO over Chandra in this wavelength range.

Note the resulting difference in the total number of counts in the IXO spectra as compared with Chandra, as well as the higher resolution of the emission lines especially in orders 2 through 4. In order to detect a spectral line shift during a stellar flare, for example, sufficient sensitivity and spectral resolution are needed (Aeff ~3000 cm2 at 1-2 keV, R=2000-3000). IXO will be able to resolve spectral features 100 times fainter than currently possible, and will have 20 times more collecting area at 1 keV than any previous X-ray telescope. Enhanced line-profile analysis of objects like active CV binaries and massive stars represents just one example of the many ways in which IXO/CAT grating spectra will advance high resolution X-ray spectroscopy.

Scripting It

The file fit.py is a Python script which performs the primary commands used above; it can be executed by typing execfile("fit.py") 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)

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.5 syntax to you.

History

14 Sep 2009	original version
22 Mar 2010	new fitting script
13 Jul 2010	updated for CIAO 4.2 Sherpa v2: removal of S-Lang version of thread.
15 Dec 2011	reviewed for CIAO 4.4: a work-around for a `save_pha` bug was added