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-03-06 15:50:38,747 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-03-06 15:50:38,757 Making observation of source in evt.fits.
soxs : [INFO     ] 2024-03-06 15:50:38,858 Detecting events from source beta_model.
soxs : [INFO     ] 2024-03-06 15:50:38,858 Applying energy-dependent effective area from xrs_hdxi_3x10.arf.
soxs : [INFO     ] 2024-03-06 15:50:38,936 Pixeling events.
soxs : [INFO     ] 2024-03-06 15:50:38,963 Scattering events with a image-based PSF.
soxs : [INFO     ] 2024-03-06 15:50:38,987 184402 events were detected from the source.
soxs : [INFO     ] 2024-03-06 15:50:38,993 Scattering energies with RMF xrs_hdxi.rmf.
soxs : [INFO     ] 2024-03-06 15:50:39,698 Adding background events.
soxs : [INFO     ] 2024-03-06 15:50:39,769 Adding in astrophysical foreground.
soxs : [INFO     ] 2024-03-06 15:50:51,085 Adding in instrumental background.
soxs : [INFO     ] 2024-03-06 15:50:51,571 Making 17079625 events from the galactic foreground.
soxs : [INFO     ] 2024-03-06 15:50:51,572 Making 203708 events from the instrumental background.
soxs : [INFO     ] 2024-03-06 15:50:52,496 Writing events to file evt.fits.
soxs : [INFO     ] 2024-03-06 15:50:54,483 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)')
../_images/cookbook_Radial_Profile_19_1.png

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)')
../_images/cookbook_Radial_Profile_21_1.png