More Advanced Thermal Emission¶
In this example, we’ll look at the emission from a disk galaxy from the Illustris TNG simulations. This dataset has metallicity information for several species in it. We’ll make a cut in phase space like we did in the previous example. The dataset we want to use for this example is available for download here.
First, import our necessary modules:
[1]:
import yt
import pyxsim
import soxs
We will make phase space cuts on the gas cells using density, temperature, and star formation rate:
[2]:
# Note that the units of all numbers in this function are CGS
# define hot gas filter
def hot_gas(pfilter, data):
pfilter1 = data[pfilter.filtered_type, "temperature"] > 3.0e5
pfilter2 = data["PartType0", "StarFormationRate"] == 0.0
pfilter3 = data[pfilter.filtered_type, "density"] < 5e-25
return pfilter1 & pfilter2 & pfilter3
yt.add_particle_filter(
"hot_gas",
function=hot_gas,
filtered_type="gas",
requires=["temperature", "density"],
)
Next, we load
the dataset with yt, and add the "hot_gas"
filter to the dataset:
[3]:
ds = yt.load(
"cutout_37.hdf5", bounding_box=[[-1000.0, 1000], [-1000.0, 1000], [-1000.0, 1000]]
)
ds.add_particle_filter("hot_gas")
yt : [INFO ] 2023-04-14 13:08:39,374 Calculating time from 1.000e+00 to be 4.356e+17 seconds
yt : [INFO ] 2023-04-14 13:08:39,410 Parameters: current_time = 4.355810528213311e+17 s
yt : [INFO ] 2023-04-14 13:08:39,410 Parameters: domain_dimensions = [1 1 1]
yt : [INFO ] 2023-04-14 13:08:39,410 Parameters: domain_left_edge = [-1000. -1000. -1000.]
yt : [INFO ] 2023-04-14 13:08:39,410 Parameters: domain_right_edge = [1000. 1000. 1000.]
yt : [INFO ] 2023-04-14 13:08:39,411 Parameters: cosmological_simulation = 1
yt : [INFO ] 2023-04-14 13:08:39,411 Parameters: current_redshift = 2.220446049250313e-16
yt : [INFO ] 2023-04-14 13:08:39,411 Parameters: omega_lambda = 0.6911
yt : [INFO ] 2023-04-14 13:08:39,411 Parameters: omega_matter = 0.3089
yt : [INFO ] 2023-04-14 13:08:39,411 Parameters: omega_radiation = 0.0
yt : [INFO ] 2023-04-14 13:08:39,411 Parameters: hubble_constant = 0.6774
yt : [WARNING ] 2023-04-14 13:08:39,412 A bounding box was explicitly specified, so we are disabling periodicity.
yt : [INFO ] 2023-04-14 13:08:39,421 Allocating for 2.761e+06 particles
Loading particle index: 100%|█████████████████| 11/11 [00:00<00:00, 1577.08it/s]
[3]:
True
We also need to tell pyXSIM which elements have fields in the dataset that should be used. To do this we create a var_elem
dictionary of (key, value) pairs corresponding to the element name and the yt field name (assuming the "hot_gas"
type).
[4]:
# metal fields to use
metals = [
"C_fraction",
"N_fraction",
"O_fraction",
"Ne_fraction",
"Mg_fraction",
"Si_fraction",
"Fe_fraction",
]
var_elem = {elem.split("_")[0]: ("hot_gas", elem) for elem in metals}
Now that we have everything we need, we’ll set up the IGMSourceModel
, which is based on Cloudy and includes resonant scattering off of the CXB (see here for more details). Because we created a hot gas filter, we will use the "hot_gas"
field type for the emission measure, temperature, and metallicity fields.
[5]:
source_model = pyxsim.IGMSourceModel(
0.1,
4.0,
5000,
("hot_gas", "metallicity"),
binscale="log",
resonant_scattering=True,
temperature_field=("hot_gas", "temperature"),
emission_measure_field=("hot_gas", "emission_measure"),
nh_field=("hot_gas", "H_nuclei_density"),
var_elem=var_elem,
)
pyxsim : [INFO ] 2023-04-14 13:08:39,874 kT_min = 0.00431 keV
pyxsim : [INFO ] 2023-04-14 13:08:39,874 kT_max = 64 keV
As in other examples, we choose big numbers for the collecting area and exposure time, and a redshift:
[6]:
exp_time = (1.0, "Ms") # exposure time
area = (5000.0, "cm**2") # collecting area
redshift = 0.01
Next, we’ll create a box object to serve as a source for the photons. The dataset consists of only the galaxy at a specific location, which we use below, and pick a width of 1 Mpc:
[7]:
c = ds.arr([0.0, 0.0, 0.0], "code_length")
width = ds.quan(1.0, "Mpc")
le = c - 0.5 * width
re = c + 0.5 * width
box = ds.box(le, re)
So, that’s everything–let’s create the photons! We use the make_photons
function for this:
[8]:
n_photons, n_cells = pyxsim.make_photons(
"cutout_37_photons", box, redshift, area, exp_time, source_model
)
pyxsim : [INFO ] 2023-04-14 13:08:39,886 Cosmology: h = 0.6774, omega_matter = 0.3089, omega_lambda = 0.6911
pyxsim : [INFO ] 2023-04-14 13:08:39,887 Using emission measure field '('hot_gas', 'emission_measure')'.
pyxsim : [INFO ] 2023-04-14 13:08:39,887 Using temperature field '('hot_gas', 'temperature')'.
pyxsim : [INFO ] 2023-04-14 13:08:39,887 Using nH field '('hot_gas', 'H_nuclei_density')'.
---------------------------------------------------------------------------
BlockingIOError Traceback (most recent call last)
Cell In[8], line 1
----> 1 n_photons, n_cells = pyxsim.make_photons(
2 "cutout_37_photons", box, redshift, area, exp_time, source_model
3 )
File ~/Source/pyxsim/pyxsim/photon_list.py:338, in make_photons(photon_prefix, data_source, redshift, area, exp_time, source_model, point_sources, parameters, center, dist, cosmology, velocity_fields, bulk_velocity, observer, fields_to_keep)
334 parameters["data_type"] = "particles"
336 source_model.set_pv(p_fields, v_fields, le, re, dw, c, ds.periodicity, observer)
--> 338 f = h5py.File(photon_file, "w")
340 # Info
341 info = f.create_group("info")
File ~/mambaforge/envs/py311/lib/python3.11/site-packages/h5py/_hl/files.py:567, in File.__init__(self, name, mode, driver, libver, userblock_size, swmr, rdcc_nslots, rdcc_nbytes, rdcc_w0, track_order, fs_strategy, fs_persist, fs_threshold, fs_page_size, page_buf_size, min_meta_keep, min_raw_keep, locking, alignment_threshold, alignment_interval, meta_block_size, **kwds)
558 fapl = make_fapl(driver, libver, rdcc_nslots, rdcc_nbytes, rdcc_w0,
559 locking, page_buf_size, min_meta_keep, min_raw_keep,
560 alignment_threshold=alignment_threshold,
561 alignment_interval=alignment_interval,
562 meta_block_size=meta_block_size,
563 **kwds)
564 fcpl = make_fcpl(track_order=track_order, fs_strategy=fs_strategy,
565 fs_persist=fs_persist, fs_threshold=fs_threshold,
566 fs_page_size=fs_page_size)
--> 567 fid = make_fid(name, mode, userblock_size, fapl, fcpl, swmr=swmr)
569 if isinstance(libver, tuple):
570 self._libver = libver
File ~/mambaforge/envs/py311/lib/python3.11/site-packages/h5py/_hl/files.py:237, in make_fid(name, mode, userblock_size, fapl, fcpl, swmr)
235 fid = h5f.create(name, h5f.ACC_EXCL, fapl=fapl, fcpl=fcpl)
236 elif mode == 'w':
--> 237 fid = h5f.create(name, h5f.ACC_TRUNC, fapl=fapl, fcpl=fcpl)
238 elif mode == 'a':
239 # Open in append mode (read/write).
240 # If that fails, create a new file only if it won't clobber an
241 # existing one (ACC_EXCL)
242 try:
File h5py/_objects.pyx:54, in h5py._objects.with_phil.wrapper()
File h5py/_objects.pyx:55, in h5py._objects.with_phil.wrapper()
File h5py/h5f.pyx:126, in h5py.h5f.create()
BlockingIOError: [Errno 35] Unable to synchronously create file (unable to lock file, errno = 35, error message = 'Resource temporarily unavailable')
And now we create events using the project_photons
function. Let’s project along the "z"
axis. We’ll use the "tbabs"
foreground absorption model this time, with a neutral hydrogen column of \(N_H = 2 \times 10^{20}~{\rm cm}^{-2}\):
[9]:
n_events = pyxsim.project_photons(
"cutout_37_photons",
"cutout_37_events",
"x",
(30.0, 45.0),
absorb_model="tbabs",
nH=0.02,
)
pyxsim : [INFO ] 2023-04-14 13:08:55,114 Foreground galactic absorption: using the tbabs model and nH = 0.02.
---------------------------------------------------------------------------
OSError Traceback (most recent call last)
Cell In[9], line 1
----> 1 n_events = pyxsim.project_photons(
2 "cutout_37_photons",
3 "cutout_37_events",
4 "x",
5 (30.0, 45.0),
6 absorb_model="tbabs",
7 nH=0.02,
8 )
File ~/Source/pyxsim/pyxsim/photon_list.py:873, in project_photons(photon_prefix, event_prefix, normal, sky_center, absorb_model, nH, abund_table, no_shifting, north_vector, sigma_pos, flat_sky, kernel, save_los, prng)
773 def project_photons(
774 photon_prefix,
775 event_prefix,
(...)
787 prng=None,
788 ):
789 r"""
790 Projects photons onto an image plane given a line of sight, and
791 stores them in an HDF5 dataset which contains an event list.
(...)
871 ... nH=0.04)
872 """
--> 873 return _project_photons(
874 "external",
875 photon_prefix,
876 event_prefix,
877 normal,
878 sky_center,
879 absorb_model=absorb_model,
880 nH=nH,
881 abund_table=abund_table,
882 no_shifting=no_shifting,
883 north_vector=north_vector,
884 sigma_pos=sigma_pos,
885 flat_sky=flat_sky,
886 kernel=kernel,
887 save_los=save_los,
888 prng=prng,
889 )
File ~/Source/pyxsim/pyxsim/photon_list.py:544, in _project_photons(obs, photon_prefix, event_prefix, normal, sky_center, absorb_model, nH, abund_table, no_shifting, north_vector, flat_sky, sigma_pos, kernel, save_los, prng)
541 if nH is None:
542 nH = 0.0
--> 544 f = h5py.File(photon_file, "r")
546 p = f["parameters"]
548 data_type = p["data_type"].asstr()[()]
File ~/mambaforge/envs/py311/lib/python3.11/site-packages/h5py/_hl/files.py:567, in File.__init__(self, name, mode, driver, libver, userblock_size, swmr, rdcc_nslots, rdcc_nbytes, rdcc_w0, track_order, fs_strategy, fs_persist, fs_threshold, fs_page_size, page_buf_size, min_meta_keep, min_raw_keep, locking, alignment_threshold, alignment_interval, meta_block_size, **kwds)
558 fapl = make_fapl(driver, libver, rdcc_nslots, rdcc_nbytes, rdcc_w0,
559 locking, page_buf_size, min_meta_keep, min_raw_keep,
560 alignment_threshold=alignment_threshold,
561 alignment_interval=alignment_interval,
562 meta_block_size=meta_block_size,
563 **kwds)
564 fcpl = make_fcpl(track_order=track_order, fs_strategy=fs_strategy,
565 fs_persist=fs_persist, fs_threshold=fs_threshold,
566 fs_page_size=fs_page_size)
--> 567 fid = make_fid(name, mode, userblock_size, fapl, fcpl, swmr=swmr)
569 if isinstance(libver, tuple):
570 self._libver = libver
File ~/mambaforge/envs/py311/lib/python3.11/site-packages/h5py/_hl/files.py:231, in make_fid(name, mode, userblock_size, fapl, fcpl, swmr)
229 if swmr and swmr_support:
230 flags |= h5f.ACC_SWMR_READ
--> 231 fid = h5f.open(name, flags, fapl=fapl)
232 elif mode == 'r+':
233 fid = h5f.open(name, h5f.ACC_RDWR, fapl=fapl)
File h5py/_objects.pyx:54, in h5py._objects.with_phil.wrapper()
File h5py/_objects.pyx:55, in h5py._objects.with_phil.wrapper()
File h5py/h5f.pyx:106, in h5py.h5f.open()
OSError: Unable to synchronously open file (file signature not found)
Now that we have a set of “events” on the sky, we can read them in and write them to a SIMPUT file:
[10]:
events = pyxsim.EventList("cutout_37_events.h5")
events.write_to_simput("cutout_37", overwrite=True)
---------------------------------------------------------------------------
BlockingIOError Traceback (most recent call last)
Cell In[10], line 1
----> 1 events = pyxsim.EventList("cutout_37_events.h5")
2 events.write_to_simput("cutout_37", overwrite=True)
File ~/Source/pyxsim/pyxsim/event_list.py:43, in EventList.__init__(self, filespec)
41 self.num_events = []
42 for i, fn in enumerate(self.filenames):
---> 43 with h5py.File(fn, "r") as f:
44 p = f["parameters"]
45 info = f["info"]
File ~/mambaforge/envs/py311/lib/python3.11/site-packages/h5py/_hl/files.py:567, in File.__init__(self, name, mode, driver, libver, userblock_size, swmr, rdcc_nslots, rdcc_nbytes, rdcc_w0, track_order, fs_strategy, fs_persist, fs_threshold, fs_page_size, page_buf_size, min_meta_keep, min_raw_keep, locking, alignment_threshold, alignment_interval, meta_block_size, **kwds)
558 fapl = make_fapl(driver, libver, rdcc_nslots, rdcc_nbytes, rdcc_w0,
559 locking, page_buf_size, min_meta_keep, min_raw_keep,
560 alignment_threshold=alignment_threshold,
561 alignment_interval=alignment_interval,
562 meta_block_size=meta_block_size,
563 **kwds)
564 fcpl = make_fcpl(track_order=track_order, fs_strategy=fs_strategy,
565 fs_persist=fs_persist, fs_threshold=fs_threshold,
566 fs_page_size=fs_page_size)
--> 567 fid = make_fid(name, mode, userblock_size, fapl, fcpl, swmr=swmr)
569 if isinstance(libver, tuple):
570 self._libver = libver
File ~/mambaforge/envs/py311/lib/python3.11/site-packages/h5py/_hl/files.py:231, in make_fid(name, mode, userblock_size, fapl, fcpl, swmr)
229 if swmr and swmr_support:
230 flags |= h5f.ACC_SWMR_READ
--> 231 fid = h5f.open(name, flags, fapl=fapl)
232 elif mode == 'r+':
233 fid = h5f.open(name, h5f.ACC_RDWR, fapl=fapl)
File h5py/_objects.pyx:54, in h5py._objects.with_phil.wrapper()
File h5py/_objects.pyx:55, in h5py._objects.with_phil.wrapper()
File h5py/h5f.pyx:106, in h5py.h5f.open()
BlockingIOError: [Errno 35] Unable to synchronously open file (unable to lock file, errno = 35, error message = 'Resource temporarily unavailable')
We can then use this SIMPUT file as an input to the instrument simulator in SOXS. We’ll observe it with the 2eV LEM model for 1 Ms. First, we’ll create a background file that we’ll use for the background:
[11]:
soxs.make_background_file(
"bkgnd_evt_37.fits", (1000.0, "ks"), "lem_2eV", [30.0, 45.0], overwrite=True
)
soxs : [INFO ] 2023-04-14 13:08:55,394 Adding in point-source background.
soxs : [INFO ] 2023-04-14 13:09:02,363 Detecting events from source ptsrc_bkgnd.
soxs : [INFO ] 2023-04-14 13:09:02,363 Applying energy-dependent effective area from lem_300522.arf.
soxs : [INFO ] 2023-04-14 13:09:03,158 Pixeling events.
soxs : [INFO ] 2023-04-14 13:09:03,635 Scattering events with a gaussian-based PSF.
soxs : [INFO ] 2023-04-14 13:09:03,833 3458745 events were detected from the source.
soxs : [INFO ] 2023-04-14 13:09:03,930 Scattering energies with RMF lem_2ev_110422.rmf.
soxs : [INFO ] 2023-04-14 13:09:19,656 Generated 3458745 photons from the point-source background.
soxs : [INFO ] 2023-04-14 13:09:19,656 Adding in astrophysical foreground.
soxs : [INFO ] 2023-04-14 13:09:25,327 Adding in instrumental background.
soxs : [INFO ] 2023-04-14 13:09:25,610 Making 6080180 events from the galactic foreground.
soxs : [INFO ] 2023-04-14 13:09:25,611 Making 3414820 events from the instrumental background.
soxs : [INFO ] 2023-04-14 13:09:26,059 Writing events to file bkgnd_evt_37.fits.
Now we simulate the source itself, adding in the background:
[12]:
soxs.instrument_simulator(
"cutout_37_simput.fits",
"evt_37.fits",
(1000.0, "ks"),
"lem_2eV",
[30.0, 45.0],
overwrite=True,
bkgnd_file="bkgnd_evt_37.fits",
)
soxs : [INFO ] 2023-04-14 13:09:27,826 Making observation of source in evt_37.fits.
soxs : [INFO ] 2023-04-14 13:09:28,021 Detecting events from source cutout_37.
soxs : [INFO ] 2023-04-14 13:09:28,022 Applying energy-dependent effective area from lem_300522.arf.
soxs : [INFO ] 2023-04-14 13:09:28,050 Pixeling events.
soxs : [INFO ] 2023-04-14 13:09:28,064 Scattering events with a gaussian-based PSF.
soxs : [INFO ] 2023-04-14 13:09:28,068 79382 events were detected from the source.
soxs : [INFO ] 2023-04-14 13:09:28,072 Scattering energies with RMF lem_2ev_110422.rmf.
soxs : [INFO ] 2023-04-14 13:09:28,774 Adding background events from the file bkgnd_evt_37.fits.
soxs : [INFO ] 2023-04-14 13:09:29,583 Adding 12953745 background events from bkgnd_evt_37.fits.
soxs : [INFO ] 2023-04-14 13:09:29,982 Writing events to file evt_37.fits.
soxs : [INFO ] 2023-04-14 13:09:31,785 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.644 and 0.65 keV, which focuses on the redshifted OVIII line:
[13]:
soxs.write_image("evt_37.fits", "img_37.fits", emin=0.644, emax=0.65, overwrite=True)
Now we can take a quick look at the image:
[14]:
soxs.plot_image("img_37.fits", stretch="log", cmap="arbre", width=0.4, vmin=0.5)
[14]:
(<Figure size 1000x1000 with 2 Axes>, <WCSAxes: >)

Now we will make spectra to look at. First, filter the events of both the combined source and background files and the background-only files within 0.15 degree of the center:
[15]:
soxs.filter_events(
"evt_37.fits",
"evt_37_filter.fits",
overwrite=True,
region='fk5\ncircle(30.0000000,45.0000000,540.000")',
)
soxs.filter_events(
"bkgnd_evt_37.fits",
"bkgnd_evt_37_filter.fits",
overwrite=True,
region='fk5\ncircle(30.0000000,45.0000000,540.000")',
)
Now bin up spectra for these new event files:
[16]:
soxs.write_spectrum("evt_37_filter.fits", "evt_37.pi", overwrite=True)
soxs.write_spectrum("bkgnd_evt_37_filter.fits", "bkgnd_evt_37.pi", overwrite=True)
Finally, we can plot the spectra. Below, the total spectrum is in blue and the background/foreground spectrum is in orange. The lines from the emission of the distant galaxy are redshifted away from the foreground Milky Way lines.
[17]:
fig, ax = soxs.plot_spectrum("evt_37.pi", xmin=0.5, xmax=0.7, xscale="linear", ymin=0.5)
soxs.plot_spectrum("bkgnd_evt_37.pi", xmin=0.5, xmax=0.7, fig=fig, ax=ax, ymin=0.5)
[17]:
(<Figure size 1000x1000 with 1 Axes>,
<Axes: xlabel='Energy (keV)', ylabel='Count Rate (counts/s/keV)'>)
