# How can I plot an unfolded source spectrum in flux units, independent of a model spectrum?

The primary function of software packages like Sherpa is to
solve the integral equation (instrument response) which
relates the total counts produced in a given pulse-height
bin in a spectrum to the incident source
spectrum. Therefore, once you have a pulse-height spectrum
and the appropriate RMF and ARF instrument response, you
can use Sherpa to fold the source data through an assumed
model to determine the "true" source flux. This means that
*unfolding* the source spectrum to plot it in flux
units, for example (Sherpa Plotting FAQ #2), is dependent upon a
model spectrum.

In order to unfold a spectrum in a model-independent way, we assume a constant source model set to unity, as shown in the following script:

# The electron charge, ~1.6e-9 ergs/1 keV photon (nist.gov) # from sherpa.astro.utils import charge_e as q # Assuming source and background data and corresponding responses have # been loaded into Sherpa: # set_analysis('energy') subtract() set_source(polynom1d.truespec) # Assume the units are in photon / s / cm^2 / keV rplot = get_ratio_plot() photon = rplot.y photerr = rplot.yerr # Assume the units are keV energy = rplot.x conv = q * energy * energy # convert the residuals spec = conv * photon specerr = conv * photerr # plot nu F_nu plt.figure() plt.errorbar(energy, spec, specerr, linestyle='none') plt.plot(energy, spec, '.') plt.xlabel("Energy (keV)") plt.ylabel(r"$\nu$ F$_{\nu}$ (ergs s$^{-1}$ cm$^{-2}$)") plt.yscale('log') # plot ph/s/cm**2/keV spectrum plt.figure() plt.errorbar(energy, photon, photerr, linestyle='none') plt.plot(energy, photon, '.') plt.xlabel("Energy (keV)") plt.ylabel("photons s$^{-1}$ cm$^{-2}$ keV$^{-1}$") plt.yscale('log') # write to disk colnames = ["keV", "FLUX", "FLUXerr"] save_arrays("plot_flux.energy.dat", [energy, spec, specerr], colnames, ascii=True) save_arrays("plot_flux.photon.dat", [energy, photon, photerr], colnames, ascii=True)

The unfolded spectrum is calculated with the plot_ratio() command, and obtained via get_ratio_plot().y, in units of photons/s/cm2/keV. The plot_ratio() command plots the ratio of a data set to its assigned source model; in this case, the observed spectrum is divided by the summation of the product of the responses, since we have assumed a constant source model set to unity.

To convert photons flux units to energy flux units, we define a conversion factor q ~ 1.6e-9 ergs/1 keV photon, based on the electron charge value obtained from the National Institute of Standards and Technology (NIST).