Chandra X-Ray Observatory (CXC)
Skip to the navigation links
Search the Sherpa website or contact the CXC HelpDesk
Last modified: 13 Dec 2012

URL: http://cxc.harvard.edu/sherpa/threads/aciss_grating_sim/thread.html
Return to Threads Page: Top | All | Fitting | Simulations

Simulating Chandra ACIS-S LETG Spectra with Sherpa

Sherpa Threads (CIAO 4.5 Sherpa v1)

Overview

Synopsis:

This thread illustrates the use of the Sherpa fake_pha command to simulate a spectrum of a point source as would be observed with a different Chandra instrument configuration than that used to produce the user's existing data and source model. We simulate ACIS-S/LETG +1/-1 LEG spectral orders using a source model expression tailored to the +1/-1 HEG and MEG orders of a real ACIS-S/HETG data set.

If you do not have experience with simulating X-ray spectral data in Sherpa, you may wish to follow the "step-by-step" example in the introductory simulation thread, before attempting the sample analysis in this thread.

Last Update: 13 Dec 2012 - reviewed for CIAO 4.5: group commands no longer clear the existing data filter

Contents

Getting started - downloading calibration response files for simulations

In order to simulate a Chandra ACIS-S/LETG spectrum with Sherpa, we must define an instrument response with the appropriate ARF (auxiliary response function) and RMF (redistribution matrix function) files. These files may be downloaded from the Chandra Proposal Planning page of the CalDB (Calibration Database) website, where sample RMFs and corresponding ARFs positioned at the aimpoint of the ACIS-S array (and at selected off-axis points) are available.

In this thread, we use the files aciss_leg1_cy13.grmf, aciss_leg-1_cy13.grmf, aciss_leg1_cy13.garf, and aciss_leg-1_cy13.garf, positioned at the default pointing for ACIS-S.

The Sherpa fake_pha command calculates a synthetic 1-D spectral data set based on a defined instrument response and source model expression. For details on the functionality of this command, with examples of its usage, see the threads "Simulating 1-D Data: the Sherpa FAKE_PHA Command" and "Simulating Chandra ACIS-S Spectra with Sherpa"

Sample ObsID used: 459 (HETG/ACIS-S, 3C 273)

The ACIS-S/HETG data files used to establish the appropriate source model expression for the simulated data were created by following several of the CIAO Grating threads:

The following files were loaded into Sherpa and analyzed according to the detailed steps shown here, in order to establish the source model expression to be used with fake_pha: 

spectra:
459_heg_m1_bin10.pha
459_heg_p1_bin10.pha
459_meg_m1_bin10.pha
459_meg_p1_bin10.pha

gARFs:
459_heg_m1.arf
459_heg_p1.arf
459_meg_m1.arf
459_meg_p1.arf

The spectra used in the fitting session were binned by a factor of 10. The data files are available in sherpa.tar.gz, as explained in the Sherpa Getting Started thread.

Defining the Instrument Response for the Simulation

Armed with a source model expression resulting from a satisfactory fit to the ACIS-S/HETG data, we may begin the simulation process by establishing the instrument response corresponding to the default pointing of the ACIS-S detector, using the ACIS-S/LETG grating ARF and RMF files retrieved from the CalDB website: 

sherpa> arf1=unpack_arf("aciss_leg1_cy13.garf")
sherpa> arf2=unpack_arf("aciss_leg-1_cy13.garf")

sherpa> print(arf1)
name     = aciss_leg1_cy13.garf
energ_lo = Float64[8192]
energ_hi = Float64[8192]
specresp = Float64[8192]
bin_lo   = Float64[8192]
bin_hi   = Float64[8192]
exposure = 28398.6297106


sherpa> print(arf2)
name     = aciss_leg-1_cy13.garf
energ_lo = Float64[8192]
energ_hi = Float64[8192]
specresp = Float64[8192]
bin_lo   = Float64[8192]
bin_hi   = Float64[8192]
exposure = 28398.4547152

sherpa> rmf1=unpack_rmf("aciss_leg1_cy13.grmf")
sherpa> rmf2=unpack_rmf("aciss_leg-1_cy13.grmf")

sherpa> print(rmf1)
name     = aciss_leg1_cy13.grmf
detchans = 8192
energ_lo = Float64[8192]
energ_hi = Float64[8192]
n_grp    = Int16[8192]
f_chan   = UInt32[8192]
n_chan   = UInt32[8192]
matrix   = Float64[827327]
offset   = 1
e_min    = Float64[8192]
e_max    = Float64[8192]

sherpa> print(rmf2)
name     = aciss_leg-1_cy13.grmf
detchans = 8192
energ_lo = Float64[8192]
energ_hi = Float64[8192]
n_grp    = Int16[8192]
f_chan   = UInt32[8192]
n_chan   = UInt32[8192]
matrix   = Float64[827353]
offset   = 1
e_min    = Float64[8192]
e_max    = Float64[8192]

Here, the ARF and RMF data are loaded and assigned to a set of variables with the unpack_arf and unpack_rmf commands. These variables will be used to assign the instrument response to the faked data sets we will create in the next section, "Running the Simulation with fake_pha."

Establishing the Source Model Expression for the Simulation

The fake_pha command has four required arguments: data set ID, ARF, RMF, and exposure time. In this thread, we simulate a LETG spectrum with the source model parameter values resulting from a fit to a HETG grating spectrum. The detailed steps taken to establish this model can be found here, and include:

  • Defining the Instrument Response for the Real Data
  • Filtering the Data
  • Defining the Source and Background Models
  • Examining Method & Statistic Settings
  • Fitting
  • Examining Fit Results

This saved Sherpa session, including all defined data sets, models, and instrument responses, may be restored with the Sherpa restore command, as follows: 

sherpa> restore("acis_hetg_fit_results.save")

Since we want to run the simulation using the source model tailored to the ACIS-S/HETG data, we use the set_source command to assign this source model expression to each of the LEG orders, as well as name the corresponding faked data sets which will be produced in our simulation.

The following steps assume that the ACIS-S/HETG spectral orders correspond to data set IDs 1 through 4. 

sherpa> show_source(1)
Model: 1
(atten.abs1 * bpl1d.bpow1)
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   abs1.hcol    frozen     1.81e+20        1e+17        1e+24           
   abs1.heiRatio frozen          0.1            0            1           
   abs1.heiiRatio frozen         0.01            0            1           
   bpow1.gamma1 thawed     0.424711          -10           10           
   bpow1.gamma2 thawed    0.0496577          -10           10           
   bpow1.eb     thawed        10.62            0         1000           
   bpow1.ref    frozen            1 -3.40282e+38  3.40282e+38           
   bpow1.ampl   thawed   0.00548927        1e-20  3.40282e+38    

sherpa> set_source("faked_leg1", get_source(1))
sherpa> set_source("faked_leg-1", get_source(1))

Running the Simulation with fake_pha

Simulating the Chandra ACIS-S/LETG spectra means taking the model expression defined for the ACIS-S/HETG data, folding it through the Chandra ACIS-S/LETG response, and applying Poisson noise to the counts predicted by the model. The simulation is run with fake_pha, which has four required arguments: data set ID, ARF, RMF, and exposure time. In this thread, we will simulate ACIS-S/LETG spectra using the gARF and gRMF files downloaded from the CalDB website, a source model tailored to the ACIS-S/HETG data, and an exposure time of 50 ks.

We are now ready to run the simulations with fake_pha: 

sherpa> fake_pha("faked_leg1", arf1, rmf1, exposure=50000)
sherpa> fake_pha("faked_leg-1", arf2, rmf2, exposure=50000)

These commands associate the data IDs "faked_leg1" and "faked_leg-1" with simulated source data sets for the +1 and -l ACIS-S/LETG orders, respectively, based on the assumed exposure time, instrument responses, and source model expression; Poisson noise is added to the modeled data.

Note that as of Sherpa in CIAO 4.2, the 'arf' and 'rmf' arguments of the fake_pha command can accept filenames directly; e.g., we could have done the following: 

sherpa> fake_pha("faked_leg1",  arf="aciss_leg1_cy13.garf", rmf="aciss_leg1_cy13.grmf", exposure=50000)
sherpa> fake_pha("faked_leg-1", arf="aciss_leg-1_cy13.garf", rmf="aciss_leg-1_cy13.grmf", exposure=50000)

For detailed information on the available options for setting the 'arf' and 'rmf' arguments of fake_pha, refer to the fake_pha ahelp file.

We may inspect some basic properties of the new data sets with the show_data command; e.g.,: 

sherpa> show_data("faked_leg1")
Data Set: faked_leg1
Filter: 0.1199-12.3219 Energy (keV)
Noticed Channels: 1-8192
header         = None
name           = faked
channel        = Int32[8192]
counts         = Float64[8192]
staterror      = None
syserror       = None
bin_lo         = None
bin_hi         = None
grouping       = None
quality        = None
exposure       = 50000
backscal       = None
areascal       = None
grouped        = False
subtracted     = False
units          = energy
rate           = True
plot_fac       = 0
response_ids   = [1]
background_ids = []

RMF Data Set: faked_leg1:1
name     = aciss_leg1_cy13.grmf
detchans = 8192
energ_lo = Float64[8192]
energ_hi = Float64[8192]
n_grp    = Int16[8192]
f_chan   = UInt32[8192]
n_chan   = UInt32[8192]
matrix   = Float64[827327]
offset   = 1
e_min    = Float64[8192]
e_max    = Float64[8192]

ARF Data Set: faked_leg1:1
name     = aciss_leg1_cy13.garf
energ_lo = Float64[8192]
energ_hi = Float64[8192]
specresp = Float64[8192]
bin_lo   = Float64[8192]
bin_hi   = Float64[8192]
exposure = 28398.6297106

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

Writing the Simulated Data to Output Files

We may use the save_pha command to write the simulated data as a PHA file, with a header containing the exposure time value and paths to the ARF and RMF files: 

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

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

sherpa> get_data("faked_leg-1").header = {}               # Make header an empty dictionary
sherpa> get_data("faked_leg-1").name = "simulation-1.pha"  # Update name for save_pha
sherpa> save_pha("faked_leg-1", "simulation-1.pha")

We also have the option to save the data to an ASCII table: 

sherpa>save_arrays("my_sim_data.fits", [get_model_plot("faked_leg1").xlo, get_model_plot("faked_leg1").y], ascii=False)

sherpa>save_arrays("my_sim_data.txt", [get_model_plot("faked_leg1").xlo, get_model_plot("faked_leg1").y], ascii=True)

Fitting the Simulated Data

The simulated data set may be filtered and fit as any other data set in Sherpa. For example, we can filter the simulated data in energy space to include only the counts in the restricted energy range, such as 5.0 keV - 11.0 keV.: 

sherpa> set_analysis("energy")

sherpa> print calc_energy_flux(5.,11., id="faked_leg1")  # ergs cm-2 s-1
1.84657469546e-12


sherpa> set_analysis("wave")
sherpa> ignore()
sherpa> notice(5.,11.)

sherpa> show_filter()
Data Set Filter: 1
4.9925-10.9925 Wavelength (Angstrom)

Data Set Filter: 2
4.9925-10.9925 Wavelength (Angstrom)

Data Set Filter: 3
4.9849-10.9849 Wavelength (Angstrom)

Data Set Filter: 4
4.9849-10.9849 Wavelength (Angstrom)

Data Set Filter: faked_leg-1
4.9937-10.9937 Wavelength (Angstrom)

Data Set Filter: faked_leg1
4.9937-10.9937 Wavelength (Angstrom)

Then, we can fit the simulated data with the source model expression we used to create it: 

sherpa> set_method("neldermead")
sherpa> set_stat("cstat")

sherpa> fit("faked_leg1", "faked_leg-1")
Datasets              = 'faked_leg1', 'faked_leg-1'
Method                = neldermead
Statistic             = cstat
Initial fit statistic = 27813.5
Final fit statistic   = 955.658 at function evaluation 852
Data points           = 962
Degrees of freedom    = 958
Probability [Q-value] = 0.515284
Reduced statistic     = 0.997555
Change in statistic   = 26857.9
   bpow1.gamma1   1.5427      
   bpow1.gamma2   10          
   bpow1.eb       10.7629     
   bpow1.ampl     0.022249    


sherpa> plot("fit","faked_leg1","fit","faked_leg-1")

The resulting plot is shown in Figure 1.

Figure 1: Plot of fit to simulated ACIS-S/LETG source spectra

[Plot of fit to simulated spectra]
[Print media version: Plot of fit to simulated spectra]

Figure 1: Plot of fit to simulated ACIS-S/LETG source spectra

Plot of fit to simulated ACIS-S/LETG +/-1 spectral orders, using a source model tailored to a real ACIS-S/HETG data set.

Next, we examine the quality of the fit with the confidence method (conf), and print the fit and confidence results with the show_fit and get_conf_results commands, respectively. 

sherpa> conf("faked_leg1", "faked_leg-1")
bpow1.gamma2 lower bound:       -4.66522
bpow1.eb lower bound:   -0.0920633
bpow1.gamma1 lower bound:       -0.0321875
bpow1.gamma2 upper bound:       14.7113
bpow1.ampl lower bound: -0.00137792
bpow1.ampl upper bound: 0.00146978
bpow1.gamma1 upper bound:       0.0322129
bpow1.eb upper bound:   0.0678082
Datasets              = 'faked_leg1', 'faked_leg-1'
Confidence Method     = confidence
Fitting Method        = neldermead
Statistic             = cstat
confidence 1-sigma (68.2689%) bounds:
   Param            Best-Fit  Lower Bound  Upper Bound
   -----            --------  -----------  -----------
   bpow1.gamma1       1.5427   -0.0321875    0.0322129
   bpow1.gamma2           10     -4.66522      14.7113
   bpow1.eb          10.7629   -0.0920633    0.0678082
   bpow1.ampl       0.022249  -0.00137792   0.00146978

sherpa> show_fit()
Optimization Method: NelderMead
name         = simplex
ftol         = 1.19209289551e-07
maxfev       = None
initsimplex  = 0
finalsimplex = 9
step         = None
iquad        = 1
verbose      = 0

Statistic: CStat
Maximum likelihood function (XSPEC style)

Fit:Datasets              = 'faked_leg1', 'faked_leg-1'
Method                = neldermead
Statistic             = cstat
Initial fit statistic = 27813.5
Final fit statistic   = 955.658 at function evaluation 852
Data points           = 962
Degrees of freedom    = 958
Probability [Q-value] = 0.515284
Reduced statistic     = 0.997555
Change in statistic   = 26857.9
   bpow1.gamma1   1.5427      
   bpow1.gamma2   10          
   bpow1.eb       10.7629     
   bpow1.ampl     0.022249    

sherpa> print get_conf_results()
datasets   = ('faked_leg1', 'faked_leg-1')
methodname = confidence
fitname    = neldermead
statname   = cstat
sigma      = 1
percent    = 68.2689492137
parnames   = ('bpow1.gamma1', 'bpow1.gamma2', 'bpow1.eb', 'bpow1.ampl')
parvals    = (1.5426960872549176, 9.9999998155233776, 10.762852399012568, 0.022249049856703392)
parmins    = (-0.032187515403694134, -4.6652191458098962, -0.092063289778602453, -0.001377917523524206)
parmaxes   = (0.032212924793423303, 14.71133964094448, 0.067808230777638201, 0.0014697786917591514)
nfits      = 39

Note that the cstat statistic is appropriate for simultaneously fitting source and background data, as well as low-counts data, but it does not calculate errors for the data points. We can group the data so that each bin contains a specified minimum number of counts, and then change the fit statistic to something more suitable to calculate the errors. Finally, we can view the results of the new fit with the plot_fit_delchi command: 

sherpa> group_counts("faked_leg1", 15)

sherpa> group_counts("faked_leg-1", 15)

sherpa> set_stat("chi2xspecvar")

sherpa> plot_fit_delchi("faked_leg1")
sherpa> plot_fit_delchi("faked_leg-1")

The new fits to the grouped, simulated source spectra, along with the residuals divided by the uncertainties, are shown in Figure 2 and Figure 3.

Figure 2: Plot of fit to simulated ACIS-S/LETG +1 order spectrum, with residuals

[]
[Print media version: ]

Figure 2: Plot of fit to simulated ACIS-S/LETG +1 order spectrum, with residuals

Figure 3: Plot of fit to simulated ACIS-S/LETG -1 order spectrum, with residuals

[]
[Print media version: ]

Figure 3: Plot of fit to simulated ACIS-S/LETG -1 order spectrum, with residuals

The plots may be saved as PostScript files accordingly: 

sherpa> plot_fit_delchi("faked_leg1")
sherpa> print_window("sim_fit_p1")

sherpa> plot_fit_delchi("faked_leg-1")
sherpa> print_window("sim_fit_m1")

Scripting It

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

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

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

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

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

History

02 Apr 2009 Created for CIAO/Sherpa 4.1
28 Apr 2009 replaced use of atten model with Sherpa user model "atten_wave"; new script command is available with CIAO 4.1.2
13 Jan 2010 updated for CIAO 4.2; new calibration response files used in simulation
13 Jul 2010 updated for CIAO 4.2 Sherpa v2: removal of S-Lang version of thread.
15 Dec 2011 reviewed for CIAO 4.4: a work-around for a save_pha bug was added; response files used in simulations were updated for Chandra proposal cycle 13
13 Dec 2012 reviewed for CIAO 4.5: group commands no longer clear the existing data filter
Return to Threads Page: Top | All | Fitting | Simulations

Last modified: 13 Dec 2012
CXC logo

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-2012. All rights reserved.