Simulating 1-D Data: the Sherpa FAKE_PHA Command

Sherpa Threads (CIAO 4.1)

Overview

Last Update: 8 Dec 2008 - Created for CIAO/Sherpa4.1, replacing CIAO/Sherpa3.x FAKEIT command

Synopsis:

Sherpa can be used to simulate 1-D data, using the FAKE_PHA (S-Lang or Python help) (formally FAKEIT) command. This thread provides examples of its usage.

Proceed to the HTML or hardcopy (PDF: A4 | letter) version of the thread.

The FAKE_PHA (S-Lang or Python help) command creates a simulated 1-D dataset calculated using a previously defined source model expression and a previously defined 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 is added.

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 (S-Lang or Python help), although an instrument defined solely by a RMF is possible.

There are several steps involved in creating simulations:

Defining an instrument model using a redistribution matrix file and an ancillary response file.
Defining a source model expression.
Defining a background by reading a background data file or specifying background model.
Running the simulation using FAKE_PHA.
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 (S-Lang or Python help) 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 (S-Lang or Python help) can be plotted and fit. The thread concludes with an additional overall example.

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 automatically loads instrument models when PHA data is called on using LOAD_PHA (S-Lang or Python help) or LOAD_DATA (S-Lang or Python help). 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 RMF and ARF may be manually loaded with the following commands:

sherpa> arf1=UNPACK_ARF (S-Lang or Python help)("dataA_arf.fits")
sherpa> rmf1=UNPACK_RMF (S-Lang or Python help)("dataA_rmf.fits")

sherpa> arf1=UNPACK_ARF (S-Lang or Python help)("dataA_arf.fits");
sherpa> rmf1=UNPACK_RMF (S-Lang or Python help)("dataA_rmf.fits");

Defining a Source Model Expression

The source model expression is defined in the standard way; however, before defining the model, it is necessary to first produce an empty PHA data set to associate the model with:

sherpa> DATASPACE1D (S-Lang or Python help)(0.,512.,id=1,dstype=DataPHA)
sherpa> GET_DATA (S-Lang or Python help)(1).channel=get_data(1).channel.astype(float)
sherpa> get_data(1).units="energy"
sherpa> SHOW_DATA (S-Lang or Python help)()
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 = []

sherpa> DATASPACE1D (S-Lang or Python help)(0.,512.;id=1,dstype=DataPHA);
sherpa> GET_DATA (S-Lang or Python help)(1).channel=get_data(1).channel.astype(Float_Type);
sherpa> get_data(1).units="energy";
sherpa> SHOW_DATA (S-Lang or Python help)();
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 created an empty, 512 element, 64-bit floating point PHA array using DATASPACE1D (S-Lang or Python help) and changing the default units of "channels" to "energy."

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

sherpa> SET_SOURCE (S-Lang or Python help)(1,"POWLAW1D (S-Lang or Python help).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

sherpa> SET_SOURCE (S-Lang or Python help)(1,"POWLAW1D (S-Lang or Python help).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:

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:

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.
If no background model stack is defined, and the data are background-subtacted, then the source model stack is evaluated directly, and the new, faked data are background-subtracted. Note that subsequently issuing an UNSUBTRACT (S-Lang or Python help) command is unwise, because the unsubtracted data will not be integer counts data.
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.

Reading a background data file

A PHA-type background file is read so that the exposure time and backscale information may be taken into account in the simulations. The exposure time and backscale from the background file is only used for apropriate scaling of the input background data.

sherpa> LOAD_BKG (S-Lang or Python help)(1,"dataA_bkg.pha")

sherpa> LOAD_BKG (S-Lang or Python help)(1,"dataA_bkg.pha");

This loads the background PHA dataset, associating it with the dataset. Unlike FAKEIT, an instrument model does not need to be explicitly set for the background, so long as a RMF and ARF are associated with the dataset ID.

Define a Background Model Expression

If a background model expression is defined then the background model amplitude will be added to the source amplitudes (taking into account differences in exposure time and backscale). Faked data are then sampled given the sum. A background model is defined in a standard way with the model name given in the bracket. In this example, the background source model is a combination of a 1-D polynomial and Gaussian models:

sherpa> SET_BKG_MODEL (S-Lang or Python help)(1,polynom1d.bkgA+gauss1dbkgB)
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

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-5;
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

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 (S-Lang or Python help), where a set of four parameters--dataset ID, ARF, RMF, and exposure time--are required, and several options may be called on:

sherpa> FAKE_PHA (S-Lang or Python help)(1,arf=arf1,rmf=rmf1,exposure=get_exposure(1),grouped=False,bkg=get_bkg(1))

sherpa> FAKE_PHA (S-Lang or Python help)(1;arf=arf1,rmf=rmf1,exposure=get_exposure(1),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 (S-Lang or Python help). The options grouped=False implies that the dataset is not grouped, and bkg calls on the background dataset associated with this synthesized data.

Note that when you're working in S-lang and you use keyword arguments, the delimiter before the first one must be a semi-colon and not a comma.

We may inspect some basic properties of the new dataset:

sherpa> sum(GET_COUNTS (S-Lang or Python help)(1))
           4939599.0
sherpa> CALC_DATA_SUM (S-Lang or Python help)(id=1,lo=None,hi=None)
           4939599.0

sherpa> sum(GET_COUNTS (S-Lang or Python help)(1));
4.93777e+06
sherpa> CALC_DATA_SUM (S-Lang or Python help)(NULL,NULL,1);
4.93777e+06

These two commands are equivalent, giving the total number of counts simulated in the dataset. GET_COUNTS (S-Lang or Python help) provides an array of the unfiltered counts per channel, and the sum function adds together each element of the array. The CALC_DATA_SUM (S-Lang or Python help) 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) (i.e. lo=1.0,hi=NULL), but to include all counts without filtering, we choose (lo=None,hi=None)(lo=NULL,hi=NULL).

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

sherpa> sum(get_counts(1))/get_exposure(1);
147.47

The simulated dataset is not saved to disk until the command WRITE_PHA (S-Lang or Python help) 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

sherpa> fake_pha(1;arf=NULL,rmf=rmf1,exposure=get_exposure(1),grouped=NULL,bkg=get_bkg(1));
sherpa> sum(get_counts(1));
28214

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 PythonS-lang 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 (S-Lang or Python help)(lo=2.,hi=10.)
            2.5706384381817505e-09

sherpa> fake_pha(1,arf=arf1,rmf=rmf1,exposure=get_exposure(1),grouped=NULL,bkg=get_bkg(1))
sherpa> CALC_ENERGY_FLUX (S-Lang or Python help)(lo=2.,hi=10.);
2.57064e-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

sherpa>1.0e-13/calc_energy_flux(;lo=2.,hi=10.);
3.89008e-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)

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

and FAKE_PHA (S-Lang or Python help) is run again:

sherpa> FAKE_PHA (S-Lang or Python help)(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

sherpa> fake_pha(1;arf=arf1,rmf=rmf1,exposure=get_exposure(1),grouped=NULL,bkg=get_bkg(1));
sherpa> calc_photon_flux(lo=2.,hi=10.);
1.54636e-05
sherpa> calc_energy_flux(lo=2.,hi=10.);
1e-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 (S-Lang or Python help) 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 (S-Lang or Python help)(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 (S-Lang or Python help)(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

sherpa> modelA.ampl=1.0;
sherpa> FAKE_PHA (S-Lang or Python help)(1;arf=arf1,rmf=rmf1,exposure=get_exposure(1),grouped=NULL,bkg=get_bkg(1));
sherpa> rate=sum(get_counts(1))/get_exposure(1);
sherpa> print(rate)
147.396
sherpa> modelA.ampl=0.4/rate;
sherpa> FAKE_PHA (S-Lang or Python help)(1;arf=arf1,rmf=rmf1,exposure=get_exposure(1),grouped=NULL,bkg=get_bkg(1));
sherpa> sum(get_counts(1));
33525
sherpa> sum(get_counts(1))/get_exposure(1);
1.00125

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.

Defining Model Normalization using Sherpa Module Functions

All the steps required to set model normalization can be efficiently scripted using Sherpa Module functions calc_energy_flux (S-Lang or Python help), get_par (S-Lang or Python help), and set_par (S-Lang or Python help):

sherpa> modflux=calc_energy_flux(id=1,lo=2,hi=10)
sherpa> print modflux
6.9728917106e-12

sherpa> obsflux=1e-13
sherpa> foo=get_par("modelA.ampl")
sherpa> print foo
val         = 0.00271251
min         = 0
max         = 3.4e+38
units       =
frozen      = False
link        = None
default_val = 0.00271251
default_min = 0
default_max = 3.4e+38

sherpa> foo.val=foo.val*(obsflux/modflux)
sherpa> print foo
val         = 3.89008e-05
min         = 0
max         = 3.4e+38
units       =
frozen      = False
link        = None
default_val = 3.89008e-05
default_min = 0
default_max = 3.4e+38

sherpa> set_par(foo)
sherpa> show_model()
Model: 1
apply_rmf(apply_arf((33483.2 * (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  3.89008e-05            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

sherpa> modflux=calc_energy_flux(;lo=2,hi=10);
sherpa> print(modflux);
6.97613e-12

sherpa> obsflux=1e-13;
sherpa> foo=get_par("modelA.ampl");
sherpa> print(foo);
val         = 0.00271377
min         = 0
max         = 3.4e+38
units       =
frozen      = False
link        = None
default_val = 0.00271377
default_min = 0
default_max = 3.4e+38

sherpa> foo.val=foo.val*(obsflux/modflux);
sherpa> print(foo);
val         = 3.89008e-05
min         = 0
max         = 3.4e+38
units       =
frozen      = False
link        = None
default_val = 3.89008e-05
default_min = 0
default_max = 3.4e+38

sherpa> set_par(foo);
sherpa> show_model();
Model: 1
apply_rmf(apply_arf((33483.2 * (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  3.89008e-05            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

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 (S-Lang or Python help) and command:

sherpa> SAVE_PHA (S-Lang or Python help)(1,"dataA_sim1.pha")

sherpa> SAVE_PHA (S-Lang or Python help)(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.

Using FAKE_PHA with a PHA File

The FAKE_PHA (S-Lang or Python help) 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 (S-Lang or Python help).

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

sherpa> CLEAN (S-Lang or Python help)()
sherpa> SHOW_ALL (S-Lang or Python help)()

sherpa> LOAD_PHA (S-Lang or Python help)("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

sherpa> CLEAN (S-Lang or Python help)();
sherpa> SHOW_ALL (S-Lang or Python help)();

sherpa> LOAD_PHA (S-Lang or Python help)("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.2
sherpa> get_backscal(1);
0.0441895

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 (S-Lang or Python help) and GET_BACKSCAL (S-Lang or Python help) 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:

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 (S-Lang or Python help).) 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 (S-Lang or Python help)(powlaw1d.modelA)
sherpa> modelA.gamma=2
sherpa> FAKE_PHA (S-Lang or Python help)(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

sherpa> SET_SOURCE (S-Lang or Python help)(powlaw1d.modelA);
sherpa> modelA.gamma=2;
sherpa> FAKE_PHA (S-Lang or Python help)(1;arf=get_arf(),rmf=get_rmf(),exposure=get_exposure(),bkg=get_bkg(),grouped=NULL);
sherpa> sum(get_counts(1));
4.91804e+06
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=NULL,bkg=get_bkg(1));
sherpa> sum(get_counts(1));
13388

Note that the FAKE_PHA (S-Lang or Python help) 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 (S-Lang or Python help)(1,"dataA_sim2.pha")

sherpa> SAVE_PHA (S-Lang or Python help)(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 (S-Lang or Python help)(1)

sherpa> PLOT_DATA (S-Lang or Python help)(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 (S-Lang or Python help)(0.5,8)
sherpa> PLOT_DATA (S-Lang or Python help)(1)

sherpa> NOTICE (S-Lang or Python help)(0.5,8);
sherpa> PLOT_DATA (S-Lang or Python help)(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 (S-Lang or Python help)(1)
sherpa> SET_METHOD (S-Lang or Python help)("NELDERMEAD (S-Lang or Python help)") 
sherpa> FIT (S-Lang or Python help)(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 (S-Lang or Python help)(1)

sherpa> SUBTRACT (S-Lang or Python help)(1);
sherpa> SET_METHOD (S-Lang or Python help)("NELDERMEAD (S-Lang or Python help)"); 
sherpa> FIT (S-Lang or Python help)(1);
Dataset               = 1
Method                = neldermead
Statistic             = chi2gehrels
Initial fit statistic = 167.62
Final fit statistic   = 153.692 at function evaluation 280
Data points           = 263
Degrees of freedom    = 261
Probability [Q-value] = 1
Reduced statistic     = 0.58886
Change in statistic   = 13.9281
   modelA.gamma   2.06324
   modelA.ampl    0.00265019

sherpa> PLOT_FIT (S-Lang or Python help)(1);

Figure 3 shows the resulting plot.

We here examine the fit using the Sherpa command COVARIANCE (S-Lang or Python help):

sherpa> covariance (S-Lang or Python help)()
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

sherpa> covariance (S-Lang or Python help)();
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.06324   -0.0196924    0.0196924
   modelA.ampl    0.00265019 -3.09154e-05  3.09154e-05

This quick method to calculate one-sigma parameter ranges, using COVAR (S-Lang or Python help), may underestimate the errors for a complex parameter space. The slower, but more appropriate, PROJECTION (S-Lang or Python help) 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

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*1.0;
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=NULL);
 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));
3.09376e+07

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

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=NULL);
sherpa> asum(get_counts(2));
10919

And finally, the second simulated dataset is plotted:

sherpa> PLOT_DATA (S-Lang or Python help)(2)

sherpa> PLOT_DATA (S-Lang or Python help)(2);

The second simulated dataset may also be written to disk:

sherpa> SAVE_PHA (S-Lang or Python help)(2,"dataB_sim1.pha")

sherpa> SAVE_PHA (S-Lang or Python help)(2,"dataB_sim1.pha");

To plot both simulated datasets at once:

sherpa> notice(0.5,8)      
sherpa>   # Plotting Both Simulated Datasets
sherpa> PLOT (S-Lang or Python help)("data",1,"data",2)

sherpa> notice(0.5,8)      
sherpa>   % Plotting Both Simulated Datasets
sherpa> PLOT (S-Lang or Python help)("data",1,"data",2);

Figure 4 shows the resulting plot.

To exit Sherpa:

sherpa> EXIT
Goodbye.
unix%

sherpa> QUIT

Using a Sherpa Script to Run Simulations

Multiple simulations with count-rate normalization

The following script does several FAKE_PHA simulations using the count rate for normalization; this is different than the previous examples which all use the flux instead. Note that this type of a script can be included in a more complex script with fake_time, energy and cnt_rate as parameters.

The response files used in the script - aciss_aimpt_cy10.rmf and aciss_aimpt_cy10.arf - can be downloaded from the ACIS Cycle 10 RMFs/ARFs CALDB page.

sherpa> !cat multi_sim.py

sherpa> !cat multi_sim.sl

Running this script from Sherpa:

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

sherpa> evalfile("multi_sim.sl");
"Count rate in specified energy range is"
0.01914
1

Where the "1" at the end indicates the script has run successfully to completion.

creates Figure 5 . Scripts may also be launched from the Unix command line as follows

unix% sherpa multi_sim.py

unix% sherpa multi_sim.sl

If your script ends with a EXITQUIT command, then the script may be run from the command line, with the screen outputs written to a file "multi_sim.out" as follows:

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

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

Simulating hot plasma emission

In this script, emission from a hot plasma is simulated with two different models. Note that the normalization is obtain within the script. The results of both simulations are then compared:

sherpa> $more thread.overallexample
RSP[acis](dataB_rmf.fits, dataB_arf.fits)
INSTRUMENT = acis
FAKEIT TIME = 5000
PARAMPROMPT OFF
SOURCE = xswabs[abs1]*xsmekal[m1]
abs1.1 = 0.03
m1.1 = 4.0
m1.norm = 0.025
FAKEIT
# Renormalize to the 2-10 keV flux of 1.e-13 erg/s/cm2
modflux1=get_eflux(1,[2,10])
obsflux1=1.e-13
foo1=get_par("m1.ampl")
foo1.value=foo1.value*(obsflux1/modflux1.value)
()=set_par(foo1)
FAKEIT
WRITE DATA sim_meka.pha PHA

INSTRUMENT 2 = acis
FAKEIT 2 TIME = 5000
SOURCE 2 = abs1 * xsapec[ap1]
ap1.kT = 4.0
ap1.norm = 0.025
FAKEIT 2
# SHOW
# Renormalize to the 2-10 keV flux of 1.e-13 erg/s/cm2
modflux2=get_eflux(2,[2,10])
obsflux2=1.e-13
foo2=get_par("ap1.ampl")
foo2.value=foo2.value*(obsflux2/modflux2.value)
()=set_par(foo2)
WRITE DATA 2 sim_apec.pha PHA
EXIT

This script can be run from the Unix command line as follows:

unix% sherpa thread.overallexample >&! thread.overallexample.out

The thread.overallexample.out file will then contain the screen output generated by the commands in the script.

History

08 Dec 2008

Created for CIAO/Sherpa4.1, replacing CIAO/Sherpa3.x FAKEIT command