pyXSIM Example

To show how to make a set of photons from a 3D dataset using pyXSIM and yt for reading into SOXS, we’ll look at is that of thermal emission from a galaxy cluster. In this case, the gas in the core of the cluster is “sloshing” in the center, producing spiral-shaped cold fronts. The dataset we want to use for this example is available for download from the yt Project at this link.

First, import our necessary modules:

[1]:
import yt
import pyxsim
import soxs

Next, we load the dataset with yt (this dataset does not have species fields, so we specify that the gas is fully ionized in this case so that the emission measure field can be computed correctly):

[2]:
ds = yt.load(
    "GasSloshing/sloshing_nomag2_hdf5_plt_cnt_0150", default_species_fields="ionized"
)
yt : [INFO     ] 2024-03-06 15:53:44,519 Parameters: current_time              = 1.1835090993823291e+17
yt : [INFO     ] 2024-03-06 15:53:44,519 Parameters: domain_dimensions         = [16 16 16]
yt : [INFO     ] 2024-03-06 15:53:44,520 Parameters: domain_left_edge          = [-3.70272e+24 -3.70272e+24 -3.70272e+24]
yt : [INFO     ] 2024-03-06 15:53:44,520 Parameters: domain_right_edge         = [3.70272e+24 3.70272e+24 3.70272e+24]
yt : [INFO     ] 2024-03-06 15:53:44,521 Parameters: cosmological_simulation   = 0

Let’s use yt to take a slice of density and temperature through the center of the dataset so we can see what we’re looking at:

[3]:
slc = yt.SlicePlot(
    ds, "z", [("gas", "density"), ("gas", "temperature")], width=(1.0, "Mpc")
)
slc.show()
yt : [INFO     ] 2024-03-06 15:53:44,972 xlim = -1542838790481162406985728.000000 1542838790481162406985728.000000
yt : [INFO     ] 2024-03-06 15:53:44,972 ylim = -1542838790481162406985728.000000 1542838790481162406985728.000000
yt : [INFO     ] 2024-03-06 15:53:44,974 xlim = -1542838790481162406985728.000000 1542838790481162406985728.000000
yt : [INFO     ] 2024-03-06 15:53:44,974 ylim = -1542838790481162406985728.000000 1542838790481162406985728.000000
yt : [INFO     ] 2024-03-06 15:53:44,978 Making a fixed resolution buffer of (('gas', 'density')) 800 by 800
yt : [INFO     ] 2024-03-06 15:53:45,235 Making a fixed resolution buffer of (('gas', 'temperature')) 800 by 800


Ok, sloshing gas as advertised. Next, we’ll create a sphere object to serve as a source for the photons. Place it at the center of the domain with "c", and use a radius of 500 kpc:

[4]:
sp = ds.sphere("c", (0.5, "Mpc"))

Now, we need to set up the emission model for our source. We said we were going to look at the thermal emission from the hot plasma, so we’ll do that here by using ThermalSourceModel. The first four arguments are the name of the underlying spectral model, the maximum and minimum energies, and the number of bins in the spectrum. We’ve chosen these numbers so that the spectrum has an energy resolution of about 1 eV. Setting thermal_broad=True turns on thermal broadening. This simulation does not include metallicity, so we’ll do something simple and say that it uses the above spectral model and the metallicity is a constant \(Z = 0.3~Z_\odot\):

[5]:
source_model = pyxsim.CIESourceModel(
    "apec", 0.5, 9.0, 9000, thermal_broad=True, Zmet=0.3
)
pyxsim : [INFO     ] 2024-03-06 15:53:46,022 kT_min = 0.025 keV
pyxsim : [INFO     ] 2024-03-06 15:53:46,022 kT_max = 64 keV

We’re almost ready to go to generate the photons from this source, but first we should decide what our redshift, collecting area, and exposure time should be. Let’s pick big numbers, because remember the point of this first step is to create a Monte-Carlo sample from which to draw smaller sub-samples for mock observations. Note these are all (value, unit) tuples:

[6]:
exp_time = (300.0, "ks")  # exposure time
area = (3.0, "m**2")  # collecting area
redshift = 0.2

So, that’s everything–let’s create the photons!

[7]:
n_photons, n_cells = pyxsim.make_photons(
    "my_photons", sp, redshift, area, exp_time, source_model
)
pyxsim : [INFO     ] 2024-03-06 15:53:46,034 Cosmology: h = 0.71, omega_matter = 0.27, omega_lambda = 0.73
pyxsim : [INFO     ] 2024-03-06 15:53:46,034 Using emission measure field '('gas', 'emission_measure')'.
pyxsim : [INFO     ] 2024-03-06 15:53:46,034 Using temperature field '('gas', 'temperature')'.
pyxsim : [INFO     ] 2024-03-06 16:00:39,974 Finished generating photons.
pyxsim : [INFO     ] 2024-03-06 16:00:39,978 Number of photons generated: 21427563
pyxsim : [INFO     ] 2024-03-06 16:00:39,980 Number of cells with photons: 4054508

Ok, that was easy. Now we have a photon list that we can use to create events, using the project_photons() function. To be realistic, we’re going to want to assume foreground Galactic absorption, using the “TBabs” absorption model and assuming a foreground absorption column of \(N_H = 4 \times 10^{20}~{\rm cm}^{-2}\). Here we’ll just do a simple projection along the z-axis, reducing the exposure time, and centering the photons at RA, Dec = (30, 45) degrees:

[8]:
n_events = pyxsim.project_photons(
    "my_photons", "my_events", "z", (30.0, 45.0), absorb_model="tbabs", nH=0.04
)
pyxsim : [INFO     ] 2024-03-06 16:00:40,002 Foreground galactic absorption: using the tbabs model and nH = 0.04.
pyxsim : [INFO     ] 2024-03-06 16:00:48,470 Detected 19261921 events.

We can then use this event list that we wrote as an input to the instrument simulator in SOXS. We’ll use a smaller exposure time (100 ks instead of 500 ks), and observe it with the Lynx calorimeter:

[9]:
soxs.instrument_simulator(
    "my_events.h5",
    "evt.fits",
    (100.0, "ks"),
    "lynx_lxm",
    [30.0, 45.0],
    overwrite=True,
)
soxs : [INFO     ] 2024-03-06 16:00:48,479 Making observation of source in evt.fits.
soxs : [INFO     ] 2024-03-06 16:00:49,143 Detecting events from source my_events.
soxs : [INFO     ] 2024-03-06 16:00:49,144 Applying energy-dependent effective area from xrs_mucal_3x10_3.0eV.arf.
soxs : [INFO     ] 2024-03-06 16:00:51,347 Pixeling events.
soxs : [INFO     ] 2024-03-06 16:00:51,761 Scattering events with a image-based PSF.
soxs : [INFO     ] 2024-03-06 16:00:52,067 2632075 events were detected from the source.
soxs : [INFO     ] 2024-03-06 16:00:52,151 Scattering energies with RMF xrs_mucal_3.0eV.rmf.
soxs : [INFO     ] 2024-03-06 16:00:57,598 Adding background events.
soxs : [INFO     ] 2024-03-06 16:00:57,689 Adding in point-source background.
soxs : [INFO     ] 2024-03-06 16:00:57,910 Detecting events from source ptsrc_bkgnd.
soxs : [INFO     ] 2024-03-06 16:00:57,910 Applying energy-dependent effective area from xrs_mucal_3x10_3.0eV.arf.
soxs : [INFO     ] 2024-03-06 16:00:57,921 Pixeling events.
soxs : [INFO     ] 2024-03-06 16:00:57,927 Scattering events with a image-based PSF.
soxs : [INFO     ] 2024-03-06 16:00:57,939 40602 events were detected from the source.
soxs : [INFO     ] 2024-03-06 16:00:57,940 Scattering energies with RMF xrs_mucal_3.0eV.rmf.
soxs : [INFO     ] 2024-03-06 16:00:58,567 Generated 40602 photons from the point-source background.
soxs : [INFO     ] 2024-03-06 16:00:58,567 Adding in astrophysical foreground.
soxs : [INFO     ] 2024-03-06 16:01:15,514 Adding in instrumental background.
soxs : [INFO     ] 2024-03-06 16:01:15,565 Making 71868 events from the galactic foreground.
soxs : [INFO     ] 2024-03-06 16:01:15,566 Making 0 events from the instrumental background.
soxs : [INFO     ] 2024-03-06 16:01:15,598 Writing events to file evt.fits.
soxs : [INFO     ] 2024-03-06 16:01:16,070 Observation complete.

We can use the write_image() function in SOXS to bin the events into an image and write them to a file, restricting the energies between 0.5 and 2.0 keV:

[10]:
soxs.write_image("evt.fits", "img.fits", emin=0.5, emax=2.0, overwrite=True)

We can show the resulting image:

[11]:
fig, ax = soxs.plot_image(
    "img.fits", stretch="sqrt", cmap="afmhot", vmax=1000.0, width=0.05
)
../_images/cookbook_pyXSIM_Example_22_0.png

We can also bin the events into a spectrum using write_spectrum() and write the spectrum to disk:

[12]:
soxs.write_spectrum("evt.fits", "evt.pha", overwrite=True)

and plot the spectrum using plot_spectrum():

[13]:
fig, ax = soxs.plot_spectrum("evt.pha", xmin=0.5, xmax=7.0)
../_images/cookbook_pyXSIM_Example_26_0.png

Let’s zoom into the region of the spectrum around the iron line to look at the detailed structure afforded by the resolution of the calorimeter:

[14]:
ax.set_xlim(5.4, 5.7)
fig
[14]:
../_images/cookbook_pyXSIM_Example_28_0.png