Simulating 1-D Data: the Sherpa FAKE_PHA Command

Getting started

Please follow the "Obtaining data used in Sherpa threads" thread.

Introduction to FAKE_PHA

The fake_pha command creates a simulated 1-D dataset using a previously defined source model expression and instrument model; Poisson deviates are then added to the modeled data. If a background file is supplied, then the backscale correction is used, and background is added to the simulated data before Poisson noise.

Note that all Sherpa library models can be used in the simulations. In addition, models created and implemented by a user may be included. Individual models can be combined into composite models, and parameters can be linked across the model components.

Note also that both the source and instrument models must be predefined. Currently, only instrument models that are defined using redistribution matrix (RMF) and ancillary (ARF) 1-D spectral files may be used with fake_pha, although an instrument defined solely by a RMF is possible.

There are several steps involved in creating simulations:

1. Defining an instrument model using a redistribution matrix file and an ancillary response file.
2. Defining a source model expression.
3. Defining a background by reading a background data file or specifying background model.
4. Running the simulation using FAKE_PHA.
5. Defining the model normalization for the simulated data, if desired.
6. Writing the simulated data to output files.

This thread illustrates each of these steps with an example. This thread will also illustrate use of the fake_pha command when a PHA data file has been previously read. In addition, the thread will show how a dataset placed into memory by fake_pha can be plotted and fit. The thread concludes with an additional overall example.

NOTE, if you already have a PHA dataset: Sherpa automatically loads instrument models and background data when PHA data is called on using load_pha or load_data. If this is the case, then you may proceed to Using FAKE_PHA with a PHA File. If there is no PHA data available, the background data and source and instrument models are manually defined and loaded with the following commands:

Defining a Source Model Expression without a Dataset

The source model expression is defined in the standard way. NOTE: The next step is required if the user wishes to run fake_pha with an instrument + model, only. Before defining the model, it is necessary to first produce an empty PHA dataset to associate the model with dataspace1d:

sherpa> dataspace1d(0.,512.,id=1,dstype=DataPHA)
sherpa> get_data(1).channel=get_data(1).channel.astype(float)
sherpa> set_analysis("energy",type="rate")
sherpa> show_data()
Data Set: 1
Filter: 0.5000-511.5000 Energy (keV)
Noticed Channels: 1.0-512.0
name           = dataspace1d
channel        = Float64[512]
counts         = Float64[512]
staterror      = None
syserror       = None
bin_lo         = Float64[512]
bin_hi         = Float64[512]
grouping       = None
quality        = None
exposure       = None
backscal       = None
areascal       = None
grouped        = False
subtracted     = False
units          = energy
response_ids   = []
background_ids = []

This creates an empty, 512 element, 64-bit floating point PHA array using dataspace1d and changing the default units of "channels" to "energy." NOTE: The use of dataspace1d to load a blank PHA data set is strictly a temporary work-around that will be addressed in the CIAO 4.1.1 update.

The source model is defined in the standard manner, where we assume an 1-D power-law for the source.

sherpa> set_source(1,"powlaw1d.modelA")
sherpa> modelA.gamma=2
sherpa> show_model(1)
Model: 1
powlaw1d.modelA
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   modelA.gamma thawed            2          -10           10
   modelA.ref   frozen            1 -3.40282e+38  3.40282e+38
   modelA.ampl  thawed            1            0  3.40282e+38

The Sherpa model library and model syntax are described in detail in the Sherpa Reference Manual, and the model parameters can then be set as desired:

Defining an Instrument Model

The default instrument response files for 1-D simulations of Chandra data include the redistribution matrix (RMF) and the ancillary response (ARF) files in standard FITS format. These files can be created with the mkrmf and mkarf tools, or with a script as described in the Using specextract to Extract ACIS Spectra and Response Files thread.

sherpa> arf1=unpack_arf("dataA_arf.fits")
sherpa> print arf1 
name     = dataA_arf.fits
energ_lo = Float64[1180]
energ_hi = Float64[1180]
specresp = Float64[1180]
bin_lo   = None
bin_hi   = None
exposure = None

sherpa> rmf1=unpack_rmf("dataA_rmf.fits")
sherpa> print rmf1
name     = dataA_rmf.fits
detchans = 512
energ_lo = Float64[1180]
energ_hi = Float64[1180]
n_grp    = Int16[1180]
f_chan   = UInt32[1189]
n_chan   = UInt32[1189]
matrix   = Float64[260775]
offset   = 0
e_min    = Float64[512]
e_max    = Float64[512]

sherpa> set_arf(arf1)
sherpa> set_rmf(rmf1)

Define a Background

A background model or background data can be used in the simulations. Note that if PHA data are input as a background file, TIME and BACKSCALE parameters will have default settings corresponding to the values of the header keywords EXPTIME and BACKSCAL. These may of course be changed.

Background data and/or models are treated as follows in FAKE_PHA:

1. If a background model stack is defined, it is evaluated on the source data grid, and the resulting background amplitudes are added to the source amplitudes (taking into account differences in exposure time and backscale). Faked data are then sampled given the sum. If background data exist, they are not altered. If the source dataset was background-subtracted prior to the command FAKE_PHA being issued, it will not be background-subtracted afterwards.

2. If no background model stack is defined, and the data are background-subtracted, then the source model stack is evaluated directly, and the new, faked data are background-subtracted. Note that subsequently issuing an unsubtract command is unwise, because the unsubtracted data will not be integer counts data.

3. If no background model stack is defined, and the data are not background-subtracted, then the source model stack is evaluated directly, and the (properly scaled) background data are added to the faked data.

sherpa> bkg1=unpack_pha(1,"dataA_bkg.pha")[0]
sherpa> set_bkg(1,bkg1)
sherpa> print bkg1
name           = dataA_bkg.pha
channel        = Float64[512]
counts         = Float64[512]
staterror      = None
syserror       = None
bin_lo         = None
bin_hi         = None
grouping       = None
quality        = Int16[512]
exposure       = 108675.66732
backscal       = 0.044189453125
areascal       = 1.0
grouped        = False
subtracted     = False
units          = channel
response_ids   = []
background_ids = []

sherpa> set_bkg_model(1,polynom1d.bkgA+gauss1d.bkgB)
sherpa> bkgA.c0=6.4e-5
sherpa> bkgA.c1=2.3e-5
sherpa> bkgA.c3=1.4e-05
sherpa> bkgB.fwhm=5.6254
sherpa> bkgB.pos=0.057
sherpa> bkgB.ampl=0.003
sherpa> show_model(1)
Model: 1
(powlaw1d.modelA + scale_factor * ((polynom1d.bkgA + gauss1d.bkgB)))
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   modelA.gamma thawed            2          -10           10
   modelA.ref   frozen            1 -3.40282e+38  3.40282e+38
   modelA.ampl  thawed            1            0  3.40282e+38
   bkgA.c0      thawed      6.4e-05 -3.40282e+38  3.40282e+38
   bkgA.c1      frozen      2.3e-05 -3.40282e+38  3.40282e+38
   bkgA.c2      frozen            0 -3.40282e+38  3.40282e+38
   bkgA.c3      frozen      1.4e-05 -3.40282e+38  3.40282e+38
   bkgA.c4      frozen            0 -3.40282e+38  3.40282e+38
   bkgA.c5      frozen            0 -3.40282e+38  3.40282e+38
   bkgA.c6      frozen            0 -3.40282e+38  3.40282e+38
   bkgA.c7      frozen            0 -3.40282e+38  3.40282e+38
   bkgA.c8      frozen            0 -3.40282e+38  3.40282e+38
   bkgA.offset  frozen            0 -3.40282e+38  3.40282e+38
   bkgB.fwhm    thawed       5.6254  1.17549e-38  3.40282e+38
   bkgB.pos     thawed        0.057 -3.40282e+38  3.40282e+38
   bkgB.ampl    thawed        0.003 -3.40282e+38  3.40282e+38

Setting Parameters for the Synthesized Dataset

Since the current PHA dataset is completely empty, we must setup the desired properties for the synthesized dataset. In this case, we would like to simulate a dataset from a 33 ks exposure, and we set the backscale and areascale to be the same as that of the background dataset.

sherpa> get_data(1).exposure=33483.2
sherpa> get_data(1).backscal=get_bkg(1).backscal
sherpa> get_data(1).areascal=get_bkg(1).areascal

Since the simulated data and the background data will have the same backscale value in this example, this means that the area of the simulated source will be the same as the supplied background file. If the backscale is set to 1 (the default value), then the scaling of the area is assumed based on backscale from the header of the background file.

Running the Simulation Using FAKE_PHA

A simulation is run with fake_pha, where a set of four parameters--dataset ID, ARF, RMF, and exposure time--are required, and several options may be called on where the convolved source model with scaled background with Poisson noise is calculated:

sherpa> fake_pha(1,arf=arf1,rmf=rmf1,exposure=get_exposure(1),grouped=False,bkg=get_bkg(1))

This creates a synthetic dataset at dataset ID 1 using an instrument model based on the ARF and RMF unpacked to variables arf1 and rmf1, respectively, for an exposure time 33483.2 seconds, as pre-defined in the empty PHA dataset, called by get_exposure. The options grouped=False implies that the dataset is not grouped, and bkg calls on the background dataset associated with this synthesized data.

We may inspect some basic properties of the new dataset:

sherpa> sum(get_counts(1))
           4939599.0
sherpa> calc_data_sum(id=1,lo=None,hi=None)
           4939599.0

These two commands are equivalent, giving the total number of counts simulated in the dataset. get_counts provides an array of the unfiltered counts per channel, and the sum function adds together each element of the array. The calc_data_sum function sums up the counts per array element for dataset 1, but includes an option for filtering, which in this case is energy in units of keV. By default, CALC_DATA_SUM introduces an energy filter that excludes everything <1.0 keV (i.e. lo=1.0,hi=None) , but to include all counts without filtering, we choose (lo=None,hi=None).

Similarly, we may determine the count rate of this dataset in units of counts/second by:

sherpa> sum(get_counts(1))/get_exposure(1)
           147.52469895350507

The simulated dataset is not saved to disk until the command write_pha is given. (see the Writing The Simulated Data To Output Files section).

If we are absent an ARF for the instrument model, we may still simulate a dataset using just the RMF by setting arf=None as follows:

sherpa> fake_pha(1,arf=None,rmf=rmf1,exposure=get_exposure(1),grouped=False,bkg=get_bkg(1))
sherpa> sum(get_counts(1))
            28301.0

Defining the Model Normalization for the Simulated Data

Simulated data should have a normalization which agrees, for example, with the observed source flux. If the flux of the simulated model is known, it can be easily converted to the model parameter values. This requires simulating the data with the default parameters, calculating the flux, and then rescaling the model parameters. We first describe the steps and then use the Python functions to demonstrate how efficient it is to renormalize the model through scripting.

Assuming that we have just simulated data using default parameters, including instrument model using both ARF and RMF, and modelA.ampl=1, we calculate the flux:

sherpa> fake_pha(1,arf=arf1,rmf=rmf1,exposure=get_exposure(1),grouped=False,bkg=get_bkg(1))
sherpa> calc_energy_flux(lo=2.,hi=10.)
            2.5706384381817505e-09

This provides the unconvolved energy flux of the model, between 2.0-10.0 keV in units of ergs cm^-2 s^-1.

To find the normalization, given a source flux of 1.0e-13 ergs cm^-2 s^-1, we divide the true source flux (1.0e-13) by the calculated flux for the above default source model (2.5706384381817505e-09); the normalization (modelA.ampl) was 1 photon keV^-1 cm^-2 s^-1:

sherpa>1.0e-13/calc_energy_flux(lo=2.,hi=10.)
            3.8900842107819502e-05

So, the normalization of the source model may be adjusted as follows:

sherpa> modelA.ampl=1.e-13/calc_energy_flux(lo=2.0,hi=10.0)

and fake_pha is run again:

sherpa> fake_pha(1,arf=arf1,rmf=rmf1,exposure=get_exposure(1),grouped=False,bkg=get_bkg(1))sherpa> calc_photon_flux(lo=2.,hi=10.)
            1.5463634382489872e-05
sherpa> calc_energy_flux(lo=2.,hi=10.)
            1.0000000000000012e-13

Which also provides the photon flux between 2.0-10.0 keV in units of photons cm^-2 s^-1. If the count rate of the source in a given instrument is known in advance (e.g. it can be obtained from PIMMS, or from previous observations), then a simple scaling can be used to adjust the normalization of the source model.

For example, say the known count rate is 0.4 cts/sec. From the first fake_pha run (or a repeat run, with modelA.ampl returned to 1.0), we note the non-normalized count rate:

sherpa> modelA.ampl=1.0
sherpa> fake_pha(1,arf=arf1,rmf=rmf1,exposure=get_exposure(1),grouped=False,bkg=get_bkg(1))
sherpa> rate=sum(get_counts(1))/get_exposure(1)
sherpa> print rate
147.392572992
sherpa> modelA.ampl=0.4/rate
sherpa> fake_pha(1,arf=arf1,rmf=rmf1,exposure=get_exposure(1),grouped=False,bkg=get_bkg(1))
sherpa> sum(get_counts(1))
            33516.0
sherpa> sum(get_counts(1))/get_exposure(1)
            1.000979595737564

Using this non-normalized count rate (147.4 cts/sec), the dataset could be normalized by scaling the model amplitude by setting the value to the known count rate and dividing it by the non-normalized count rate.

NOTE for FAKEIT users: The sum of the simulated counts is less than Sherpa 3.4. This is because the scaled background model counts are now correctly handled in Sherpa 4.1.

Writing the Simulated Data to Output Files

The simulated data are not automatically saved to disk. The user may write the data as a PHA file where the PHA header will contain exposure time and backscale values, as well as paths to the RMF, ARF, and background files.

This is done with the save_pha and command:

sherpa> save_pha(1,"dataA_sim1.pha")

This creates a PHA file with columns of channel and counts. This file may then be grouped using dmgroup if necessary.

To save the simulated data to an ASCII table:

sherpa> set_analysis("energy",type="rate")
sherpa> load_arrays("flux",get_model_plot(1).xlo,get_model_plot(1).y)
sherpa> save_data("flux","dataA_sim1.dat")

Using FAKE_PHA with a PHA File

The fake_pha command will operate whether or not the user has previously read data from a PHA file. If the user has previously read data from a PHA file:

• An instrument model may have been automatically loaded if the response files are properly referenced in the header of the data file.
• If the data file's backfile keyword supplies a background file, then it will be used in the simulations.
• Since input PHA files also usually contain the exposure time and backscale keywords pertinent to an observation, this information can be easily called on, especially useful for scripting.
• The input dataset will be overwritten by the simulated dataset created by fake_pha.

For example, erase all current Sherpa settings (or EXIT) and begin a new session), and then input a PHA data file:

sherpa> clean()
sherpa> show_all()

sherpa> load_pha("dataA.pha")
systematic errors were found in file 'dataA.pha'
but not used; to use them, re-read with use_errors=True
read ARF file dataA_arf.fits
read RMF file dataA_rmf.fits
read background file dataA_bkg.pha
sherpa> get_exposure(1)
            33483.25
sherpa> get_backscal(1)
            0.044189453125

Loading this PHA data file resulted in an instrument model being automatically defined (using RMF and ARF files), and a background data file being automatically input. And, from the get_exposure and get_backscal commands we see that the exposure time and backscale are readily called upon, automatically.

That is, input of this PHA data file has automatically completed two steps involved in creating simulations:

• Defining an instrument model using a redistribution matrix file and an ancillary response file.
• Reading a background data file.

The following commands will complete the remaining steps:

• Defining a source model expression.
• Running the simulation using FAKE_PHA.
• Defining the model normalization for the simulated data, if desired.
• Writing the simulated data to output files.

Important note: If your input PHA data are grouped, the simulated data produced by FAKE_PHA will also be grouped. (The grouping cannot be removed with ungroup.) When you write the simulated data to a file, the channels in the output file represent the grouped channels. However, no grouping information is written to the file, so when you read it back into Sherpa, the bin energies are incorrect. Currently, there is no workaround for this problem, so WRITE_PHA should be used only with ungrouped PHA data.

sherpa> set_source(powlaw1d.modelA)
sherpa> modelA.gamma=2
sherpa> fake_pha(1,arf=get_arf(),rmf=get_rmf(),exposure=get_exposure(),bkg=get_bkg(),grouped=False)
sherpa> sum(get_counts(1))
            4919567.0
sherpa> modelA.ampl=0.4/(sum(get_counts(1))/get_exposure(1))
sherpa> fake_pha(1,arf=get_arf(1),rmf=get_rmf(1),exposure=get_exposure(1),grouped=False,bkg=get_bkg(1))
sherpa> sum(get_counts(1))
           13465.0

Note that the fake_pha command overwrites the input dataset with the simulated data, and we have futher normalized the simulated dataset.

The results may be written out to a PHA file:

sherpa> save_pha(1,"dataA_sim2.pha")

Plotting and Fitting Simulated Data

The simulated dataset may be plotted in the same way as any other dataset:

sherpa> plot_data(1)

Figure 1  shows the resulting plot.

The simulated dataset may be also filtered and fit in the same way as any other dataset. For example, we next filter the simulated dataset, to include only data between 0.5 and 8.0 keV:

sherpa> notice(0.5,8)
sherpa> plot_data(1)

Figure 2  shows the resulting plot. Then, we fit the simulated dataset using the same source model expression that had been used to create it:

sherpa> subtract(1)
sherpa> set_method("neldermead") 
sherpa> fit(1)
Dataset               = 1
Method                = neldermead
Statistic             = chi2gehrels
Initial fit statistic = 192.902
Final fit statistic   = 162.514 at function evaluation 286
Data points           = 263
Degrees of freedom    = 261
Probability [Q-value] = 1
Reduced statistic     = 0.62266
Change in statistic   = 30.3872
   modelA.gamma   2.04951
   modelA.ampl    0.00263902
sherpa> plot_fit(1)

Figure 3  shows the resulting plot.

We here examine the fit using the Sherpa command covariance:

sherpa> covariance()
Dataset               = 1
Confidence Method     = covariance
Fitting Method        = neldermead
Statistic             = chi2gehrels
covariance 1-sigma (68.2689%) bounds:
   Param            Best-Fit  Lower Bound  Upper Bound
   -----            --------  -----------  -----------
   modelA.gamma      2.04951   -0.0193856    0.0193856
   modelA.ampl    0.00263902 -3.07702e-05  3.07702e-05

This quick method to calculate one-sigma parameter ranges, using COVAR, may underestimate the errors for a complex parameter space. The slower, but more appropriate, projection algorithm should be used in such cases.

Simulating Multiple Datasets

One may simulate multiple datasets within the same Sherpa session. Each dataset is assigned a number.

For example, we next simulate a second dataset, as we did above by defining the instrument with the ARF and RMF, and creating an empty PHA dataset, except matching the exposure, backscale, and area scale to the second background dataset:

sherpa> arf2=unpack_arf("dataB_arf.fits")
sherpa> rmf2=unpack_rmf("dataB_rmf.fits")
sherpa> dataspace1d(0.0,1024.0,id=2,dstype=DataPHA)
sherpa> get_data(id=2).channel = get_data(id=2).channel.astype(float)
sherpa> get_data(id=2).units="energy"

sherpa> set_source(2,xsphabs.modelB1*powlaw1d.modelB2)
sherpa> modelB1.nH=0.025
sherpa> modelB2.gamma=1.7
sherpa> modelB2.ampl=1.0

sherpa> load_bkg(2,"dataB_bkg.pha")
sherpa> get_data(2).exposure=get_bkg(2).exposure
sherpa> get_data(2).backscal=get_bkg(2).backscal
sherpa> get_data(2).areascal=get_bkg(2).areascal
sherpa> fake_pha(2,arf=arf2,rmf=rmf2,exposure=get_exposure(2),bkg=get_bkg(2),grouped=False)
 Solar Abundance Vector set to angr:  Anders E. & Grevesse N. Geochimica et Cosmochimica Acta 53, 197 (1989)
 Cross Section Table set to bcmc:  Balucinska-Church and McCammon, 1998
sherpa> sum(get_counts(2))
           30933676.0

Note that in order to simulate dataset number 2, there must be a ARF and RMF for a new instrument model, a background dataset number 2.

We will now normalize this second dataset, in the same manner as the first simulated dataset (dividing the known count rate of say 0.4, with the non-normalized count rate, to obtain a scaling factor; and re-simulate the second dataset:

sherpa> modelB2.ampl=0.4/(sum(get_counts(2))/get_exposure(2))
sherpa> fake_pha(2,arf=arf2,rmf=rmf2,exposure=get_exposure(2),bkg=get_bkg(2),grouped=False)
sherpa> sum(get_counts(2))
            10852.0

And finally, the second simulated dataset is plotted:

sherpa> plot_data(2)

The second simulated dataset may also be written to disk:

sherpa> save_pha(2,"dataB_sim1.pha")

To plot both simulated datasets at once:

sherpa> notice(0.5,8)      
sherpa> plot("data",1,"data",2)

Figure 4  shows the resulting plot.

To exit Sherpa:

sherpa> EXIT
Goodbye.
unix%

Using a Sherpa Script to Run Simulations

sherpa> !cat multi_sim.py
#!/usr/bin/env python
from sherpa.astro.ui import *

# Create a fake spectrum with a count rate of <cnt_rate> in the energy 
# range <energy>, with a total exposure time of <fake_time>.

fake_time = 100000
elo = 0.5
ehi = 2.0
cnt_rate = 0.02

# Define instrument response and source model
rmf = unpack_rmf("aciss_aimpt_cy10.rmf")
arf = unpack_arf("aciss_aimpt_cy10.arf")
set_source(powlaw1d.pow1)
pow1.gamma = 1.9

# Define structure to allow manupulation of the pow1.ampl parameter
# value in Python, and give an initial guess for fake_pha
pow_ampl = get_par("pow1.ampl")
pow_ampl.val = 1e-5

# We now create an empty PHA dataset once to define an initial dataset.  We can
# then set the exposure time set_exposure() command to function
dataspace1d(0,1024,dstype=DataPHA)
get_data().channel = get_data().channel.astype(float)
get_data().units="energy"
set_exposure(fake_time)

# Run fake_pha once to get model counts for initial amplitude, then
# normalize to get desired count rate
set_par(pow_ampl)
fake_pha(1,arf=arf,rmf=rmf,exposure=fake_time,grouped=False)
pow_ampl.val *= cnt_rate / (calc_model_sum(elo,ehi)/fake_time)
set_par(pow_ampl)
fake_pha(1,arf=arf,rmf=rmf,exposure=fake_time,grouped=False)

# Verify correct counts in fake dataset and plot data
print("Count rate in specified energy range is")
print(calc_data_sum(lo=elo,hi=ehi)/fake_time)
plot_data()

sherpa> execfile("multi_sim.py") 
Count rate in specified energy range is
0.01989

unix% sherpa multi_sim.py

unix% sherpa multi_sim.py >&! multi_sim.out

sherpa> !cat plasma.py
#!/usr/bin/env python
from sherpa.astro.ui import *

# this script simulates emission from a hot plasma with two different models.

arf=unpack_arf("dataB_arf.fits")
rmf=unpack_rmf("dataB_rmf.fits")

fake_time=5000
set_source(1,xswabs.abs1*xsmekal.m1)

abs1.nH=0.03
m1.kT=4.0
m1.norm=0.025

dataspace1d(0,1024,id=1,dstype=DataPHA)
get_data(1).channel = get_data(1).channel.astype(float)
get_data(1).units="energy"

fake_pha(1,arf=arf,rmf=rmf,exposure=fake_time,grouped=False)

# renormalize to 2-10 keV flux of 1e-13 erg/s/cm**2
modflux1=calc_energy_flux(id=1,lo=2,hi=10)
obsflux1=1.e-13

foo1=get_par("m1.norm")
foo1.val=foo1.val*(obsflux1/modflux1)
set_par(foo1)

fake_pha(1,arf=arf,rmf=rmf,exposure=fake_time,grouped=False)
save_pha(1,"sim_meka.pha")


set_source(2,xswabs.abs1*xsapec.ap1)

ap1.kT=4.0
ap1.norm=0.025

dataspace1d(0,1024,id=2,dstype=DataPHA)
get_data(2).channel = get_data(2).channel.astype(float)
get_data(2).units="energy"

fake_pha(2,arf=arf,rmf=rmf,exposure=fake_time,grouped=False)

# renormalize to the 2-10 keV flux of 1.e-13 erg/s/cm**2
modflux2=calc_energy_flux(id=2,lo=2,hi=10)
obsflux2=1.e-13

foo2=get_par("ap1.norm")
foo2.val=foo2.val*(obsflux2/modflux2)
set_par(foo2)

fake_pha(2,arf=arf,rmf=rmf,exposure=fake_time,grouped=False)
save_pha(2,"sim_apec.pha")

exit()

unix% sherpa plasma.py >&! plasma.out

unix% sherpa plasma.sl >&! plasma.out

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)

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.1 syntax to you.