Simulating Chandra ACIS-S LETG Spectra with Sherpa
Sherpa Threads (CIAO 4.9 Sherpa v1)
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.
- Getting started - downloading calibration response files for simulations
- Defining the Instrument Response for the Simulation
- Establishing the Source Model Expression for the Simulation
- Running the Simulation with fake_pha
- Writing the Simulated Data to Output Files
- Fitting the Simulated Data
- Scripting It
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:
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.
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 energ_hi = Float64 specresp = Float64 bin_lo = Float64 bin_hi = Float64 exposure = 28398.6296948 sherpa> print(arf2) name = aciss_leg-1_cy19.garf energ_lo = Float64 energ_hi = Float64 specresp = Float64 bin_lo = Float64 bin_hi = Float64 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 energ_hi = Float64 n_grp = UInt64 f_chan = UInt32 n_chan = UInt32 matrix = Float64 offset = 1 e_min = Float64 e_max = Float64 sherpa> print(rmf2) name = aciss_leg-1_cy19.grmf detchans = 8192 energ_lo = Float64 energ_hi = Float64 n_grp = UInt64 f_chan = UInt32 n_chan = UInt32 matrix = Float64 offset = 1 e_min = Float64 e_max = Float64
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."
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
- 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:
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))
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 counts = Float64 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 =  background_ids =  RMF Data Set: faked_leg1:1 name = aciss_leg1_cy18.grmf detchans = 8192 energ_lo = Float64 energ_hi = Float64 n_grp = UInt64 f_chan = UInt32 n_chan = UInt32 matrix = Float64 offset = 1 e_min = Float64 e_max = Float64 ARF Data Set: faked_leg1:1 name = aciss_leg1_cy18.garf energ_lo = Float64 energ_hi = Float64 specresp = Float64 bin_lo = Float64 bin_hi = Float64 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."
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:
There is also an option to save the data to an ASCII file:
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.
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 _ 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 ---------- ..  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 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")
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.
|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.|