Simultaneously Fitting Source and Background Spectra
Sherpa Threads (CIAO 4.12 Sherpa v1)
Overview
Synopsis:
If you would like to fit a backgroundsubtracted 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: 11 Dec 2019  Updated to use Matplotlib in CIAO 4.12 and to take advantage of changes to the plot support (e.g. plot_fit(ylog=True)); the text and commands have been revised slightly.
Contents
 Getting Started
 Statistical Issues: Background Subtraction
 Fit a BackgroundSubtracted Source
 Simultaneously Fit Source and Background with the Same Responses
 Simultaneously Fit Source and Background with Independent Responses
 Scripting It
 Summary
 History

Images
 Figure 1: The backgroundsubtracted source spectrum: all the data
 Figure 2: The backgroundsubtracted source spectrum: 0.3  3 keV
 Figure 3: Fit of the backgroundsubtracted data
 Figure 4: Simultaneous fit of the source and background data (same responses)
 Figure 5: Plot of the ungrouped source and background data
 Figure 6: Adding responses
 Figure 7: Plot of the ungrouped source and background data with default fits
 Figure 8: Plot of the ungrouped source data with custom fits
 Figure 9: Simultaneous fit of the ungrouped source data (different responses)
 Figure 10: Simultaneous fit of the ungrouped background data (different responses)
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 Xray background, the particle background, and so on, but we ignore this complication.)
The proper way to treat background data is to model them. However, many Xray 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 backgroundsubtracted data are not Poissondistributed; one cannot fit them with the Poisson likelihood or the Cash statistic, even in the lowcounts 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 BackgroundSubtracted Source
Reading FITS data for source and background
By default, specextract 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
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")
WARNING: systematic errors were not found in file '3c273.pi'
statistical errors were found in file '3c273.pi'
but not used; to use them, reread with use_errors=True
read ARF file 3c273.arf
read RMF file 3c273.rmf
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, reread 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, reread with use_errors=True
Subtracting the background
The user can then subtract the background from the source spectrum and fit the backgroundsubtracted 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.33 keV backgroundsubtracted source data in Sherpa:
sherpa> subtract() sherpa> set_xlog() sherpa> set_ylog() sherpa> plot_data()
Figure 1: The backgroundsubtracted source spectrum: all the data
sherpa> notice(0.3, 3.)
sherpa> plot_data(xlog=False, ylog=False)
Which produces Figure 2.
Figure 2: The backgroundsubtracted source spectrum: 0.3  3 keV
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 complete scaling factor used to scale the background counts in these cases consists of the product of the sourcetobackground exposure and backscale (extraction region area).
The backgroundsubtracted spectrum is then:
Backgroundsubtracted spectrum (cts/s) = SOURCE_COUNTS/SOURCE_EXPOSURE  BGD_COUNTS/BGD_EXPOSURE/BGD_BACKSCAL*SOURCE_BACKSCAL.and the background scaling factor is:
Background scale factor = (SOURCE_EXPOSURE/BACKGROUND_EXPOSURE)*(SOURCE_BACKSCAL/BACKGROUND_BACKSCAL)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 background scaling factor is also returned by the show_data and show_all functions. The value returned by the get_bkg_scale function 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 ethresh = 1e10 sherpa> print(get_rmf()) name = 3c273.rmf detchans = 1024 energ_lo = Float64[1090] energ_hi = Float64[1090] n_grp = UInt64[1090] f_chan = UInt32[2004] n_chan = UInt32[2004] matrix = Float64[58564] offset = 1 e_min = Float64[1024] e_max = Float64[1024] ethresh = 1e10 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 ethresh = 1e10 sherpa> print(get_bkg_rmf()) name = 3c273.rmf detchans = 1024 energ_lo = Float64[1090] energ_hi = Float64[1090] n_grp = UInt64[1090] f_chan = UInt32[2004] n_chan = UInt32[2004] matrix = Float64[58564] offset = 1 e_min = Float64[1024] e_max = Float64[1024] ethresh = 1e10
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 ethresh = 1e10 sherpa> print(get_bkg_rmf()) name = 3c273.rmf detchans = 1024 energ_lo = Float64[1090] energ_hi = Float64[1090] n_grp = UInt64[1090] f_chan = UInt32[2004] n_chan = UInt32[2004] matrix = Float64[58564] offset = 1 e_min = Float64[1024] e_max = Float64[1024] ethresh = 1e10
Defining and fitting a source model
The following commands may be used as an example script for fitting the backgroundsubtracted data with a source model consisting of an absorbed powerlaw:
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> print(p1) powlaw1d.p1 Param Type Value Min Max Units       p1.gamma thawed 1 10 10 p1.ref frozen 1 3.40282e+38 3.40282e+38 p1.ampl thawed 0.000144301 1.44301e07 0.144301 sherpa> fit() Dataset = 1 Method = levmar Statistic = chi2gehrels Initial fit statistic = 89.8323 Final fit statistic = 15.4612 at function evaluation 184 Data points = 33 Degrees of freedom = 30 Probability [Qvalue] = 0.98686 Reduced statistic = 0.515375 Change in statistic = 74.3711 abs1.nH 0.000216916 +/ 0.0149809 p1.gamma 1.70584 +/ 0.153887 p1.ampl 0.000161164 +/ 1.53492e05 sherpa> plot_fit_resid() sherpa> plt.yscale('linear')
The results are shown Figure 3.
Figure 3: Fit of the backgroundsubtracted data
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σ errors:
sherpa> get_conf_opt() {'sigma': 1, 'eps': 0.01, 'maxiters': 200, 'soft_limits': False, 'remin': 0.01, 'fast': False, 'parallel': True, 'numcores': 4, 'maxfits': 5, 'max_rstat': 3, 'tol': 0.2, 'verbose': False, 'openinterval': False} sherpa> set_conf_opt("sigma", 1) sherpa> conf() abs1.nH lower bound:  abs1.nH upper bound: 0.0113339 p1.gamma lower bound: 0.106925 p1.gamma upper bound: 0.123277 p1.ampl lower bound: 8.44533e06 p1.ampl upper bound: 1.43411e05 Dataset = 1 Confidence Method = confidence Iterative Fit Method = None Fitting Method = levmar Statistic = chi2gehrels confidence 1sigma (68.2689%) bounds: Param BestFit Lower Bound Upper Bound     abs1.nH 0.000216916  0.0113339 p1.gamma 1.70584 0.106925 0.123277 p1.ampl 0.000161164 8.44533e06 1.43411e05
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") WARNING: systematic errors were not found in file '3c273.pi' statistical errors were found in file '3c273.pi' but not used; to use them, reread with use_errors=True read ARF file 3c273.arf read RMF file 3c273.rmf 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, reread with use_errors=True read background file 3c273_bg.pi
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 = UInt64[1090] f_chan = UInt32[2004] n_chan = UInt32[2004] matrix = Float64[58564] offset = 1 e_min = Float64[1024] e_max = Float64[1024] ethresh = 1e10 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 ethresh = 1e10 sherpa> print(get_bkg_rmf()) name = 3c273.rmf detchans = 1024 energ_lo = Float64[1090] energ_hi = Float64[1090] n_grp = UInt64[1090] f_chan = UInt32[2004] n_chan = UInt32[2004] matrix = Float64[58564] offset = 1 e_min = Float64[1024] e_max = Float64[1024] ethresh = 1e10
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.
As with the Defining the source instrument response section, the load_bkg_arf and load_bkg_rmf routines can be used to set the responses.
Defining and fitting source and background models
The fit_bkg function is available (and the equivalent 'bkg_only' option of the fit command) allows you to fit models to just the background components. It 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 powerlaw model for the source and an unabsorbed powerlaw 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(abs1 * powlaw1d.bkgp1) sherpa> abs1.nh = 0.1 sherpa> guess(srcp1) sherpa> guess(bkgp1) sherpa> show_model() Model: 1 apply_rmf(apply_arf((38564.608926889 * ((xswabs.abs1 * powlaw1d.srcp1) + 0.134921 * ((xswabs.abs1 * 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.44301e07 0.144301 bkgp1.gamma thawed 1 10 10 bkgp1.ref frozen 1 3.40282e+38 3.40282e+38 bkgp1.ampl thawed 0.000144301 1.44301e07 0.144301 sherpa> show_bkg_model() Background Model: 1:1 apply_rmf(apply_arf((38564.608926889 * (xswabs.abs1 * powlaw1d.bkgp1)))) Param Type Value Min Max Units       abs1.nH thawed 0.1 0 100000 10^22 atoms / cm^2 bkgp1.gamma thawed 1 10 10 bkgp1.ref frozen 1 3.40282e+38 3.40282e+38 bkgp1.ampl thawed 0.000144301 1.44301e07 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 [Qvalue] = 0.000478885 Reduced statistic = 1.57766 Change in statistic = 4466.75 abs1.nH 0.0178333 +/ 0.00594283 srcp1.gamma 1.92493 +/ 0.0776745 srcp1.ampl 0.000174607 +/ 1.27987e05 bkgp1.gamma 0.167143 +/ 0.517068 bkgp1.ampl 2.46653e06 +/ 2.08853e06 sherpa> plot_fit_resid(xlog=True, ylog=True) sherpa> plt.yscale('linear') sherpa> plt.xlim(0.2, 3.5) sherpa> plt.sca(plt.gcf().axes[0]) sherpa> plt.ylim(7e4, 2e2)
The results are shown Figure 4.
Figure 4: Simultaneous fit of the source and background data (same responses)
The calc_stat_info command (see also get_stat_info) may be used to access the goodnessoffit statistics for the individual source and background data set included in the simultaneous fit, plus the simultaneous fit, without having to rerun the fit.
sherpa> calc_stat_info() Dataset = 1 Statistic = chi2gehrels Fit statistic value = 36.8354 Data points = 46 Degrees of freedom = 41 Probability [Qvalue] = 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 [Qvalue] = 1.71508e06 Reduced statistic = 2.33538 Dataset = 1 Statistic = chi2gehrels Fit statistic value = 137.257 Data points = 92 Degrees of freedom = 87 Probability [Qvalue] = 0.000478885 Reduced statistic = 1.57766
The confidence 1σ error estimates on the individual parameters of the fit are:
sherpa> conf() bkgp1.ampl lower bound: 2.0436e06 abs1.nH lower bound: 0.0052732 bkgp1.gamma lower bound: 1.03539 srcp1.ampl lower bound: 1.24576e05 bkgp1.gamma upper bound: 0.659018 bkgp1.ampl upper bound: 4.28518e06 abs1.nH upper bound: 0.00672868 srcp1.gamma lower bound: 0.0750823 srcp1.ampl upper bound: 1.28531e05 srcp1.gamma upper bound: 0.0763325 Dataset = 1 Confidence Method = confidence Iterative Fit Method = None Fitting Method = levmar Statistic = chi2gehrels confidence 1sigma (68.2689%) bounds: Param BestFit Lower Bound Upper Bound     abs1.nH 0.0178333 0.0052732 0.00672868 srcp1.gamma 1.92493 0.0750823 0.0763325 srcp1.ampl 0.000174607 1.24576e05 1.28531e05 bkgp1.gamma 0.167143 1.03539 0.659018 bkgp1.ampl 2.46653e06 2.0436e06 4.28518e06
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, reread with use_errors=True sherpa> load_bkg("back.pi") statistical errors were found in file 'source.pi' but not used; to use them, reread 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", "bkg")
The plots are displayed in Figure 5. The data are plotted in bin (i.e. channel) space, since we have not yet defined the instrument response (RMFs and ARFs).
Figure 5: Plot of the ungrouped source and background data
Defining source and background instrument responses
Here we establish 1D 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")
Now when we plot the data the X axis is in energy, rather than channel (Figure 6).
Figure 6: Adding responses
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 ethresh = 1e10 sherpa> print(get_rmf()) name = source.rmf detchans = 1024 energ_lo = Float64[1090] energ_hi = Float64[1090] n_grp = UInt64[1090] f_chan = UInt32[1090] n_chan = UInt32[1090] matrix = Float64[572598] offset = 1 e_min = Float64[1024] e_max = Float64[1024] ethresh = 1e10 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 ethresh = 1e10 sherpa> print(get_bkg_rmf()) name = back.rmf detchans = 1024 energ_lo = Float64[1090] energ_hi = Float64[1090] n_grp = UInt64[1090] f_chan = UInt32[1090] n_chan = UInt32[1090] matrix = Float64[552166] offset = 1 e_min = Float64[1024] e_max = Float64[1024] ethresh = 1e10
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 and model (with the default parameter values) produced by plot commands are shown in Figure 7.
sherpa> notice(0.3, 10.) sherpa> plot_fit(xlog=True, ylog=True, color='green') sherpa> plot_bkg_fit(overplot=True, color='orange')
Figure 7: Plot of the ungrouped source and background data with default fits
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 fitbyeye, 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 = 1e02 sherpa> b2.kT = 0.5 sherpa> b2.norm.min = 1e6 sherpa> b2.norm = 1e05
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> plt.subplot(2, 1, 1) sherpa> plot_fit(xlog=True, ylog=True, clearwindow=False) sherpa> plt.subplot(2, 1, 2) sherpa> plot_bkg_fit(xlog=True, ylog=True, clearwindow=False) sherpa> plt.subplots_adjust(hspace=0.5)
The plot is visible Figure 8.
Figure 8: Plot of the ungrouped source 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> Dataset = 1 Method = levmar Statistic = chi2gehrels Initial fit statistic = 828.712 Final fit statistic = 731.835 at function evaluation 196 Data points = 1056 Degrees of freedom = 1051 Probability [Qvalue] = 1 Reduced statistic = 0.696322 Change in statistic = 96.877 a1.nH 0.0352536 +/ 0.0122578 b1.kT 4.23927 +/ 0.636896 b1.norm 0.000146661 +/ 5.15993e05 b2.kT 0.556463 +/ 0.0103887 b2.norm 1.12408e05 +/ 2.44213e07
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.29937.9935 Energy (keV) Bkg Scale: 0.694444 Noticed Channels: 21548 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.233467924 backscal = 1.872535141462e05 areascal = 1.0 grouped = False subtracted = False units = energy rate = True plot_fac = 0 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 = UInt64[1090] f_chan = UInt32[1090] n_chan = UInt32[1090] matrix = Float64[572598] offset = 1 e_min = Float64[1024] e_max = Float64[1024] ethresh = 1e10 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.94775939 ethresh = 1e10 Background Data Set: 1:1 Filter: 0.29937.9935 Energy (keV) Noticed Channels: 21548 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.233467924 backscal = 2.6964506037052e05 areascal = 1.0 grouped = False subtracted = False units = energy rate = True plot_fac = 0 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 = UInt64[1090] f_chan = UInt32[1090] n_chan = UInt32[1090] matrix = Float64[552166] offset = 1 e_min = Float64[1024] e_max = Float64[1024] ethresh = 1e10 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.94775939 ethresh = 1e10 Model: 1 apply_rmf(apply_arf((49429.233467924 * ((xswabs.a1 * xsbbody.b1) + 0.694444 * ((xswabs.a1 * xsbbody.b2)))))) Param Type Value Min Max Units       a1.nH thawed 0.0352536 0 100000 10^22 atoms / cm^2 b1.kT thawed 4.23927 0.01 100 keV b1.norm thawed 0.000146661 0 1e+24 L39 / (D10)**2 b2.kT thawed 0.556463 0.01 100 keV b2.norm thawed 1.12408e05 1e06 1e+24 L39 / (D10)**2 Optimization Method: LevMar name = levmar ftol = 1.1920928955078125e07 xtol = 1.1920928955078125e07 gtol = 1.1920928955078125e07 maxfev = None epsfcn = 1.1920928955078125e07 factor = 100.0 numcores = 1 verbose = 0 Statistic: Chi2Gehrels Chi Squared with Gehrels variance. The variance is estimated from the number of counts in each bin, but unlike `Chi2DataVar`, the Gaussian approximation is not used. This makes it moresuitable for use with lowcount data. The standard deviation for each bin is calculated using the approximation from [1]_: sigma(i,S) = 1 + sqrt(N(i,s) + 0.75) where the higherorder terms have been dropped. This is accurate to approximately one percent. For data where the background has not been subtracted then the error term is: sigma(i) = sigma(i,S) whereas with background subtraction, sigma(i)^2 = sigma(i,S)^2 + [A(S)/A(B)]^2 sigma(i,B)^2 See Also  Chi2DataVar, Chi2ModVar, Chi2XspecVar Notes  The accuracy of the error term when the background has been subtracted has not been determined. A preferable approach to background subtraction is to model the background as well as the source signal. References  .. [1] "Confidence limits for small numbers of events in astrophysical data", Gehrels, N. 1986, ApJ, vol 303, p. 336346. http://adsabs.harvard.edu/abs/1986ApJ...303..336G Fit:Dataset = 1 Method = levmar Statistic = chi2gehrels Initial fit statistic = 828.712 Final fit statistic = 731.835 at function evaluation 196 Data points = 1056 Degrees of freedom = 1051 Probability [Qvalue] = 1 Reduced statistic = 0.696322 Change in statistic = 96.877 a1.nH 0.0352536 +/ 0.0122578 b1.kT 4.23927 +/ 0.636896 b1.norm 0.000146661 +/ 5.15993e05 b2.kT 0.556463 +/ 0.0103887 b2.norm 1.12408e05 +/ 2.44213e07
One may visually examine the fit with the Sherpa plot functions:
sherpa> set_xlog() sherpa> set_ylog() sherpa> plot_fit_delchi() sherpa> plt.yscale('linear')
Figure 9: Simultaneous fit of the ungrouped source data (different responses)
sherpa> plot_bkg_fit_resid()
sherpa> plt.yscale('linear')
Figure 10: Simultaneous fit of the ungrouped background data (different responses)
We may also calculate the confidence 1σ error estimates on the individual parameters of the fit:
sherpa> conf() b2.norm lower bound: 2.55525e07 b1.norm lower bound: 4.66583e05 b2.kT lower bound: 0.0107563 b2.kT upper bound: 0.0106946 a1.nH lower bound: 0.0108897 b1.norm upper bound: 9.71841e05 b2.norm upper bound: 2.57521e07 a1.nH upper bound: 0.0114366 b1.kT lower bound: 0.656143 b1.kT upper bound: 1.00158 Dataset = 1 Confidence Method = confidence Iterative Fit Method = None Fitting Method = levmar Statistic = chi2gehrels confidence 1sigma (68.2689%) bounds: Param BestFit Lower Bound Upper Bound     a1.nH 0.0352536 0.0108897 0.0114366 b1.kT 4.23927 0.656143 1.00158 b1.norm 0.000146661 4.66583e05 9.71841e05 b2.kT 0.556463 0.0107563 0.0106946 b2.norm 1.12408e05 2.55525e07 2.57521e07
Manually defining source and background models
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> src_model = rsp(xswabs.a1 * xsbbody.b1) sherpa> bkg_model = bkg_rsp(a1 * xsbbody.b2) sherpa> set_full_model(src_model + bkg_scale * bkg_model) sherpa> set_bkg_full_model(bkg_model)
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> src_model = rmf(arf(xswabs.a1 * xsbbody.b1)) sherpa> bkg_model = bgrmf(bgarf(a1 * xsbbody.b2)) sherpa> set_full_model(src_model + bkg_scale * bkg_model) sherpa> set_bkg_full_model(bkg_model)
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 %run i 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, et cetera.)
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 1sigma 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. SLang 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 
13 Dec 2012  reviewed for CIAO 4.5: removed an outdated warning about filtering/grouping both source and background data, as the associated bug was fixed 
04 Dec 2013  reviewed for CIAO 4.6: removed references to deprecated tools, updated fit results 
03 Dec 2015  updated for CIAO 4.8: no content change 
09 Nov 2016  reviewed for CIAO 4.9: updated fits, no content change. 
25 May 2018  reviewed for CIAO 4.10: no content change; typeset equations with LaTeX. 
10 Dec 2018  reviewed for CIAO 4.11, screen output revised 
11 Dec 2019  Updated to use Matplotlib in CIAO 4.12 and to take advantage of changes to the plot support (e.g. plot_fit(ylog=True)); the text and commands have been revised slightly. 