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

Proposer Threads (Cycle 16)

## Contents

**1. Thread Overview****2. Preliminary Considerations**-
**3. Simulate Source Spectrum** -
**4. Assessing the Simulations** **5. Assessing Pile-up****6. Complete Target Form****7. Further Resources****8. Thread Summary**-
**Images**

## 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_cy16.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_cy16.garf") sherpa> hpa = unpack_arf("aciss_heg1_cy16.garf") sherpa> mma = unpack_arf("aciss_meg-1_cy16.garf") sherpa> mpa = unpack_arf("aciss_meg1_cy16.garf") sherpa> hmr = unpack_rmf("aciss_heg-1_cy16.grmf") sherpa> hpr = unpack_rmf("aciss_heg1_cy16.grmf") sherpa> mmr = unpack_rmf("aciss_meg-1_cy16.grmf") sherpa> mpr = unpack_rmf("aciss_meg1_cy16.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_cy16.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_cy16.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_cy16.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_cy16.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_cy16.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_cy16.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_cy16.grmf detchans = 8192 RMF Data Set: 4:1 name = aciss_meg1_cy16.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_cy16.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_cy16.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_cy16.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_cy16.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_cy16.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_cy16.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_cy16.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_cy16.grmf detchans = 8192 RMF Data Set: 4:1 name = aciss_meg1_cy16.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_cy16.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.

## 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 16 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.