Fitting a PHA Data Set with Multiple Responses
![[CXC Logo]](/ciao/imgs/cxc-logo.gif){kind=logo}
Sherpa Threads (CIAO 4.5 Sherpa v1)
Overview
Synopsis:
This thread demonstrates the use of the Sherpa functions set_full_model, set_bkg_full_model, and get_response to explicitly define complex source and background model expressions in which some of the model components are convolved with the instrument response, while others are not. It is not possible to define such a model expression in the usual way with set_source and set_bkg_source, as these functions do not allow for applying separate responses to individual model components within a single model expression.
We consider the scenario in which one would choose to model the X-ray background together with the source emission, instead of subtracting it. For example, this would be useful for modeling extended, diffuse source emission which covers the entire field of view of an observation, thereby leaving no source-free region from which to extract a background spectrum. As a result, the various components contributing to the background level would be modeled together with the source emission (assuming that using an average background level from a blank-sky data set is inappropriate).
If you would like to fit a background-subtracted source spectrum using a common RMF and ARF for the source and background, or fit source and background spectra simultaneously with proper and distinct RMFs and ARFs (but with one response per model expression, unlike in this thread), refer to the Sherpa thread Simultaneously Fitting Source and Background Spectra.
The sample data files used in this thread are available in sherpa.tar.gz; they can be generated by following the CIAO thread Using specextract to Extract ACIS Spectra and Response Files for Extended Sources .
Last Update: 13 Dec 2012 - updated for CIAO 4.5: group commands no longer clear the existing data filter
Contents
- Getting Started
- Reading the PHA Data into Sherpa
- Defining the Instrument Responses
- Defining the Model Expressions
- Modifying Method & Statistic Settings
- Fitting
- Plotting individual fitted model components
- Examining Fit Results
- Saving the Sherpa Session
- Scripting It
- Summary
- History
- Images
Getting Started
Please follow the "Obtaining data used in Sherpa threads" thread.
Reading the PHA Data into Sherpa
In this thread, we wish to fit the the PHA spectrum contained in the FITS file source.pi. This data set is read into Sherpa with the load_pha command, and assigned to default data set ID "1":
sherpa> load_pha("source.pi") read ARF file source.warf read RMF file source.wrmf read ARF (background) file bkg.warf read RMF (background) file bkg.wrmf read background file bkg.pi
The associated background PHA file and ARF and RMF instrument response files are automatically read into the session along with the source data by the load_pha function (watch out for the set_full_model bug when using this option, however). If the names of these files had not been recorded in the header of source.pi, they could have been loaded with load_bkg, load_arf, and load_rmf; however, watch out for the Sherpa 4.3 set_full_model
Information on all loaded data sets may be viewed with the show_all command, or individually with show_data, show_bkg, print(get_arf())/print(get_rmf()), and print(get_bkg_arf())/print(get_bkg_rmf()). The source, background, and source ARF data may be visualized with plot_data, plot_bkg and plot_arf (or simply the plot command, as shown below).
sherpa> set_xlog()
sherpa> set_ylog()
sherpa> ignore(":1.0, 7.0:")
sherpa> group_counts(1,30,bkg_id=1)
sherpa> group_counts(1,30)
sherpa> plot("data", 1, "bkg",1)
Before plotting, we set our preferences for the x- and y-axis scale of all data, model, and fit plots created in this session to logarithmic by using the set_xlog and set_ylog commands with no arguments (to learn how to change the default axis scale from linear to logarithmic so that these commands do not have to be run in each Sherpa session, see this Sherpa FAQ.) Then, we filter the data with the ignore command in order to include only the 1.0-7.0 keV range in our analysis, and group both spectra so that each bin contains a minumum of 30 counts to ensure that the spectral features are clearly visible in the plot displayed in Figure 1. Finally, the source and background data sets are plotted in separate frames within one window using the plot command.
![[Plot of the source and background spectrum]](1.png)
Figure 1: Plot of PHA source and background spectra
The grouped source spectrum (top) and associated background spectrum (bottom)
After plotting the spectra we ungroup the data sets, since we wish to fit the data ungrouped.
sherpa> ungroup() sherpa> ungroup(1,1)
Defining the Instrument Responses
Source Instrument Response
The ARF and RMF response files for source.pi were automatically loaded with load_pha, therefore we do not have to load them separately with load_arf and load_rmf. We assign the instrument response associated with the source data to the variable 'rsp' using the get_response function, which returns the instrument response model (ARF*RMF) associated with a data set.
sherpa> rsp = get_response()
The source data set ID is "1", which is the default, therefore it is not necessary to explicitly enter a data set ID as an argument to get_response. The current definition of the instrument response may be examined using show_all or the get_arf/get_bkg_arf and get_rmf/get_bkg_rmf commands:
sherpa> print(get_arf()) name = source.warf energ_lo = Float64[1070] energ_hi = Float64[1070] specresp = Float64[1070] bin_lo = None bin_hi = None exposure = 110163.239678 sherpa> print(get_rmf()) name = source.wrmf detchans = 1024 energ_lo = Float64[1070] energ_hi = Float64[1070] n_grp = UInt64[1070] f_chan = UInt32[1346] n_chan = UInt32[1346] matrix = Float64[382838] offset = 1 e_min = Float64[1024] e_max = Float64[1024] sherpa> print(get_bkg_arf()) name = bkg.warf energ_lo = Float64[1070] energ_hi = Float64[1070] specresp = Float64[1070] bin_lo = None bin_hi = None exposure = 110163.239678 sherpa> print(get_bkg_rmf()) name = bkg.wrmf detchans = 1024 energ_lo = Float64[1070] energ_hi = Float64[1070] n_grp = UInt64[1070] f_chan = UInt32[1364] n_chan = UInt32[1364] matrix = Float64[382698] offset = 1 e_min = Float64[1024] e_max = Float64[1024]
The output shows that source.warf, source.wrmf, bkg.warf, and bkg.wrmf currently define the source and background instrument responses.
Background Instrument Responses
In this example, we choose to model the background emission together with the source emission, which in this case involves defining a complex background model consisting of multiple components. We assume that the same ARF and RMF associated with the source data set apply to the background spectrum (i.e., there is not an independent response for the background since the background spectrum was extracted from a source-free free region of the source observation).
The background level of an X-ray observation consists of contributions from both astrophysical and instrumental sources, therefore some components of the background model need to be convolved with the source instrument response, whereas others do not. This requires that the instrumental components of the background model be unconvolved, i.e. convolved with a unit effective area, defined as follows:
sherpa> copy_data(1,2) sherpa> unit_arf = get_bkg_arf(2) sherpa> unit_arf.specresp = 0.*unit_arf.specresp + 1.0 sherpa> bunitrsp = get_response(2,bkg_id=1)
Here, we simply copy source data set 1 and its associated background and responses (id=1, bkg_id=1) to data set 2 (id=2, bkg_id=1) to avoid overwriting the background ARF data contained in data set 1, which we will need for the background model components which should be folded through the instrument effective area. Then, we assign the background ARF model in data set 2 (which is exactly the same as the background ARF in data set 1) returned by get_arf to the variable 'unit_arf' and set the spectral response of 'unit_arf' to contain all 1s. This is the unit ARF response model which is appropriate for convolving the instrumental background components of the model expression. The ARF*RMF instrument response model stored in variable 'bunitrsp' contains the unit background ARF and the untouched background RMF.
The background ARF*RMF instrument response model returned by 'get_response(bkg_id=1)' should be used to convolve the cosmic background model component, as shown in the next section. We define another response variable for the background, 'brsp', which includes the untouched ARF and untouched RMF for convolving the astrophysical background component:
sherpa> brsp = get_response(bkg_id=1)
Defining the Model Expressions
Background Model
The set_full_model and set_bkg_full_model functions are used to explicitly define complex source and background model expressions for simultaneous fitting. Since the background model in this example is complex, in that it contains many model components each convolved by a separate response, we consider the background model components first.
We assume that the background contributions to the emission spectrum consist of an absorbed thermal plasma model to represent the astrophysical diffuse Cosmic X-ray Background (CXB), and several Gaussian models plus an exponential cutoff power law model to fit the instrumental quiescent background. (Note that this example does not fully account for all possible background contributions to a Chandra ACIS data set - e.g., we are ignoring the contribution from resolved point sources to the CXB - it is meant only for demonstration purposes.)
We define the explicit background model expression using the set_bkg_full_model function, which is used to model the background associated with a source modeled by the set_full_model function. Note that it is not possible to define such a complex background model expression in the usual way with set_bkg_source, as this function does not allow for applying separate responses to individual model components within a single model expression.
sherpa> set_bkg_full_model(bunitrsp(gauss1d.bg1 + gauss1d.bg2 + xscutoffpl.bpl) + brsp(xsphabs.gal*xsapec.bth))
Multiplication by the background exposure time is implicit, and the necessary scaling resulting from differences between the source and background exposure times and spectral extraction areas is applied to the background model within the set_full_model expression, as shown in the next section (note in this example, the exposure ratio is 1 since the source and background spectra were extracted from the same observation).
Source Model
To describe the source contributions to the total emission in the spectrum, we define a source model consisting of an absorbed power law multiplied by an absorbed thermal plasma model, to describe the bright central point source and the surrounding diffuse thermal gas, respectively. We specify that the source model be convolved by the instrument response contained in the variable 'rsp', defined above. Multiplication by the source exposure time is done implicitly. We define the full model expression to include the source and background components, including the scaling of the background to account for the different spectral extraction areas of the source and background.
sherpa> my_bkg_model = (bunitrsp(bg1+bg2+bpl) + brsp(gal*bth)) sherpa> scale = get_bkg_scale() sherpa> set_full_model(rsp((xsphabs.gal+xsphabs.intrin)*xspowerlaw.spl+gal*xsapec.sth) + scale*my_bkg_model)
Note that the source component in the set_full_model expression above would normally defined by doing 'set_source((xsphabs.gal+xsphabs.intrin)*xspowerlaw.spl + gal*xsapec.sth)', which automatically and implicitly convolves the source model expression with the appropriate response. We do not define the source model in this usual way in this example, since the associated background spectrum requires a complex model expression consisting of multiple components convolved by different responses. Otherwise, set_bkg_source would be used to set the background source model in the usual way, for simultaneous source and background fitting.
(It is important to keep in mind 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, source and background models should be defined in the usual way with the automatic functions set_source and set_bkg_source.)
Modifying Method & Statistic Settings
The show_method and show_stat functions may be used to view the current method and statistics settings:
sherpa> show_method() 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 sherpa> show_stat() Statistic: Chi2Gehrels Chi Squared with Gehrels variance
We decide to use the Cstat fit statistic and Nelder-Mead optimization method for fitting this data, which we specify using the set_stat and set_method commands.
sherpa> set_method("neldermead")
sherpa> set_stat("cstat")
Details about the optimization methods and fit statistics available in Sherpa may be accessed by typing 'ahelp [method name]' at the Sherpa prompt, or by visiting the Sherpa Statistics and Optimization Methods pages.
sherpa> ahelp "levmar" sherpa> ahelp "chi2datavar"
Note that it is usually advisable to switch the optimization method and/or fit statistic at least once during a fitting trial to get a feel for the model parameter space and try to determine the 'true' best-fit model parameters. Observe an example of this in the Fitting section below, where we first fit using the chosen Nelder-Mead optimization method, and then switch to the Levenberg-Marquardt method to see how this affects the fit results. The Sherpa Optimization Methods page contains a detailed comparison of the available methods in Sherpa.
Fitting
The current definition of the Sherpa source model expression, including initial parameter values for each model component, may be examined using show_model:
sherpa> show_model() Model: 1 (apply_rmf(apply_arf((110163.239678 * (((xsphabs.gal + xsphabs.intrin) * xspowerlaw.spl) + (xsphabs.gal * xsapec.sth))))) + ((1.83660825173 * const2d.c0) * (apply_rmf(apply_arf((110163.239678 * ((gauss1d.bg1 + gauss1d.bg2) + xscutoffpl.bpl)))) + apply_rmf(apply_arf((110163.239678 * (xsphabs.gal * xsapec.bth))))))) Param Type Value Min Max Units ----- ---- ----- --- --- ----- gal.nH thawed 1 0 100000 10^22 atoms / cm^2 intrin.nH thawed 1 0 100000 10^22 atoms / cm^2 spl.PhoIndex thawed 1 -2 9 spl.norm thawed 1 0 1e+24 sth.kT thawed 1 0.008 64 keV sth.Abundanc frozen 1 0 5 sth.redshift frozen 0 0 10 sth.norm thawed 1 0 1e+24 c0.c0 thawed 1 0 3.40282e+38 bg1.fwhm thawed 10 1.17549e-38 3.40282e+38 bg1.pos thawed 0 -3.40282e+38 3.40282e+38 bg1.ampl thawed 1 -3.40282e+38 3.40282e+38 bg2.fwhm thawed 10 1.17549e-38 3.40282e+38 bg2.pos thawed 0 -3.40282e+38 3.40282e+38 bg2.ampl thawed 1 -3.40282e+38 3.40282e+38 bpl.PhoIndex thawed 1 -2 9 bpl.HighECut thawed 15 1 500 keV bpl.norm thawed 1 0 1e+24 bth.kT thawed 1 0.008 64 keV bth.Abundanc frozen 1 0 5 bth.redshift frozen 0 0 10 bth.norm thawed 1 0 1e+24
Since we are fitting a complex model expression to the data containing multiple components and responses, we approach the final fit to the data in a series of smaller trials in which we fit each model component separately, before fitting the entire model.
Background
The Sherpa guess function may be used to estimate initial parameter values for a model, as well as the minima and maxima for their ranges. To have Sherpa automatically query for the initial parameter values when a model is established, set 'paramprompt(True)' (it is 'False' by default).
We guess the initial parameter values for some background model components and set other values ourselves. For the thermal plasma model, we assume the solar abundances contained in Grevesse, N. & Sauval, A.J. (1998, Space Science Reviews 85, 161).
sherpa> set_xsabund("grsa")
sherpa> guess(bpl)
sherpa> guess(bg1)
sherpa> guess(bg2)
sherpa> guess(bth)
sherpa> set_par(gal.nH, 0.041, frozen=True)
sherpa> bpl.phoindex= 0.15
sherpa> bpl.HighECut = 5.6
sherpa> bg1.pos = 1.7
sherpq> bg2.pos= 2.1
sherpa> bth.kT = 0.15
First, we fit only the instrumental power law component using the fit_bkg function, which fits only the background data set(s) current in the session; this requires temporarily re-defining the background model expression to include only the model component we are interested in fitting:
sherpa> set_bkg_full_model(bunitrsp(bpl)) sherpa> fit_bkg(1) Dataset = 1 Method = neldermead Statistic = cstat Initial fit statistic = 912.151 Final fit statistic = 912.151 at function evaluation 341 Data points = 512 Degrees of freedom = 509 Probability [Q-value] = 2.68299e-25 Reduced statistic = 1.79204 Change in statistic = 5.05873e-06 bpl.PhoIndex 0.442125 bpl.HighECut 500 bpl.norm 0.0071846
Next, we successively add the two instrumental Gaussian model components to the background model expression, and finally the cosmic thermal plasma component, re-fitting at each step:
sherpa> set_bkg_full_model(bunitrsp(bpl+bg1))
sherpa> fit_bkg(1)
Dataset = 1
Method = neldermead
Statistic = cstat
Initial fit statistic = 789.642
Final fit statistic = 674.536 at function evaluation 3124
Data points = 512
Degrees of freedom = 506
Probability [Q-value] = 7.06269e-07
Reduced statistic = 1.33307
Change in statistic = 115.106
bpl.PhoIndex 0.683139
bpl.HighECut 499.988
bpl.norm 0.00791478
bg1.fwhm 2.0365
bg1.pos 7.9862
bg1.ampl 0.00375058
sherpa> set_bkg_full_model(bunitrsp(bpl+bg1+bg2))
sherpa> fit_bkg(1)
Dataset = 1
Method = neldermead
Statistic = cstat
Initial fit statistic = 39628.8
Final fit statistic = 696.579 at function evaluation 7041
Data points = 512
Degrees of freedom = 503
Probability [Q-value] = 2.12696e-08
Reduced statistic = 1.38485
Change in statistic = 38932.2
bpl.PhoIndex 2.32043
bpl.HighECut 37.9477
bpl.norm 0.00353292
bg1.fwhm 15.3817
bg1.pos 7.9862
bg1.ampl 0.00349173
bg2.fwhm 0.919855
bg2.pos 1.96145
bg2.ampl 0.00385232
sherpa> set_bkg_full_model(bunitrsp(bpl+bg1+bg2) + brsp(gal*bth))
sherpa> fit_bkg(1)
Dataset = 1
Method = neldermead
Statistic = cstat
Initial fit statistic = 820679
Final fit statistic = 712.799 at function evaluation 6200
Data points = 512
Degrees of freedom = 500
Probability [Q-value] = 1.12331e-09
Reduced statistic = 1.4256
Change in statistic = 819966
bpl.PhoIndex 1.38392
bpl.HighECut 42.595
bpl.norm 0.00496462
bg1.fwhm 14.0695
bg1.pos 7.9862
bg1.ampl 0.00342836
bg2.fwhm 0.0151716
bg2.pos 2.1474
bg2.ampl 0.0547313
gal.nH 3.08825e-06
bth.kT 0.068824
bth.norm 0.000181551
sherpa> bpl.PhoIndex = 0.1
sherpa> bpl.HighECut = 5.6
sherpa> set_method("levmar")
sherpa> fit_bkg(1)
Dataset = 1
Method = levmar
Statistic = cstat
Initial fit statistic = 803.149
Final fit statistic = 487.134 at function evaluation 841
Data points = 410
Degrees of freedom = 399
Probability [Q-value] = 0.00164593
Reduced statistic = 1.22089
Change in statistic = 316.015
bpl.PhoIndex 0.250295
bpl.HighECut 46.4581
bpl.norm 0.004966
bg1.fwhm 0.0646018
bg1.pos 2.14831
bg1.ampl 0.0127971
bg2.fwhm 0.00263481
bg2.pos 1.75962
bg2.ampl 0.246075
bth.kT 0.272472
bth.norm 9.19381e-06
We can check that the fit was successful by examining the fit residuals:
sherpa> plot_bkg_fit_resid()
sherpa> set_current_plot("plot2")
sherpa> linear_scale(Y_AXIS)
Figure 2 shows the resulting plot.
![[bitmap image of fit and residuals]](2.png)
Source plus background
Now that we have reasonable background model parameter values resulting from the background fitting trials, we can fit the entire source-plus-background model expression, first with the background model parameters frozen to constrain the source model, and finally with them thawed.
We fit the source and background model components together in one fitting trial after setting initial model parameter values for the source:
sherpa> show_model(1) Model: 1 (apply_rmf(apply_arf((110163.239678 * (((xsphabs.gal + xsphabs.intrin) * xspowerlaw.spl) + (xsphabs.gal * xsapec.sth))))) + (1.83660825173 * (apply_rmf(apply_arf((110163.239678 * ((gauss1d.bg1 + gauss1d.bg2) + xscutoffpl.bpl)))) + apply_rmf(apply_arf((110163.239678 * (xsphabs.gal * xsapec.bth))))))) Param Type Value Min Max Units ----- ---- ----- --- --- ----- gal.nH frozen 0.041 0 100000 10^22 atoms / cm^2 intrin.nH thawed 1 0 100000 10^22 atoms / cm^2 spl.PhoIndex thawed 1 -2 9 spl.norm thawed 1 0 1e+24 sth.kT thawed 1 0.008 64 keV sth.Abundanc frozen 1 0 5 sth.redshift frozen 0 0 10 sth.norm thawed 1 0 1e+24 bg1.fwhm thawed 0.0646018 0.0022776 2277.6 bg1.pos thawed 2.14831 1.0074 6.9934 bg1.ampl thawed 0.0127971 6.37784e-06 6.37784 bg2.fwhm thawed 0.00263481 0.0022776 2277.6 bg2.pos thawed 1.75962 1.0074 6.9934 bg2.ampl thawed 0.246075 6.37784e-06 6.37784 bpl.PhoIndex thawed 0.250295 -2 9 bpl.HighECut thawed 46.4581 1 500 keV bpl.norm thawed 0.004966 6.37784e-06 6.37784 bth.kT thawed 0.272472 0.008 64 keV bth.Abundanc frozen 1 0 5 bth.redshift frozen 0 0 10 bth.norm thawed 9.19381e-06 6.37784e-06 6.37784 sherpa> freeze(bpl, bg1, bg2, bth.norm, bth.kT) sherpa> guess(spl) sherpa> guess(sth) sherpa> intrin.nH = 0.038 sherpa> spl.phoindex= 2.0 sherpa> set_par(sth.redshift, 2.0, frozen=True) sherpa> sth.kT=1.5
The source and background data assigned to data set 1 may be simultaneously fit using the fit function, as follows:
sherpa> fit(1)
Dataset = 1
Method = levmar
Statistic = cstat
Initial fit statistic = 1.14228e+06
Final fit statistic = 928.394 at function evaluation 114
Data points = 820
Degrees of freedom = 815
Probability [Q-value] = 0.00342189
Reduced statistic = 1.13913
Change in statistic = 1.14135e+06
intrin.nH 6.18161e-05
spl.PhoIndex 0.97105
spl.norm 6.37784e-06
sth.kT 26.1037
sth.norm 0.00106604
sherpa> thaw(bpl, bg1, bg2, bth.norm, bth.kT)
sherpa> fit(1)
Dataset = 1
Method = levmar
Statistic = cstat
Initial fit statistic = 928.394
Final fit statistic = 894.235 at function evaluation 4101
Data points = 820
Degrees of freedom = 804
Probability [Q-value] = 0.0143468
Reduced statistic = 1.11223
Change in statistic = 34.1583
intrin.nH 0.701925
spl.PhoIndex 1.76154
spl.norm 2.90006e-05
sth.kT 17.1906
sth.norm 0.000550728
bg1.fwhm 0.0291892
bg1.pos 2.14697
bg1.ampl 0.0232272
bg2.fwhm 0.00304577
bg2.pos 1.75931
bg2.ampl 0.286318
bpl.PhoIndex 0.237348
bpl.HighECut 269.221
bpl.norm 0.00464436
bth.kT 0.287499
bth.norm 6.37784e-06
sherpa> set_method("neldermead")
sherpa> fit(1)
Dataset = 1
Method = neldermead
Statistic = cstat
Initial fit statistic = 894.235
Final fit statistic = 893.885 at function evaluation 6211
Data points = 820
Degrees of freedom = 804
Probability [Q-value] = 0.014646
Reduced statistic = 1.1118
Change in statistic = 0.350213
intrin.nH 0.688133
spl.PhoIndex 1.75548
spl.norm 2.89329e-05
sth.kT 17.2815
sth.norm 0.00054637
bg1.fwhm 0.0315643
bg1.pos 2.14706
bg1.ampl 0.0215552
bg2.fwhm 0.00320581
bg2.pos 1.75926
bg2.ampl 0.272416
bpl.PhoIndex 0.242769
bpl.HighECut 499.643
bpl.norm 0.00464142
bth.kT 0.28771
bth.norm 6.37789e-06
The plot_fit_resid and plot_fit_delchi functions may be used to visualize the quality of the fit, and the Chips function print_window to save the plot as a PostScript file:
sherpa> plot_fit_resid()
sherpa> print_window("sherpa_fit.ps")
Figure 3 shows the resulting plot.
![[bitmap image of fit and residuals]](3.png)
Plotting individual fitted model components
To visualize the contributions to the fit of the emission from the individual source and background components, we can use the Sherpa plot_model_component command together with the ChIPS set_histogram command for adjusting the color of the line describing each of the model components. Note that due to a bug in Sherpa 4.3, the set_xlog and set_ylog commands issued earlier in the session will not produce model component plots in log scale, as they should, therefore until the bug is resolved it is necessary to use the get_model_component_plot command to do so, as shown below.
sherpa> get_model_component_plot(gal*spl).histo_prefs['xlog'] = True
sherpa> get_model_component_plot(gal*spl).histo_prefs['ylog'] = True
sherpa> plot_model_component(gal*spl)
sherpa> set_histogram("line.color=yellow")
sherpa> get_model_component_plot(bg1+bg2).histo_prefs['xlog'] = True
sherpa> get_model_component_plot(bg1+bg2).histo_prefs['ylog'] = True
sherpa> plot_model_component(bg1+bg2, overplot=True)
sherpa> set_histogram("line.color=cyan")
sherpa> get_model_component_plot(bpl).histo_prefs['xlog'] = True
sherpa> get_model_component_plot(bpl).histo_prefs['ylog'] = True
sherpa> plot_model_component(bpl, overplot=True)
sherpa> set_histogram("line.color=blue")
sherpa> get_model_component_plot(gal*bth).histo_prefs['xlog'] = True
sherpa> get_model_component_plot(gal*bth).histo_prefs['ylog'] = True
sherpa> plot_model_component(gal*bth, overplot=True)
sherpa> set_histogram("line.color=darkred")
sherpa> set_plot_title("Source and background model components")
The resulting plot is shown in Figure 4
![[bitmap image of confidence interval]](4.png)
Figure 4: Individual components of best-fit model to source and background spectra
Plot of the best-fit model to the combined source and background spectra in red, with the individual model components overlaid in different colors (yellow=source absorbed power law, blue=background power law, cyan=background Gaussians, darkred=background absorbed thermal plasma).
As demonstrated above, the plot_model_component command may be used to plot one, or a combination of, individual source model components to allow the user to quickly visualize the contribution to the full model being used to fit the data. The model components plotted with this command are convolved with any assigned convolution models, e.g. PSF or PHA responses. The plot_source_component command is available for visualizing the unconvolved model components, and get_model_component_plot/get_source_component_plot for accessing the data arrays and preferences which define the corresponding plot.
Examining Fit Results
Several algorithms are available in Sherpa for examining fit results, such as confidence, covariance, interval-projection, and region-projection. Here, we use the preferred confidence method to estimate the 1-sigma confidence intervals for the thawed model parameters.
sherpa> set_conf_opt("fast","True")
sherpa> conf(1)
...
{verbose output omitted for brevity}
...
Dataset = 1
Confidence Method = confidence
Iterative Fit Method = None
Fitting Method = levmar
Statistic = cstat
confidence 1-sigma (68.2689%) bounds:
Param Best-Fit Lower Bound Upper Bound
----- -------- ----------- -----------
intrin.nH 0.701925 -0.299223 0.551248
spl.PhoIndex 1.76154 -0.210291 0.278333
spl.norm 2.90006e-05 -6.27555e-06 3.39516e-05
sth.kT 17.1906 ----- -----
sth.norm 0.000550728 -0.000307154 0.000194673
bg1.fwhm 0.0291892 -0.0278343 0.0622292
bg1.pos 2.14697 -0.0124926 0.0214192
bg1.ampl 0.0232272 -0.0171875 0.169657
bg2.fwhm 0.00304577 ----- 0.00967904
bg2.pos 1.75931 -0.00936602 0.0099398
bg2.ampl 0.286318 -0.195799 0.144778
bpl.PhoIndex 0.237348 -0.0581616 0.0664979
bpl.HighECut 269.221 -220.068 -----
bpl.norm 0.00464436 -0.000427341 0.000653372
bth.kT 0.287499 ----- 0.0710831
bth.norm 6.37784e-06 ----- 0.000503317
Note that the lines of dashes "----" in the confidence results indicate that the hard upper or lower limit was reached for one or more parameters. This occurs when the parameter bound found by the confidence method lies outside the hard limit boundary for a model parameter - this could result from an issue with the signal-to-noise of the data, the applicability of the model to the data, systematic errors in the data, among other things. The interval-projection and region-projection routines are useful for exploring these issues further.
The quality of a fit is also available in the chi-square goodness-of-fit information reported after each fit, as well as in the output of the functions get_fit_results and show_fit .
sherpa> show_fit() 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: CStat Maximum likelihood function (XSPEC style) Fit:Dataset = 1 Method = levmar Statistic = cstat Initial fit statistic = 928.394 Final fit statistic = 894.235 at function evaluation 4101 Data points = 820 Degrees of freedom = 804 Probability [Q-value] = 0.0143468 Reduced statistic = 1.11223 Change in statistic = 34.1583 intrin.nH 0.701925 spl.PhoIndex 1.76154 spl.norm 2.90006e-05 sth.kT 17.1906 sth.norm 0.000550728 bg1.fwhm 0.0291892 bg1.pos 2.14697 bg1.ampl 0.0232272 bg2.fwhm 0.00304577 bg2.pos 1.75931 bg2.ampl 0.286318 bpl.PhoIndex 0.237348 bpl.HighECut 269.221 bpl.norm 0.00464436 bth.kT 0.287499 bth.norm 6.37784e-06
Saving the Sherpa Session
Before exiting Sherpa, you may wish to save the session in order to return to the analysis at a later point:
sherpa> save("session1.save") sherpa> save_all("session1.ascii")
The save function records all the information about the current session to the binary file session1.save, and the save_all function records the session settings to an editable ASCII file.
To restore the session that was saved to the binary file session1.save or ASCII file session1.ascii:
sherpa> restore("session1.save") sherpa> execfile("session1.ascii")
One may verify that the session has been properly restored by examining the information returned by the show_all() command.
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.5 syntax to you.
Summary
This thread is complete, so we can exit the Sherpa session:
sherpa> quit
History
| 05 Jul 2010 | original version |
| 15 Dec 2010 | updated for Sherpa in CIAO 4.3: new functions plot_model_component, get_bkg_scale, and set_xlog/set_ylog. |
| 29 Jun 2011 | title of a referenced thread was changed from "Independent Background Responses" to "Simultaneously Fitting Source and Background S pectra" |
| 06 Jan 2012 | reviewed for CIAO 4.4: reordered the ignore and group_counts commands in the section "Reading the PHA Data into Sherpa" to work around a grouping/filtering bug |
| 13 Dec 2012 | updated for CIAO 4.5: group commands no longer clear the existing data filter |
