An example combining heasp and pyxspec

The following example uses pyxspec to calculate a set of fluxes given the input response, multiplies the result by the response to give predicted counts and writes the result out as a spectrum.

import os
from heasp import *
import xspec as xsp
import nump as np

responseName = "response.rsp"
inputRSP = rmf(responseName)

# energies needs to be the standard internal xspec energy array
numEnergies = inputRSP.NumberEnergyBins() + 1
energies = np.empty((numEnergies))
energies[0] = inputRSP.getLowEnergyElement(0)
for j in range(1,numEnergies): energies[j] =
                inputRSP.getHighEnergyElement(j-1)

# set up the parameters array as a list
params = [ 1.7 ]

# calculate the model and output flux as a list
fluxlist = []
xsp.callModelFunction("powerlaw", energies.tolist(), params, fluxlist)

# convert the flux to a numpy array and multiply by the response
flux = np.array(fluxlist)
phaValues = inputRSP.multiplyByModel(flux)

# set up the spectrum channel array
channel = np.arange(phaValues.size).astype(np.int32)

# pick an arbitrary exposure time of 10000
exposure = 10000.0
phaValues *= exposure

phaOut = pha()
phaOut.setFirstChannel(0)
phaOut.setPha(phaValues)
phaOut.setChannel(channel)
phaOut.setExposure(exposure)
phaOut.setDetChans(phaValues.size)
phaOut.setPoisserr(True)
phaOut.setDatatype("COUNT")
phaOut.setSpectrumtype("TOTAL")
phaOut.setResponseFile(responseName)
phaOut.setTelescope(inputRSP.getTelescope())
phaOut.setInstrument(inputRSP.getInstrument())
phaOut.setFilter(inputRSP.getFilter())

phafile = "outputspectrum.pha"
if (os.path.exists(phafile)): os.remove(phafile)
status = phaOut.write(phafile)
if status != 0: print("Failed to write outputspectrum.pha: status = ", status)

HEASARC Home | Observatories | Archive | Calibration | Software | Tools | Students/Teachers/Public

Last modified: Wednesday, 28-Aug-2024 16:49:27 EDT