Estimate Flux Errors

A demo for ADASS 2010 to estimate the integrated flux from an X-ray spectral fit using simulations.

#!/usr/bin/env python
import sherpa.astro.ui as sherpa
import numpy

# load the data and fit with an absorbed power-law
sherpa.load_pha("3c273.pi")
sherpa.notice(0.5,7.0)
sherpa.subtract()

sherpa.set_stat("chi2datavar")
sherpa.set_model(sherpa.xsphabs.abs1*sherpa.xspowerlaw.p1)
abs1.nh=0.07
sherpa.freeze(abs1.nh)
sherpa.fit()

# <demo> stop

# simulate fluxes by sampling model parameters for 100 iterations.

flux100=sherpa.sample_energy_flux(0.5,7.,num=100)

# slice the flux column from simulation table.

fluxes = flux100[:,0]

print "Mean:", numpy.mean(fluxes)

print "Median:", numpy.median(fluxes)

print "Std Dev:", numpy.std(fluxes)

# determine 95% and 50% quantiles using numpy array sorting

sf=numpy.sort(fluxes)

print "95% quantile:", sf[0.95*(len(fluxes)-1)]
print "50% quantile:", sf[0.5*(len(fluxes)-1)]

# <demo> stop

# run 7.5e3 flux simulations and bin a histogram using the result

hist=sherpa.get_energy_flux_hist(0.5,7.,num=7500)
hist.histo_prefs['linecolor']=None
sherpa.plot_energy_flux(recalc=False)

# <demo> stop

# fit a Gaussian to the histogram
sherpa.load_arrays(2,hist.xlo,hist.xhi,hist.y,sherpa.Data1DInt)

sherpa.set_model(2,sherpa.gauss1d.g0)
g0.integrate=False
sherpa.guess(2,g0)
g0.ampl.min=0
g0.ampl.max=100
g0.ampl=1

sherpa.set_stat("leastsq")
sherpa.fit(2)

sherpa.plot_model(2, overplot=True)

# <demo> stop

# calculate sigma instead of FWHM
# FWHM = sigma * sqrt(8 ln 2)

fac0=numpy.sqrt(8*numpy.log(2))
sig=g0.fwhm.val/fac0

print "Sigma: ",sig,"Position: ", g0.pos.val

# write binned histogram arrays to file

sherpa.save_arrays("hist.dat",[hist.xlo,hist.xhi,hist.y],
                   fields=["xlo","xhi","y"], clobber=True)

# <demo> stop