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: >)
../_images/cookbook_More_Advanced_Thermal_Emission_28_1.png

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