Simulating Chandra ACIS-S Spectra with Sherpa

Sherpa Threads (CIAO 4.1)

[Python Syntax]

Overview

Last Update: 29 Apr 2009 - new script command is available with CIAO 4.1.2

Synopsis:

This thread illustrates the use of the Sherpa fake_pha command to simulate a spectrum of a point source obtained with the ACIS-S detector aboard Chandra, with and without consideration of a background component.

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

Getting started - downloading calibration response files for simulations

In order to simulate a Chandra ACIS-S spectrum with Sherpa, we must define an instrument response with the appropriate ARF (auxiliary response function) and RMF (redistribution matrix function) files. These files may be downloaded from the Chandra Proposal Planning page of the CalDB (Calibration Database) website, where sample RMFs and corresponding ARFs positioned at the aimpoint of the ACIS-S array (and at selected off-axis points) are available.

In this thread, we use the files aciss_aimpt_cy11.rmf and aciss_aimpt_cy11.arf, positioned at the default pointing for ACIS-S.

The Sherpa fake_pha command calculates a synthetic 1-D spectral data set based on a defined instrument response and source model expression. For details on the functionality of this command, with examples of its usage, see the thread Simulating 1-D Data: the Sherpa FAKE_PHA Command.

Defining the Instrument Response

without a background component

with a background component

We begin by establishing the instrument response corresponding to the default pointing of the ACIS-S detector:

sherpa> arf1=unpack_arf("aciss_aimpt_cy11.arf")

sherpa> print(arf1)
name     = aciss_aimpt_cy11.arf
energ_lo = Float64[1024]
energ_hi = Float64[1024]
specresp = Float64[1024]
bin_lo   = None
bin_hi   = None
exposure = 8034.23745459

sherpa> rmf1=unpack_rmf("aciss_aimpt_cy11.rmf")

sherpa> print(rmf1)
name     = aciss_aimpt_cy11.rmf
detchans = 1024
energ_lo = Float64[1024]
energ_hi = Float64[1024]
n_grp    = Int16[1024]
f_chan   = UInt32[1504]
n_chan   = UInt32[1504]
matrix   = Float64[387227]
offset   = 1
e_min    = Float64[1024]
e_max    = Float64[1024]

Here, the ARF and RMF data are loaded and assigned to a set of variables with the unpack_* commands. These variables will be used to assign the instrument response to the faked data set we will create in the next section, "Defining a Source Model Expression --> without a background component".

with a background component

Here we establish the instrument response corresponding to the default pointing of the ACIS-S detector, for both a source and background component:

sherpa> arf1=unpack_arf("aciss_aimpt_cy11.arf")

sherpa> bkg1_arf = arf1

sherpa> print(arf1)
name     = aciss_aimpt_cy11.arf
energ_lo = Float64[1024]
energ_hi = Float64[1024]
specresp = Float64[1024]
bin_lo   = None
bin_hi   = None
exposure = 8034.23745459

sherpa> rmf1=unpack_rmf("aciss_aimpt_cy11.rmf")

sherpa> bkg1_rmf = rmf1

sherpa> print(rmf1)
name     = aciss_aimpt_cy11.rmf
detchans = 1024
energ_lo = Float64[1024]
energ_hi = Float64[1024]
n_grp    = Int16[1024]
f_chan   = UInt32[1504]
n_chan   = UInt32[1504]
matrix   = Float64[387227]
offset   = 1
e_min    = Float64[1024]
e_max    = Float64[1024]

The source ARF and RMF data are loaded and assigned to a set of variables with the unpack_* commands. These variables will be used to assign the instrument response to both the source and background components of the faked data set we will create in the next section, "Defining Source Model Expressions --> with a background component". If the background response is different than the source response, we load the appropriate background ARF and RMF files accordingly:

sherpa> bkg1_rmf = unpack_rmf("background.rmf"); # separate background response
sherpa> bkg1_arf = unpack_arf("background.arf");

Defining a Source Model Expression

without a background component

with a background component

Now that we have defined the ARF and RMF instrument responses, we establish a source model expression for our simulation. To do this, we must first use the Sherpa dataspace1d command to create a blank 1-D PHA data array with as many elements as our instrument response (1024). Then, we assign the ARF, RMF, and source model to the blank data set, which we have given an ID of "1":


sherpa> dataspace1d(0.0, 1024.0, id=1, dstype=DataPHA) 

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

sherpa> set_source(1, xsphabs.abs1*powlaw1d.m1)
sherpa> m1.gamma=2
sherpa> abs1.nh=0.2

sherpa> show_model()
Model: 1
apply_rmf(apply_arf((xsphabs.abs1 * powlaw1d.m1)))
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   abs1.nh      thawed          0.2            0       100000 10^22 atoms / cm^2
   m1.gamma     thawed            2          -10           10           
   m1.ref       frozen            1 -3.40282e+38  3.40282e+38           
   m1.ampl      thawed            1            0  3.40282e+38

We have defined a source model expression with the set_source command, using for this simulation an absorbed 1-D power law model with a Galactic neutral hydrogen column density of 2e21 cm^-2 and a power law photon index of 2.

with a background component

Here we use the dataspace1d command to create an additional data set (id=2) to which to assign the background model expression for our simulation:


sherpa> dataspace1d(0.0, 1024.0, id=1, dstype=DataPHA) 

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

sherpa> set_source(1, xsphabs.abs1*powlaw1d.m1)
sherpa> m1.gamma=2
sherpa> abs1.nh=0.2

sherpa> dataspace1d(0.0,1024.0,id=2,dstype=DataPHA) 

sherpa> set_arf(2, bkg1_arf)
sherpa> set_rmf(2, bkg1_rmf)

sherpa> set_model(2, polynom1d.bkgA)
sherpa> bkgA.c0=1.

sherpa> show_model()
Model: 1
apply_rmf(apply_arf((xsphabs.abs1 * powlaw1d.m1)))
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   abs1.nh      thawed          0.2            0       100000 10^22 atoms / cm^2
   m1.gamma     thawed            2          -10           10           
   m1.ref       frozen            1 -3.40282e+38  3.40282e+38           
   m1.ampl      thawed            1            0  3.40282e+38           

Model: 2
apply_rmf(apply_arf(polynom1d.bkgA))
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   bkgA.c0      thawed            1 -3.40282e+38  3.40282e+38           
   bkgA.c1      frozen            0 -3.40282e+38  3.40282e+38           
   bkgA.c2      frozen            0 -3.40282e+38  3.40282e+38           
   bkgA.c3      frozen            0 -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

We have defined source and background model expressions with the set_source and set_model commands. For the source simulation data set, we use an absorbed 1-D power law model with a Galactic neutral hydrogen column density of 2e21 cm^-2 and a photon index of 2. For the background simulation data set, we assume a flat profile with a 1-D polynomial function.

Running the Simulation with fake_pha

without a background component

with a background component

Simulating the Chandra spectrum means taking the defined model expression, folding it through the Chandra ACIS-S response, and applying Poisson noise to the counts predicted by the model. The simulation is run with fake_pha, which has four required arguments: dataset ID, ARF, RMF, and exposure time. We decide to simulate an ACIS-S spectrum resulting from a 50 ks exposure of a point source.

sherpa> fake_pha(1, arf1, rmf1, exposure=50000, grouped=False)

This command replaces the empty PHA data set id=1 with a simulated data set based on the assumed exposure time, instrument response, and source model expression; Poisson noise is added to the modeled data.

We may inspect some basic properties of the new data set with the show_data command:

sherpa> show_data()
Data Set: 1
Filter: 0.0110-14.9431 Energy (keV)
Noticed Channels: 1-1024
name           = faked
channel        = Int32[1024]
counts         = Float64[1024]
staterror      = None
syserror       = None
bin_lo         = Float64[1024]
bin_hi         = Float64[1024]
grouping       = None
quality        = None
exposure       = 50000
backscal       = None
areascal       = None
grouped        = False
subtracted     = False
units          = energy
response_ids   = [1]
background_ids = []

RMF Data Set: 1:1
name     = aciss_aimpt_cy11.rmf
detchans = 1024
energ_lo = Float64[1024]
energ_hi = Float64[1024]
n_grp    = Int16[1024]
f_chan   = UInt32[1504]
n_chan   = UInt32[1504]
matrix   = Float64[387227]
offset   = 1
e_min    = Float64[1024]
e_max    = Float64[1024]

ARF Data Set: 1:1
name     = aciss_aimpt_cy11.arf
energ_lo = Float64[1024]
energ_hi = Float64[1024]
specresp = Float64[1024]
bin_lo   = None
bin_hi   = None
exposure = 8034.23745459

However, the simulated data set currently does not have the correct normalization - the flux of the simulated data is incorrect because the default power law normalization is arbitrarily set to 1.0, as shown with the show_model command:

sherpa> show_model()
Model: 1
apply_rmf(apply_arf((50000 * (xsphabs.abs1 * powlaw1d.m1))))
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   abs1.nh      thawed          0.2            0       100000 10^22 atoms / cm^2
   m1.gamma     thawed            2          -10           10
   m1.ref       frozen            1 -3.40282e+38  3.40282e+38
   m1.ampl      thawed            1            0  3.40282e+38

To correct the flux we need to adjust the normalization, as demonstrated in the section "Defining the Model Normalization for the Simulated Data --> without a background component".

with a background component

Here we run an additional fake_pha simulation for the background data set:

sherpa> fake_pha(1, arf1, rmf1, exposure=50000, grouped=False)

sherpa> fake_pha(2, bkg1_arf, bkg1_rmf, 50000, grouped=False)

These commands replace the empty PHA data sets id=1 and id=2 with simulated source and background data sets, respectively, based on the assumed exposure times, instrument responses, and model expressions; Poisson noise is added to the modeled data.

We may inspect some basic properties of the new simulated data sets with the show_data command:

sherpa> show_data()
Data Set: 1
Filter: 0.0110-14.9431 Energy (keV)
Noticed Channels: 1-1024
name           = faked
channel        = Int32[1024]
counts         = Float64[1024]
staterror      = None
syserror       = None
bin_lo         = Float64[1024]
bin_hi         = Float64[1024]
grouping       = None
quality        = None
exposure       = 50000
backscal       = None
areascal       = None
grouped        = False
subtracted     = False
units          = energy
response_ids   = [1]
background_ids = []

RMF Data Set: 1:1
name     = aciss_aimpt_cy11.rmf
detchans = 1024
energ_lo = Float64[1024]
energ_hi = Float64[1024]
n_grp    = Int16[1024]
f_chan   = UInt32[1504]
n_chan   = UInt32[1504]
matrix   = Float64[387227]
offset   = 1
e_min    = Float64[1024]
e_max    = Float64[1024]

ARF Data Set: 1:1
name     = aciss_aimpt_cy11.arf
energ_lo = Float64[1024]
energ_hi = Float64[1024]
specresp = Float64[1024]
bin_lo   = None
bin_hi   = None
exposure = 8034.23745459

Data Set: 2
Filter: 0.0110-14.9431 Energy (keV)
Noticed Channels: 1-1024
name           = faked
channel        = Int32[1024]
counts         = Float64[1024]
staterror      = None
syserror       = None
bin_lo         = Float64[1024]
bin_hi         = Float64[1024]
grouping       = None
quality        = None
exposure       = 50000
backscal       = None
areascal       = None
grouped        = False
subtracted     = False
units          = energy
response_ids   = [1]
background_ids = []

RMF Data Set: 2:1
name     = aciss_aimpt_cy11.rmf
detchans = 1024
energ_lo = Float64[1024]
energ_hi = Float64[1024]
n_grp    = Int16[1024]
f_chan   = UInt32[1504]
n_chan   = UInt32[1504]
matrix   = Float64[387227]
offset   = 1
e_min    = Float64[1024]
e_max    = Float64[1024]

ARF Data Set: 2:1
name     = aciss_aimpt_cy11.arf
energ_lo = Float64[1024]
energ_hi = Float64[1024]
specresp = Float64[1024]
bin_lo   = None
bin_hi   = None
exposure = 8034.23745459

In the next section, we will correct the normalization of the simulated source and background data sets.

Defining the Model Normalization for the Simulated Data

without a background component

with a background component

Before we can use the simulated data set for scientific analysis, it must be re-normalized to match the flux (or total counts) required by our selected source.

The .2-10 keV flux in the simulated spectrum is 3.86e-9 ergs cm^-2 s^-1:

sherpa> calc_energy_flux(0.2,10)
3.8623970060321874e-09

The .2-10 keV flux of a source in our Chandra proposal, for example, has been measured at 1.e-12 ergs cm^-2 s^-1. Therefore, the correct normalization is (1.e-12)/(3.8623970060321874e-09)=0.000258906580147:

sherpa> my_flux=1.e-12    

sherpa> norm=my_flux/calc_energy_flux(0.2,10)  

sherpa> print(norm)
0.000258906580147

sherpa> set_par(m1.ampl,norm)

sherpa> show_model()
Model: 1
apply_rmf(apply_arf((50000 * (xsphabs.abs1 * powlaw1d.m1))))
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   abs1.nh      thawed          0.2            0       100000 10^22 atoms / cm^2
   m1.gamma     thawed            2          -10           10
   m1.ref       frozen            1 -3.40282e+38  3.40282e+38
   m1.ampl      thawed  0.000258907            0  3.40282e+38

sherpa> fake_pha(1,arf1,rmf1,exposure=50000)

sherpa> prefs = get_data_plot_prefs()

sherpa> prefs["yerrorbars"] = 0      # remove y-error bars from plot

sherpa> plot_data()

With the new normalization, the simulated flux is correctly set at the measured flux of 1.e-12 ergs cm^-2 s^-1. A plot of the data is shown in Figure 1.

[Print media version: Plot of simulated source spectrum]

Figure 1: Plot of simulated source spectrum

Note that we could have chosen to re-normalize the simulated data set to match the required total counts instead of flux. For example:

sherpa> my_counts=10000

sherpa> norm_counts=my_counts/calc_data_sum(0.5,8.)  

sherpa> print(norm_counts)
0.000432559029601

with a background component

We determine the normalization for the background data set in the same way as with the source data set, except we use a measure of total counts instead of flux to specify that we want 200 counts in the background simulation:

sherpa> my_flux=1.e-12    
sherpa> norm=my_flux/calc_energy_flux(0.2,10,id=1)  
sherpa> print(norm)
0.000258906580147

sherpa> bkg_counts = 200
sherpa> bkg_norm = bkg_counts/calc_data_sum(0.5,8.,id=2)
sherpa> print(bkg_norm)
1.63598181011e-06

Now we apply the calculated values to the amplitude parameters of each model, and re-evaluate the simulated data sets with the desired normaliztion using fake_pha:

sherpa> set_par(m1.ampl,norm)
sherpa> set_par(bkgA.c0,bkg_norm)

sherpa> show_model()
Model: 1
apply_rmf(apply_arf((50000 * (xsphabs.abs1 * powlaw1d.m1))))
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   abs1.nh      thawed          0.2            0       100000 10^22 atoms / cm^2
   m1.gamma     thawed            2          -10           10           
   m1.ref       frozen            1 -3.40282e+38  3.40282e+38           
   m1.ampl      thawed  0.000258907            0  3.40282e+38           

Model: 2
apply_rmf(apply_arf((50000 * polynom1d.bkgA)))
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   bkgA.c0      thawed  1.63598e-06 -3.40282e+38  3.40282e+38           
   bkgA.c1      frozen            0 -3.40282e+38  3.40282e+38           
   bkgA.c2      frozen            0 -3.40282e+38  3.40282e+38           
   bkgA.c3      frozen            0 -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         

sherpa> fake_pha(1,arf1,rmf1,exposure=50000)

sherpa> fake_pha(2,bkg1_arf,bkg1_rmf,50000)

Finally, we assign background data set id=2 as the background component of source data set id=1 with fake_pha, to produce a single source-plus-background simulated data set:

sherpa> mybkg = get_data(2)

sherpa> fake_pha(1,arf1,rmf1,exposure=50000,bkg=mybkg)

sherpa> prefs = get_data_plot_prefs()

sherpa> prefs[yerrorbars] = 0      # remove y-error bars from plot

sherpa> plot_data()

The resulting plot is shown in Figure 2.

[Plot of simulated source+background spectrum]

[Print media version: Plot of simulated source+background spectrum]

Figure 2: Plot of simulated source-plus-background spectrum

Writing the Simulated Data to Output Files

We may use the save_pha command to write the simulated data as a PHA file, with a header containing the exposure time value and paths to the ARF and RMF files:

sherpa> save_pha("simulation1.pha")

We also have the option to save the data to an ASCII table:

sherpa> load_arrays("my_sim_data",get_model_plot(1).xlo,get_model_plot(1).y)
sherpa> save_data("my_sim_data","simulation1.dat")

Fitting the Simulated Data

without a background component

with a background component

The simulated data set may be filtered and fit as any other data set in Sherpa. For example, we can choose to filter the simulated data to include only the counts in a restricted energy range, such as 0.5 keV - 7.0 keV.:

sherpa> print calc_energy_flux(.2,10)  # ergs cm^-2 s^-1
1e-12

sherpa>  print calc_energy_flux(.5,7)
8.3483818796e-13

sherpa> print calc_data_sum(.5,7)  # counts
5949.0 

sherpa> notice(0.5,7)

sherpa> show_filter()
Data Set Filter: 1
0.5037-7.0007 Energy (keV)

Then, we can fit the simulated data set with the source model expression we used to create it:

sherpa> set_method("neldermead")
sherpa> set_stat("cstat")

sherpa> fit()
Dataset               = 1
Method                = neldermead
Statistic             = cstat
Initial fit statistic = 494.794
Final fit statistic   = 494.794 at function evaluation 440
Data points           = 446
Degrees of freedom    = 443
Probability [Q-value] = 0.0446884
Reduced statistic     = 1.11692
Change in statistic   = 0
   abs1.nh        0.204026    
   m1.gamma       1.99511     
   m1.ampl        0.000259167 

sherpa> plot_fit()
WARNING: unable to calculate errors using current statistic: cstat

The resulting plot is shown in Figure 3.

[Print media version: Plot of fit to simulated source spectrum]

Figure 3: Plot of fit to simulated source spectrum

Next, we examine the quality of the fit with projection, and print the results of the fit with get_fit_results().

sherpa> projection()
Dataset               = 1
Confidence Method     = projection
Fitting Method        = neldermead
Statistic             = cstat
projection 1-sigma (68.2689%) bounds:
   Param            Best-Fit  Lower Bound  Upper Bound
   -----            --------  -----------  -----------
   abs1.nh          0.204026   -0.0133087    0.0134127
   m1.gamma          1.99511   -0.0432112    0.0434606
   m1.ampl       0.000259167 -1.05413e-05  1.10058e-05


sherpa> print get_fit_results()
datasets   = (1,)
methodname = neldermead
statname   = cstat
succeeded  = True
parnames   = ('abs1.nh', 'm1.gamma', 'm1.ampl')
parvals    = (0.20402568146670663, 1.9951139821543675, 0.00025916707008903498)
covarerr   = None
statval    = 494.79355811
istatval   = 494.79355811
dstatval   = 0.0
numpoints  = 446
dof        = 443
qval       = 0.0446883744616
rstat      = 1.11691548106
message    = Optimization terminated successfully
nfev       = 440

sherpa> print get_proj_results()
datasets   = (1,)
methodname = projection
fitname    = neldermead
statname   = cstat
sigma      = 1
percent    = 68.2689492137
parnames   = ('abs1.nh', 'm1.gamma', 'm1.ampl')
parvals    = (0.20402568146670663, 1.9951139821543675, 0.00025916707008903498)
parmins    = (-0.013308665108096251, -0.043211177751672381, -1.0541331230766168e-05)
parmaxes   = (0.01341271582614093, 0.043460608451337548, 1.1005800676672011e-05)
nfits      = 50

Note that the cstat statistic is appropriate for fitting low-counts data, but it does not calculate errors for the data points. We can group the data so that each bin contains a specified minimum number of counts, and then change the fit statistic to something more suitable to calculate the errors. Finally, we can view the results of the new fit with the plot_fit_delchi command:

sherpa> group_counts(15) 
WARNING: grouping flags have changed, noticing all bins

sherpa> notice(0.5,7)  #re-establish energy filter

sherpa> set_stat("chi2xspecvar")

sherpa> plot_fit_delchi()

The new fit to the grouped simulated spectrum, along with the residuals divided by the uncertainties, is shown in Figure 4.

[Print media version: Plot of fit to simulated source spectrum, with residuals]

Figure 4: Plot of fit to simulated source spectrum, with residuals

The plot may be saved as a PostScript file with the ChIPS print_window command:

sherpa> print_window("simulation_fit")

with a background component

Here, we filter the simulated source-plus-background data set (id=1) to include only the counts in the energy range 0.5 keV - 7.0 keV (recalling that data set id=1 now contains the background information stored in data set id=2; id=2 is no longer needed in the context of this thread).

sherpa> notice(0.5,7)  

sherpa> show_filter()
Data Set Filter: 1
0.4964-7.0518 Energy (keV)

Data Set Filter: 2
0.5037-7.0007 Energy (keV)

Next, we fit the simulated source data with the source model expression we used to create it, and use the set_bkg_model command to incorporate the background model into the fit:

sherpa> set_bkg_model(1, bkgA, 1)  # set model for data set id=1, bkg_id=1

sherpa> set_method("neldermead")
sherpa> set_stat("cstat")

sherpa> fit(1)  

Dataset               = 1
Method                = neldermead
Statistic             = cstat
Initial fit statistic = 616.563
Final fit statistic   = 596.406 at function evaluation 4087
Data points           = 636
Degrees of freedom    = 632
Probability [Q-value] = 0.84177
Reduced statistic     = 0.94368
Change in statistic   = 20.1574
   abs1.nh        0.218007    
   m1.gamma       2.14133     
   m1.ampl        0.000272952 
   bkgA.c0        1.73323e-06 

sherpa> plot_fit()
WARNING: unable to calculate errors using current statistic: cstat

The resulting plot is shown in Figure 5.

[Print media version: Plot of fit to simulated source-plus-background spectrum]

Figure 5: Plot of fit to simulated source-plus-background spectrum

Now we can examine the quality of the fit with projection, and return the results of the fit with get_fit_results().

sherpa> projection(1)
Dataset               = 1
Confidence Method     = projection
Fitting Method        = neldermead
Statistic             = cstat
projection 1-sigma (68.2689%) bounds:
   Param            Best-Fit  Lower Bound  Upper Bound
   -----            --------  -----------  -----------
   abs1.nh          0.218007   -0.0138014    0.0138177
   m1.gamma          2.14133    -0.048666    0.0483401
   m1.ampl       0.000272952 -1.16871e-05  1.21085e-05
   bkgA.c0       1.73323e-06 -1.17295e-07  1.21749e-07

sherpa> print get_fit_results()
datasets   = (1,)
methodname = neldermead
statname   = cstat
succeeded  = True
parnames   = ('abs1.nh', 'm1.gamma', 'm1.ampl', 'bkgA.c0')
parvals    = (0.21800727108863005, 2.1413252822977542, 0.00027295212454511546, 1.7332287723348897e-06)
covarerr   = None
statval    = 596.405796689
istatval   = 616.563237505
dstatval   = 20.1574408158
numpoints  = 636
dof        = 632
qval       = 0.841770402233
rstat      = 0.943680058053
message    = Optimization terminated successfully
nfev       = 4087

sherpa> print get_proj_results()
datasets   = (1,)
methodname = projection
fitname    = neldermead
statname   = cstat
sigma      = 1
percent    = 68.2689492137
parnames   = ('abs1.nh', 'm1.gamma', 'm1.ampl', 'bkgA.c0')
parvals    = (0.21800727108863005, 2.1413252822977542, 0.00027295212454511546, 1.7332287723348897e-06)
parmins    = (-0.013801448801628813, -0.048666029086867076, -1.1687130049730994e-05, -1.1729479581444341e-07)
parmaxes   = (0.013817744117379816, 0.048340119118419476, 1.2108523120806421e-05, 1.2174866808871212e-07)
nfits      = 64

Since the cstat fit statistic does not calculate errors for the data points, we group the data and change the fit statistic to chi2xspecvar to do so. Finally, we view the results of the new fit with the plot_fit_delchi command:

sherpa> group_counts(15) 
WARNING: grouping flags have changed, noticing all bins

sherpa> notice(0.5,7)  #re-establish energy filter

sherpa> set_stat("chi2xspecvar")

sherpa> plot_fit_delchi()

The new fit to the grouped simulated source-plus-background spectrum, along with the residuals divided by the uncertainties, is shown in Figure 6.

[Print media version: Plot of fit to simulated source-plus-background spectrum, with residuals]

Figure 6: Plot of fit to simulated source-plus-background spectrum, with residuals

The plot may now be saved as a PostScript file:

sherpa> print_window("simulation_fit_w_bkg")

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.

History

02 Feb 2009	Created for CIAO/Sherpa 4.1
29 Apr 2009	new `script` command is available with CIAO 4.1.2