Simulating X-ray Spectral Data (PHA): the fake_pha command
Sherpa Threads (CIAO 4.9 Sherpa v1)
The Sherpa fake_pha command can be used to simulate 1-D PHA data, e.g., for Chandra proposal planning. This thread provides both "quick start" and detailed "step-by-step" examples of the fake_pha functionlity.
If you have some experience with simulating X-ray spectral data in Sherpa, you may wish to skip the introductory material in this thread and jump to one of the more succint simulation threads tailored to a specific type of analysis, such as Simulating Chandra ACIS-S Spectra with Sherpa.
Last Update: 14 Dec 2015 - updated for CIAO 4.9 and Cycle 19, no content change.
- Quick Start to Simulating Data
- Step-by-Step: Getting started
- Introduction to fake_pha
- Defining a Source Model Expression for the Simulation
- Defining an Instrument Response
- Defining a Background (optional)
- Running the Simulation Using fake_pha
- Defining the Model Normalization for the Simulated Data
- Writing the Simulated Data to Output Files
- Using fake_pha with a PHA File
- Plotting and Fitting Simulated Data
- Simulating Multiple Data Sets
- Using a Sherpa Script to Run Simulations
- Scripting It
All that is required to simulate a 1-D PHA data set in Sherpa is the following:
- a source model expression defined with set_source
- a .rsp or ARF & RMF instrument response files
- exposure time for the simulation in seconds
fake_pha will do the rest:
unix% ciao ciao% sherpa # Search available Sherpa models, define a source model expression and assign it an ID such as "1". sherpa> list_models() ['atten', 'bbody', 'bbodyfreq', 'beta1d', 'beta2d', ... 'xszvarabs', 'xszvfeabs', 'xszvphabs', 'xszwabs', 'xszwndabs', 'xszxipcf'] sherpa> set_source(1, "xsphabs.abs1*xspowerlaw.p1") # Set initial model parameter values. sherpa> set_par(abs1.nH, val=0.07, frozen=True) sherpa> p1.PhoIndex = 0.8 # Fake a PHA data set over the grid defined by the input responses. sherpa> fake_pha(1, "source.arf", "source.rmf", 50000) # or sherpa> fake_pha(1, None, "source.rsp", 50000) # Filter the simulated data on energy or wavelength and plot. sherpa> notice(0.3, 7.) sherpa> plot_data() # Return information on the faked data set and associated responses. sherpa> show_data(1) Data Set: 1 Filter: 0.2993-7.0007 Energy (keV) Noticed Channels: 21-480 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: 1:1 name = source.rmf detchans = 1024 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: 1:1 name = source.arf energ_lo = Float64 energ_hi = Float64 specresp = Float64 bin_lo = None bin_hi = None exposure = 8034.23745459 # Save simulated data. sherpa> save_pha(1, "simulation1.pha")
This represents the quickest and simplest way to simulate a PHA data set in Sherpa: start CIAO and Sherpa; define a source model expression in Sherpa and assign it a data set ID such as "1" or "sim"; set some initial model parameter values; and then run fake_pha with your chosen exposure time and instrument response file(s).
Read on to learn the details of all options available for faking PHA data in Sherpa, with step-by-step instructions; or, skip to one of the more-succint simulation threads tailored to a specific type of analysis:
- Simulating Chandra ACIS-S Spectra with Sherpa
- Simulating Chandra ACIS-S LETG Spectra with Sherpa
- Simulating IXO/CAT Spectra with Sherpa
- Creating an Input Spectrum for Running MARX Simulations
- Simulating NuSTAR X-Ray Spectra with Sherpa
The fake_pha command creates a simulated 1-D PHA data set using a previously defined source model expression and instrument response; Poisson deviates are then added to the modeled data. If a background file is supplied via the 'bkg' argument of fake_pha, then the backscale correction is used, and background is added to the simulated data before Poisson noise. Note that while it is not necessary to supply a PHA data set to the 'id' argument of fake_pha—e.g., a real data set loaded into the session or an empty one created with the dataspace1d command—the 'bkg' argument does require such a PHA data set. All that is necessary for the id argument is a source model ID, i.e. the ID defined with the set_source command in defining the source model expression for the simuation (demonstrated in the section Defining a Source Model Expression without a Data Set).
All Sherpa models can be used in the simulations, as well as models created and implemented by the user. Individual models can be combined into composite models, and parameters can be linked across the model components.
Both source models and instrument responses must be pre-defined. It is possible to use both a redistribution matrix file (RMF) and an ancillary response file (ARF) with fake_pha, or just an RMF. The 'rmf' argument accepts both .rmf and .rsp files.
There are several steps involved in creating simulations:
- Loading an appropriate ancillary response file and/or redistribution matrix file.
- Defining a source model expression.
- (optional) Defining a background by reading a background data file or specifying a background model.
- Running the simulation using fake_pha.
- Defining the model normalization for the simulated data, if desired.
- Writing the simulated data to output files.
This thread illustrates each of these steps with an example. It also illustrates use of the fake_pha command when a source PHA data file has been previously read into the session, as well as how to plot and fit simulated data. The thread concludes with an example script for simulating multiple data sets within a single Sherpa session, as well as one for simulating hot plasma emission.
If you already have a source PHA data set loaded into the session, Sherpa automatically loads instrument response and background data files recorded in the header of a PHA data set when the PHA file is loaded 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 source PHA data set available, then proceed to the next section of this thread.
The source model expression is defined in the standard way with the set_source command. Before defining the model, one may use the dataspace1d command to create an empty PHA data set with which to associate the model—which will eventually store the faked data values—but this step is no longer necessary on account of enhancements to the fake_pha functionality. One may simply use the set_source command to assign a model to an ID, which will later be interpeted by fake_pha as the ID of the simulated data set (i.e, the ID defined with set_source should be supplied as the fake_pha 'id' value, at which point a data set will be created and associated with that ID).
sherpa> set_source(1, "powlaw1d.modelA") sherpa> modelA.gamma=2
In this example, we assume a 1-D power-law model for the source.
sherpa> print(get_source()) 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 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 specextract script, as described in the "Extract Spectrum and Response Files for" point-like, extended, or multiple source threads. Furthermore, calibration files for simulations can be downloaded from the Chandra Proposal Planning page of the CalDB (Calibration Database) website.
The ARF and RMF files may be loaded into the Sherpa session via the unpack_arf/unpack_rmf and load_arf/load_rmf commands. Here, we choose to use the unpack commands so that the response data can be stored in the variables arf1 and rmf1:
sherpa> arf1=unpack_arf("dataA_arf.fits") sherpa> print(arf1) name = dataA_arf.fits energ_lo = Float64 energ_hi = Float64 specresp = Float64 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 energ_hi = Float64 n_grp = UInt64 f_chan = UInt32 n_chan = UInt32 matrix = Float64 offset = 0 e_min = Float64 e_max = Float64
Background data and/or models are treated as follows with fake_pha:
If a background model 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 data set was background-subtracted prior to the command fake_pha being issued, it will not be background-subtracted afterwards.
If no background model is defined, and the data are background-subtracted, then the source model 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.
If no background model is defined, and the data are not background-subtracted, then the source model is evaluated directly, and the (properly scaled) background data are added to the faked data.
A type I PHA 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 appropriate scaling of the input background data.
sherpa> bkg1=unpack_pha("dataA_bkg.pha") sherpa> print(bkg1) name = dataA_bkg.pha channel = Float64 counts = Float64 staterror = None syserror = None bin_lo = None bin_hi = None grouping = None quality = Int16 exposure = 108675.66732 backscal = 0.044189453125 areascal = 1.0 grouped = False subtracted = False units = energy rate = True plot_fac = 0 response_ids =  background_ids = 
We have used the unpack_pha command to load the background PHA data set into the Sherpa session, and print to list the details of the data set. Note that we will not use the set_bkg command (or set_bkg_model in the next section) to associate the background with source data ID=1, as is usually done, since technically, data set 1 does not yet exist, only an ID defined with the set_source command. When fake_pha is run with 'id=1', a source data set will be created at that point, associated with ID=1. The set_bkg* commands cannot be used unless the source data set already exists. To work around this fact, we will set up the background component of source data ID=1 as a separate PHA data set, already named bkg1. It will be assigned as the background component of source data set 1 when the fake_pha command is run with the 'bkg' argument set to 'get_data(bkg1)'.
An instrument response does not need to be explicitly set for the background, as long as an ARF and RMF are associated with the source data set ID; this is because the scaled expected counts will be included before the source instrument model is convolved.
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). The simulated data is then sampled given the sum. In this instance, the background model is defined with the set_source command, as opposed to the set_bkg_model command, since source data ID=1 is not yet associated with a data set. We set as the background model a combination of a 1-D polynomial and Gaussian model:
sherpa> set_source("bkg1",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> print(get_source("bkg1")) (polynom1d.bkgA + gauss1d.bkgB) Param Type Value Min Max Units ----- ---- ----- --- --- ----- 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
Simulating a spectrum means taking the defined model expression, folding it through the instrument 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. To include the background component in the simulation and set the backscale and area scale values to those associated with the background PHA data set, we also use the optional 'bkg', 'backscal', and 'areascal' arguments. We choose to simulate a spectrum resulting from a 33 ks exposure.
sherpa> fake_pha(1, arf=arf1, rmf=rmf1, exposure=33483.2, bkg=bkg1, backscal=bkg1.backscal, areascal=bkg1.areascal)
This command associates data set ID "1" with a simulated, ungrouped source data set based on the assumed exposure time, instrument response, source model expression, and background data; Poisson noise is added to the modeled data. If an ARF is not to be included in the instrument response for the simulation, the 'arf' argument may be set to None.
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(1, arf="dataA_arf.fits", rmf="dataA_rmf.fits", exposure=33483.2, bkg=bkg1, backscal=bkg1.backscal, areascal=bkg1.areascal)
Furthermore, had we decided to load the ARF and RMF into the Sherpa session in association with a data set using the load_arf and load_rmf commands (or load_pha, as shown in the section Using fake_pha with a PHA File), then we could have set the 'arf' and 'rmf' arguments using a third option, shown below.
sherpa> fake_pha(1, arf=get_arf(), rmf=get_rmf(), exposure=33483.2, bkg=bkg1, backscal=bkg1.backscal, areascal=bkg1.areascal)
For detailed information on the available options for setting the 'arf' and 'rmf' arguments of fake_pha, refer to the fake_pha ahelp file.
Data Set: 1 Filter: 0.0276-14.6640 Energy (keV) Bkg Scale: 0.308102 Noticed Channels: 1-512 name = faked channel = Float64 counts = Float64 staterror = None syserror = None bin_lo = None bin_hi = None grouping = None quality = None exposure = 33483.2 backscal = 0.044189453125 areascal = 1.0 grouped = False subtracted = False units = energy rate = True plot_fac = 0 response_ids =  background_ids =  RMF Data Set: 1:1 name = dataA_rmf.fits detchans = 512 energ_lo = Float64 energ_hi = Float64 n_grp = UInt64 f_chan = UInt32 n_chan = UInt32 matrix = Float64 offset = 0 e_min = Float64 e_max = Float64 ARF Data Set: 1:1 name = dataA_arf.fits energ_lo = Float64 energ_hi = Float64 specresp = Float64 bin_lo = None bin_hi = None exposure = None sherpa> sum(get_counts(1)) 4917112.0 sherpa> calc_data_sum(id=1) 4917112.0
The latter two commands are equivalent, giving the total number of counts simulated in the data set. 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 data set 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 data set in units of counts/second by doing the following:
sherpa> calc_data_sum(id=1)/get_exposure(1) 146.85310842452336
In order to use simulated data for scientific analysis, the data should be re-normalized to match the flux (or total counts) of the observed source. If the flux of the simulated data is known, the parameter values of the model used to create the simulated data can easily be converted to the appropriate values. This requires first simulating the data with the default model parameters, calculating the resulting simulated flux, and then rescaling the model parameters appropriately. Then, the simulation is run again with the updated model.
Assuming that we have just simulated data using default model parameters in the source model expression, with modelA.ampl=1, we calculate the flux accordingly:
sherpa> reset(modelA) sherpa> fake_pha(1, arf=arf1, rmf=rmf1, exposure=33483.2, bkg=bkg1, backscal=bkg1.backscal, areascal=bkg1.areascal) sherpa> calc_energy_flux(lo=2.,hi=10.) 2.5706384381817513e-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, we divide a given true source flux of 1.0e-13 ergs cm-2 s-1, by the calculated flux for the above default source model (2.5706384381817513e-09); the resulting value will replace the current normalization (modelA.ampl) value of 1 photon keV-1 cm-2 s-1:
sherpa> 1.0e-13/calc_energy_flux(lo=2.,hi=10.) 3.8900842107819488e-05
The normalization of the source model may be adjusted as follows:
and fake_pha is run again with the updated source model:
sherpa> fake_pha(1, arf=arf1, rmf=rmf1, exposure=33483.2, bkg=bkg1, backscal=bkg1.backscal, areascal=bkg1.areascal) sherpa> calc_photon_flux(lo=2.,hi=10.) 1.5463634382489859e-05 sherpa> calc_energy_flux(lo=2.,hi=10.) 9.999999999999999e-14
In addition to the energy flux, the photon flux between 2.0-10.0 keV (in units of photons cm-2 s-1) is also available, via the calc_photon_flux command.
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=33483.2, bkg=bkg1, backscal=bkg1.backscal, areascal=bkg1.areascal) sherpa> rate=sum(get_counts(1))/get_exposure(1) sherpa> print(rate) 146.733794858 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> calc_data_sum(1) 285.0 sherpa> calc_data_sum(id=1)/get_exposure(1) 0.42654226597218908
The data set could then be normalized by scaling the model amplitude with this non-normalized count rate (146.8 cts/sec); this would be achieved by setting the model amplitude value equal to the known count rate divided by the non-normalized count rate.
By default, fake_pha adds Poisson noise to the synthesized data, but it is possible to suppress the noise using two possible approaches:
- artifically increase the time so Poisson noise becomes irrelevant, or
- evaluate the given model, and take the predicted model counts as the data, as shown below.
sherpa> set_source("faked","powlaw1d.modelA") sherpa> fake_pha("faked",arf="arf1.fits",rmf="rmf1.fits",exposure=33483.2,grouped=False) sherpa> mdl=get_model("faked") sherpa> spectral_coords=get_indep("faked") sherpa> counts=mdl(*spectral_coords) sherpa> set_counts("faked",counts)
While model evaluation will eliminate the noise, the counts will be non-integer valued. If integer counts are needed, either artificially increase the counts or round the model array to the nearest integer.
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 command:
This creates a PHA file with columns of channel and counts. This file may then be grouped using dmgroup if necessary.
There is also the option to save the data to a FITS or ASCII table file with the save_arrays command:
sherpa> save_arrays("dataA_sim1.fits", [get_model_plot(1).xlo, get_model_plot(1).y], ascii=False) sherpa> save_arrays("dataA_sim1.txt", [get_model_plot(1).xlo, get_model_plot(1).y], ascii=True)
If the user has previously read data from a PHA file:
- An instrument response 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 data set will be overwritten by the simulated data set created by fake_pha.
For example, erase all current Sherpa settings with the clean command (or by quitting the session and starting a new one), 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 response being automatically defined (using RMF and ARF files), and a background data file being automatically input. Furthermore, from the get_exposure and get_backscal commands, we see that the exposure time and backscale are readily accessed.
That is, input of this PHA data file has automatically completed two steps involved in creating simulations:
- Defining an instrument response using a redistribution matrix file and an ancillary response file.
- Reading a background data file.
The remaining steps are:
- 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.
The following commands complete the remaining steps listed above:
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()) sherpa> calc_data_sum(id=1) 4919163.0 sherpa> modelA.ampl=0.4/(calc_data_sum(id=1)/get_exposure(1)) sherpa> fake_pha(1,arf=get_arf(1),rmf=get_rmf(1),exposure=get_exposure(1),bkg=get_bkg(1)) sherpa> calc_data_sum(id=1) 14030.0
Note that the fake_pha command overwrites the input data set with the simulated data, and we have further normalized the simulated data set.
The results may be written to a PHA file:
The simulated data set may be plotted as any other data set in Sherpa:
Figure 2 shows the resulting plot.
The simulated data set may also be filtered and fit as any other data set in Sherpa. For example, we can filter the the simulated data set to include only data between 0.5 and 8.0 keV:
Figure 3 shows the resulting plot.
Then, we fit the simulated data set using the same source model expression that was used to create it:
sherpa> subtract(1) sherpa> set_method("neldermead") sherpa> set_stat("chi2xspecvar") sherpa> fit(1) Dataset = 1 Method = neldermead Statistic = chi2xspecvar Initial fit statistic = 329.997 Final fit statistic = 312.549 at function evaluation 295 Data points = 263 Degrees of freedom = 261 Probability [Q-value] = 0.0157183 Reduced statistic = 1.19751 Change in statistic = 17.448 modelA.gamma 2.04458 modelA.ampl 0.00269349 sherpa> plot_fit()
Figure 4 shows the resulting plot.
Here, we examine the fit using the Sherpa covariance command:
sherpa> covariance() Dataset = 1 Confidence Method = covariance Iterative Fit Method = None Fitting Method = neldermead Statistic = chi2xspecvar covariance 1-sigma (68.2689%) bounds: Param Best-Fit Lower Bound Upper Bound ----- -------- ----------- ----------- modelA.gamma 2.04458 -0.015838 0.015838 modelA.ampl 0.00269349 -2.87148e-05 2.87148e-05
This quick method of calculating 1σ parameter ranges using covar() may underestimate the errors for a complex parameter space. The slower, but more appropriate, confidence algorithm should be used in such cases.
sherpa> conf() modelA.ampl lower bound: -2.87148e-05 modelA.gamma lower bound: -0.015838 modelA.ampl upper bound: 2.87148e-05 modelA.gamma upper bound: 0.015838 Dataset = 1 Confidence Method = confidence Iterative Fit Method = None Fitting Method = neldermead Statistic = chi2xspecvar confidence 1-sigma (68.2689%) bounds: Param Best-Fit Lower Bound Upper Bound ----- -------- ----------- ----------- modelA.gamma 2.04458 -0.015838 0.015838 modelA.ampl 0.00269349 -2.87148e-05 2.87148e-05
One may simulate multiple data sets within the same Sherpa session; each data set is assigned an individual ID.
For example, we can simulate a second data set as we did the first, above, by defining the instrument response with the ARF and RMF, defining a source model expression, and matching the exposure, backscale, and area scale to the background data set:
sherpa> arf2=unpack_arf("dataB_arf.fits") sherpa> rmf2=unpack_rmf("dataB_rmf.fits") sherpa> set_source(2,xsphabs.modelB1*powlaw1d.modelB2) sherpa> modelB1.nH=0.025 sherpa> modelB2.gamma=1.7 sherpa> modelB2.ampl=1.0 sherpa> bkg2=unpack_pha("dataB_bkg.pha") WARNING: systematic errors were not found in file 'dataB_bkg.pha' statistical errors were found in file 'dataB_bkg.pha' but not used; to use them, re-read with use_errors=True sherpa> fake_pha(2,arf=arf2,rmf=rmf2,exposure=bkg2.exposure, bkg=bkg2, backscal=bkg2.backscal, areascal=bkg2.areascal) sherpa> calc_data_sum(id=2) 30935237.0
We will now normalize the second simulated data set in the same manner as we did the first: dividing the known count rate, e.g. 0.4, by the non-normalized count rate to obtain a scaling factor, and then re-simulating the data set.
sherpa> modelB2.ampl=0.4/(calc_data_sum(id=2)/get_exposure(2)) sherpa> fake_pha(2, arf=arf2,rmf=rmf2,exposure=bkg2.exposure, bkg=bkg2, backscal=bkg2.backscal, areascal=bkg2.areascal) sherpa> calc_data_sum(id=2) 11542.0
Finally, the second simulated data set can be plotted as follows:
and written to file:
sherpa> save_pha(2,"dataB_sim1.pha") sherpa> save_arrays("dataA_sim2.fits", [get_model_plot(2).xlo, get_model_plot(2).y], ascii=False) sherpa> save_arrays("dataA_sim2.txt", [get_model_plot(2).xlo, get_model_plot(2).y], ascii=True)
To plot both simulated data sets at once:
sherpa> notice(0.5,8) sherpa> plot("data",1,"data",2)
Figure 5 shows the resulting plot.
To exit Sherpa:
sherpa> quit() unix%
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. It also uses the dataspace1d command to create an empty source PHA data set before the simulaton, also different than the previous examples (since this step is not required to simulate a source data set). This script defines the variables fake_time, energy, and cnt_rate, which means it can be included in a more complex script where these variables are used as parameters.
The response files used in the script—aciss_aimpt_cy19.rmf and aciss_aimpt_cy19.arf—can be downloaded from the ACIS Cycle 19 RMFs/ARFs CALDB page.
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_cy19.rmf") arf = unpack_arf("aciss_aimpt_cy19.arf") set_source(1, powlaw1d.pow1) pow1.gamma = 1.9 # Define structure to allow manipulation 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 # Run fake_pha once to get model counts for initial amplitude, then # normalize to get desired count rate fake_pha(1,arf=arf,rmf=rmf,exposure=fake_time) 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) # Verify correct counts in fake data set and plot data print("Count rate in specified energy range is") print(calc_data_sum(lo=elo,hi=ehi)/fake_time) plot_data()
Running this script from Sherpa:
sherpa> exec(open("multi_sim.py").read()) Count rate in specified energy range is 0.01976
Scripts may also be launched from the Unix command line as follows
unix% sherpa multi_sim.py
If a script ends with quit(), then running the script from the command line as follows will write the screen output generated by the script to a file, such as multi_sim.out:
unix% sherpa multi_sim.py >&! multi_sim.out
In this script, emission from a hot plasma is simulated using two different models. Note that the appropriate normalization is calculated within the script. The results of both simulations are then compared:
sherpa> !cat plasma.py#!/usr/bin/env python from sherpa.astro.ui import * # This script simulates emission from a hot plasma using 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 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) save_pha(1,"sim_meka.pha") set_source(2,xswabs.abs1*xsapec.ap1) ap1.kT=4.0 ap1.norm=0.025 fake_pha(2,arf=arf,rmf=rmf,exposure=fake_time) # 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) save_pha(2,"sim_apec.pha") exit()
This script can be run from the Unix command line as follows:
unix% sherpa plasma.py >&! plasma.out
The plasma.out file will contain the screen output generated by the commands in the script.
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.
|08 Dec 2008||Created for CIAO/Sherpa 4.1, replacing CIAO/Sherpa 3.x FAKEIT command|
|29 Apr 2009||new script command is available with CIAO 4.1.2|
|18 Jan 2010||Updated for CIAO 4.2: simulated data is now ungrouped by default, and the save_arrays and confidence commands are available; examples modified to show that it is not necessary to run dataspace1d to fake a source data set.|
|13 Jul 2010||updated for CIAO 4.2 Sherpa v2: removal of S-Lang version of thread.|
|08 Sep 2010||figures moved inline with text|
|15 Dec 2010||updated with "quick start" section|
|15 Dec 2011||reviewed for CIAO 4.4: a work-around for a save_pha bug was added; response files used in examples were updated for Chandra proposal cycle 13|
|14 Mar 2012||updated with links to the CALDB Chandra Proposal Planning page, and the other simulation threads|
|13 Dec 2012||reviewed for CIAO 4.5: updated introductory information in the 'Define a Background' section|
|12 Dec 2013||reviewed for CIAO 4.6: updated screen output|
|15 Apr 2014||simulating without Poisson noise.|
|06 Feb 2015||updated for CIAO 4.7 and Cycle 17, no content change.|
|14 Dec 2015||updated for CIAO 4.8 and Cycle 18, no content change.|
|14 Dec 2015||updated for CIAO 4.9 and Cycle 19, no content change.|