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.
# Load marx in local environment and set up pfiles for tests
from utils import setup_marx
setup_marx()
Absorbed powerlaw on ACIS-S¶
This test checks the internal consistency of a marx spectral simulation by simulating a source placed on a back-illuminated chip of ACIS-S. The input spectrum is an absorbed powerlaw.
'''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:
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.We use the specextract tool without any special qualifiers. While
specextracthas 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.
'''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)
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()
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')
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¶
my_src
Model
| Component | Parameter | Thawed | Value | Min | Max | Units |
|---|---|---|---|---|---|---|
| xsphabs.a | nH | 1.0 | 0.0 | 1000000.0 | 10^22 atoms / cm^2 | |
| xspowerlaw.p | PhoIndex | 1.8 | -3.0 | 10.0 | ||
| norm | 0.001 | 0.0 | 1e+24 |
Fit with special extraction¶
conf_special
confidence 1σ (68.2689%) bounds
| Parameter | Best-fit value | Lower Bound | Upper Bound |
|---|---|---|---|
| a1.nH | 0.96852 | -0.0488344 | 0.0513345 |
| p1.PhoIndex | 1.80204 | -0.0573022 | 0.0585523 |
| p1.norm | 0.000932026 | -6.62985e-05 | 7.28197e-05 |
Summary (3)
Fit with CIAO default extraction¶
conf_specextract
confidence 1σ (68.2689%) bounds
| Parameter | Best-fit value | Lower Bound | Upper Bound |
|---|---|---|---|
| a2.nH | 0.897425 | -0.0465747 | 0.0487622 |
| p2.PhoIndex | 1.66154 | -0.0552411 | 0.0564912 |
| p2.norm | 0.000833138 | -5.76134e-05 | 6.27658e-05 |
Summary (3)
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.
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)
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()
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')
Input parameters¶
(some as above, but printed here again)
my_src
Model
| Component | Parameter | Thawed | Value | Min | Max | Units |
|---|---|---|---|---|---|---|
| xsphabs.a | nH | 1.0 | 0.0 | 1000000.0 | 10^22 atoms / cm^2 | |
| xspowerlaw.p | PhoIndex | 1.8 | -3.0 | 10.0 | ||
| norm | 0.001 | 0.0 | 1e+24 |
Fit with special extraction¶
conf_special
confidence 1σ (68.2689%) bounds
| Parameter | Best-fit value | Lower Bound | Upper Bound |
|---|---|---|---|
| a1.nH | 0.901507 | -0.0527481 | 0.0543107 |
| p1.PhoIndex | 1.73311 | -0.0528 | 0.0540501 |
| p1.norm | 0.000852955 | -5.91562e-05 | 6.44466e-05 |
Summary (3)
Fit with CIAO default extraction¶
conf_specextract
confidence 1σ (68.2689%) bounds
| Parameter | Best-fit value | Lower Bound | Upper Bound |
|---|---|---|---|
| a2.nH | 0.900009 | -0.0526273 | 0.0545024 |
| p2.PhoIndex | 1.73017 | -0.0530052 | 0.0539428 |
| p2.norm | 0.000870168 | -6.02136e-05 | 6.55985e-05 |
Summary (3)
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.
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)
Reading APEC data from 3.0.9
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-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)
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()
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')
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¶
my_src
Model
| Component | Parameter | Thawed | Value | Min | Max | Units |
|---|---|---|---|---|---|---|
| xsvapec.apec1 | kT | 0.7 | 0.0808 | 68.447 | keV | |
| He | 1.0 | 0.0 | 1000.0 | |||
| C | 1.0 | 0.0 | 1000.0 | |||
| N | 1.0 | 0.0 | 1000.0 | |||
| O | 1.0 | 0.0 | 1000.0 | |||
| Ne | 2.0 | 0.0 | 1000.0 | |||
| Mg | 1.0 | 0.0 | 1000.0 | |||
| Al | 1.0 | 0.0 | 1000.0 | |||
| Si | 1.0 | 0.0 | 1000.0 | |||
| S | 1.0 | 0.0 | 1000.0 | |||
| Ar | 1.0 | 0.0 | 1000.0 | |||
| Ca | 1.0 | 0.0 | 1000.0 | |||
| Fe | 1.0 | 0.0 | 1000.0 | |||
| Ni | 1.0 | 0.0 | 1000.0 | |||
| Redshift | 0.0 | -0.999 | 10.0 | |||
| norm | 0.0002 | 0.0 | 1e+24 | |||
| xsvapec.apec2 | kT | 2.0 | 0.0808 | 68.447 | keV | |
| He | 1.0 | 0.0 | 1000.0 | |||
| C | 1.0 | 0.0 | 1000.0 | |||
| N | 1.0 | 0.0 | 1000.0 | |||
| O | 1.0 | 0.0 | 1000.0 | |||
| Ne | linked | 2.0 | ⇐ apec1.Ne | |||
| Mg | 1.0 | 0.0 | 1000.0 | |||
| Al | 1.0 | 0.0 | 1000.0 | |||
| Si | 1.0 | 0.0 | 1000.0 | |||
| S | 1.0 | 0.0 | 1000.0 | |||
| Ar | 1.0 | 0.0 | 1000.0 | |||
| Ca | 1.0 | 0.0 | 1000.0 | |||
| Fe | 1.0 | 0.0 | 1000.0 | |||
| Ni | 1.0 | 0.0 | 1000.0 | |||
| Redshift | 0.0 | -0.999 | 10.0 | |||
| norm | 0.0004 | 0.0 | 1e+24 | |||
Fit with special extraction¶
conf_special
confidence 1σ (68.2689%) bounds
| Parameter | Best-fit value | Lower Bound | Upper Bound |
|---|---|---|---|
| apec1_1.kT | 0.728632 | -0.0292126 | 0.0304627 |
| apec1_1.Ne | 1.07123 | ----- | 1.16581 |
| apec1_1.norm | 0.000219788 | -1.22787e-05 | 1.34332e-05 |
| apec2_1.kT | 1.96796 | -0.126638 | 0.147016 |
| apec2_1.norm | 0.000381227 | -2.44487e-05 | 2.50262e-05 |
Summary (3)
Fit with CIAO default extraction¶
conf_specextract
confidence 1σ (68.2689%) bounds
| Parameter | Best-fit value | Lower Bound | Upper Bound |
|---|---|---|---|
| apec1_2.kT | 0.727378 | -0.0282397 | 0.0298023 |
| apec1_2.Ne | 1.18932 | ----- | 1.13253 |
| apec1_2.norm | 0.000228965 | -1.2408e-05 | 1.3536e-05 |
| apec2_2.kT | 1.98416 | -0.132099 | 0.152265 |
| apec2_2.norm | 0.000399334 | -2.6079e-05 | 2.65234e-05 |
Summary (3)
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.
