Chandra X-Ray Observatory
	(CXC)
Skip to the navigation links
Last modified: 13 December 2012

URL: http://cxc-dmz-prev.cfa.harvard.edu/proposer/threads/binary_sherpa/thread.html

Simulate a High Energy Transmission Grating Spectrum of an X-ray Binary with Sherpa

Proposer Threads (Cycle 15)


Contents


1. Thread Overview

This thread shows how to simulate a simple high energy transmission gratings (HETG) spectrum of an X-ray binary, and to assess whether it is significantly affected by pile-up. The simulations are based upon a successful Chandra proposal to study a low mass X-ray binary (LMXB). This source was detected with XMM-Newton, with the CCD spectra revealing flourescent Fe Kalpha emission at 6.4 keV, absorption from highly ionized Fe, as well as strong absorption from the interstellar medium. The high energies of these features are well-beyond the upper cutoff energies of the Reflection Gratings Spectrometer (RGS) on-board XMM-Newton, but are well-suited to the capabilities of HETG. The goal of the feasibility study is to determine the constraints that can be placed on the Fe features, as well as on the neutral column, using the high spectral resolution of the HETG instrument, when integrating over the X-ray bright phases of the binary orbit. Additionally, we wish to find out whether or not we must make any instrumental modifications to mitigate the possibility of photon pile-up.

During the bright phases of the binary orbit, which lasts approximately 65 ksec, the source was observed to have a relatively stable spectrum. We therefore only need to simulate a single spectrum. Similar procedures, however, could be followed for simulating the X-ray faint phases of the binary orbit. XMM-Newton observations show that the unabsorbed luminosity, presuming a distance of 15 kpc, is 4x10^37 ergs/s, yielding an unabsorbed flux of 1.5e-9 ergs/cm^2/s. This is sufficiently bright that the question of photon pile-up must be considered; however, observations also show the NH=7.5x10^22 atoms/cm^2. This high column, coupled with the high energy features, implies that the High Energy Gratings (HEG - covering 0.7-8 keV) will be the prime Chandra instrument, although the Medium Energy Gratings (MEG - covering 0.4-7 keV) will be important for studies of the absorption column, and may even directly measure absorption edges (e.g., Si at 6.67 Angstrom, S at 5 Angstrom).

This thread uses a Sherpa simulation. It you prefer to work with ISIS please refer to the companion ISIS thread.


2. Preliminary Considerations

Here are some things to consider when following this thread.
  • You need to have Sherpa installed. Sherpa, the CIAO plotting and fitting application, comes bundled with the CIAO software package; it may also be downloaded independently from CIAO as a standalone application (click here for details).
    This thread is saved as a Python script which a user can simply edit and customize, and rerun for their own source.
  • Download the latest instrument responses. This thread uses Cycle 13 responses and CIAO4.4 for illustration purposes: there might be some minor numerical differences when using the latest responses and some slight cosmetic differences in the Sherpa prompts.
    Note: before running it, the script needs to be edited so that the proper responses are used. For example the command hma = unpack_arf("aciss_heg-1_cy13.garf") will become hma = unpack_arf("aciss_heg-1_cy15.garf") when using cycle 15 responses. As a reminder pre-made responses are provided for proposal purposes only and THEY SHOULD NOT BE USED FOR DATA ANALYSIS.
  • Check for previous observations of this source. You will need to justify additional observations if any are found. For detailed instructions, please refer to the "Proposing to Observe an On-Axis Point Source with ACIS" thread.
  • Check for other bright sources in the field of view. The Observation Visualizer, ObsVis, available as part of CIAO, is used for this task. For further information, please refer to the the "Proposing to Observe an On-Axis Point Source with ACIS " thread. Note, however, that due to photon order sorting - i.e., the concept that for a dispersed gratings spectrum a specific location along the dispersion arms must correspond to a narrow range of energies in order to be associated with the spectrum from a known source location - gratings observations (especially of bright point sources) can be fairly tolerant of background and other sources within the field of view. Complicated regions with many sources, might require the use of MARX simulations (http://space.mit.edu/ASC/MARX/).
  • An estimate of the neutral hydrogen column to the source can be obtained from colden if it is not already known. For this source, however, previous observations indicate an NH of 7.5e22 cm^-2.

3. Simulate Source Spectrum

Simulation Script

The steps in this thread have been incorporated into a script, lmxb.py (which imports the "save_pars()" function defined in the separate script save_pars.py). Both scripts need to be downloaded for running in Sherpa. The easist way to follow this thread is to cut and paste from the script, or edit the script and run at the Sherpa or CIAO prompt:

unix% ciao
CIAO configuration is complete... 
CIAO 4.4 Friday, December  2, 2011
  bindir      : /soft/ciao-4.4/bin

unix% sherpa lmxb.py

- OR - 

unix% sherpa
-----------------------------------------------------
Welcome to Sherpa: CXC's Modeling and Fitting Package
-----------------------------------------------------
CIAO 4.4 Sherpa version 1 Friday, December 2, 2011

sherpa-1> execfile("lmxb.py")


Run Sherpa and read in responses

Sherpa includes the "additive" and "multiplicative" models of XSpec version 12.6.01. The XSpec model names have the prefix "xs", e.g. XSpec phabs model is called "xsphabs." In this thread, we will only use standard XSPEC spectral models.

unix% sherpa
-----------------------------------------------------
Welcome to Sherpa: CXC's Modeling and Fitting Package
-----------------------------------------------------
CIAO 4.4 Sherpa version 1 Thursday, December 2, 2011

sherpa-1> list_models("xspec")
     
['xsabsori',
 'xsacisabs',
 'xsapec',
 'xsbapec',
 'xsbbody',
 'xsbbodyrad',    
    . . .

 'xszvphabs',
 'xszwabs',
 'xszwndabs',
 'xszxipcf']

The CXC provides a set of instrument responses for use in proposal planning. For this proposal, we use a set of responses appropriate for the (plus and minus) first order HEG and MEG instruments. This consititutes eight files in all: two ancillary response files (ARF) each for HEG and MEG, and two redistribution matrix files (RMF) each for HEG and MEG. The ARF files are the most crucial for estimating count rates, but the RMF files are necessary to accurately estimate line widths.

Load the ARFs and RMFs and assign them to variables using the Sherpa "unpack_arf" and "unpack_rmf" commands.

sherpa> hma = unpack_arf("aciss_heg-1_cy13.garf")
sherpa> hpa = unpack_arf("aciss_heg1_cy13.garf")
sherpa> mma = unpack_arf("aciss_meg-1_cy13.garf")
sherpa> mpa = unpack_arf("aciss_meg1_cy13.garf")

sherpa> hmr = unpack_rmf("aciss_heg-1_cy13.grmf")
sherpa> hpr = unpack_rmf("aciss_heg1_cy13.grmf")
sherpa> mmr = unpack_rmf("aciss_meg-1_cy13.grmf")
sherpa> mpr = unpack_rmf("aciss_meg1_cy13.grmf")

Define the Model

For this spectrum, we choose the sum of a blackbody plus a power law plus a gaussian (6.4 keV fluorescence) minus two gaussian absorption components (ionized Fe absorption), all multiplied by neutral column absorption. This is done with the Sherpa "set_source" command:

sherpa> for i in [1,2,3,4]:
            set_source(i, "xsphabs.abs1*(xsbbodyrad.bb1+xspowerlaw.p1+xsgaussian.g1-xsgaussian.g2-xsgaussian.g3)")


# Dataset 1 = HEG-1
# Dataset 2 = HEG+1
# Dataset 3 = MEG-1
# Dataset 4 = MEG+1

Sherpa model parameter values may be set on the command line with the "set_par" command. We take values from the fits to the XMM-Newton data:

sherpa> set_par(abs1.nH,     val = 0.)
sherpa> set_par(bb1.norm,    val = 41.)
sherpa> set_par(bb1.kT,      val = 1.18)
sherpa> set_par(p1.norm,     val = 0.15)
sherpa> set_par(p1.PhoIndex, val = 2.02)
sherpa> set_par(g1.norm,     val = 3.e-5)
sherpa> set_par(g1.linee,    val = 6.4)
sherpa> set_par(g1.sigma,    val = 0.1)
sherpa> set_par(g2.norm,     val = 3.e-5)
sherpa> set_par(g2.linee,    val = 6.72)
sherpa> set_par(g2.sigma,    val = 0.02)
sherpa> set_par(g3.norm,     val = 3.e-5)
sherpa> set_par(g3.linee,    val = 7.00)
sherpa> set_par(g3.sigma,    val = 0.02)

Note that the paper describing the XMM-Newton spectra only give the line equivalent widths, thus we initially set the line normalizations to low values, and discuss below how we set the normalizations based upon the published equivalent widths. Additionally, the XMM-Newton paper quotes the unabsorbed flux. Thus, we begin by setting the absorption to 0, to verify that we have properly understood the model normalization described in the XMM-Newton paper. The XMM-Newton spectra only give 90% confidence value upper limits to the absorption line physical widths; we choose physical line widths 40% of those values.


Run fake_pha, and Set Model Normalizations

Having defined a model (which, by default, is the same for all four gratings arms) and instrument response, we now create a fake dataset (complete with Poisson noise) with the Sherpa "fake_pha" command. The exposure time for the faked data files is set via the 'exposure' argument. We choose to set the exposure times for the faked datasets to 650 ksec - a factor of 10 higher integration time than desired, at first, so that we have sufficiently good statistics to verify our settings of the model normalizations. The "fake_pha" command should be followed by "show_data" or "show_all" to verify that the datasets have indeed been created.

sherpa> expos = 650000.

sherpa> for i in [hma,hpa,mma,mpa]:
            i.exposure = expos 

sherpa> fake_pha(1, arf = hma, rmf = hmr, exposure = expos)
sherpa> fake_pha(2, arf = hpa, rmf = hpr, exposure = expos)
sherpa> fake_pha(3, arf = mma, rmf = mmr, exposure = expos)
sherpa> fake_pha(4, arf = mpa, rmf = mpr, exposure = expos)

sherpa> set_analysis("wavelength")

sherpa> show_data()
Data Set: 1
Filter: 1.0012-21.4788 Wavelength (Angstrom)
Noticed Channels: 1-8192
name           = faked
channel        = Float64[8192]
counts         = Float64[8192]
staterror      = None
syserror       = None
bin_lo         = None
bin_hi         = None
grouping       = None
quality        = None
exposure       = 650000.0
backscal       = None
areascal       = None
grouped        = False
subtracted     = False
units          = wavelength
rate           = True
plot_fac       = 0
response_ids   = [1]
background_ids = []

RMF Data Set: 1:1
name     = aciss_heg-1_cy13.grmf
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]

ARF Data Set: 1:1
name     = aciss_heg-1_cy13.garf
energ_lo = Float64[8192]
energ_hi = Float64[8192]
specresp = Float64[8192]
bin_lo   = Float64[8192]
bin_hi   = Float64[8192]
exposure = 650000.0

Data Set: 2
Filter: 1.0012-21.4788 Wavelength (Angstrom)
Noticed Channels: 1-8192
name           = faked
channel        = Float64[8192]
counts         = Float64[8192]
staterror      = None
syserror       = None
bin_lo         = None
bin_hi         = None
grouping       = None
quality        = None
exposure       = 650000.0
backscal       = None
areascal       = None
grouped        = False
subtracted     = False
units          = wavelength
rate           = True
plot_fac       = 0
response_ids   = [1]
background_ids = []

RMF Data Set: 2:1
name     = aciss_heg1_cy13.grmf
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]

ARF Data Set: 2:1
name     = aciss_heg1_cy13.garf
energ_lo = Float64[8192]
energ_hi = Float64[8192]
specresp = Float64[8192]
bin_lo   = Float64[8192]
bin_hi   = Float64[8192]
exposure = 650000.0

Data Set: 3
Filter: 1.0025-41.9575 Wavelength (Angstrom)
Noticed Channels: 1-8192
name           = faked
channel        = Float64[8192]
counts         = Float64[8192]
staterror      = None
syserror       = None
bin_lo         = None
bin_hi         = None
grouping       = None
quality        = None
exposure       = 650000.0
backscal       = None
areascal       = None
grouped        = False
subtracted     = False
units          = wavelength
rate           = True
plot_fac       = 0
response_ids   = [1]
background_ids = []

RMF Data Set: 3:1
name     = aciss_meg-1_cy13.grmf
detchans = 8192
energ_lo = Float64[8192]
energ_hi = Float64[8192]
n_grp    = UInt64[8192]
f_chan   = UInt32[8192]
n_chan   = UInt32[8192]
matrix   = Float64[2367624]
offset   = 1
e_min    = Float64[8192]
e_max    = Float64[8192]

ARF Data Set: 3:1
name     = aciss_meg-1_cy13.garf
energ_lo = Float64[8192]
energ_hi = Float64[8192]
specresp = Float64[8192]
bin_lo   = Float64[8192]
bin_hi   = Float64[8192]
exposure = 650000.0

Data Set: 4
Filter: 1.0025-41.9575 Wavelength (Angstrom)
Noticed Channels: 1-8192
name           = faked
channel        = Float64[8192]
counts         = Float64[8192]
staterror      = None
syserror       = None
bin_lo         = None
bin_hi         = None
grouping       = None
quality        = None
exposure       = 650000.0
backscal       = None
areascal       = None
grouped        = False
subtracted     = False
units          = wavelength
rate           = True
plot_fac       = 0
response_ids   = [1]
background_ids = []

RMF Data Set: 4:1
name     = aciss_meg1_cy13.grmf
detchans = 8192
RMF Data Set: 4:1
name     = aciss_meg1_cy13.grmf
detchans = 8192
energ_lo = Float64[8192]
energ_hi = Float64[8192]
n_grp    = UInt64[8192]
f_chan   = UInt32[8192]
n_chan   = UInt32[8192]
matrix   = Float64[2470385]
offset   = 1
e_min    = Float64[8192]
e_max    = Float64[8192]

ARF Data Set: 4:1
name     = aciss_meg1_cy13.garf
energ_lo = Float64[8192]
energ_hi = Float64[8192]
specresp = Float64[8192]
bin_lo   = Float64[8192]
bin_hi   = Float64[8192]
exposure = 650000.0

The HEG detectors (datasets 1 and 2) have less effective area, and hence fewer counts, than the MEG detectors (datasets 3 and 4), as shown by the output of the Sherpa "calc_data_sum" command.

sherpa> for i in [1,2,3,4]:
            calc_data_sum(id=i)

        2752752.0
        2587280.0
        4974808.0
        5356748.0

To set the normalizations, we use the "calc_energy_flux" and "calc_photon_flux" commands to define variables "pflux" and "eflux", that will return the flux in units of ergs and photons per cm2/sec. We also use the "eqwidth" function to returns line equivalent widths (for a given detector and model component normalization) in units of Angstroms and keV, which we convert to mAngstrom and eV.

sherpa> set_analysis("energy")

sherpa> eflux = calc_energy_flux(0.1, 8, id=3)   # Meg -1 flux, from 0.1 to 8 keV
sherpa> pflux = calc_photon_flux(0.1, 8, id=3)
sherpa> print("pflux="+str(pflux)+", eflux="+str(eflux))

pflux=0.653083725594, eflux=1.56797250026e-09


sherpa> g1_e = eqwidth(abs1*(bb1+p1),abs1*(bb1+p1+g1), id=1) # HEG -1 EW in keV
sherpa> g2_e = eqwidth(abs1*(bb1+p1),abs1*(bb1+p1-g2), id=1) # HEG -1 EW in keV
sherpa> g3_e = eqwidth(abs1*(bb1+p1),abs1*(bb1+p1-g3), id=1) # HEG -1 EW in keV

sherpa> set_analysis("wavelength")

sherpa> g1_a = eqwidth(abs1*(bb1+p1),abs1*(bb1+p1+g1), id=1) # HEG -1 EW in Ang
sherpa> g2_a = eqwidth(abs1*(bb1+p1),abs1*(bb1+p1-g2), id=1) # HEG -1 EW in Ang
sherpa> g3_a = eqwidth(abs1*(bb1+p1),abs1*(bb1+p1-g3), id=1) # HEG -1 EW in Ang


# Print EWs in mA and eV
sherpa> print(g1_a*1e3, g1_e*1e3)
(0.80933931459992092, 2.6757312831787838)

sherpa> print(g2_a*1e3, g2_e*1e3)
(-0.85266245548596598, -3.1057260876663078)

sherpa> print(g3_a*1e3, g3_e*1e3)
(-0.8960585465056421, -3.5414413630000472)

We see that the unabsorbed flux is very close to the quoted value, therefore we do not alter the normalizations of the continuum components. The line equivalent widths, however, are all too small. The three features should have equivalent widths of 95 eV, -14 eV, and -23 eV (from low to high energy). We therefore scale the line normalizations by the ratio of desired to measured equivalent widths, run fake_pha again, and then verify that the new equivalent widths are those that we desire.

sherpa> set_analysis("energy")

sherpa> set_par(g1.norm, g1.norm.val*(95./(g1_e*1e3)))
sherpa> set_par(g2.norm, g2.norm.val*((-14.)/(g2_e*1e3)))
sherpa> set_par(g3.norm, g3.norm.val*((-23.)/(g3_e*1e3)))

sherpa> fake_pha(1, hma, hmr, exposure = expos)
sherpa> fake_pha(2, hpa, hpr, exposure = expos)
sherpa> fake_pha(3, mma, mmr, exposure = expos)
sherpa> fake_pha(4, mpa, mpr, exposure = expos)

sherpa> g1_e = eqwidth(abs1*(bb1+p1),abs1*(bb1+p1+g1), id=1)
sherpa> g2_e = eqwidth(abs1*(bb1+p1),abs1*(bb1+p1-g2), id=1)
sherpa> g3_e = eqwidth(abs1*(bb1+p1),abs1*(bb1+p1-g3), id=1)

sherpa> set_analysis("wavelength")
sherpa> g1_a = eqwidth(abs1*(bb1+p1),abs1*(bb1+p1+g1), id=1) 
sherpa> g2_a = eqwidth(abs1*(bb1+p1),abs1*(bb1+p1-g2), id=1)
sherpa> g3_a = eqwidth(abs1*(bb1+p1),abs1*(bb1+p1-g3), id=1)

# Print EWs in mA and eV
sherpa> print(g1_a*1e3, g1_e*1e3)
(28.735037635426384, 95.000005220072708)

sherpa> print(g2_a*1e3, g2_e*1e3)
(-3.843634151870051, -14.000000563344146)

sherpa> print(g3_a*1e3, g3_e*1e3)
(-5.8194800607770567, -23.000000924861549)

Although not exact matches to the desired equivalent widths, these values are well-within the XMM-Newton error bars, so we deem the normalizations sufficiently accurate for simulation purposes.

Finally, we set the neutral column to the value actually measured by XMM-Newton, we reduce the exposure time to 65 ksec, and we run fake_pha once more. Furthermore, we save the final model parameters to variable "lmxb_pars" using the imported "save_pars" function, in order to restore the settings later in the session.

sherpa> from save_pars import *

sherpa> set_par(abs1.nH, val=7.5)  

sherpa> comps = ["abs1", "bb1", "p1", "g1", "g2", "g3"]  # or list_model_components()

sherpa> lmxb_pars = save_pars(comps)
    
sherpa> expos = 65000.

sherpa> for i in [hma,hpa,mma,mpa]:
            i.exposure = expos 

sherpa> fake_pha(1, hma, hmr, exposure = expos)
sherpa> fake_pha(2, hpa, hpr, exposure = expos)
sherpa> fake_pha(3, mma, mmr, exposure = expos)
sherpa> fake_pha(4, mpa, mpr, exposure = expos)

sherpa> show_data()
Data Set: 1
Filter: 1.0012-21.4788 Wavelength (Angstrom)
Noticed Channels: 1-8192
name           = faked
channel        = Float64[8192]
counts         = Float64[8192]
staterror      = None
syserror       = None
bin_lo         = None
bin_hi         = None
grouping       = None
quality        = None
exposure       = 65000.0
backscal       = None
areascal       = None
grouped        = False
subtracted     = False
units          = wavelength
rate           = True
plot_fac       = 0
response_ids   = [1]
background_ids = []

RMF Data Set: 1:1
name     = aciss_heg-1_cy13.grmf
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]

ARF Data Set: 1:1
name     = aciss_heg-1_cy13.garf
energ_lo = Float64[8192]
energ_hi = Float64[8192]
specresp = Float64[8192]
bin_lo   = Float64[8192]
bin_hi   = Float64[8192]
exposure = 65000.0

Data Set: 2
Filter: 1.0012-21.4788 Wavelength (Angstrom)
Noticed Channels: 1-8192
name           = faked
channel        = Float64[8192]
counts         = Float64[8192]
staterror      = None
syserror       = None
bin_lo         = None
bin_hi         = None
grouping       = None
quality        = None
exposure       = 65000.0
backscal       = None
areascal       = None
grouped        = False
subtracted     = False
units          = wavelength
rate           = True
plot_fac       = 0
response_ids   = [1]
background_ids = []

RMF Data Set: 2:1
name     = aciss_heg1_cy13.grmf
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]

ARF Data Set: 2:1
name     = aciss_heg1_cy13.garf
energ_lo = Float64[8192]
energ_hi = Float64[8192]
specresp = Float64[8192]
bin_lo   = Float64[8192]
bin_hi   = Float64[8192]
exposure = 65000.0

Data Set: 3
Filter: 1.0025-41.9575 Wavelength (Angstrom)
Noticed Channels: 1-8192
name           = faked
channel        = Float64[8192]
counts         = Float64[8192]
staterror      = None
syserror       = None
bin_lo         = None
bin_hi         = None
grouping       = None
quality        = None
exposure       = 65000.0
backscal       = None
areascal       = None
grouped        = False
subtracted     = False
units          = wavelength
rate           = True
plot_fac       = 0
response_ids   = [1]
background_ids = []

RMF Data Set: 3:1
name     = aciss_meg-1_cy13.grmf
detchans = 8192
energ_lo = Float64[8192]
energ_hi = Float64[8192]
n_grp    = UInt64[8192]
f_chan   = UInt32[8192]
n_chan   = UInt32[8192]
matrix   = Float64[2367624]
offset   = 1
e_min    = Float64[8192]
e_max    = Float64[8192]

ARF Data Set: 3:1
name     = aciss_meg-1_cy13.garf
energ_lo = Float64[8192]
energ_hi = Float64[8192]
specresp = Float64[8192]
bin_lo   = Float64[8192]
bin_hi   = Float64[8192]
exposure = 65000.0

Data Set: 4
Filter: 1.0025-41.9575 Wavelength (Angstrom)
Noticed Channels: 1-8192
name           = faked
channel        = Float64[8192]
counts         = Float64[8192]
staterror      = None
syserror       = None
bin_lo         = None
bin_hi         = None
grouping       = None
quality        = None
exposure       = 65000.0
backscal       = None
areascal       = None
grouped        = False
subtracted     = False
units          = wavelength
rate           = True
plot_fac       = 0
response_ids   = [1]
background_ids = []

RMF Data Set: 4:1
name     = aciss_meg1_cy13.grmf
detchans = 8192
RMF Data Set: 4:1
name     = aciss_meg1_cy13.grmf
detchans = 8192
energ_lo = Float64[8192]
energ_hi = Float64[8192]
n_grp    = UInt64[8192]
f_chan   = UInt32[8192]
n_chan   = UInt32[8192]
matrix   = Float64[2470385]
offset   = 1
e_min    = Float64[8192]
e_max    = Float64[8192]

ARF Data Set: 4:1
name     = aciss_meg1_cy13.garf
energ_lo = Float64[8192]
energ_hi = Float64[8192]
specresp = Float64[8192]
bin_lo   = Float64[8192]
bin_hi   = Float64[8192]
exposure = 65000.0

sherpa> for i in [1,2,3,4]:
            calc_data_sum(id=i)

        70330.0
        70711.0
        58189.0
        68961.0

4. Assessing the Simulations

Fe lines viewed with HEG data

We begin by considering the Fe line features in the 6-7.5 keV region of the HEG spectrum. Since we are mostly concerned with just the line features, we perform a local fit for the continuum model in just the 1.65-2.3 Angstrom (5.4-7.5 keV) range. Given the restricted wavelength/energy range, we freeze the values of the neutral column, the blackbody temperature, and the power law photon index (i.e., we only allow continuum normalizations and line parameters to vary in the fit.) Furthermore, we at first will only consider fits to the higher spectral resolution HEG data.

To increase the signal-to-noise, we group the HEG data sets by a factor of 2 and the MEG data by a factor of 4 (i.e., to about twice and one times, respectively, the intrinsic resolution of the gratings detectors). Additionally, for fitting purposes, we will combine the two HEG datasets, and we will combine the two MEG datasets.

sherpa> group_bins(1, len(get_data().channel)/2) # Group HEG data by a factor of 2
sherpa> group_bins(2, len(get_data().channel)/2) 
sherpa> group_bins(3, len(get_data().channel)/4) # Group MEG data by a factor of 4
sherpa> group_bins(4, len(get_data().channel)/4)

sherpa> copy_data(1,"heg")
sherpa> set_source("heg", get_source(1))

sherpa> copy_data(3,"meg")
sherpa> set_source("meg", get_source(3))

sherpa> get_data("heg").counts = get_data(1).counts+get_data(2).counts
sherpa> get_data("meg").counts = get_data(3).counts+get_data(4).counts

sherpa> get_arf("heg").specresp = get_arf(1).specresp+get_arf(2).specresp
sherpa> get_arf("meg").specresp = get_arf(3).specresp+get_arf(4).specresp

sherpa> set_analysis("wavelength")
sherpa> notice_id([1,2,"heg"], 1.65, 2.3)        # Only notice the 1.65-2.3 A data for HEG

sherpa> freeze(abs1.nH, bb1.kT, p1.PhoIndex)     # Freeze Nh, blackbody kT, Photon Index

sherpa> fit(1,2,"heg")                           # Ignore MEG data for now
Datasets              = 1, 2, 'heg'
Method                = levmar
Statistic             = chi2gehrels
Initial fit statistic = 364.744
Final fit statistic   = 332.827 at function evaluation 193
Data points           = 396
Degrees of freedom    = 385
Probability [Q-value] = 0.974305
Reduced statistic     = 0.864486
Change in statistic   = 31.917
   bb1.norm       41.0482     
   p1.norm        0.146431    
   g1.LineE       6.39515     
   g1.Sigma       0.091802    
   g1.norm        0.00115142  
   g2.LineE       6.70795     
   g2.Sigma       0.0281879   
   g2.norm        0.000195587 
   g3.LineE       7.00995     
   g3.Sigma       0.0309599   
   g3.norm        0.000265202 


sherpa> heg_fit_pars = save_pars(comps)          # Save the HEG fit parameters

sherpa> set_analysis(1, "wavelength", "counts", factor=0)
sherpa> set_analysis(2, "wavelength", "counts", factor=0)

sherpa> plot_fit_delchi(1)
sherpa> limits(X_AXIS, 1.65, 2.3)
sherpa> current_plot("plot1")
sherpa> lin_scale(X_AXIS)
sherpa> log_scale(Y_AXIS)
sherpa> current_plot("plot2")
sherpa> lin_scale(XY_AXIS)
sherpa> current_plot("plot1")
sherpa> set_plot_title("HEG -1")
sherpa> print_window("heg_m1_plot.ps")

sherpa> plot_fit_delchi(2)
sherpa> limits(X_AXIS, 1.65, 2.3)
sherpa> current_plot("plot1")
sherpa> lin_scale(X_AXIS)
sherpa> log_scale(Y_AXIS)
sherpa> current_plot("plot1")
sherpa> set_plot_title("HEG +1")
sherpa> current_plot("plot2")
sherpa> lin_scale(XY_AXIS)
sherpa> print_window("heg_p1_plot.ps")

sherpa> set_analysis("heg", "wavelength", "counts", factor=0)
sherpa> set_analysis("meg", "wavelength", "counts", factor=0)

sherpa> plot_fit("heg")
sherpa> set_plot_title("Combined HEG Data")
sherpa> print_window("heg_combined_plot.ps")

The resulting plots are seen here (Figure 1), here (Figure 2), and here (Figure 3). The "print_window" command is used to save the current plot to a color PostScript file, e.g., to generate a hardcopy of the plot.

We now turn to determining the confidence bounds for the line parameters found above. For this, we use the Sherpa "conf" command, specifying 90% confidence limits (default is 68%). The conf command will return the upper and lower parameter values for these confidence intervals.

sherpa> set_conf_opt("sigma", 1.6)

sherpa> conf(1,2,"heg", g1.norm,g1.LineE,g1.Sigma)
g1.norm lower bound:	-0.000110968
g1.Sigma lower bound:	-0.00964446
g1.Sigma upper bound:	0.0108352
g1.norm upper bound:	0.00011528
g1.LineE lower bound:	-0.0090216
g1.LineE upper bound:	0.00921093
Datasets              = 1, 2, 'heg'
Confidence Method     = confidence
Iterative Fit Method  = None
Fitting Method        = levmar
Statistic             = chi2gehrels
confidence 1.6-sigma (89.0401%) bounds:
   Param            Best-Fit  Lower Bound  Upper Bound
   -----            --------  -----------  -----------
   g1.norm        0.00115142 -0.000110968   0.00011528
   g1.LineE          6.39515   -0.0090216   0.00921093
   g1.Sigma         0.091802  -0.00964446    0.0108352


sherpa> conf(1,2,"heg", g2.norm,g2.LineE,g2.Sigma)
g2.norm lower bound:	-4.81348e-05
g2.Sigma lower bound:	-0.00893762
g2.Sigma upper bound:	0.00940272
g2.norm upper bound:	5.0964e-05
g2.LineE lower bound:	-0.00982178
g2.LineE upper bound:	0.00979638
Datasets              = 1, 2, 'heg'
Confidence Method     = confidence
Iterative Fit Method  = None
Fitting Method        = levmar
Statistic             = chi2gehrels
confidence 1.6-sigma (89.0401%) bounds:
   Param            Best-Fit  Lower Bound  Upper Bound
   -----            --------  -----------  -----------
   g2.norm       0.000195587 -4.81348e-05   5.0964e-05
   g2.LineE          6.70795  -0.00982178   0.00979638
   g2.Sigma        0.0281879  -0.00893762   0.00940272


sherpa> set_method("neldermead") 
sherpa> conf(1,2,"heg", g2.Sigma)
g2.Sigma lower bound:	-0.00866418
g2.Sigma +: WARNING: The confidence level lies within (3.622711e-02, 3
.872711e-02)
g2.Sigma upper bound:	0.00928919
Datasets              = 1, 2, 'heg'
Confidence Method     = confidence
Iterative Fit Method  = None
Fitting Method        = neldermead
Statistic             = chi2gehrels
confidence 1.6-sigma (89.0401%) bounds:
   Param            Best-Fit  Lower Bound  Upper Bound
   -----            --------  -----------  -----------
   g2.Sigma        0.0281879  -0.00866418   0.00928919


sherpa> set_method("levmar") 
sherpa> conf(1,2,"heg", g3.norm,g3.LineE,g3.Sigma)
g3.norm lower bound:	-5.63587e-05
g3.Sigma lower bound:	-0.0089006
g3.Sigma upper bound:	0.00972749
g3.norm upper bound:	6.11435e-05
g3.LineE lower bound:	-0.00855174
g3.LineE upper bound:	0.00834949
Datasets              = 1, 2, 'heg'
Confidence Method     = confidence
Iterative Fit Method  = None
Fitting Method        = levmar
Statistic             = chi2gehrels
confidence 1.6-sigma (89.0401%) bounds:
   Param            Best-Fit  Lower Bound  Upper Bound
   -----            --------  -----------  -----------
   g3.norm       0.000265202 -5.63587e-05  6.11435e-05
   g3.LineE          7.00995  -0.00855174   0.00834949
   g3.Sigma        0.0309599   -0.0089006   0.00972749

We find that the line positions and widths are all determined with 90% confidence level uncertainties of approximately 0.01 keV ; in terms of velocities, this corresponds to approximately 500 km/sec. Note that for one of the error bars searches, the fit failed, so we switched to a different (slower) fit method, and searched for the error bar again.

Above, we also find the uncertainties in the line normalizations. To properly cast that in terms of equivalent width uncertainties, one should individually change and freeze the line normalizations to their confidence limits, re-fit the model, and then recalculate the equivalent width; we show this below.

sherpa> set_method("neldermead")

sherpa> set_analysis("wavelength")
sherpa> g1_a = eqwidth(abs1*(bb1+p1),abs1*(bb1+p1+g1), id=1)
sherpa> set_analysis("energy")
sherpa> g1_e = eqwidth(abs1*(bb1+p1),abs1*(bb1+p1+g1), id=1)

sherpa> print(g1_a*1e3, g1_e*1e3)
(31.246131570801506, 103.13365776626956)

sherpa> set_par(g1.norm, val = 0.000814, frozen=True)
sherpa> fit(1,2) 
Datasets              = 1, 2
Method                = neldermead
Statistic             = chi2gehrels
Initial fit statistic = 250.323
Final fit statistic   = 242.306 at function evaluation 2613
Data points           = 264
Degrees of freedom    = 254
Probability [Q-value] = 0.690518
Reduced statistic     = 0.953961
Change in statistic   = 8.01708
   bb1.norm       38.1959     
   p1.norm        0.175624    
   g1.LineE       6.39209     
   g1.Sigma       0.0739747   
   g2.LineE       6.70859     
   g2.Sigma       0.0297662   
   g2.norm        0.000220109 
   g3.LineE       7.00882     
   g3.Sigma       0.0335095   
   g3.norm        0.000308113 


sherpa> set_analysis("wavelength")
sherpa> g1_a = eqwidth(abs1*(bb1+p1),abs1*(bb1+p1+g1), id=1)
sherpa> set_analysis("energy")
sherpa> g1_e = eqwidth(abs1*(bb1+p1),abs1*(bb1+p1+g1), id=1)

sherpa> print(g1_a*1e3, g1_e*1e3)
(21.782190608093064, 71.81023581167176)

sherpa> set_par(g1.norm, val = 0.00112, frozen=True)
sherpa> fit(1,2) 
Datasets              = 1, 2
Method                = neldermead
Statistic             = chi2gehrels
Initial fit statistic = 238.931
Final fit statistic   = 229.128 at function evaluation 3252
Data points           = 264
Degrees of freedom    = 254
Probability [Q-value] = 0.866915
Reduced statistic     = 0.902078
Change in statistic   = 9.80279
   bb1.norm       41.1162     
   p1.norm        0.145203    
   g1.LineE       6.39461     
   g1.Sigma       0.0898398   
   g2.LineE       6.70751     
   g2.Sigma       0.0282028   
   g2.norm        0.000196481 
   g3.LineE       7.00906     
   g3.Sigma       0.03021     
   g3.norm        0.000269458 

sherpa> set_pars(comps, heg_fit_pars)

sherpa> set_analysis("wavelength")
sherpa> g2_a = eqwidth(abs1*(bb1+p1),abs1*(bb1+p1-g2), id=1)
sherpa> set_analysis("energy")
sherpa> g2_e = eqwidth(abs1*(bb1+p1),abs1*(bb1+p1-g2), id=1)

sherpa> print(g2_a*1e3, g2_e*1e3)
(-5.5871334076782784, -20.278184904584421)

sherpa> set_par(g2.norm, 0.000135569, frozen=True)
sherpa> fit(1,2)
Datasets              = 1, 2
Method                = neldermead
Statistic             = chi2gehrels
Initial fit statistic = 232.487
Final fit statistic   = 231.181 at function evaluation 2175
Data points           = 264
Degrees of freedom    = 254
Probability [Q-value] = 0.844854
Reduced statistic     = 0.910163
Change in statistic   = 1.30546
   bb1.norm       40.7684     
   p1.norm        0.146349    
   g1.LineE       6.39413     
   g1.Sigma       0.0912387   
   g1.norm        0.00115623  
   g2.LineE       6.70766     
   g2.Sigma       0.0213835   
   g3.LineE       7.00915     
   g3.Sigma       0.0298476   
   g3.norm        0.000265441 

 
sherpa> set_analysis("wavelength")
sherpa> g2_a = eqwidth(abs1*(bb1+p1),abs1*(bb1+p1-g2), id=1)
sherpa> set_analysis("energy")
sherpa> g2_e = eqwidth(abs1*(bb1+p1),abs1*(bb1+p1-g2), id=1)

sherpa> print(g2_a*1e3, g2_e*1e3)
(-3.8910199882129457, -14.120633819485851)


sherpa> set_par(g2.norm, val = 0.000306819, frozen=True)
sherpa> fit(1,2)
Datasets              = 1, 2
Method                = neldermead
Statistic             = chi2gehrels
Initial fit statistic = 247.575
Final fit statistic   = 234.832 at function evaluation 2904
Data points           = 264
Degrees of freedom    = 254
Probability [Q-value] = 0.800403
Reduced statistic     = 0.924535
Change in statistic   = 12.7429
   bb1.norm       40.9209     
   p1.norm        0.148468    
   g1.LineE       6.39683     
   g1.Sigma       0.093589    
   g1.norm        0.00114488  
   g2.LineE       6.70679     
   g2.Sigma       0.039331    
   g3.LineE       7.00898     
   g3.Sigma       0.0309477   
   g3.norm        0.000278112 


sherpa> set_analysis("wavelength")
sherpa> g2_a = eqwidth(abs1*(bb1+p1),abs1*(bb1+p1-g2), id=1)
sherpa> set_analysis("energy")
sherpa> g2_e = eqwidth(abs1*(bb1+p1),abs1*(bb1+p1-g2), id=1)

sherpa> print(g2_a*1e3, g2_e*1e3)
(-8.7418789696637482, -31.718855343259548)

sherpa> set_pars(comps, heg_pars)

sherpa> set_analysis("wavelength")
sherpa> g3_a = eqwidth(abs1*(bb1+p1),abs1*(bb1+p1-g3), id=1)
sherpa> set_analysis("energy")
sherpa> g3_e = eqwidth(abs1*(bb1+p1),abs1*(bb1+p1-g3), id=1)

sherpa> print(g3_a*1e3, g3_e*1e3)
(-7.9959207472859797, -31.693013352890716)


sherpa> set_par(g2.norm, val=0.000187434, frozen=True)
sherpa> fit(1,2)
Datasets              = 1, 2
Method                = neldermead
Statistic             = chi2gehrels
Initial fit statistic = 229.31
Final fit statistic   = 229.021 at function evaluation 2640
Data points           = 264
Degrees of freedom    = 254
Probability [Q-value] = 0.868002
Reduced statistic     = 0.901658
Change in statistic   = 0.288615
   bb1.norm       40.9951     
   p1.norm        0.145504    
   g1.LineE       6.39459     
   g1.Sigma       0.0913126   
   g1.norm        0.00115409  
   g2.LineE       6.70739     
   g2.Sigma       0.0273525   
   g3.LineE       7.00914     
   g3.Sigma       0.0302963   
   g3.norm        0.00026987  

sherpa> set_analysis("wavelength")
sherpa> g3_a = eqwidth(abs1*(bb1+p1),abs1*(bb1+p1-g3), id=1)
sherpa> set_analysis("energy")
sherpa> g3_e = eqwidth(abs1*(bb1+p1),abs1*(bb1+p1-g3), id=1)

sherpa> print(g3_a*1e3, g3_e*1e3)
(-8.1600770411226247, -32.336099034837218)

sherpa> set_par(g3.norm, val=0.000352399, frozen=True)
sherpa> fit(1,2)
Datasets              = 1, 2
Method                = neldermead
Statistic             = chi2gehrels
Initial fit statistic = 234.301
Final fit statistic   = 231.859 at function evaluation 2705
Data points           = 264
Degrees of freedom    = 255
Probability [Q-value] = 0.847958
Reduced statistic     = 0.909249
Change in statistic   = 2.44273
   bb1.norm       39.6691     
   p1.norm        0.157773    
   g1.LineE       6.39273     
   g1.Sigma       0.0906418   
   g1.norm        0.00115555  
   g2.LineE       6.70765     
   g2.Sigma       0.0272429   
   g3.LineE       7.00859     
   g3.Sigma       0.038366    

sherpa> set_analysis("wavelength")
sherpa> g3_a = eqwidth(abs1*(bb1+p1),abs1*(bb1+p1-g3), id=1)
sherpa> set_analysis("energy")
sherpa> g3_e = eqwidth(abs1*(bb1+p1),abs1*(bb1+p1-g3), id=1)

sherpa> print(g3_a*1e3, g3_e*1e3)
(-10.575598436702972, -41.903129900055376)


Spectrum Viewed with MEG data

We now consider the spectrum, specifically the fitted neutral column, by fitting the combined MEG+/-1 data. First, we exclude the HEG data from consideration, and only consider the 2.1-7.2 Angstrom (1.7-5.9 keV) range of the MEG data. This range has more than 20 counts per bin everywhere, and avoids the Fe line region (which is better viewed with the HEG data).

To facilitate the fitting process, we load the saved simulation parameters (which includes the lines), evaluate the model once (without fitting), and then rewrite the model without the line parameters. In this manner, we restore the unfrozen, starting fit values for the neutral column, blackbody, and power-law.

sherpa> set_analysis("wavelength")
sherpa> notice_id([3,4], 2.1, 7.2)    # Only notice 1.7-7.2 Angstrom (1.7-5.9 keV)

sherpa> set_pars(comps, lmxb_pars);   # Simulation parameters, including lines

sherpa> fit(3,4)
Datasets              = 3, 4
Method                = neldermead
Statistic             = chi2gehrels
Initial fit statistic = 382.335
Final fit statistic   = 370.292 at function evaluation 4552
Data points           = 514
Degrees of freedom    = 500
Probability [Q-value] = 0.999997
Reduced statistic     = 0.740584
Change in statistic   = 12.0429
   abs1.nH        7.86179     
   bb1.kT         1.14034     
   bb1.norm       42.2951     
   p1.PhoIndex    2.00106     
   p1.norm        0.193864    
   g1.LineE       6.93721     
   g1.Sigma       0.255125    
   g1.norm        2.27541     
   g2.LineE       6.25855     
   g2.Sigma       0.158773    
   g2.norm        0.00847606  
   g3.LineE       8.73998     
   g3.Sigma       0.750566    
   g3.norm        0.212904  

sherpa> set_method("levmar")
sherpa> set_source(3, abs1*(bb1+p1))
sherpa> set_source(4, abs1*(bb1+p1))

sherpa> fit(3,4)
Datasets              = 3, 4
Method                = levmar
Statistic             = chi2gehrels
Initial fit statistic = 372.664
Final fit statistic   = 370.937 at function evaluation 32
Data points           = 514
Degrees of freedom    = 509
Probability [Q-value] = 0.999999
Reduced statistic     = 0.728756
Change in statistic   = 1.72748
   abs1.nH        7.83956     
   bb1.kT         1.1233      
   bb1.norm       41.045      
   p1.PhoIndex    1.85041     
   p1.norm        0.179416    


sherpa> show_model(3)    # same model as for data set 4
Model: 3
apply_rmf(apply_arf((65000.0 * (xsphabs.abs1 * (xsbbodyrad.bb1 + xspowerlaw.p1))))
)
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   abs1.nH      thawed      7.83956            0       100000 10^22 atoms / cm^2
   bb1.kT       thawed       1.1233        0.001          100        keV
   bb1.norm     thawed       41.045            0        1e+24           
   p1.PhoIndex  thawed      1.85041           -2            9           
   p1.norm      thawed     0.179416            0        1e+24  

sherpa> meg_pars = save_pars(comps)   # Save the MEG fit parameters

sherpa> conf(3,4,abs1.nH)
abs1.nH lower bound:	-1.04912
abs1.nH upper bound:	1.21502
Datasets              = 3, 4
Confidence Method     = confidence
Iterative Fit Method  = None
Fitting Method        = levmar
Statistic             = chi2gehrels
confidence 1.6-sigma (89.0401%) bounds:
   Param            Best-Fit  Lower Bound  Upper Bound
   -----            --------  -----------  -----------
   abs1.nH           7.83956     -1.04912      1.21502

sherpa> set_analysis(3, "wavelength", "counts", factor=0)
sherpa> set_analysis(4, "wavelength", "counts", factor=0)


# Plot the MEG data            # MEG-1 data (although fit combined with MEG+1)
sherpa> plot_fit_delchi(3)
sherpa> limits(X_AXIS, 2.1, 7.2)       # Only plot 2.1-7.2 Angstroms
sherpa> current_plot("plot1")
sherpa> lin_scale(X_AXIS)
sherpa> log_scale(Y_AXIS)
sherpa> set_plot_title("MEG -1")
sherpa> current_plot("plot2")
sherpa> lin_scale(XY_AXIS)
sherpa> print_window("meg_m1_plot.ps")

sherpa> plot_fit_delchi(4)           # MEG+1 data (although fit combined with MEG-1)
sherpa> limits(X_AXIS, 2.1, 7.2)          
sherpa> current_plot("plot1")
sherpa> lin_scale(X_AXIS)
sherpa> log_scale(Y_AXIS)
sherpa> set_plot_title("MEG 1")
sherpa> current_plot("plot2")
sherpa> lin_scale(XY_AXIS)
sherpa> print_window("meg_p1_plot.ps")

sherpa> clear()                   # close any open plots

sherpa> set_analysis(3, "wavelength")
sherpa> notice_id(3, 2.1, 7.2)

# Calculated the "unfolded" spectrum for MEG-1

sherpa> xx = get_fit_plot(3).dataplot.x     
sherpa> dd = get_fit_plot(3).dataplot.y
sherpa> ee = get_fit_plot(3).dataplot.yerr
sherpa> mm = get_fit_plot(3).modelplot.y

sherpa> src = get_source(3)(xx)
sherpa> add_curve(xx, dd/mm*src, [ee/mm*src],        #Plot the MEG-1 unfolded spectrum
["line.style", "0", "err.style", "line",
"symbol.style", "4","symbol.size", "3",
"symbol.fill", "false"])

sherpa> log_scale(Y_AXIS)
sherpa> lin_scale(X_AXIS)
sherpa> limits(Y_AXIS, 1e-08, AUTO)
sherpa> limits(X_AXIS, 1.8, 7)
sherpa> set_plot_ylabel("Photons sec^{-1} cm^{-2} Ang^{-1}")
sherpa> set_plot_xlabel("Wavelength (Angstrom)")
sherpa> set_plot_title("MEG -1")
sherpa> print_window("meg_m1_unfolded_plot.ps")

Plots of the fits can be found here (Figure 4) and here (Figure 5). Note that most of the structure seen in this counts/bin spectrum is instrumental, as can be seen by comparing to the "unfolded" MEG-1 spectrum shown here (Figure 6). Note, however, that real structure is seen at the Si (6.67 A) and S (5 A) edges.

With these fits, we find that the neutral column is found to be 7.84(+1.22)(-1.05)x10^22 cm^-2. There are fair number of counts shortward of the Si edge; therefore, a proposer could additionally explore what kinds of limits can be set on higher energy lines (e.g., Sulfur lines, using the MEG and HEG data).


5. Assessing Pile-up

Photon pile-up occurs when two or more photons land in the same, or neighboring, detector pixels in the same readout frame (see the Chandra ABC Guide to Pileup and the CIAO AHELP Pileup Page.) Pile-up is a concern for point sources viewed with the Chandra CCDs; however, it can also be a concern when using the gratings to view bright objects, such as X-ray binaries. Here we come up with a simple estimate for the peak pile-up fraction in these proposed observations.

Roughly speaking, the peak pile-up fraction (which translates to the fraction by which the observed count rate will be reduced in the observed spectrum) goes as the peak counts per 3 pixel region (along the gratings arm) per frame readout time. This can be determined graphically using a combination of Sherpa plot and Python Numpy functions (note the availability of the Sherpa "jdpileup" model for non-grating data).

The MEG detectors have the highest area, and the coarsest spectral resolution, of the HETG detectors, and hence typically have the greatest pile-up. Often this peak occurs close to the peak effective area of MEG, near 2 keV. However, given the very large neutral column to this particular source, the peak pile-up occurs at somewhat higher energies/shorter wavelengths (near 3 A/4 keV, as we show below).

The simulations indicate that the MEG+1 detector (dataset 4) has the highest count rate of the two MEG detectors; therefore, we assess the pile-up first for that detector as a "worst case" scenario. We use the Sherpa "plot_data" command to view the data plotted as counts/sec/Angstrom, in linear x- and y-scale. We access the x and y data arrays defining the plot with the "get_data_plot" command, and determine the wavelength and y-value corresponding to the peak of the plot using simple Python functions. We then multiply the y-value by a typical frame time (3.2 sec), 3 pixels (events are read in 3X3 pixel detector cells), and 0.011 (the number of Angstroms per MEG pixel), to arrive at an estimate of the peak pile-up fraction.


sherpa> set_analysis(4, "wavelength")     # Create plots as counts/sec/Angstrom
sherpa> plot_data(4)
sherpa> lin_scale(XY_AXIS)
sherpa> print_window("assess_pileup.ps")

sherpa> counts_array = get_data_plot(4).y  
sherpa> max_counts_y = numpy.max(counts_array)  # 719 counts
sherpa> max_counts_index = int(np.where(counts_array==max_counts_y)[0])
sherpa> wave_array = get_data_plot(4).x
sherpa> max_counts_x = wave_array[max_counts_index]

sherpa> print "x = "+str(max_counts_x), "y = "+str(max_counts_y)
'x = 2.80999780556', 'y = 0.569996470648'

The plot used for this assessment can be seen here (Figure 7). Thus, we find that for a full readout frame time of 3.2 sec, the peak pile-up fraction is expected to be 6%, and occurs near 2.9 Angstroms. Although not very severe, one can almost always reduce the frame time for gratings observations down to 1.7 sec by only using half of the chips along the chip readout direction. This loses none of the dispersed spectra, and further reduces the peak pile-up fraction to 3.4%. Given the very large neutral column, and hence the very low flux at long wavelengths (i.e., far along the dispersed arms, away from the center of the dispersed image), one could consider using an even smaller sub-array and further reducing the readout time and pile-up fraction.

[MEG1 Counts per Angstrom]
[Print media version: MEG1 Counts per Angstrom]

Figure 7: MEG1 Counts per Angstrom

Peak pile-up shown to occur near 2.9 Angstroms


6. Complete Target Form

For general instructions on how to submit a proposal, please see the Using RPS to Prepare and Submit a Chandra Proposal thread.

The above simulations provide details to use in the science justification (e.g., expected limits on Fe line wavelengths, physical widths, and strengths; estimates of the accuracy of the fitted neutral column), and can be used to fill out the proposal target forms.

For this observation, we will use the HETG detector with the ACIS-S chips. To minimize pile-up, we will use only 1/2 of the array. The 1st order count rate can be found by summing the counts from the simulations (use the "calc_data_sum" command, as shown above) and dividing by 65 ksec. This yields 4.2 counts/sec. For the total field rate, we note that approximately an additional 20% count rate will appear in the higher orders, and for a source dominated by counts > 2 keV, the 0th order count rate will be approximately the same as the total first order rate. (The 0th order image will suffer severely from pile-up, hence reducing the detected rate. Calculation of that effect is beyond the scope of this thread. By ignoring pile-up in the 0th order image, we obtain a "worst case" estimate of the count rate. For successful proposals, the Chandra gratings contact scientists will work with the proposer to optimize settings to minimize the count rate in the 0th order image.) Thus, as a very conservative estimate, we use a total field rate of 10 counts/sec.

This count rate is sufficiently low such that we can use "Faint Mode" for the telemetry format. If the source were several times brighter, "Graded Mode" would be required. Furthermore, pile-up would increase. Sufficiently bright sources would require the use of Continuous Clocking mode. Currently, the calibration is better understood for Timed Event mode, as opposed to Continuous Clocking mode. It is the recommendation of the gratings team that so long as estimates place the peak pile-up fraction at less than 10%, Timed Event mode should be used instead of Continuous Clocking mode (unless sub-1.7 sec time resolution is required to meet the science goals of the observation). For brighter sources, one should consider using Continuous Clocking mode. Users can contact the gratings team for further advice.

Given the low effective area of the frontside chips at low energies/long wavelengths, the S0 chip, covering the low energies of HEG-1/MEG-1, can often be designated as an "optional chip". This is especially true in this particular case where the large neutral column already ensures very little flux at low energy.

The RPS form should have the following parameter values for this observation. If a parameter isn't listed here, use the default RPS value or leave the field blank. Note that the detector and SIM offsets have been chosen to place a custom 1/2 subarray on the half of the chips nearest the readout. This not only reduces the readout time, but minimizes Charge Transfer Inefficiency (CTI). The gratings contact scientists can help optimize these choices for successful proposals.

  • SIM Translation Offset -- -6.14
  • Total Observing Time (ksec) -- 65 ksec
  • Count Rate -- 4.2
  • Total Field Count Rate -- 10
  • Exposure Mode -- TE
  • Event Telemetry Format -- Faint
  • CCDs On -- S0 = O1, S1-S5 = Y
  • Subarray Type -- Use Custom Subarray
  • Custom Subarray Start Row -- 1
  • Custom Subarray End Ros -- 512

Use this link to view the completed Target Form.


7. Further Resources

All of the steps taken in steps 1 through 5 above have been placed in the Python script found here (lmxb.py). Save the script as lmxb.py in the same directory as the Cycle 13 response files, start Sherpa in that directory, and then type:

sherpa> execfile("lmxb.py")

The script can be edited to change the values of the parameters, as well as the file format of saved plots.

For further information about installing and running Sherpa, see the Sherpa webpage on the CIAO website.


8. Thread Summary

This thread shows how to simulate a simple gratings spectrum of an X-ray binary, wherein the spectrum shows narrow Fe emission and absorption lines at energies > 6 keV. Furthermore, it shows how to estimate the constraints that the MEG place upon measurements of the neutral column. Simulating a 65 ksec spectrum, we found that HEG constrains the line positions and widths to within approximately 500 km/sec. The neutral column is constrained to within 20%. We looked at the estimated level of pile-up in the dispersed spectrum, and concluded that with a 1/2 sub-array, the peak pile-up fraction will be less than 4%. Accordingly, we choose a custom 1/2 sub-array to perform the 65 ksec observation.


Return to Threads Page

Last modified: 13 December 2012
Smithsonian Institute Smithsonian Institute

The Chandra X-Ray Center (CXC) is operated for NASA by the Smithsonian Astrophysical Observatory. 60 Garden Street, Cambridge, MA 02138 USA.   Email:   cxcweb@head.cfa.harvard.edu Smithsonian Institution, Copyright © 1998-2014. All rights reserved.