Reproducing an input spectrum

marx allows users to specify either a flat input spectrum or to pass in a file. Here, we generate input spectra with the spectral modeling program Sherpa. We run marx to simulate a point source and extract the data from that simulation using CIAO, read it back into Sherpa and fit the simulated data. If everything works well and all detector effects are treated consistently, we should recover the same spectral parameters (up to Poisson noise). On the other hand, if marx and CIAO are inconsistent the fit parameters will deviate from the input parameters.

In the past, this has happened after the release of marx 5.0, which contains some files to describe the ACIS contamination. This contamination changes with time and several CalDB releases had happened before we released marx 5.1. At that time, marx always predicted too many counts at low energies.

The following tests are designed to test consistency with CIAO. Since there is always some uncertainty about the intrinsic spectrum of astrophysical sources, this test is best done with simulated input spectra.

Absorbed powerlaw on ACIS-S

This test checks the internal consistency of a simulated marx spectrum by simulating a source placed on a back-illuminated chip of ACIS-S. The input spectrum is an absorbed powerlaw.

Click to hide code

 
'''Set up a Sherpa model and use that to generate an spectrum as input for Marx.'''

import numpy as np
from sherpa.astro import ui

ui.set_source(1, ui.xsphabs.a * ui.xspowerlaw.p)
a.nH = 1.0
p.PhoIndex = 1.8
p.norm = 0.001

# get source
my_src = ui.get_source()

# set energy grid
bin_width = 0.01
energies = np.arange(0.15, 12., bin_width)

# evaluate source on energy grid
flux = my_src(energies[:-1], energies[1:])

ui.save_arrays("marx_input.tbl", [energies[1:], flux / bin_width], 
               ["keV","photons/s/cm**2/keV"], ascii=True, clobber=True)
 
import subprocess

out = subprocess.run('marx SourceFlux=-1 SpectrumType=FILE SpectrumFile=marx_input.tbl ExposureTime=30000 ' + 
                     'TStart=2015.5 OutputDir=marx GratingType=NONE DetectorType=ACIS-S DitherModel=INTERNAL ' +
                     'RA_Nom=30 Dec_Nom=40 SourceRA=30 SourceDEC=40', shell=True, capture_output=True)
out2 = subprocess.run('marx2fits --pixadj=EDSER marx marx_evt2.fits', shell=True, capture_output=True)
out3 = subprocess.run('marxasp MarxDir=marx OutputFile=marx_asol1.fits', shell=True, capture_output=True)

We generate an input spectrum and run marx. We then extract the spectrum in two different ways:

  1. We use CIAO tools, but specify special qualifiers to account for the fact that marx does not include ACIS QE maps and uses FEF rmfs (FEF’s are parameterized functions to describe the shape of the RMF. They are faster to compute, but not as accurate as the tabulated values usually used in CIAO). We do that by setting the detector name to "ACIS-7;UNIFORM;bpmask=0" in mkarf and by passing the appropriate FEF to mkrmf.

  2. We use the specextract tool without any special qualifiers. While specextract has a number of parameters, it does not allow us to set the dectector or force the use of an FEF file. It requires only a call to a single tool, but is inconsistent with the simulations done by marx, so fitted spectra will not match the input parameters as well as in method 1.

Click to hide code

 
'''Extract spectra

Use special settings to account for steps that |marx| does not include:
ACIS QE maps and current (non FEF) rmfs.

We warp this into a function because we are going to use the same steps
for more examples further down in the notebook and thus would be a fairly long
piece of code to copy and paste.
'''
import os

import ciao_contrib.runtool as rt
from coords.chandra import sky_to_chandra
from astropy.io import fits

def special_extraction(x, y, radius):
    phagrid="pi=1:1024:1"

    evtfile="marx_evt2.fits"
    asolfile="marx_asol1.fits"
    phafile="marx_pha.fits"
    rmffile="marx_rmf.fits"
    arffile="marx_arf.fits"
    asphistfile="marx_asp.fits"

    rt.asphist(infile=asolfile, outfile=asphistfile, evtfile=evtfile, clobber=True)

    rt.dmextract(infile=f"{evtfile}[sky=circle({x},{y},{radius})][bin {phagrid}]", 
                outfile=phafile, clobber=True)
    coos = sky_to_chandra(fits.getheader(evtfile, 1), x, y)
    ccdid = coos['chip_id'][0]
    chipx = coos['chipx'][0]
    chipy = coos['chipy'][0]

    detname=f"ACIS-{ccdid};UNIFORM;bpmask=0"
    # For ACIS-I, use engrid="0.3:11.0:0.003". This reflects a limitation of mkrmf.
    engrid="0.3:12.0:0.003"

    rt.mkarf(mirror="hrma", detsubsys=detname, grating="NONE",
    outfile=arffile, obsfile=evtfile, engrid=engrid, asphistfile=asphistfile,
    sourcepixelx=x, sourcepixely=y, maskfile="NONE", pbkfile="NONE", dafile="NONE", clobber=True, verbose=0)

    fef=os.environ['CALDB'] + "/data/chandra/acis/fef_pha/acisD2000-01-29fef_phaN0005.fits"

    cxfilter=f"chipx_hi>={chipx},chipx_lo<={chipx}"
    cyfilter=f"chipy_hi>={chipy},chipy_lo<={chipy}"
    rt.mkrmf(infile=f"{fef}[ccd_id={ccdid},{cxfilter},{cyfilter}]", outfile=rmffile,
            axis1=f"energy={engrid}", axis2=phagrid, thresh=1e-8, clobber=True, verbose=0)
    
special_extraction(4096.4, 4096.4, 20.0)

Click to hide code

 
from sherpa.utils.logging import SherpaVerbosity
# Only display error messages, but hide other output from specextract
with SherpaVerbosity('ERROR'):
    rt.specextract(f"marx_evt2.fits[sky=circle(4096.5,4096.5,20)]", "spec", asp="marx_asp.fits", weight=False,
                clobber=True, mskfile="NONE", badpixfile="NONE")
       

with SherpaVerbosity('ERROR'):
    ui.load_pha("special extraction", 'marx_pha.fits')
    ui.load_rmf("special extraction", 'marx_rmf.fits')
    ui.load_arf("special extraction", 'marx_arf.fits')
    ui.load_data("specextract", 'spec_grp.pi')
    ui.group_counts("special extraction", 25)
    ui.group_counts("specextract", 25)
    ui.ignore(None, 0.3)
    ui.ignore(7, None)
    # set source properties
    ui.set_source("special extraction", ui.xsphabs.a1 * ui.xspowerlaw.p1)
    ui.set_source("specextract", ui.xsphabs.a2 * ui.xspowerlaw.p2)

with SherpaVerbosity('ERROR'):
    ui.fit("special extraction")
    ui.conf("special extraction")
    conf_special = ui.get_confidence_results()
    
with SherpaVerbosity('ERROR'):
    ui.fit("specextract")
    ui.conf("specextract")
    conf_specextract = ui.get_confidence_results()

Click to hide code

 
import matplotlib.pyplot as plt

# Better setting to display labels on a log axis
import matplotlib as mpl
mpl.rcParams['axes.formatter.min_exponent'] = 2

a1 = ui.get_arf("special extraction")
a2 = ui.get_arf("specextract")
# bring on same energy grid
engrid = np.arange(0.3, 10., 0.05)
a1 = np.interp(engrid, 0.5 * (a1.energ_lo + a1.energ_hi), a1.get_y())
a2 = np.interp(engrid, 0.5 * (a2.energ_lo + a2.energ_hi) , a2.get_y())


fig, axes = plt.subplots(2, 1, sharex=True, figsize=(8, 6))
plt.sca(axes[0])
ui.plot_data("special extraction", xlog=True, ylog=True, clearwindow=False)
ui.plot_model("special extraction", overplot=True, label="special extraction")
ui.plot_model("specextract", overplot=True, label="specextract")
axes[0].legend()
axes[0].set_title("Simulated data with fits using different ARFs and RMFs")

axes[1].loglog(engrid, a1, label="special extraction", c="C1")
axes[1].loglog(engrid, a2, label="specextract", c="C2")
axes[1].set_xlabel("Energy (keV)")
axes[1].set_ylabel("Effective Area (cm$^2$)")
axes[1].legend()

# Better setting to display labels on a log axis
axes[1].get_xaxis().set_minor_formatter(mpl.ticker.LogFormatterMathtext(labelOnlyBase=False,
                                                               minor_thresholds=(2, .5)))
axes[1].tick_params(axis='x', labelsize=mpl.rcParams['xtick.labelsize'], which='both')
../_images/1afa358c21f4b216cf9a832451b368a8d1ced6b01230af9cfbb218a33f1fe64f.png

top: Simulated data and two model fits, assuming different ARFs and RMFs. “Special extraction” uses special qualifiers to match marx simulations, while “specextract” uses standard CIAO tools without any special qualifiers.

bottom: Shape of the ARFs used in the two extractions. Since the ARF is different, the model values are different when fitting to the same data (e.g. where the ARF is lower, the model flux has to be higher to match the data, leading to different slopes and normalizations in the two fits).

Input parameters

Click to hide code

 
my_src
<BinaryOpModel model instance 'xsphabs.a * xspowerlaw.p'>

Fit with special extraction

Click to hide code

 
conf_special
<confidence results instance>

Fit with CIAO default extraction

Click to hide code

 
conf_specextract
<confidence results instance>

The tables above show the input parameters and the best-fit parameters for the two extraction methods. As expected, the fit using the special extraction method recovers the input parameters better (and typically within the uncertainties), while the standard specextract method shows larger deviations.

Powerlaw spectrum of an off-axis source

Same as above, but for an off-axis source on ACIS-I.

Click to hide code

 
out = subprocess.run('marx SourceFlux=-1 SpectrumType=FILE SpectrumFile=marx_input.tbl ExposureTime=60000 ' + 
                     'TStart=2015.5 OutputDir=marx GratingType=NONE DetectorType=ACIS-I DitherModel=INTERNAL ' +
                     'RA_Nom=30 Dec_Nom=40 SourceRA=29.9814917 SourceDEC=40.1508083', shell=True, capture_output=True)
out2 = subprocess.run('marx2fits --pixadj=EDSER marx marx_evt2.fits', shell=True, capture_output=True)
out3 = subprocess.run('marxasp MarxDir=marx OutputFile=marx_asol1.fits', shell=True, capture_output=True)

Click to hide code

 
special_extraction(4200, 5200, 200)


with SherpaVerbosity('ERROR'):
    rt.specextract(f"marx_evt2.fits[sky=circle(4200,5200,200)]", "spec", asp="marx_asp.fits", weight=False,
                clobber=True, mskfile="NONE", badpixfile="NONE")

with SherpaVerbosity('ERROR'):
    ui.load_pha("special extraction", 'marx_pha.fits')
    ui.load_rmf("special extraction", 'marx_rmf.fits')
    ui.load_arf("special extraction", 'marx_arf.fits')
    ui.load_data("specextract", 'spec_grp.pi')
    ui.group_counts("special extraction", 25)
    ui.group_counts("specextract", 25)
    ui.ignore(None, 0.3)
    ui.ignore(7, None)
    # set source properties
    ui.set_source("special extraction", ui.xsphabs.a1 * ui.xspowerlaw.p1)
    ui.set_source("specextract", ui.xsphabs.a2 * ui.xspowerlaw.p2)

with SherpaVerbosity('ERROR'):
    ui.fit("special extraction")
    ui.conf("special extraction")
    conf_special = ui.get_confidence_results()
with SherpaVerbosity('ERROR'):
    ui.fit("specextract")
    ui.conf("specextract")
    conf_specextract = ui.get_confidence_results()

Click to hide code

 
a1 = ui.get_arf("special extraction")
a2 = ui.get_arf("specextract")
# bring on same energy grid
engrid = np.arange(0.3, 10., 0.05)
a1 = np.interp(engrid, 0.5 * (a1.energ_lo + a1.energ_hi), a1.get_y())
a2 = np.interp(engrid, 0.5 * (a2.energ_lo + a2.energ_hi) , a2.get_y())


fig, axes = plt.subplots(2, 1, sharex=True, figsize=(8, 6))
plt.sca(axes[0])
ui.plot_data("special extraction", xlog=True, ylog=True, clearwindow=False)
ui.plot_model("special extraction", overplot=True, label="special extraction")
ui.plot_model("specextract", overplot=True, label="specextract")
axes[0].legend()
axes[0].set_title("Simulated data with fits using different ARFs and RMFs")

axes[1].loglog(engrid, a1, label="special extraction", c="C1")
axes[1].loglog(engrid, a2, label="specextract", c="C2")
axes[1].set_xlabel("Energy (keV)")
axes[1].set_ylabel("Effective Area (cm$^2$)")
axes[1].legend()

# Better setting to display labels on a log axis
axes[1].get_xaxis().set_minor_formatter(mpl.ticker.LogFormatterMathtext(labelOnlyBase=False,
                                                               minor_thresholds=(2, .5)))
axes[1].tick_params(axis='x', labelsize=mpl.rcParams['xtick.labelsize'], which='both')
../_images/106e5fd9961c89ffb0c4f3bd9e4fdff9e08c554f2a0ea8326e9e042d62adb207.png

Input parameters

(some as above, but printed here again)

Click to hide code

 
my_src
<BinaryOpModel model instance 'xsphabs.a * xspowerlaw.p'>

Fit with special extraction

Click to hide code

 
conf_special
<confidence results instance>

Fit with CIAO default extraction

Click to hide code

 
conf_specextract
<confidence results instance>

For the off-axis simulation, the difference between the two extraction methods is smaller, probably because we are using the ACIS-I detector here.

Two thermal components on ACIS-I, on-axis

Same as above, but for a front-illuminated ACIS-I chip and using an input spectrum with two thermal components similar to a stellar corona.

Click to hide code

 
ui.set_source(1, ui.xsvapec.apec1 + ui.xsvapec.apec2)
apec1.Ne.frozen = False
apec2.Ne = apec1.Ne
apec1.kT = 0.7
apec1.Ne = 2.0
apec1.norm = 0.0002
apec2.kT = 2.0
apec2.norm = 0.0004

# get source
my_src = ui.get_source(1)

# set energy grid
bin_width = 0.01
energies = np.arange(0.15, 12., bin_width)

# evaluate source on energy grid
flux = my_src(energies[:-1], energies[1:])

ui.save_arrays("marx_input.tbl", [energies[1:], flux / bin_width], 
               ["keV","photons/s/cm**2/keV"], ascii=True, clobber=True)

out = subprocess.run('marx SourceFlux=-1 SpectrumType=FILE SpectrumFile=marx_input.tbl ExposureTime=30000 ' + 
                     'TStart=2015.5 OutputDir=marx GratingType=NONE DetectorType=ACIS-I DitherModel=INTERNAL ' +
                     'RA_Nom=30 Dec_Nom=40 SourceRA=30 SourceDEC=40', shell=True, capture_output=True)
out2 = subprocess.run('marx2fits --pixadj=EDSER marx marx_evt2.fits', shell=True, capture_output=True)
out3 = subprocess.run('marxasp MarxDir=marx OutputFile=marx_asol1.fits', shell=True, capture_output=True)

Reading APEC data from 3.0.9

Click to hide code

 
special_extraction(4096.4, 4096.4, 20.0)
with SherpaVerbosity('ERROR'):
    rt.specextract(f"marx_evt2.fits[sky=circle(4096.5,4096.5,20)]", "spec", asp="marx_asp.fits", weight=False,
                clobber=True, mskfile="NONE", badpixfile="NONE")
       

with SherpaVerbosity('ERROR'):
    ui.load_pha("special extraction", 'marx_pha.fits')
    ui.load_rmf("special extraction", 'marx_rmf.fits')
    ui.load_arf("special extraction", 'marx_arf.fits')
    ui.load_data("specextract", 'spec_grp.pi')
    ui.group_counts("special extraction", 25)
    ui.group_counts("specextract", 25)
    ui.ignore(None, 0.3)
    ui.ignore(7, None)
    # set source properties
    ui.set_source("special extraction", ui.xsvapec.apec1_1 + ui.xsvapec.apec2_1)
    ui.set_source("specextract", ui.xsvapec.apec1_2 + ui.xsvapec.apec2_2)

    apec1_1.Ne.frozen = False
    apec2_1.Ne = apec1_1.Ne
    apec1_2.Ne.frozen = False
    apec2_2.Ne = apec1_2.Ne

    # Set some starting values close to the true value
    apec1_1.kT = 0.9
    apec1_1.Ne = 1.5
    apec1_1.norm = 0.00015
    apec2_1.kT = 1.8
    apec2_1.norm = 0.00035
    apec1_2.kT = 0.9
    apec1_2.Ne = 1.5
    apec1_2.norm = 0.00015
    apec2_2.kT = 1.8
    apec2_2.norm = 0.00035  
with SherpaVerbosity('ERROR'):
    ui.fit("special extraction")
    ui.conf("special extraction")
    conf_special = ui.get_confidence_results()
with SherpaVerbosity('ERROR'):
    ui.fit("specextract")
    ui.conf("specextract")
    conf_specextract = ui.get_confidence_results()

Click to hide code

 
a1 = ui.get_arf("special extraction")
a2 = ui.get_arf("specextract")
# bring on same energy grid
engrid = np.arange(0.3, 10., 0.05)
a1 = np.interp(engrid, 0.5 * (a1.energ_lo + a1.energ_hi), a1.get_y())
a2 = np.interp(engrid, 0.5 * (a2.energ_lo + a2.energ_hi) , a2.get_y())


fig, axes = plt.subplots(2, 1, sharex=True, figsize=(8, 6))
plt.sca(axes[0])
ui.plot_data("special extraction", xlog=True, ylog=True, clearwindow=False)
ui.plot_model("special extraction", overplot=True, label="special extraction")
ui.plot_model("specextract", overplot=True, label="specextract")
axes[0].legend()
axes[0].set_title("Simulated data with fits using different ARFs and RMFs")

axes[1].loglog(engrid, a1, label="special extraction", c="C1")
axes[1].loglog(engrid, a2, label="specextract", c="C2")
axes[1].set_xlabel("Energy (keV)")
axes[1].set_ylabel("Effective Area (cm$^2$)")
axes[1].legend()

# Better setting to display labels on a log axis
axes[1].get_xaxis().set_minor_formatter(mpl.ticker.LogFormatterMathtext(labelOnlyBase=False,
                                                               minor_thresholds=(2, .5)))
axes[1].tick_params(axis='x', labelsize=mpl.rcParams['xtick.labelsize'], which='both')
../_images/3a8fe9f47d52b84d5fbfd582b3b70133d47c344bd65fa0ccbb8a95893fe8e0f8.png

top: Simulated data and two model fits, assuming different ARFs and RMFs. “Special extraction” uses special qualifiers to match marx simulations, while “specextract” uses standard CIAO tools without any special qualifiers.

bottom: Shape of the ARFs used in the two extractions. Since the ARF is different, the model values are different when fitting to the same data (e.g. where the ARF is lower, the model flux has to be higher to match the data, leading to different slopes and normalizations in the two fits).

Input parameters

Click to hide code

 
my_src
<BinaryOpModel model instance 'xsvapec.apec1 + xsvapec.apec2'>

Fit with special extraction

Click to hide code

 
conf_special
<confidence results instance>

Fit with CIAO default extraction

Click to hide code

 
conf_specextract
<confidence results instance>

Again we are looking at an ACIS-I simulation and the differences between the two extraction methods are smaller than for ACIS-S. While the Ne abundance is not well-constrained in this case, the fitted values are in general in agreement with the input values.