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).