Chandra X-Ray Observatory
Skip to the navigation links
Last modified: 11 December 2013

URL: /data/cdoweb/Evan_test2/proposer/threads/binary/

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

Proposer Threads (Cycle 16)


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 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 an ISIS simulation. It you prefer to work with Sherpa please refer to the companion Sherpa thread.

2. Preliminary Considerations

Here are some things to consider when following this thread.
  • You need to have the Interactive Spectral Interpretation System ( ISIS) installed. This can be installed from source code, or from binary. Similar simulation procedures can be followed using Sherpa or XSPEC; however, ISIS incorporates in a programmable manner access to atomic databases that could be useful in subsequent analysis, it allows for convenient in memory addition and regrouping of datasets, and should photon pile-up be deemed an issue, it allows for modeling of gratings pile-up. (See the Chandra ABC Guide to Pileup and the CIAO Pileup AHELP Page for further details.)
    ISIS can be programmed using the S-lang interpretted language ( . This thread is saved as a script which a user can simply edit and customize, and rerun for their own source
  • Download the latest instrument responses. The example and figure presented below used Cycle 13 responses and ISIS 1.6.1-52 for illustration purposes. There might be some minor numerical differences when using the latest responses and some slight cosmetic differences in the ISIS prompts.
    Note: the linked script references Cycle 15 response files - e.g., the script has a command
    as opposed to the example
    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 (
  • 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, The easist way to follow this thread is to cut and paste from the script, or edit the script and run at the ISIS prompt:

unix% isis

Welcome to ISIS Version 1.6.1-52
Copyright (C) 1998-2011 Massachusetts Institute of Technology

          Isis web page:
   Mailing list archive:
 Send questions to the mailing list: <>.
     For a summary of recent changes, type:  "help changes"

isis> require("xspec");
isis> .load

Run ISIS and read in responses

The ISIS binary contains all the current XSPEC models, and if one builds ISIS from source, ISIS also can have access to any XSPEC local models (see the ISIS web page,, for details.) In this thread, we will only use standard XSPEC spectral models.

unix% isis

Welcome to ISIS Version 1.6.1-52
Copyright (C) 1998-2011 Massachusetts Institute of Technology

          Isis web page:
   Mailing list archive:
 Send questions to the mailing list: <>.
     For a summary of recent changes, type:  "help changes"

isis> require("xspec");

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, assign them to data set numbers 1-4, and then set the ARF exposure times to 650 ksec. (We also set the data exposure values to 650 ksec, although this is not strictly necessary - the fake data file integration time are set by the ARF exposures.) We choose 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.

isis> hma = load_arf("aciss_heg-1_cy13.garf");
isis> hpa = load_arf("aciss_heg1_cy13.garf");
isis> mma = load_arf("aciss_meg-1_cy13.garf");
isis> mpa = load_arf("aciss_meg-1_cy13.garf");

isis> hmr = load_rmf("aciss_heg-1_cy13.grmf");
isis> hpr = load_rmf("aciss_heg1_cy13.grmf");
isis> mmr = load_rmf("aciss_meg-1_cy13.grmf");
isis> mpr = load_rmf("aciss_meg1_cy13.grmf");

isis> assign_arf(hma,1);  % Dataset 1 = HEG-1
isis> assign_arf(hpa,2);  % Dataset 2 = HEG+1
isis> assign_arf(mma,3);  % Dataset 3 = MEG-1
isis> assign_arf(mpa,4);  % Dataset 4 = MEG+1

isis> assign_rmf(hmr,1);  % Dataset 1 = HEG-1
isis> assign_rmf(hpr,2);  % Dataset 2 = HEG+1
isis> assign_rmf(mmr,3);  % Dataset 3 = MEG-1
isis> assign_rmf(mpr,4);  % Dataset 4 = MEG+1

isis> set_arf_exposure([hma,hpa,mma,mpa],6.5e5);
isis> set_data_exposure(1,6.5e5);
isis> set_data_exposure(2,6.5e5);
isis> set_data_exposure(3,6.5e5);
isis> set_data_exposure(4,6.5e5);

Define the Model

ISIS for the most part uses model component names identical to those found in XSPEC. ISIS, however, does not allow the names to be abbreviated, and appends to each name "(#)", where # is an integer. By using different values of #, multiple model components (e.g., lines) of the same type can be introduced. If each model component is viewed as a vector of values (i.e., counts in each wavelength bin), then nearly any "mathematically sensible" combination of model components will be properly treated by ISIS as a complete 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 ISIS "fit_fun" command:

isis> fit_fun("phabs(1)*(bbodyrad(1)+powerlaw(1)+gaussian(1)-gaussian(2)-gaussian(3))");

Invoking the ISIS "edit_par" command will open an editor window (e.g., using emacs or vi) wherein individual parameter values can be entered. Alternatively, one can enter parameter values on the command line with the ISIS "set_par" command. We take parameter values from the fits to the XMM-Newton data, and set the parameter minimum and maximum values to reasonable limits for subsequent fits (in order to prevent, for example, gaussians becoming very broad and fitting the continuum):

isis> set_par("phabs(1).nH", 0., 0, 0, 20);
isis> set_par("bbodyrad(1).norm", 41., 0, 0, 100);
isis> set_par("bbodyrad(1).kT", 1.18, 0, 0.5, 3);
isis> set_par("powerlaw(1).PhoIndex", 2.02, 0, 1, 3);
isis> set_par("gaussian(1).norm", 3.e-5, 0, 0, 0.01);
isis> set_par("gaussian(1).LineE", 6.4, 0, 6, 6.7);
isis> set_par("gaussian(1).Sigma", 0.1, 0, 0, 1);
isis> set_par("gaussian(2).norm", 3.e-5, 0, 0, 0.01);
isis> set_par("gaussian(2).LineE", 6.72, 0, 6.4, 6.9);
isis> set_par("gaussian(2).Sigma", 0.02, 0, 0, 0.1);
isis> set_par("gaussian(3).norm", 3.e-5, 0, 0, 0.01);
isis> set_par("gaussian(3).LineE", 7.00, 0, 6.7, 7.3);
isis> set_par("gaussian(3).Sigma", 0.02, 0, 0, 1);

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 fakeit, and Set Model Normalizations

Having defined a model (which, by default, is the same for all four gratings arms), we now create a fake dataset (complete with Poisson noise) with the ISIS "fakeit" command. The exposure time has already been set by assigning that exposure time to the respective ARFs. One then simply has to type "fakeit", followed by a "list_data" to verify that the datasets have indeed been created.

isis> fakeit;
isis> list_data;
Current Spectrum List:
 id    instrument   m prt src    use/nbins   A   R     totcts   exp(ksec)  target
  1                -1  0   0    8192/ 8192   1   1  2.9792e+06   650.000  
  2                 1  0   0    8192/ 8192   2   2  2.8209e+06   650.000  
  3                -1  0   0    8192/ 8192   3   3  5.5118e+06   650.000  
  4                 1  0   0    8192/ 8192   4   4  5.5037e+06   650.000  

The HEG detectors (datasets 1 and 2) have less effective area, and hence fewer counts, than the MEG detectors (datasets 3 and 4).

ISIS allows a great deal of customization. Here to help us set the normalizations, we define a function, "flux", that will return the flux (both as a print to the screen, and as a variable) in units of photons and ergs per cm2/sec. We also define a function, "eqw", that similarly returns line equivalent widths (for a given detector and model component normalization) in units of mAngstrom and eV.

isis> public define flux (id, e_lo, e_hi)
   variable l_lo, l_hi, m, e_avg, pflux, eflux;

   l_hi = _A(e_lo);        % Convert keV to A, or visa versa
   l_lo = _A(e_hi);        % Convert keV to A, or visa versa
   m = get_model_flux(id);
   pflux = rebin (l_lo, l_hi, m.bin_lo, m.bin_hi, m.value)[0];   
   e_avg = _A(1.)*(1./m.bin_hi+1./m.bin_lo)/2.*1.602e-9;
   eflux = rebin (l_lo, l_hi, m.bin_lo, m.bin_hi, m.value*e_avg)[0];   

   () = printf("\n photons/cm2/sec: %11.3e, ergs/cm2/sec: %11.3e \n \n", pflux, eflux);

   return pflux, eflux;  % Photon flux & ergs flux (per cm^2/sec) 

isis> public define eqw(indx,par)
   variable eva, eve;
   variable fv_hold = Fit_Verbose, pmin=0;

   Fit_Verbose = -1;
   () = eval_counts;

   variable a = get_model_flux(indx);
   variable phold = get_par_info(par);
   () = eval_counts;
   variable b = get_model_flux(indx);

   () = eval_counts;

   variable iw = [0:length(a.bin_lo)-1];

   variable diff = a.value - b.value;
   variable bdena = b.value/(a.bin_hi-a.bin_lo);
   variable bdene = reverse(reverse(b.value)/(_A(a.bin_lo)-_A(a.bin_hi)));
   variable fdena = diff/(a.bin_hi-a.bin_lo);
   variable fdene = reverse(reverse(diff)/(_A(a.bin_lo) - _A(a.bin_hi)));
   variable flux = sum( diff[iw] );
   variable iwa = min(where( abs(fdena[iw]) == max(abs(fdena[iw])) ));
   variable iwe = min(where( abs(fdene[iw]) == max(abs(fdene[iw])) ));
   eva = flux/bdena[iw[iwa]]*1000.;
   eve = flux/bdene[iw[iwa]]*1000.;

   Fit_Verbose = fv_hold;
   () = printf("\n EqW (mA): %11.3e, EqW (eV): %11.3e \n \n", eva, eve);

   return eva, eve;

isis> (pflux,eflux) = flux(3,0.1,8);  % Meg -1 flux, from 0.1 to 8 keV

 photons/cm2/sec:   6.530e-01, ergs/cm2/sec:   1.566e-09 
isis> (g1_a,g1_e) = eqw(1,"gaussian(1).norm"); % Heg -1 Equivalent width

 EqW (mA):   8.094e-01, EqW (eV):   2.677e+00 
isis> (g2_a,g2_e) = eqw(1,"gaussian(2).norm"); % Heg -1 Equivalent width

 EqW (mA):  -8.533e-01, EqW (eV):  -3.112e+00
isis> (g3_a,g3_e) = eqw(1,"gaussian(3).norm"); % Heg -1 Equivalent width

 EqW (mA):  -8.960e-01, EqW (eV):  -3.541e+00

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 fakeit again, and then verify that the new equivalent widths are those that we desire.

isis> set_par("gaussian(1).norm",get_par("gaussian(1).norm")*95/g1_e);
isis> set_par("gaussian(2).norm",get_par("gaussian(2).norm")*(-14)/g2_e);
isis> set_par("gaussian(3).norm",get_par("gaussian(3).norm")*(-23)/g3_e);

isis> fakeit;
isis> (g1_a,g1_e) = eqw(1,"gaussian(1).norm"); % Heg -1 Equivalent width

 EqW (mA):   2.828e+01, EqW (eV):   9.500e+01 
isis> (g2_a,g2_e) = eqw(1,"gaussian(2).norm"); % Heg -1 Equivalent width

 EqW (mA):  -3.571e+00, EqW (eV):  -1.302e+01 
isis> (g3_a,g3_e) = eqw(1,"gaussian(3).norm"); % Heg -1 Equivalent width

 EqW (mA):  -5.604e+00, EqW (eV):  -2.215e+01 

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 fakeit once more. Furthermore, we save a parameter file, lmxb.par, with the final model.

isis> set_par("phabs(1).nH",7.5);
isis> save_par("lmxb.par");

isis> set_arf_exposure([hma,hpa,mma,mpa],6.5e4);
isis> set_data_exposure(1,6.5e4);
isis> set_data_exposure(2,6.5e4);
isis> set_data_exposure(3,6.5e4);
isis> set_data_exposure(4,6.5e4);
isis> fakeit;
isis> list_data;
Current Spectrum List:
 id    instrument part/m  src    use/nbins   A   R     totcts   exp(ksec)  target
  1                   -1   0    8192/ 8192   1   1  6.9668e+04    65.000  
  2                   +1   0    8192/ 8192   2   2  7.1021e+04    65.000  
  3                   -1   0    8192/ 8192   3   3  5.7348e+04    65.000  
  4                   +1   0    8192/ 8192   4   4  5.7378e+04    65.000  

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. This is an "in memory" combination only. ISIS allows such regroupings of data and combinations of data during analysis in order to facilitate easier exploration of the optimal groupings.

isis> group_data([1,2]; min_chan=2);   % Group HEG data by a factor of 2
isis> group_data([3,4], min_chan=2);   % Group MEG data by a factor of 4
isis> heg = combine_datasets(1,2);     % In memory combination of HEG data
isis> meg = combine_datasets(3,4);     % In memory combination of MEG data

isis> notice_values([1,2],1.65,2.3; unit="a"); % Only notice the 1.65-2.3 A data for HEG
isis> exclude(3,4);                            % Ignore MEG data for now

isis> freeze([1,3,5]);                         % Freeze Nh, blackbody kT, Photon Index

isis> fit_counts;
 Parameters[Variable] = 14[11]
            Data bins = 130
           Chi-square = 118.9635
   Reduced chi-square = 0.9996931
isis> list_par;
 idx  param             tie-to  freeze         value         min         max
  1  phabs(1).nH            0     1              7.5           0          20  10^22
  2  bbodyrad(1).norm       0     0         40.88613           0         100  
  3  bbodyrad(1).kT         0     1             1.18         0.5           3  keV
  4  powerlaw(1).norm       0     0        0.1447135           0          10  
  5  powerlaw(1).PhoIndex   0     1             2.02           1           3  
  6  gaussian(1).norm       0     0     0.0009578358           0        0.01  
  7  gaussian(1).LineE      0     0         6.394854           6         6.7  keV
  8  gaussian(1).Sigma      0     0        0.1082421           0           1  keV
  9  gaussian(2).norm       0     0     0.0001736271           0        0.01  
 10  gaussian(2).LineE      0     0         6.717939         6.4         6.9  keV
 11  gaussian(2).Sigma      0     0       0.03438081           0         0.1  keV
 12  gaussian(3).norm       0     0     0.0001914165           0        0.01  
 13  gaussian(3).LineE      0     0         7.000227         6.7         7.3  keV
 14  gaussian(3).Sigma      0     0       0.02443937           0           1  keV

isis> save_par("heg_fit.par");  % Save the HEG fit parameters

isis> putenv("PGPLOT_BACKGROUND=white");
isis> putenv("PGPLOT_FOREGROUND=black");  % Make plots black on white

isis> set_line_width(3); set_frame_line_width(3); charsize(1.1); % Nicer plot defaults
isis> xrange(1.65,2.3);  % Plot x-range to only show what we've noticed
isis> ylog;              % Logarithmic y-axis

isis> rplot_counts(1);  % HEG-1 data (although fit includes HEG+1);
isis> rplot_counts(2);  % HEG+1 data (although fit includes HEG-1);

isis> hd = get_combined(heg, &get_data_counts);  % Get the summed data
isis> hm = get_combined(heg, &get_model_counts); % Get the summed model
isis> label("Angstrom","Counts/Bin","Combined HEG Data");

isis> hplot(hd.bin_lo,hd.bin_hi,hd.value);
isis> ohplot(hm.bin_lo,hm.bin_hi,hm.value);

The resulting plots are seen here (Figure 1), here (Figure 2), and here (Figure 3). Note that the ISIS commands for generating hardcopies of plots follows the format of:

isis> id = open_plot(""); % Create a color postscript file
isis> resize(21,0.85);   % 21 cm, with an aspect ratio of 0.85:1
... plot commands ...
isis> close_plot(id);

We now turn to determining the confidence bounds for the line parameters found above. For this, we use the ISIS "conf" command, which defaults to 90% confidence limits. 68% or 99% can also be specified. (Arbritrary values of the confidence level can be specified with the "fconf" command.) The conf command will return the upper and lower parameter values for these confidence intervals.

isis> conf("gaussian(1).norm");
isis> conf("gaussian(1).LineE");
isis> conf("gaussian(1).Sigma");
isis> conf("gaussian(2).norm");
isis> conf("gaussian(2).LineE");
isis> conf("gaussian(2).Sigma");
isis> conf("gaussian(3).norm");
isis> conf("gaussian(3).LineE");
isis> conf("gaussian(3).Sigma");

We find that the line positions are all determined with 90% confidence level uncertainties of 0.01-0.02 keV, and likewise the line widths are determined with confidences of 0.01-0.02 keV. In terms of velocities, this corresponds to approximately 500-1000 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.

isis> set_fit_method("subplex");  % Slow fit method, but often useful for gratings data

isis> (g1_a,g1_e) = eqw(1,"gaussian(1).norm");  

 EqW (mA):   2.616e+01, EqW (eV):   8.653e+01 
isis> set_par("gaussian(1).norm",0.000803,1);
() = fit_counts();
 Parameters[Variable] = 14[10]
            Data bins = 130
           Chi-square = 121.6589
   Reduced chi-square = 1.013824

isis> (g1_a,g1_e) = eqw(1,"gaussian(1).norm");

 EqW (mA):   2.176e+01, EqW (eV):   7.177e+01 
isis> set_par("gaussian(1).norm",0.00114,1);

isis> () = fit_counts();
 Parameters[Variable] = 14[10]
            Data bins = 130
           Chi-square = 121.6729
   Reduced chi-square = 1.013941

isis> (g1_a,g1_e) = eqw(1,"gaussian(1).norm");

 EqW (mA):   3.127e+01, EqW (eV):   1.037e+02 
isis> load_par("heg_fit.par");  % Restore the HEG fit parameters

isis> (g2_a,g2_e) = eqw(1,"gaussian(2).norm"); 

 EqW (mA):  -4.976e+00, EqW (eV):  -1.810e+01 
isis> set_par("gaussian(2).norm",0.00000983,1);

isis> () = fit_counts();
 Parameters[Variable] = 14[10]
            Data bins = 130
           Chi-square = 137.6527
   Reduced chi-square = 1.147106
isis> (g2_a,g2_e) = eqw(1,"gaussian(2).norm");

 EqW (mA):  -2.868e-01, EqW (eV):  -1.043e+00 
isis> set_par("gaussian(2).norm",0.00037,1);

isis> () = fit_counts();
 Parameters[Variable] = 14[10]
            Data bins = 130
           Chi-square = 277.286
   Reduced chi-square = 2.310717

isis> (g2_a,g2_e) = eqw(1,"gaussian(2).norm");

 EqW (mA):  -1.036e+01, EqW (eV):  -3.767e+01 
isis> (g3_a,g3_e) = eqw(1,"gaussian(3).norm");

 EqW (mA):  -5.799e+00, EqW (eV):  -2.292e+01 
isis> set_par("gaussian(3).norm",0.000130,1);

isis> () = fit_counts();
 Parameters[Variable] = 14[10]
            Data bins = 130
           Chi-square = 121.6507
   Reduced chi-square = 1.013756

isis> (g3_a,g3_e) = eqw(1,"gaussian(3).norm");

 EqW (mA):  -3.984e+00, EqW (eV):  -1.574e+01 
isis> set_par("gaussian(3).norm",0.000263,1);

isis> () = fit_counts();
 Parameters[Variable] = 14[10]
            Data bins = 130
           Chi-square = 121.6708
   Reduced chi-square = 1.013923

isis> (g3_a,g3_e) = eqw(1,"gaussian(3).norm");

 EqW (mA):  -7.861e+00, EqW (eV):  -3.107e+01 

In the above, we find the fitted line equivalent width to be: 87+17/-15 eV, -18-20/+8 eV, and -23-8/+7 eV, respectively. The original input values are encompassed within these error bars, but as expected from a Monte Carlo simulation, they are not exact matches.

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. (ISIS will always remember the last used parameter values for any model component, which makes it easy to subtract and add these components in a useful manner.)

isis> exclude([1,2]);                         % Ignore HEG
isis> include(3,4);                           % Notice MEG
isis> notice_values([3,4],2.1,7.2; unit="a"); % Only notice 1.7-7.2 Angstrom (1.7-5.9 keV)

isis> load_par("lmxb.par");  % Simulation parameters, including lines
isis> () = eval_counts;
 Parameters[Variable] = 14[14]
            Data bins = 509
           Chi-square = 444.8406
   Reduced chi-square = 0.898668

isis> fit_fun("phabs(1)*(bbodyrad(1)+powerlaw(1))"); % Don't fit lines
isis> () = fit_counts;
 Parameters[Variable] = 5[5]
            Data bins = 509
           Chi-square = 429.7703
   Reduced chi-square = 0.8527188

isis> list_par;
 idx  param             tie-to  freeze         value         min         max
  1  phabs(1).nH            0     0         8.107716           0          20  10^22
  2  bbodyrad(1).norm       0     0         34.06908           0         100  
  3  bbodyrad(1).kT         0     0         1.217122         0.5           3  keV
  4  powerlaw(1).norm       0     0        0.3114823           0          10  
  5  powerlaw(1).PhoIndex   0     0         2.341737           1           3  

isis> save_par("meg_fit.par");  % Save the MEG fit parameters

isis> conf("phabs(1).nH");

isis> title("");       % We set the plot title in "label" above, so remove it
isis> xrange(2.1,7.2); % Only plot 2.1-7.2 Angstroms
isis> rplot_counts(3); % MEG-1 data (although fit combined with MEG+1)
isis> rplot_counts(4); % MEG+1 data (although fit combined with MEG-1)
isis> flux_corr(3);      % Calculated the "unfolded" spectrum for MEG-1
isis> plot_data_flux(3); % Plot the MEG-1 unfolded spectrum 

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.0(+1.4)(-0.3)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 ISIS.

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 ISIS "plot_bin_density" plotting option to see the data plotted as counts/sec/Angstrom, we plot with a linear y-scale, and we use the ISIS "cursor" function, which draws interactive cross hairs on the plot. The cursor functions lets us click on the peak of the plot, and then prints the corresponding wavelength and y-value. (Note: the cursor function should be used on a linear-linear scale plot in this case.) 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.

isis> ylin;                % Linear y-axis
isis> plot_bin_density;    % Create plots as counts/sec/Angstrom
isis> plot_data_counts(4); % MEG+1 is the brightest arm
isis> cursor;              % Measure with the pgplot cursor
Click for cursor coordinates, type 'q' in plot window to quit

x=  2.917383e+00  y=  5.821421e-01  ch=D

isis> 0.58*3.2*3*(0.011); % cts/s/A * 3.2 sec frame * 3 pixels * 0.011 A/pixel

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 "list_data" 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 above have been placed in this S-lang script found here ( Save the script as in the same directory as the Cycle 13 response files, start ISIS in that directory, and then type:

isis> .load

The script can be edited to change the values of the parameters, and to save the plots as postscript files.

For further information about installing and running ISIS, see:

For useful scripts, specifically utilities that will aid those more familiar with XSPEC, and useful plotting scripts (including ones that can handle both combined data and unfolded data), start at:

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 500-1000 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.

Last modified: 11 December 2013
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: Smithsonian Institution, Copyright © 1998-2014. All rights reserved.