Gallery: Fitting Data

Examples

1D PHA data:

2D Image data:

Estimating fit parameter bounds (for data of any dimensionality):
- Parameter bounds with interval-projection: Confidence plot of fit statistic vs. parameter value
- Parameter bounds with region-projection: Confidence contour plot of fit statistic vs. two parameter values

1D PHA data: Power-law model fit to source counts spectrum, with residuals

Here we display an absorbed power-law fit to a PHA source counts spectrum, along with the fit delta chi residuals (data-minus-model counts residuals normalized by their errors).

Sherpa 4.13 (Python)

load_pha("3c273.pi")

load_arf("3c273.arf")
load_rmf("3c273.rmf")
load_bkg("3c273_bg.pi")

notice_id(1, 0.1, 6.0)
subtract()

set_source(xsphabs.abs1 * powlaw1d.p1)
abs1.nH = 0.07
freeze(abs1.nH)
guess(p1)

fit()

plot_fit_delchi()

Sherpa 3.4

data 3c273.pi

instrument source 1 = RSP[mydetector]("3c273.rmf", "3c273.arf")
instrument back 1 = mydetector


notice energy 0.1:6.0
subtract

paramprompt off
source = xsphabs[abs1]*powlaw1d[p1]
abs1.nH = 0.07
freeze abs1
guess p1 

fit

lplot 2 fit delchi

For detailed information about each of the steps in the script, see the Sherpa thread "Introduction to Fitting PHA Spectra", which also contains links to the ahelp files for each Sherpa function used.

1D PHA data: Simultaneous power-law fit to ACIS-S/HETG source grating spectra

Here we display a simultaneous fit of four source-plus-background ACIS-S/HETG grating spectra, with the source modeled by an absorbed broken power-law, and the background modeled by an absorbed power-law.

Sherpa 4.13 (Python)

load_pha(1, "459_heg_m1_bin10.pha")
load_pha(2, "459_heg_p1_bin10.pha")
load_pha(3, "459_meg_m1_bin10.pha")
load_pha(4, "459_meg_p1_bin10.pha")

load_arf(1, "459_heg_m1.arf")
load_arf(2, "459_heg_p1.arf")
load_arf(3, "459_meg_m1.arf")
load_arf(4, "459_meg_p1.arf")

load_arf(1, "459_heg_m1.arf", bkg_id=1)
load_arf(1, "459_heg_m1.arf", bkg_id=2)
load_arf(2, "459_heg_p1.arf", bkg_id=1)
load_arf(2, "459_heg_p1.arf", bkg_id=2)
load_arf(3, "459_meg_m1.arf", bkg_id=1)
load_arf(3, "459_meg_m1.arf", bkg_id=2)
load_arf(4, "459_meg_p1.arf", bkg_id=1)
load_arf(4, "459_meg_p1.arf", bkg_id=2)

set_analysis("wave")

ignore()
notice(1., 15.)

set_analysis("wave")
hc = 12.39841874  #in [keV-Angstrom] 
dummy = atten.dummy
dummy.integrate = False
def atten_wave(p, *energ_args, **kwargs):
            wave_args = [hc/arg for arg in energ_args[::-1]]
            return dummy.calc(p, *wave_args, **kwargs)

load_user_model(atten_wave, "abs1")
add_user_pars("abs1", ['hcol','heiRatio','heiiRatio'],
    [dummy.hcol.val,dummy.heiRatio.val,dummy.heiiRatio.val],
    [dummy.hcol.min,dummy.heiRatio.min,dummy.heiiRatio.min],
    [dummy.hcol.max,dummy.heiRatio.max,dummy.heiiRatio.max],
    [dummy.hcol.units,dummy.heiRatio.units,dummy.heiiRatio.units],
    [False,False,False]) 

abs1.hcol = 1e+20
abs1.heiRatio = 0.1 
abs1.heiiRatio = 0.01

create_model_component("bpl1d", "bpow1")
bpow1.gamma1 = 0 
bpow1.gamma2 = 0 
bpow1.eb = 7.99625 
bpow1.ref = 1
freeze(bpow1.ref)
bpow1.ampl = 0.001

create_model_component("powlaw1d","pow1d")
pow1d.gamma = 1 
pow1d.ref = 1 
pow1d.ampl = 1e-5
pow1d.ampl.min = 2.383e-10

abs1.hcol = 1.81e20 
freeze(abs1)

set_model(1, abs1*bpow1)
set_model(2, abs1*bpow1)
set_model(3, abs1*bpow1)
set_model(4, abs1*bpow1)

set_bkg_model(1, abs1*pow1d, 1)
set_bkg_model(1, abs1*pow1d, 2)
set_bkg_model(2, abs1*pow1d, 1)
set_bkg_model(2, abs1*pow1d, 2)
set_bkg_model(3, abs1*pow1d, 1)
set_bkg_model(3, abs1*pow1d, 2)
set_bkg_model(4, abs1*pow1d, 1)
set_bkg_model(4, abs1*pow1d, 2)

set_method("neldermead")
set_stat("cstat")
fit()

plot("fit", 1, "fit", 2, "fit", 3, "fit", 4)

current_plot("all")
set_plot_title("3C 273 (ObsID 459)")

current_plot("plot1")
add_label(10, 0.1, "HEG -1")
set_label(["color","green"])
current_plot("plot2")
add_label(10, 0.1, "HEG +1")
set_label(["color","green"])
current_plot("plot3")
add_label(10, 0.15, "MEG -1")
set_label(["color","green"])
current_plot("plot4")
add_label(10, 0.15, "MEG +1")
set_label(["color","green"])

Sherpa 3.4

data 1 459_heg_m1_bin10.pha
data 2 459_heg_p1_bin10.pha
data 3 459_meg_m1_bin10.pha
data 4 459_meg_p1_bin10.pha

paramprompt off

rsp[hegm1]
rsp[hegp1]
rsp[megm1]
rsp[megp1]

hegm1.arf = 459_heg_m1.arf
hegp1.arf = 459_heg_p1.arf
megm1.arf = 459_meg_m1.arf
megp1.arf = 459_meg_p1.arf

instrument 1 = hegm1
instrument 2 = hegp1
instrument 3 = megm1
instrument 4 = megp1

ignore allsets all
notice allsets wave 1:15

paramprompt on

atten[abs]

bpl1d[bpow]

bpow1.gamma1 = 0 
bpow1.gamma2 = 0 
bpow1.eb = 7.99625 
bpow1.ref = 1
freeze(bpow1.ref)
bpow1.ampl = 0.001

powlaw1d[pow1d]

pow1d.gamma = 1 
pow1d.ref = 1 
pow1d.ampl = 1e-5
pow1d.ampl.min = 2.383e-10

abs.hcol=1.81e20 
freeze abs

source 1:4 = abs*bpow

background 1:4 = abs*pow1d

fit

lplot 4 fit 1 fit 2 fit 3 fit 4

d 1,3,4 ylabel ""
title "3C 273 (ObsID 459)"

d 1 label 12 0.075 "HEG -1"
d 2 label 12 0.075 "HEG +1"
d 3 label 12 0.125 "MEG -1"
d 4 label 12 0.125 "MEG +1"
d all l all green

redraw

For detailed information about each of the steps in the script, see the Sherpa thread "Fitting Grating Data", which also contains links to the ahelp files for each Sherpa function used.

High-resolution power-law model components fit to HRC-S/LETG source grating spectrum

Here we display the unresolved -1 HRC-S/LETG spectral order fit with the high-resolution model components corresponding to the -1/-2/-3 HRC-S/LETG spectral orders. This results from a simultaneous fit of the +1/-1 background-subtracted HRC-S/LETG source grating spectra, with the source modeled by the XSpec version of an absorbed broken power-law.

Sherpa 4.13 (Python)

load_pha(1, "460_leg_m1_bin10.pha") 
load_pha(2, "460_leg_p1_bin10.pha")

load_multi_arfs(1, ["460_LEG_-1.garf","460_LEG_-2.garf","460_LEG_-3.garf"], [1,2,3])
load_multi_rmfs(1, ["460_leg_-1.grmf","460_leg_-2.grmf","460_leg_-3.grmf"], [1,2,3])
load_multi_arfs(2, ["460_LEG_1.garf","460_LEG_2.garf","460_LEG_3.garf"], [1,2,3])
load_multi_rmfs(2, ["460_leg_1.grmf","460_leg_2.grmf","460_leg_3.grmf"], [1,2,3])

set_analysis("wave")

ignore('49.:56., 59.:66.')

set_model(xswabs.abs1*xsbknpower.bpow1)
set_model(2, xswabs.abs1*xsbknpower.bpow1)

abs1.nh = 0.1 
bpow1.phoindx1 = 1 
bpow1.breake = 5
bpow1.phoindx2 = 2 
bpow1.norm = 0.0434012
bpow1.integrate = "true"

abs1.nh = 1.81e-02
freeze(abs1.nh)

bpow1.breake = 1

set_method("neldermead")
set_method_opt("finalsimplex", 0)
set_stat("chi2xspecvar")

subtract()
subtract(2)
group_counts(1, 1)
group_counts(2, 1)

fit()

plot_data()

plot_model(overplot=1)

plot_order(1,1,overplot=1)
set_histogram("line.color=darkred")

plot_order(1,2,overplot=1)
set_histogram("line.color=pink")

plot_order(1,3,overplot=1)
set_histogram("line.color=orange")

log_scale(Y_AXIS)

title=("Model Orders + ([ " +
              "[\\color{darkred}-1, ]" +
              "[\\color{pink}-2, ]" +
              "[\\color{orange}-3, ]")

set_plot_title(title)

set_plot_ylabel("normalized counts sec^{-1} \\AA^{-1}")

set_plot_xlabel("m\\lambda [\\AA]")

Sherpa 3.4

data 1 460_leg_m1_bin10.pha 
data 2 460_leg_p1_bin10.pha 

data 3 460_leg_m1_bin10.pha 
data 4 460_leg_m1_bin10.pha 
data 5 460_leg_m1_bin10.pha 

frmf1d[rmfm1](460_leg_-1.grmf)
frmf1d[rmfm2](460_leg_-2.grmf)
frmf1d[rmfm3](460_leg_-3.grmf)
frmf1d[rmfp1](460_leg_1.grmf)
frmf1d[rmfp2](460_leg_2.grmf)
frmf1d[rmfp3](460_leg_3.grmf)

farf1d[arfm1](460_LEG_-1.garf)
farf1d[arfm2](460_LEG_-2.garf)
farf1d[arfm3](460_LEG_-3.garf)
farf1d[arfp1](460_LEG_1.garf)
farf1d[arfp2](460_LEG_2.garf)
farf1d[arfp3](460_LEG_3.garf)

instrument 1 = arfm1*rmfm1 + arfm2*rmfm2 + arfm3*rmfm3
instrument 2 = arfp1*rmfp1 + arfp2*rmfp2 + arfp3*rmfp3

instrument 3 = arfm1*rmfm1
instrument 4 = arfm2*rmfm2
instrument 5 = arfm3*rmfm3

d all limits x 45 70       

ignore 1 wave 49:56
ignore 2 wave 59:66

xswabs[abs]
 
xsbknpower[bpow]

abs.nH=1.81E-02
freeze abs.nh

bpow.2=1

source 1,2,3,4,5 = (abs*bpow)

subtract 1,2,3,4,5

fit 1,2,3,4,5

oplot fit 1 model 1 model 3 model 4 model 5 

c 2 red
c 3 green
c 4 blue
c 5 yellow

d 1 label 125 0.04 "Model orders -1, -2, -3"
l 1 size 1.5

redraw

For detailed information about each of the steps in the script, see the Sherpa thread "Fitting Multiple Orders of HRC-S/LETG Data", which also contains links to the ahelp files for each Sherpa function used.

2D Image data: DS9 data-to-model counts ratio image

Here we display the data-to-model ratio image of a 2D Image source data set in DS9, Sherpa's default imager. The multi-component, source-plus-background model expression used to fit the data is the sum of two 2D Gaussians (bulk and core emission) and a constant (background).

Sherpa 4.13 (Python)

load_image("image2.fits")

set_coord("physical")
set_stat("cash")
set_method("simplex")

notice2d("circle(4072.46,4249.34,108)")

set_model(gauss2d.g1)
g1.ampl = 20
g1.fwhm = 20
g1.xpos = 4065.5
g1.ypos = 4250.5
g1.fwhm.max = 4300
g1.xpos.max = 4300
g1.ypos.max = 4300
g1.ampl.min = 1
g1.ampl.max = 1000

set_model(g1+const2d.bgnd)
bgnd.c0 = 0.2
bgnd.c0.max = 100

fit()

freeze(g1)
set_model(gauss2d.g2+g1+bgnd)
g2.fwhm = 10
g2.ampl = 100
g2.xpos = 4065.5
g2.ypos = 4250.5
g2.fwhm.max = 4300
g2.xpos.max = 4300
g2.ypos.max = 4300
g2.ampl.min = 1
g2.ampl.max = 1000

fit()

thaw(g2.ellip)
thaw(g2.theta)
g2.ellip = 0.1
g2.theta = 1

fit()

image_ratio()

Sherpa 3.4

data image2.fits

coord physical
statistic cash
method simplex

notice physical "circle(4072.46,4249.34,108)"

source = gauss2d.g1
g1.ampl = 20
g1.fwhm = 20
g1.xpos = 4065.5
g1.ypos = 4250.5
g1.fwhm.max = 4300
g1.xpos.max = 4300
g1.ypos.max = 4300
g1.ampl.min = 1
g1.ampl.max = 1000

source = g1+const2d.bgnd 
bgnd.c0 = 0.2
bgnd.c0.max = 100

fit

freeze g1
source =  gauss2d.g2+g1+bgnd
g2.fwhm = 10
g2.ampl = 100
g2.xpos = 4065.5
g2.ypos = 4250.5
g2.fwhm.max = 4300
g2.xpos.max = 4300
g2.ypos.max = 4300
g2.ampl.min = 1
g2.ampl.max = 1000

fit

thaw g2.ellip
thaw g2.theta
g2.ellip = 0.1
g2.theta = 1

fit

image ratio

For detailed information about the steps in the script, see the Sherpa thread "Fitting FITS Image Data", which also contains links to the ahelp files for each Sherpa function used.

2D Image data: 2D Gaussian model radial profile fit with residuals

Here we display the radial profile of a 2D image data set, for which the profile center, ellipticity, and position angle are determined by source model values, and the bin width and radial extent in pixels are taken from the data. The model fit to the radial profile with delta chi residuals is also shown; the source-plus-background emission is modeled by the sum of two 2D Gaussians (bulk and core emission) and a constant (background).

Sherpa 4.13 (Python)

from sherpa_contrib import *
load_image("image2.fits")

set_coord("physical")
set_stat("cash")
set_method("simplex")

notice2d("circle(4072.46,4249.34,108)")
 
set_model(gauss2d.src)

guess(src)

set_par(src.fwhm, 20, 0.1, 300)
set_par(src.ampl, 20, 0.1, 1000)

set_model(src + const2d.bgnd)
bgnd.c0 = 0.2
bgnd.c0.max = 100

src.fwhm = src.fwhm.val * 2

fit()

set_source(bgnd + src + gauss2d.core)

guess(core)
set_par(core.fwhm, 10, 0.1, 100)
set_par(core.ampl, 100, 0.1, 1000)
set_par(core.ellip, 0.1)
set_par(core.theta, 1)


fit()

get_data_prof_prefs()["ylog"] = True
get_delchi_prof_prefs()["xlog"] = True

prof_fit_delchi(model=src, group_counts=200)

limits(X_AXIS, 0.5, AUTO)

Sherpa 3.4

data image2.fits

coord physical
statistic cash
method simplex

ignore all
notice physical "circle(4072.46,4249.34,108)"

paramprompt off

source = gauss2d[src1]

src1.fwhm = 20
src1.fwhm.min = 0.1
src1.fwhm.max = 300

src1.ampl = 20
src1.ampl.min = 0.1
src1.ampl.max = 1000

source = src1 + const2d[bgnd1]
bgnd1 integrate off
bgnd1.c0 = 0.2
bgnd1.c0.max = 100

src1.fwhm = src1.fwhm.val * 2

fit

source = gauss2d[core1] + src1 + bgnd1
core1.fwhm = 10
core1.fwhm.min = 0.1
core1.fwhm.max = 100

core1.ampl = 100
core1.ampl.min = 0.1
core1.ampl.max = 1000

thaw core1.ellip core1.theta
core1.ellip = 0.1
core1.theta = 1
simplex.iters = 3000

fit

set_log

() = evalfile("sherpa_plotfns.sl");

plot_rprofr("src1",0,150,5)

For detailed information about each of the steps in the script, see the Sherpa thread "Radial and elliptical profiles of Image Data."

2D Image data: 1D Beta model radial profile fit with residuals

Here we display a 1D Beta model fit to a background-subtracted radial profile extracted from a 2D Image data set (either spatial table or image file) with the CIAO tool dmextract.

Sherpa 4.13 (Python)

load_data(1, "1838_rprofile_rmid.fits", 3, \ 
["RMID","SUR_BRI","SUR_BRI_ERR"])

set_source("beta1d.sbr1")
sbr1.r0 = 105
sbr1.beta = 4
sbr1.ampl = 0.00993448

freeze(sbr1.xpos)

fit() 

plot_fit()
log_scale(XY_AXIS)

limits(X_AXIS, 10, 200)
limits(Y_AXIS, 0.0001, 10)

Sherpa 3.4

read data 1 "1838_rprofile_rmid.fits[columns rmid,sur_bri]" FITSBIN
read errors 1 "1838_rprofile_rmid.fits[columns rmid,sur_bri_err]" FITSBIN

beta1d[sbr1]

sbr1.ampl.max=10
show sbr1                     

source=sbr1
fit

lplot fit
log
limits y 0.0001 10  
limits x 10 200
redraw

For detailed information about the input data in the script, see the CIAO thread "Obtain and Fit a Radial Profile."

2D Image data: DS9 counts images of data, PSF-convolved model, and fit residuals

Here we display three images in one DS9 window: the counts image of a 2D image source data set, the corresponding 2D PSF-convolved model to be fitted to the data, and the counts residuals (data-model) of the fit. The multi-component source-plus-background model expression is defined by convolving a PSF kernel read from an image file with a 2D Gaussian plus constant model.

Sherpa 4.13 (Python)

load_image("center_box_0.25pix.fits")
image_data()
image_getregion()
notice2d("rotbox(88.16875,75.8625,70.416667,68.508334,0)")

load_psf("psf0", "psf_f1_norm_0.25pix.fits")

set_psf(psf0)
psf0.size=[32,32]
psf0.center=[128,129]
print get_psf()

set_psf(psf0)
psf0.center=[128,129]
psf0.size=[72,72]

set_source(const2d.c1+gauss2d.g1)
c1.c0=1
g1.fwhm = 6
g1.xpos = 90
g1.ypos = 80
g1.ampl = 100


set_stat("cash")
freeze(c1.c0)
set_method("neldermead")
set_method_opt("iquad",0)
set_method_opt("finalsimplex",0)
fit()

thaw(c1.c0)
fit()

image_fit()

Sherpa 3.4

data center_box_0.25pix.fits

ignore all
notice filter "box(88.16875,75.8625,70.416667,68.508334)"

paramprompt off

fpsf[psf0]
psf0.file = psf_f1_norm_0.25pix.fits
                        
instrument = psf0

psf0.ysize=72
psf0.xsize=72
psf0.xoff=0
psf0.yoff=0

source = const2d[c1] + gauss2d[g1]
g1 integrate on
c1.c0=1
freeze c1.c0


statistic cash
method simplex
fit

thaw c1.c0
fit

image fit

For detailed information about the steps in the script, see the Sherpa thread "Accounting for PSF Effects in 2D Image Fitting", which also contains links to the ahelp files for each Sherpa function used.

2D Image data: Contour plot of residuals of PSF-convolved model fit

Here we display a contour plot of the counts residuals (data-model) resulting from the fit of a PSF-convolved Gaussian model fit to a 2D image source data set (see the fitting example DS9 counts images of data, PSF-convolved model, and fit residuals for model details).

Sherpa 4.13 (Python)

load_image("center_box_0.25pix.fits")

notice2d("rotbox(88.16875,75.8625,70.416667,68.508334,0)")

load_psf("psf0", "psf_f1_norm_0.25pix.fits")

set_psf(psf0)
psf0.size=[32,32]
psf0.center=[128,129]

set_psf(psf0)
psf0.center=[128,129]
psf0.size=[72,72]

set_source(const2d.c1+gauss2d.g1)
c1.c0=1
g1.fwhm = 6
g1.xpos = 90
g1.ypos = 80
g1.ampl = 100

show_model()

set_stat("cash")
freeze(c1.c0)
set_method("neldermead")
set_method_opt("iquad",0)
set_method_opt("finalsimplex",0)
fit()

thaw(c1.c0)
fit()

image_fit()

contour_resid()
limits(X_AXIS, 55, 120)
limits(Y_AXIS, 45, 110)

Sherpa 3.4

data center_box_0.25pix.fits

ignore all
notice filter "box(88.16875,75.8625,70.416667,68.508334)"

paramprompt off

fpsf[psf0]
psf0.file = psf_f1_norm_0.25pix.fits

instrument = psf0
psf0.ysize=72
psf0.xsize=72
psf0.xoff=0
psf0.yoff=0

source = const2d[cc1] + gauss2d[g2]
g2 integrate on
cc1.c0=1
freeze cc1.c0

statistic cash
method simplex 
fit

thaw cc1.c0
fit

image fit

cplot residuals

For detailed information about the steps in the script, see the Sherpa thread "Accounting for PSF Effects in 2D Image Fitting."

Parameter bounds with interval-projection: Confidence plot of fit stastic vs. parameter value

Here we display the confidence plot for a particular model parameter after finding the best-fit of the model to a 1D PHA source data set. The minimum of the parabola corresponds to the best-fit value for this particular model parameter (lowest chi-square fit statistic), while the arms of the parabola demonstrate how the chi-square fit statistic varies within the 90% (1.6 sigma) confidence bounds of the parameter value.

Sherpa 4.13 (Python)

load_data("source_grouped_pi.fits")

notice(0.5,8.0)

set_source(xswabs.abs1*powlaw1d.p1)

fit()

set_proj_opt("sigma",1.6)

int_proj(p1.gamma,min=1,max=2)

Sherpa 3.4

data source_grouped_pi.fits

ignore energy :0.5,8:

source = xswabs[abs1] * powlaw1d[p1]

fit

restore_intproj

sherpa.intproj.sigma = [1,1.6]
sherpa.intproj.arange = 0
sherpa.intproj.min = 1
sherpa.intproj.max = 2

intproj p1.gamma

For detailed information about the steps in the script, see the Sherpa thread "Step-by-Step Guide to Estimating Errors and Confidence Levels."

Parameter bounds with region-projection: Confidence contour of fit statistic vs. two parameter values

Here we display the confidence contour plot for two different model parameters after finding the best-fit of the model to a 1D PHA source data set. The center point of the ellipses corresponds to the set of best-fit values for the two model parameters (where the chi square fit statistic is at a minimum), while the inner ellipse represents the 68.3% (1 sigma) confidence bounds on the correlated values, and the outer ellipse the 90% (1.6 sigma) bounds. The chi square fit statistic varies along the axis perpendicular to the plot (i.e. out of the screen).

Sherpa 4.13 (Python)

load_data("source_grouped_pi.fits")

notice(0.5,8.0)

set_source(xswabs.abs1*powlaw1d.p1)

fit()

set_proj_opt("sigma",1.6)

reg_proj(p1.gamma,abs1.nH,min=[1.2,2],max=[1.9,2.8], \ 
nloop=[100,100],sigma=[1,1.6])

Sherpa 3.4

data source_grouped_pi.fits

ignore energy :0.5,8:

source = xswabs[abs1] * powlaw1d[p1]

fit

restore_regproj

sherpa.regproj.sigma = [1,1.6]
sherpa.regproj.arange = 0
sherpa.regproj.min = [1.2,2]
sherpa.regproj.max = [1.9,2.8]
sherpa.regproj.nloop = [100,100]

regproj p1.gamma abs1.nh

For detailed information about the steps in the script, see the Sherpa thread "Step-by-Step Guide to Estimating Errors and Confidence Levels."

/ sherpa / gallery / fit.html