Radial Profile¶
This example shows how to create a radial profile from a SOXS event file, including using an exposure map to get flux-based quantities. We’ll simulate a simple isothermal cluster.
[1]:
import matplotlib
matplotlib.rc("font", size=18)
import matplotlib.pyplot as plt
import soxs
import astropy.io.fits as pyfits
First, create the spectrum for the cluster using an absorbed thermal APEC model:
[2]:
emin = 0.05 # keV
emax = 20.0 # keV
nbins = 20000
agen = soxs.ApecGenerator(emin, emax, nbins)
[3]:
kT = 6.0
abund = 0.3
redshift = 0.05
norm = 1.0
spec = agen.get_spectrum(kT, abund, redshift, norm)
spec.rescale_flux(1.0e-13, emin=0.5, emax=2.0, flux_type="energy")
spec.apply_foreground_absorption(0.02)
And a spatial distribution based on a \(\beta\)-model:
[4]:
pos = soxs.BetaModel(30.0, 45.0, 50.0, 0.67)
Generate a SIMPUT catalog from these two models, and write it to a file:
[5]:
width = 10.0 # arcmin by default
nx = 1024 # resolution of image
cluster = soxs.SimputSpectrum.from_models("beta_model", spec, pos, width, nx)
cluster_cat = soxs.SimputCatalog.from_source(
"beta_model.simput", cluster, overwrite=True
)
soxs : [INFO ] 2024-06-19 16:01:58,897 Appending source 'beta_model' to beta_model.simput.
and run the instrument simulation (for simplicity we’ll turn off the point-source background):
[6]:
soxs.instrument_simulator(
"beta_model.simput",
"evt.fits",
(100.0, "ks"),
"lynx_hdxi",
[30.0, 45.0],
overwrite=True,
ptsrc_bkgnd=False,
)
soxs : [INFO ] 2024-06-19 16:01:58,907 Making observation of source in evt.fits.
soxs : [INFO ] 2024-06-19 16:01:59,001 Detecting events from source beta_model.
soxs : [INFO ] 2024-06-19 16:01:59,002 Applying energy-dependent effective area from xrs_hdxi_3x10.arf.
soxs : [INFO ] 2024-06-19 16:01:59,082 Pixeling events.
soxs : [INFO ] 2024-06-19 16:01:59,111 Scattering events with a image-based PSF.
soxs : [INFO ] 2024-06-19 16:01:59,136 184616 events were detected from the source.
soxs : [INFO ] 2024-06-19 16:01:59,142 Scattering energies with RMF xrs_hdxi.rmf.
soxs : [INFO ] 2024-06-19 16:01:59,845 Adding background events.
soxs : [INFO ] 2024-06-19 16:01:59,911 Adding in astrophysical foreground.
soxs : [INFO ] 2024-06-19 16:02:11,411 Adding in instrumental background.
soxs : [INFO ] 2024-06-19 16:02:11,997 Making 17083337 events from the galactic foreground.
soxs : [INFO ] 2024-06-19 16:02:11,998 Making 261668 events from the instrumental background.
soxs : [INFO ] 2024-06-19 16:02:13,057 Writing events to file evt.fits.
soxs : [INFO ] 2024-06-19 16:02:14,560 Observation complete.
Make an exposure map so that we can obtain flux-based quantities:
[7]:
soxs.make_exposure_map("evt.fits", "expmap.fits", 2.3, overwrite=True)
Make the radial profile, using energies between 0.5 and 5.0 keV, between radii of 0 and 200 arcseconds, with 50 bins:
[8]:
soxs.write_radial_profile(
"evt.fits",
"profile.fits",
[30.0, 45.0],
0,
200,
50,
emin=0.5,
emax=5.0,
expmap_file="expmap.fits",
overwrite=True,
)
Now we can use AstroPy’s FITS reader to open the profile and have a look at the columns that are inside:
[9]:
f = pyfits.open("profile.fits")
f["PROFILE"].columns
[9]:
ColDefs(
name = 'RLO'; format = 'D'; unit = 'arcsec'
name = 'RHI'; format = 'D'; unit = 'arcsec'
name = 'RMID'; format = 'D'; unit = 'arcsec'
name = 'AREA'; format = 'D'; unit = 'arcsec**2'
name = 'NET_COUNTS'; format = 'D'; unit = 'count'
name = 'NET_ERR'; format = 'D'; unit = 'count'
name = 'NET_RATE'; format = 'D'; unit = 'count/s'
name = 'ERR_RATE'; format = 'D'; unit = 'count/s'
name = 'SUR_BRI'; format = 'D'; unit = 'count/s/arcsec**2'
name = 'SUR_BRI_ERR'; format = '1D'; unit = 'count/s/arcsec**2'
name = 'MEAN_SRC_EXP'; format = 'D'; unit = 'cm**2'
name = 'NET_FLUX'; format = 'D'; unit = 'count/s/cm**2'
name = 'NET_FLUX_ERR'; format = 'D'; unit = 'count/s/cm**2'
name = 'SUR_FLUX'; format = 'D'; unit = 'count/s/cm**2/arcsec**2'
name = 'SUR_FLUX_ERR'; format = 'D'; unit = 'count/s/cm**2/arcsec**2'
)
and use Matplotlib to plot some quantities. We can plot the surface brightness:
[10]:
plt.figure(figsize=(8, 8))
plt.errorbar(
f["profile"].data["rmid"],
f["profile"].data["sur_bri"],
lw=2,
yerr=f["profile"].data["sur_bri_err"],
)
plt.xscale("log")
plt.yscale("log")
plt.xlabel("r (arcsec)")
plt.ylabel("S (cts/s/arcsec**2)")
[10]:
Text(0, 0.5, 'S (cts/s/arcsec**2)')
and, since we used an exposure map, the surface flux:
[11]:
plt.figure(figsize=(8, 8))
plt.errorbar(
f["profile"].data["rmid"],
f["profile"].data["sur_flux"],
lw=2,
yerr=f["profile"].data["sur_flux_err"],
)
plt.xscale("log")
plt.yscale("log")
plt.xlabel("r (arcsec)")
plt.ylabel("S (cts/s/cm**2/arcsec**2)")
[11]:
Text(0, 0.5, 'S (cts/s/cm**2/arcsec**2)')