Chandra X-Ray Observatory
	(CXC)
Skip to the navigation links
Last modified: 6 Dec 2016

URL: http://cxc.harvard.edu/sherpa/threads/aciss_grating_sim/

Simulating Chandra ACIS-S LETG Spectra with Sherpa

Sherpa Threads (CIAO 4.9 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: 6 Dec 2016 - updated for Chandra proposal cycle 19; include defining source model found in the restored state file, that would fail to load in newer versions of Sherpa.


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 (ancillary 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_cy19.grmf, aciss_leg-1_cy19.grmf, aciss_leg1_cy19.garf, and aciss_leg-1_cy19.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_cy19.garf")
sherpa> arf2=unpack_arf("aciss_leg-1_cy19.garf")

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

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

sherpa> rmf1=unpack_rmf("aciss_leg1_cy19.grmf")
sherpa> rmf2=unpack_rmf("aciss_leg-1_cy19.grmf")

sherpa> print(rmf1)
detchans = 8192
energ_lo = Float64[8192]
energ_hi = Float64[8192]
n_grp    = UInt64[8192]
f_chan   = UInt32[8192]
n_chan   = UInt32[8192]
matrix   = Float64[826552]
offset   = 1
e_min    = Float64[8192]
e_max    = Float64[8192]

sherpa> print(rmf2)
name     = aciss_leg-1_cy19.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[826506]
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")

If restore fails to load the state file, due to version incompatibility, the source model defined in the file can be loaded into Sherpa:

sherpa> set_source(1,atten.abs1*bpl1d.bpow1)

sherpa> abs1.hcol=1.81e+20
sherpa> abs1.heiRatio=0.1
sherpa> abs1.heiiRatio=0.01
sherpa> bpow1.gamma1=0.424711
sherpa> bpow1.gamma2=0.0496577
sherpa> bpow1.eb=10.62      
sherpa> bpow1.ref=1
sherpa> bpow1.ampl=0.00548927

sherpa> freeze(abs1)

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 in the saved state file.

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.7 the WARNING message is printed out when the previously saved sherpa session is being used.

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_cy19.garf", rmf="aciss_leg1_cy19.grmf", exposure=50000)
sherpa> fake_pha("faked_leg-1", arf="aciss_leg-1_cy19.garf", rmf="aciss_leg-1_cy19.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
name           = faked
channel        = Float64[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_cy18.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[826552]
offset   = 1
e_min    = Float64[8192]
e_max    = Float64[8192]

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

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:

sherpa>save_pha("faked_leg1","simulation1.pha")
sherpa>save_pha("faked_leg-1","simulation-1.pha")

There is also an option to save the data to an ASCII file:

sherpa>save_pha("faked_leg1","simulation1.dat", ascii=True)
sherpa>save_pha("faked_leg-1","simulation-1.dat", ascii=True)

Fitting the Simulated Data

The simulated data set may be filtered and fit as any other data set in Sherpa.


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

sherpa> show_filter()
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 = 18980.6
Final fit statistic   = 957.705 at function evaluation 591
Data points           = 962
Degrees of freedom    = 958
Probability [Q-value] = 0.496615
Reduced statistic     = 0.999692
Change in statistic   = 18022.8
   bpow1.gamma1   1.54257     
   bpow1.gamma2   1.32829     
   bpow1.eb       14.4672 
   bpow1.ampl     0.0223159      

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

The resulting plot is shown in Figure 1.

[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.ampl lower bound:	-0.00165232
bpow1.gamma2 lower bound:	-----
bpow1.gamma1 lower bound:	-0.0391858
bpow1.gamma2 upper bound:	-----
bpow1.gamma1 upper bound:	0.0390244
bpow1.ampl upper bound:	0.00178702
bpow1.eb -: WARNING: The confidence level lies within (4.482115e-05, 0.000000e+00)
bpow1.eb lower bound:	-14.4672
bpow1.eb upper bound:	-----
Datasets              = 'faked_leg1', 'faked_leg-1'
Confidence Method     = confidence
Iterative Fit Method  = None
Fitting Method        = neldermead
Statistic             = cstat
confidence 1-sigma (68.2689%) bounds:
   Param            Best-Fit  Lower Bound  Upper Bound
   -----            --------  -----------  -----------
   bpow1.gamma1      1.54257   -0.0391858    0.0390244
   bpow1.gamma2      1.32829        -----        -----
   bpow1.eb          14.4672     -14.4672        -----
   bpow1.ampl      0.0223159  -0.00165232   0.00178702

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).

    This is equivalent to the XSpec implementation of the
    Cash statistic [1]_ except that it requires a model to be fit
    to the background. To handle the background in the same manner
    as XSpec, use the WStat statistic.

    Counts are sampled from the Poisson distribution, and so the best
    way to assess the quality of model fits is to use the product of
    individual Poisson probabilities computed in each bin i, or the
    likelihood L:

    L = (product)_i [ M(i)^D(i)/D(i)! ] * exp[-M(i)]

    where M(i) = S(i) + B(i) is the sum of source and background model
    amplitudes, and D(i) is the number of observed counts, in bin i.

    The cstat statistic is derived by (1) taking the logarithm of the
    likelihood function, (2) changing its sign, (3) dropping the
    factorial term (which remains constant during fits to the same
    dataset), (4) adding an extra data-dependent term (this is what
    makes it different to `Cash`, and (5) multiplying by two:

    C = 2 * (sum)_i [ M(i) - D(i) + D(i)*[log D(i) - log M(i)] ]

    The factor of two exists so that the change in the cstat statistic
    from one model fit to the next, (Delta)C, is distributed
    approximately as (Delta)chi-square when the number of counts in
    each bin is high.  One can then in principle use (Delta)C instead
    of (Delta)chi-square in certain model comparison tests. However,
    unlike chi-square, the cstat statistic may be used regardless of
    the number of counts in each bin.

    The inclusion of the data term in the expression means that,
    unlike the Cash statistic, one can assign an approximate
    goodness-of-fit measure to a given value of the cstat statistic,
    i.e. the observed statistic, divided by the number of degrees of
    freedom, should be of order 1 for good fits.

    Notes
    -----
    The background should not be subtracted from the data when this
    statistic is used.  It should be modeled simultaneously with the
    source.

    The cstat statistic function evaluates the logarithm of each data
    point. If the number of counts is zero or negative, it's not
    possible to take the log of that number. The behavior in this case
    is controlled by the `truncate` and `trunc_value` settings in the
    .sherpa.rc file:

    - if `truncate` is `True` (the default value), then
      `log(trunc_value)` is used whenever the data value is <= 0.  The
      default is `trunc_value=1.0e-25`.

    - when `truncate` is `False` an error is raised.

    References
    ----------

    .. [1] The description of the Cash statistic (`cstat`) in
           https://heasarc.gsfc.nasa.gov/xanadu/xspec/manual/XSappendixStatistics.html

    

Fit:Datasets              = 'faked_leg1', 'faked_leg-1'
Method                = neldermead
Statistic             = cstat
Initial fit statistic = 18980.6
Final fit statistic   = 957.705 at function evaluation 591
Data points           = 962
Degrees of freedom    = 958
Probability [Q-value] = 0.496615
Reduced statistic     = 0.999692
Change in statistic   = 18022.8
   bpow1.gamma1   1.54257     
   bpow1.gamma2   1.32829     
   bpow1.eb       14.4672     
   bpow1.ampl     0.0223159   

sherpa> print get_conf_results()
datasets    = ('faked_leg1', 'faked_leg-1')
methodname  = confidence
iterfitname = none
fitname     = neldermead
statname    = cstat
sigma       = 1
percent     = 68.2689492137
parnames    = ('bpow1.gamma1', 'bpow1.gamma2', 'bpow1.eb', 'bpow1.ampl')
parvals     = (1.5425700646177956, 1.3282897030064635, 14.467198800903374, 0.022315909860366474)
parmins     = (-0.039185829887720569, None, -14.467176390327081, -0.0016523183493867485)
parmaxes    = (0.03902443229449637, None, None, 0.0017870162192871603)
nfits       = 97

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.

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 exec(open("fit.py").read()) 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, et cetera.)

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.9 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
12 Dec 2013 reviewed for CIAO 4.6: fixed typos, updated calibration files, show restore warning messages.
03 Mar 2014 updated for Chandra proposal cycle 16
30 Jan 2015 updated for Chandra proposal cycle 17
15 Dec 2015 updated for Chandra proposal cycle 18
06 Dec 2016 updated for Chandra proposal cycle 19; include defining source model found in the restored state file, that would fail to load in newer versions of Sherpa.


Last modified: 6 Dec 2016
Smithsonian Institute Smithsonian Institute

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