Simulating spectra with PyXSPEC

Contact: honghui.liu<at>uni-tuebingen.de; menglei.zhou<at>astro.uni-tuebingen.de

Why doing simulation:

• To test performance of a model with a certain instrument, e.g., Liu et al. 2025.

• Monte Carlo simulations to compare models (e.g., line detection) on real data.

• To test capability of future instruments.

• And more…

import xspec
import os
import glob

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[1], line 1
----> 1 import xspec
      2 import os
      3 import glob

ModuleNotFoundError: No module named 'xspec'

def fake_pha(infile):

    xspec.AllModels.setEnergies("0.01 500. 10000 log")
    
    # Define the model for simulation
    # Here `infile` is a additive model stored in FITS format
    model = xspec.Model(f"const*tbabs*(cflux*cutoffpl + cflux*atable{{{infile}}})")
    
    # Define the parameters and the flux for each component.
    # Parameters must be resonable
    model.setPars({2:"0.6 -0.1", 3:"20 -0.1", 4:"40 -0.1", 6:"2.0", 
                   7:"300. -0.1", 9:"=p3", 10:"=p4"})


    # Set const_1 to get the right overall flux
    # Set the total constant to ensure the flux between 2.0 and 10.0 keV 
    # is 2e-8 erg/s/cm^2

    xspec.AllModels.calcFlux("2. 10.")
    const1 = (2.0e-8/xspec.AllModels(1).flux[0])*model(1).values[0]
    model(1).values=const1

    xspec.AllModels.calcFlux('')

    # Save the model

    label=infile.split(".fits")[0]
    xspec.Xset.save(f"{label}_ref.mdl", info="m")

    # Clear the spectrum
    # xspec.AllData.clear()

    ## NuSTAR/FPM

    fs1 = xspec.FakeitSettings(response="nustar.rmf",
                         arf="point_60arcsecRad_1arcminOA.arf",
                         background="bgd_60arcsec.pha",
                         exposure=30000.0, correction=1.0, backExposure=30000.0,
                         fileName=f"fpm_{label}_30ks.pha")
    ## NICER/XTI

    fs2 = xspec.FakeitSettings(response="nixtiref20170601v003.rmf",
                         arf="nixtiaveonaxis20170601v005.arf",
                         background="nixtiback20190807.pi",
                         exposure=5000.0, correction=1.0, backExposure=5000.0,
                         fileName=f"xti_{label}_5ks.pha")

    xspec.AllData.fakeit(2, [fs1, fs2])

    xspec.AllModels.clear()
    xspec.AllData.clear()

if __name__ == "__main__":
    os.environ['HEADASNOQUERY'] = ''
    os.environ['HEADASPROMPT'] = '/dev/null'

    # Simulate spectra for all models in the dir
    for infile in glob.glob("*.fits"):
        #infile = sys.argv[1]

        fake_pha(infile)

After faking the spectra, analyze them as if they are real data in XSPEC.