Simultaneously Fitting Source and Background Spectra

Sherpa Threads (CIAO 4.4 Sherpa v1)

Overview

Synopsis:

If you would like to fit a background-subtracted source spectrum using a common RMF and ARF for source and background, simply read the source spectrum fits file into Sherpa, subtract the background, and fit it. To fit source and background spectra simultaneously with proper and distinct RMFs and ARFs, load the source and background as different data sets. This thread illustrates both cases.

The sample data files used in this thread are available in sherpa.tar.gz; they can be generated by following the CIAO thread Extract ACIS Spectra and Response Files for Pointlike Sources.

Last Update: 15 Dec 2011 - reviewed for CIAO 4.4: added a warning about filtering/grouping source and background

Getting Started

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

Statistical Issues: Background Subtraction

A typical dataset may contain multiple spectra, one of which contains contributions from the source and the background, and one or more which contain background counts alone. (The background itself may contain contributions from the cosmic X-ray background, the particle background, an so on, but we ignore this complication.)

The proper way to treat background data is to model them. However, many X-ray astronomers subtract background data from the raw data.

Why should one not subtract background?

It reduces the amount of statistical information in the analysis - the final fit parameter values will be a less accurate estimate of the true values.
The background subtracted data are not Poisson-distributed; one cannot fit them with the Poisson likelihood or the Cash statistic, even in the low-count limit. For example, subtracting a background can give negative counts; this is definitely not Poissonian!
Fluctuations, particularly in the vicinity of localized features, can adversely affect analysis.

Fit a Background-Subtracted Source

Reading FITS data for source and background

By default, psextract updates the header keywords BACKFILE, RESPFILE, and ANCRFILE in the source spectrum file:

unix% dmkeypar 3c273.pi BACKFILE echo+
3c273_bg.pi

unix% dmkeypar 3c273.pi RESPFILE echo+
3c273.rmf

unix% dmkeypar 3c273.pi ANCRFILE echo+
3c273.arf

The RESPFILE and ANCRFILE keywords are not updated in the ungrouped background file, however, and should be manually added with dmhedit.

On account of these keywords, Sherpa automatically reads in the (ungrouped) background spectrum "3c273_bg.pi" and the source response matrices "3c273.rmf" and "3c273.arf" when the source spectrum is read:

sherpa> load_pha("3c273.pi")
statistical errors were found in file '3c273.pi' 
but not used; to use them, re-read with use_errors=True
read ARF file 3c273.arf
read RMF file 3c273.rmf
statistical errors were found in file '3c273_bg.pi' 
but not used; to use them, re-read with use_errors=True
read background file 3c273_bg.pi

If Sherpa does not automatically read in the background data, then it can be input as follows:

sherpa> load_bkg("3c273_bg.pi")
WARNING: systematic errors were not found in file '3c273_bg.pi'          
statistical errors were found in file '3c273_bg.pi' 
but not used; to use them, re-read with use_errors=True

Subtracting the background

The user can then subtract the background from the source spectrum and fit the background-subtracted source spectrum using response matrices for the source only. The following commands may be used as an example script for subtracting the background and then visualizing the 0.3-3 keV background-subtracted source data in Sherpa:

sherpa> subtract()

sherpa> set_xlog
sherpa> set_ylog
sherpa> plot_data()

sherpa> ignore()
sherpa> notice(0.3, 3.)

sherpa> plot_data()

sherpa> limits(Y_AXIS, 1.e-3, 0.2)
sherpa> limits(X_AXIS, 0.2, 3.5)

Which produces Figure 1.

Figure 1: The background-subtracted source spectrum

[Plot of the background-subtracted source spectrum]

[Print media version: Plot of the background-subtracted source spectrum]

Figure 1: The background-subtracted source spectrum

Plot of background-subtracted spectrum of 3c273

The get_bkg_scale function - not to be confused with get_backscal - returns the value of the coefficient which is used to scale background counts during background subtraction of a source spectrum, or a simultaneous fit of source and background spectra. (The get_backscal function returns the value associated with the OGIP PHA header keyword BACKSCAL, which is only one component of the complete scaling factor.) The complete scaling factor used to scale the background counts in these cases consists of the product of the source-to-background exposure and backscal (extraction region area) ratios:

 Background-subtracted spectrum (cts/s) = SOURCE_COUNTS/SOURCE_EXPOSURE - BGD_COUNTS/BGD_EXPOSURE/BGD_BACKSCAL*SOURCE_BACKSCAL.

      Background scale factor = (SOURCE_EXPOSURE/BACKGROUND_EXPOSURE)*(SOURCE_BACKSCAL/BACKGROUND_BACKSCAL)

The background scaling factor is also returned by the show_data and show_all functions. The value returned by the get_bkg_scale funcion can be checked by doing the following in Sherpa:

sherpa> get_bkg_scale()
        0.13492064388809599


sherpa> manual_bkg_scale_check = get_exposure()*get_backscal()/get_exposure(bkg_id=1)/get_backscal(bkg_id=1)

sherpa> print(manual_bkg_scale_check) 
0.134920643888

Defining the source instrument response

Since the header of the input source spectrum file referenced the instrument response files, an instrument response was automatically defined:

sherpa> print(get_arf())

name     = 3c273.arf
energ_lo = Float64[1090]
energ_hi = Float64[1090]
specresp = Float64[1090]
bin_lo   = None
bin_hi   = None
exposure = 38570.6924812

sherpa> print(get_rmf())

name     = 3c273.rmf
detchans = 1024
energ_lo = Float64[1090]
energ_hi = Float64[1090]
n_grp    = Int16[1090]
f_chan   = UInt32[2004]
n_chan   = UInt32[2004]
matrix   = Float64[58564]
offset   = 1
e_min    = Float64[1024]
e_max    = Float64[1024]

sherpa> print(get_bkg_arf())

name     = 3c273.arf
energ_lo = Float64[1090]
energ_hi = Float64[1090]
specresp = Float64[1090]
bin_lo   = None
bin_hi   = None
exposure = 38570.6924812

sherpa> print(get_bkg_rmf())

name     = 3c273.rmf
detchans = 1024
energ_lo = Float64[1090]
energ_hi = Float64[1090]
n_grp    = Int16[1090]
f_chan   = UInt32[2004]
n_chan   = UInt32[2004]
matrix   = Float64[58564]
offset   = 1
e_min    = Float64[1024]
e_max    = Float64[1024]

Note that the instrument response was automatically set up for both the source and background. However, unless a background model is defined for fitting, using the set_bkg_model command, then the background data will be ignored. If you manually read in a background file after reading the source data, the background instrument response will be reset. To set it to be identical to the source instrument response, use load_bkg_arf and load_bkg_rmf with the source response files:

sherpa> load_bkg_rmf("3c273.rmf")
sherpa> load_bkg_arf("3c273.arf")

sherpa> print(get_bkg_arf())
name     = 3c273.arf
energ_lo = Float64[1090]
energ_hi = Float64[1090]
specresp = Float64[1090]
bin_lo   = None
bin_hi   = None
exposure = 38570.6924812

sherpa> print(get_bkg_rmf())
name     = 3c273.rmf
detchans = 1024
energ_lo = Float64[1090]
energ_hi = Float64[1090]
n_grp    = Int16[1090]
f_chan   = UInt32[2004]
n_chan   = UInt32[2004]
matrix   = Float64[58564]
offset   = 1
e_min    = Float64[1024]
e_max    = Float64[1024]

Defining and fitting a source model

The following commands may be used as an example script for fitting the background-subtracted data with a source model consisting of an absorbed power law:

sherpa> set_source(xswabs.abs1*powlaw1d.p1)
sherpa> show_source()
Model: 1
(xswabs.abs1 * powlaw1d.p1)
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   abs1.nh      thawed            1            0       100000 10^22 atoms / cm^2
   p1.gamma     thawed            1          -10           10           
   p1.ref       frozen            1 -3.40282e+38  3.40282e+38           
   p1.ampl      thawed            1            0  3.40282e+38  

sherpa> abs1.nh = 0.1
sherpa> guess(p1) 

sherpa> fit()
 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
Data Set               = 1
Method                = levmar
Statistic             = chi2gehrels
Initial fit statistic = 89.8323
Final fit statistic   = 15.478 at function evaluation 94
Data points           = 33
Degrees of freedom    = 30
Probability [Q-value] = 0.986745
Reduced statistic     = 0.515934
Change in statistic   = 74.3543
   abs1.nh        0           
   p1.gamma       1.68306     
   p1.ampl        0.000160777 

sherpa> plot_fit_resid()


sherpa> limits(Y_AXIS, 1.e-3, 0.1)
sherpa> limits(X_AXIS, 0.2, 3.5)

The results are shown Figure 2.

Figure 2: Fit of the background-subtracted data

[Plot of the fit to the background-subtracted source spectrum]

[Print media version: Plot of the fit to the background-subtracted source spectrum]

Figure 2: Fit of the background-subtracted data

Plot of the source model fit to the background-subtracted spectrum of 3c273

The confidence function may be used to estimate the errors on the individual parameters of the fit; in this example, we choose to calculate the 1-sigma errors:

sherpa> print(get_conf_opt())
{'verbose': False, 'openinterval': False, 'max_rstat': 3, 'maxiters': 200, 'soft_limits': False, 'eps': 0.01, 'fast': False, 'maxfits': 5, 'remin': 0.01, 'tol': 0.20000000000000001, 'sigma': 3, 'parallel': True}

sherpa> set_conf_opt("sigma", 1)

sherpa> conf()
abs1.nH lower bound:    -----
abs1.nH upper bound:    0.0117646
p1.gamma lower bound:   -0.0849583
p1.gamma upper bound:   0.147068
p1.ampl lower bound:    -8.16938e-06
p1.ampl upper bound:    1.49111e-05
Dataset               = 1
Confidence Method     = confidence
Fitting Method        = levmar
Statistic             = chi2gehrels
confidence 1-sigma (68.2689%) bounds:
   Param            Best-Fit  Lower Bound  Upper Bound
   -----            --------  -----------  -----------
   abs1.nH                 0        -----    0.0117646
   p1.gamma          1.68306   -0.0849583     0.147068
   p1.ampl       0.000160777 -8.16938e-06  1.49111e-05

Simultaneously Fit Source and Background with the Same Responses

Instead of subtracting the background, the user may choose to simultaneously fit the source and background spectra, using the same (source) RMF and ARF.

Reading FITS data for source and background

Again, Sherpa automatically reads in the (ungrouped) background spectrum and the source response matrices when the source spectrum is read:

sherpa> clean()
sherpa> load_pha("3c273.pi")
statistical errors were found in file '3c273.pi' 
but not used; to use them, re-read with use_errors=True
read ARF file 3c273.arf
read RMF file 3c273.rmf
statistical errors were found in file '3c273_bg.pi' 
but not used; to use them, re-read with use_errors=True
read background file 3c273_bg.pi

If Sherpa does not automatically read in the background data, then it can be input as follows:

sherpa> load_bkg("3c273_bg.pi")
WARNING: systematic errors were not found in file '3c273_bg.pi'
statistical errors were found in file '3c273_bg.pi' 
but not used; to use them, re-read with use_errors=True

Warning about filtering and grouping source and background data

Note that if the source and background data sets are to be grouped (not considered in this thread), then due to a bug in CIAO 4.4, you must ensure that the background data is grouped before the source, as shown in the example below; otherwise, any data filter in place for the background will be removed after the grouping command is issued.

sherpa> notice(0.5,7.)             # source and background filtered

sherpa> group_counts(20, bkg_id=1) # background data must be grouped before
                                   # source data to retain background filter                            

sherpa> group_counts(20)           # both source and background filtered and gro

Defining the source and background instrument response

Since the header of the input source spectrum file referenced the instrument response files, an instrument response was automatically defined:

sherpa> print(get_arf())

name     = 3c273.arf
energ_lo = Float64[1090]
energ_hi = Float64[1090]
specresp = Float64[1090]
bin_lo   = None
bin_hi   = None
exposure = 38570.6924812

sherpa> print(get_rmf())

name     = 3c273.rmf
detchans = 1024
energ_lo = Float64[1090]
energ_hi = Float64[1090]
n_grp    = Int16[1090]
f_chan   = UInt32[2004]
n_chan   = UInt32[2004]
matrix   = Float64[58564]
offset   = 1
e_min    = Float64[1024]
e_max    = Float64[1024]

sherpa> print(get_bkg_arf())

name     = 3c273.arf
energ_lo = Float64[1090]
energ_hi = Float64[1090]
specresp = Float64[1090]
bin_lo   = None
bin_hi   = None
exposure = 38570.6924812

sherpa> print(get_bkg_rmf())

name     = 3c273.rmf
detchans = 1024
energ_lo = Float64[1090]
energ_hi = Float64[1090]
n_grp    = Int16[1090]
f_chan   = UInt32[2004]
n_chan   = UInt32[2004]
matrix   = Float64[58564]
offset   = 1
e_min    = Float64[1024]
e_max    = Float64[1024]

The instrument response was automatically set up for both the source and background. However, unless a background model is defined for fitting, using the set_bkg_model() command, then the background data will be ignored. In order to fit the background simultaneously with the source, the set_bkg_model() command will be issued in the next step.

If you manually read in a background file after reading the source data, the background instrument response will be reset. To set it to be identical to the source instrument response, use load_bkg_arf() and load_bkg_rmf() with the source response files:

sherpa> load_bkg_rmf("3c273.rmf")
sherpa> load_bkg_arf("3c273.arf")

sherpa> print(get_bkg_arf());

name     = 3c273.arf
energ_lo = Float64[1090]
energ_hi = Float64[1090]
specresp = Float64[1090]
bin_lo   = None
bin_hi   = None
exposure = 38570.6924812

sherpa> print(get_bkg_rmf())

name     = 3c273.rmf
detchans = 1024
energ_lo = Float64[1090]
energ_hi = Float64[1090]
n_grp    = Int16[1090]
f_chan   = UInt32[2004]
n_chan   = UInt32[2004]
matrix   = Float64[58564]
offset   = 1
e_min    = Float64[1024]
e_max    = Float64[1024]

Defining and fitting source and background models

As of Sherpa in CIAO 4.2, the fit_bkg function is available (and the equivalent 'bkg_only' option of the fit command). This function fits the defined background model(s) to the background data set(s) by ID. It may be called with no arguments, in which case a fit is done simultaneously on all background data sets for which the user has defined a model to be fit; source data sets ready for fitting will be ignored in this case.

Here, we choose to simultaneously fit the source and background with the fit function, using an absorbed power law model for the source and an unabsorbed power law model for the background. The source and background models are defined using the set_source and set_bkg_model() commands, respectively.

sherpa> set_source(xswabs.abs1*powlaw1d.srcp1)
sherpa> set_bkg_model(powlaw1d.bkgp1)
sherpa> abs1.nh = 0.1
sherpa> guess(srcp1)
sherpa> guess(bkgp1)

sherpa> show_model()
Model: 1
apply_rmf(apply_arf((38564.6089269 * ((xswabs.abs1 * powlaw1d.srcp1) + scale_factor * (powlaw1d.bkgp1)))))
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   abs1.nH      thawed          0.1            0       100000 10^22 atoms / cm^2
   srcp1.gamma  thawed            1          -10           10           
   srcp1.ref    frozen            1 -3.40282e+38  3.40282e+38           
   srcp1.ampl   thawed  0.000144301  1.44301e-07     0.144301           
   bkgp1.gamma  thawed            1          -10           10           
   bkgp1.ref    frozen            1 -3.40282e+38  3.40282e+38           
   bkgp1.ampl   thawed  0.000144301  1.44301e-07     0.144301      

sherpa> print(get_bkg_model())
apply_rmf(apply_arf((38564.6089269 * powlaw1d.bkgp1)))
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   bkgp1.gamma  thawed            1          -10           10           
   bkgp1.ref    frozen            1 -3.40282e+38  3.40282e+38           
   bkgp1.ampl   thawed  0.000144301  1.44301e-07     0.144301          

sherpa> fit()
Dataset               = 1
Method                = levmar
Statistic             = chi2gehrels
Initial fit statistic = 4604.01
Final fit statistic   = 137.257 at function evaluation 86
Data points           = 92
Degrees of freedom    = 87
Probability [Q-value] = 0.000478885
Reduced statistic     = 1.57766
Change in statistic   = 4466.75
   abs1.nH        0.0178333   
   srcp1.gamma    1.92493     
   srcp1.ampl     0.000174607 
   bkgp1.gamma    -0.167137   
   bkgp1.ampl     2.46655e-06 

sherpa> plot_fit_resid()


sherpa> limits(Y_AXIS, 1.e-3, 0.2)
sherpa> limits(X_AXIS, 0.2, 3.5)

The results are shown Figure 3.

Figure 3: Simultaneous fit of the source and background data (same responses)

[Plot of the simultaneous fit of the 3c273 source and background spectra]

Figure 3: Simultaneous fit of the source and background data (same responses)

Plot of the simultaneous fit of the 3c273 source and background spectra.

The calc_stat_info command (see also get_stat_info) - similar to the CIAO 3.4 GOODNESS command - may be used to access the goodness-of-fit statistics for the individua source and background data set included in the simultaneous fit, plus the simultaneous fit, without having to re-run the fit.

sherpa> calc_stat_info()
Dataset               = 1
Statistic             = chi2gehrels
Fit statistic value   = 36.8354
Data points           = 46
Degrees of freedom    = 41
Probability [Q-value] = 0.656205
Reduced statistic     = 0.898423

Background 1 in Dataset = 1
Statistic             = chi2gehrels
Fit statistic value   = 100.421
Data points           = 46
Degrees of freedom    = 43
Probability [Q-value] = 1.71508e-06
Reduced statistic     = 2.33538

Dataset               = 1
Statistic             = chi2gehrels
Fit statistic value   = 137.257
Data points           = 92
Degrees of freedom    = 87
Probability [Q-value] = 0.000478885
Reduced statistic     = 1.57766

The confidence 1-sigma error estimates on the individual parameters of the fit are:

sherpa> conf()
bkgp1.gamma lower bound:        -1.03473
bkgp1.gamma upper bound:        0.653029
srcp1.ampl lower bound: -1.24582e-05
abs1.nH lower bound:    -0.00528628
srcp1.gamma lower bound:        -0.0750981
bkgp1.ampl lower bound: -2.12015e-06
srcp1.ampl upper bound: 1.28537e-05
srcp1.gamma upper bound:        0.0763482
bkgp1.ampl upper bound: 4.32302e-06
abs1.nH upper bound:    0.00682136
Dataset               = 1
Confidence Method     = confidence
Fitting Method        = levmar
Statistic             = chi2gehrels
confidence 1-sigma (68.2689%) bounds:
   Param            Best-Fit  Lower Bound  Upper Bound
   -----            --------  -----------  -----------
   abs1.nH          0.017875  -0.00528628   0.00682136
   srcp1.gamma       1.92517   -0.0750981    0.0763482
   srcp1.ampl    0.000174634 -1.24582e-05  1.28537e-05
   bkgp1.gamma     -0.144094     -1.03473     0.653029
   bkgp1.ampl    2.55807e-06 -2.12015e-06  4.32302e-06

Simultaneously Fit Source and Background with Independent Responses

Instead of subtracting the background, the user may choose to simultaneously fit the source and background spectra, each with its own RMF and ARF.

Reading FITS data for source and background

In this thread, we wish to fit spectral data from the FITS (ungrouped) data files source.pi and back.pi. These data sets are input into Sherpa with load commands:

sherpa> clean()

sherpa> load_pha("source.pi")
statistical errors were found in file 'source.pi' 
but not used; to use them, re-read with use_errors=True

sherpa> load_bkg("back.pi")
statistical errors were found in file 'source.pi' 
but not used; to use them, re-read with use_errors=True

The source and background data sets can be plotted in separate windows with plot_data() and plot_bkg(), or together in one window with the following syntax:

sherpa> plot("data",1,"bkg",1)

The plots are displayed in Figure 4. The data are plotted in bin (i.e. channel) space, since we have not yet defined the instrument response (RMFs and ARFs).

Figure 4: Plot of the ungrouped source and background data

[Plot of the ungrouped source and background spectra]

[Print media version: Plot of the ungrouped source and background spectra]

Figure 4: Plot of the ungrouped source and background data

Defining source and background instrument responses

Here we establish 1-D instrument responses by loading the source and background response files with load commands:

sherpa> load_arf("source.arf")
sherpa> load_rmf("source.rmf")

sherpa> load_bkg_arf("back.arf")
sherpa> load_bkg_rmf("back.rmf")

The current definition of the instrument response may be examined using show_all() or get commands:

sherpa> print(get_arf())

name     = source.arf
energ_lo = Float64[1090]
energ_hi = Float64[1090]
specresp = Float64[1090]
bin_lo   = None
bin_hi   = None
exposure = 49470.9477594

sherpa> print(get_rmf())

name     = source.rmf
detchans = 1024
energ_lo = Float64[1090]
energ_hi = Float64[1090]
n_grp    = Int16[1090]
f_chan   = UInt32[1090]
n_chan   = UInt32[1090]
matrix   = Float64[572598]
offset   = 1
e_min    = Float64[1024]
e_max    = Float64[1024]

sherpa> print(get_bkg_arf())

name     = back.arf
energ_lo = Float64[1090]
energ_hi = Float64[1090]
specresp = Float64[1090]
bin_lo   = None
bin_hi   = None
exposure = 49470.9477594

sherpa> print(get_bkg_rmf())

name     = back.rmf
detchans = 1024
energ_lo = Float64[1090]
energ_hi = Float64[1090]
n_grp    = Int16[1090]
f_chan   = UInt32[1090]
n_chan   = UInt32[1090]
matrix   = Float64[552166]
offset   = 1
e_min    = Float64[1024]
e_max    = Float64[1024]

The output shows that source.arf and source.rmf currently define the source instrument response.

Defining and fitting source and background models

We now need to define a source model for our source data set and a background model for our background data set. For these data, we find that absorbed blackbody spectra are appropriate. We define a model which includes an absorbing column and a blackbody for both the source and background:

sherpa> set_source(xswabs.a1*xsbbody.b1)
sherpa> set_bkg_model(a1*xsbbody.b2)

Note that we also ignore the high and low energy regions of the spectrum, as they generally have lower quality data. The plots of the data with (default value) fits produced by plot commands are shown in Figure 5 and Figure 6.

sherpa> ignore()
sherpa> notice(0.3, 10.)
sherpa> plot_fit()
sherpa> plot_bkg_fit()

Figure 5: Plot of the ungrouped source data with default fit

[Plot of the ungrouped source spectrum with default fit]

[Print media version: Plot of the ungrouped source spectrum with default fit]

Figure 5: Plot of the ungrouped source data with default fit

Figure 6: Plot of the ungrouped background data with default fit

[Plot of the ungrouped background spectrum with default fit]

[Print media version: Plot of the ungrouped background spectrum with default fit]

Figure 6: Plot of the ungrouped background data with default fit

Clearly, the default values are not a very good fit to the data. We may choose to begin the fit now, using the fit() command, or we can refine the values of the fit by eye, prior to fitting. By varying parameter values, while plotting and replotting the data, we can help the fitting algorithm find the best minimum. By setting the parameters to values near what we expect, we also help avoid local minima in the parameter space. (Fitting spectra is something of an art; one generally gets better fits when they have a priori knowledge of the source.) We set the values of the parameters as follows:

sherpa> a1.nH=0.0336676
sherpa> b1.kT=20
sherpa> b1.norm=1e-02
sherpa> b2.kT=0.5
sherpa> b2.norm.min=1e-6
sherpa> b2.norm=1e-05

Note that the normalization of b2, "b2.norm", appears to have a better fit at values less than the default minimum. We set the b2.norm.min to a new, smaller value to give us the optimal fit. Now we have a reasonable start to the fit:

sherpa> plot_fit()
sherpa> plot_bkg_fit()

The plot is visible Figure 7.

Figure 7: Plot of the ungrouped source data with custom fits

[Plots of the ungrouped source and data with custom fits]

[Print media version: Plots of the ungrouped source and data with custom fits]

Figure 7: Plot of the ungrouped source data with custom fits

Figure 8: Plot of the ungrouped background data with custom fits

[Print media version: Plot of the ungrouped background data with custom fits]

Figure 8: Plot of the ungrouped background data with custom fits

In general, we may also want to change our fit statistic or optimization method with set_stat and set_method; here, we choose to run the fit with the current fit statisic and method:

sherpa> fit()
Data Set               = 1
Method                = levmar
Statistic             = chi2gehrels
Initial fit statistic = 1024.04
Final fit statistic   = 945.974 at function evaluation 393
Data points           = 1330
Degrees of freedom    = 1325
Probability [Q-value] = 1
Reduced statistic     = 0.713942
Change in statistic   = 78.0655
   a1.nh          0.0336454   
   b1.kt          100         
   b1.norm        1.07253     
   b2.kt          0.56523     
   b2.norm        1.16969e-05

Once the fit is complete, Sherpa will display the fit values to the screen. You may also display the overall status of the Sherpa session with the show_all() and show_bkg() command:

sherpa> show_all()
Data Set: 1
Filter: 0.2993-9.9937 Energy (keV)
Noticed Channels: 21.0-685.0

name           = source.pi
channel        = Float64[1024]
counts         = Float64[1024]
staterror      = None
syserror       = None
bin_lo         = None
bin_hi         = None
grouping       = None
quality        = None
exposure       = 49429.2334679
backscal       = 1.87253514146e-05
areascal       = 1.0
grouped        = False
subtracted     = False
units          = energy
response_ids   = [1]
background_ids = [1]

RMF Data Set: 1:1

name     = source.rmf
detchans = 1024
energ_lo = Float64[1090]
energ_hi = Float64[1090]
n_grp    = Int16[1090]
f_chan   = UInt32[1090]
n_chan   = UInt32[1090]
matrix   = Float64[572598]
offset   = 1
e_min    = Float64[1024]
e_max    = Float64[1024]

ARF Data Set: 1:1

name     = source.arf
energ_lo = Float64[1090]
energ_hi = Float64[1090]
specresp = Float64[1090]
bin_lo   = None
bin_hi   = None
exposure = 49470.9477594

sherpa> show_bkg()
Background Data Set: 1:1

name           = back.pi
channel        = Float64[1024]
counts         = Float64[1024]
staterror      = None
syserror       = None
bin_lo         = None
bin_hi         = None
grouping       = None
quality        = None
exposure       = 49429.2334679
backscal       = 2.69645060371e-05
areascal       = 1.0
grouped        = False
subtracted     = False
units          = energy
response_ids   = [1]
background_ids = []

Background RMF Data Set: 1:1

name     = back.rmf
detchans = 1024
energ_lo = Float64[1090]
energ_hi = Float64[1090]
n_grp    = Int16[1090]
f_chan   = UInt32[1090]
n_chan   = UInt32[1090]
matrix   = Float64[552166]
offset   = 1
e_min    = Float64[1024]
e_max    = Float64[1024]

Background ARF Data Set: 1:1

name     = back.arf
energ_lo = Float64[1090]
energ_hi = Float64[1090]
specresp = Float64[1090]
bin_lo   = None
bin_hi   = None
exposure = 49470.9477594

Model: 1
apply_rmf(apply_arf((49429.2334679 * ((xswabs.a1 * xsbbody.b1) + scale_factor * ((xswabs.a1 * xsbbody.b2))))))
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   a1.nh        thawed    0.0336454            0       100000 10^22 atoms / cm^2
   b1.kt        thawed          100         0.01          100        keV
   b1.norm      thawed      1.07253            0        1e+24 L39 / (D10)**2
   b2.kt        thawed      0.56523         0.01          100        keV
   b2.norm      thawed  1.16969e-05        1e-06        1e+24 L39 / (D10)**2

Optimization Method: LevMar
name    = levmar
ftol    = 1.19209289551e-07
xtol    = 1.19209289551e-07
gtol    = 1.19209289551e-07
maxfev  = None
epsfcn  = 1.19209289551e-07
factor  = 100.0
verbose = 0

Statistic: Chi2Gehrels

Fit:Data Set               = 1
Method                = levmar
Statistic             = chi2gehrels
Initial fit statistic = 1024.04
Final fit statistic   = 945.974 at function evaluation 393
Data points           = 1330
Degrees of freedom    = 1325
Probability [Q-value] = 1
Reduced statistic     = 0.713942
Change in statistic   = 78.0655
   a1.nh          0.0336454   
   b1.kt          100         
   b1.norm        1.07253     
   b2.kt          0.56523     
   b2.norm        1.16969e-05

One may visually examine the fit with the Sherpa plot functions:

sherpa> plot_fit()
sherpa> plot_bkg_fit()

The plots are visible in Figure 9 and Figure 10.

Figure 9: Simultaneous fit of the ungrouped source data (different responses)

[Print media version: Simultaneous fit of the ungrouped source data (different responses)]

Figure 9: Simultaneous fit of the ungrouped source data (different responses)

Figure 10: Simultaneous fit of the ungrouped background data (different responses)

[Print media version: Simultaneous fit of the ungrouped background data (different responses)]

Figure 10: Simultaneous fit of the ungrouped background data (different responses)

We may also calculate the confidence 1-sigma error estimates on the individual parameters of the fit:

sherpa> conf()
b2.kT lower bound:      -0.00992703
b2.kT upper bound:      0.00992703
b2.norm lower bound:    -2.24597e-07
b1.kT lower bound:      -74.2716
b1.kT upper bound:      -----
b1.norm lower bound:    -1.05292
a1.nH lower bound:      -0.0107257
b2.norm upper bound:    2.26351e-07
b1.norm upper bound:    0.0300121
a1.nH upper bound:      0.0110382
Dataset               = 1
Confidence Method     = confidence
Fitting Method        = levmar
Statistic             = chi2gehrels
confidence 1-sigma (68.2689%) bounds:
   Param            Best-Fit  Lower Bound  Upper Bound
   -----            --------  -----------  -----------
   a1.nH           0.0336454   -0.0107257    0.0110382
   b1.kT                 100     -74.2716        -----
   b1.norm           1.07253     -1.05292    0.0300121
   b2.kT             0.56523  -0.00992703   0.00992703
   b2.norm       1.16969e-05 -2.24597e-07  2.26351e-07

Manually defining source and background models

As of the release of CIAO 4.2 Sherpa version 2, we may choose to explicitly set the complete convolved model expression to be used for fitting a source data set and its associated background, using the new functions set_full_model and set_bkg_full_model, together with get_response or get_arf/get_rmf. These functions offer an alternative to using the functions set_source and set_bkg_model, which automatically convolve the appropriate instrument response with a defined source and background model expression.

The following set of commands represents the manual definition of the chosen source and background models, which were automatically (and equivalently) defined above with set_source and set_bkg_model.

sherpa> rsp = get_response()
sherpa> bkg_rsp = get_response(bkg_id=1)
sherpa> bkg_scale = get_bkg_scale()

sherpa> set_full_model(rsp(xswabs.a1*xsbbody.b1) + bkg_scale*bkg_rsp(a1*xsbbody.b2))
sherpa> set_bkg_full_model(bkg_rsp(a1*xsbbody.b2))

You could also use get_arf and get_rmf in place of get_response, as follows:

sherpa> arf = get_arf()
sherpa> rmf = get_rmf()
sherpa> bgarf = get_arf(bkg_id=1)
sherpa> bgrmf = get_rmf(bkg_id=1)

sherpa> bkg_scale = get_bkg_scale()

sherpa> set_full_model(rmf(arf(xswabs.a1*xsbbody.b1)) + bkg_scale*bgrmf(bgarf(a1*xsbbody.b2)))
sherpa> set_bkg_full_model(bgrmf(bgarf(a1*xsbbody.b2)))

Note that the Sherpa functions which are related to a source or background model defined with set_source or set_bkg_source, such as plot_source/plot_bkg_source or calc_energy_flux, are not compatible with the complete model expression defined by the set_full_model or set_bkg_full_model functions. In order to use these Sherpa functions, users should define source and background models in the usual way with the automatic functions set_source and set_bkg_source.

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)

(Note that restoring a Sherpa session from such a file could be problematic since it may include syntax errors, unwanted fitting trials, etcetera.)

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

Summary

This thread is complete, so we can exit the Sherpa session:

sherpa> quit

History

14 Jan 2005	updated for CIAO 3.2: minor change in fit results
21 Dec 2005	reviewed for CIAO 3.3: no changes
01 Dec 2006	reviewed for CIAO 3.4: no changes
08 Dec 2008	updated for Sherpa 4.1
29 Apr 2009	new `script` command is available with CIAO 4.1.2
17 Dec 2009	updated for CIAO 4.2: the `fit_bkg` function is available
14 Jun 2010	updated to include `confidence` 1-sigma error estimates on individual parameters of fits
14 Jun 2010	updated with an example using the CIAO 4.2 Sherpa v2 functions `set_full_model`, `set_bkg_full_model`, `get_response`, and a new usage of the functions `get_arf` and `get_rmf`. S-Lang version of thread removed.
15 Dec 2010	updated for CIAO 4.3: new functions `get_bkg_scale`, `calc_stat_info`, and `set_xlog`/`set_ylog` are available
22 Jun 2011	renamed thread from "Independent Background Responses" to "Simultaneously Fitting Source and Background Spectra"
15 Dec 2011	reviewed for CIAO 4.4: added a warning about filtering/grouping source and background