MCRT Analysis

The thor.mcrt module (part of thor-rt) provides Python classes for loading and analyzing Monte Carlo Radiative Transfer output.

PhotonMCRT

Top-level entry point that loads a THOR config file and discovers all available photon datasets.

from thor.mcrt import PhotonMCRT

mcrt = PhotonMCRT("config.yaml")
print(mcrt.datasets.keys())  # e.g., dict_keys(['input', 'peel', 'original'])

Constructor

PhotonMCRT(outputpath: str, overwrite: bool = False)

Parameter	Description
`outputpath`	Path to the THOR configuration YAML file
`overwrite`	Overwrite cached dataset files

The constructor reads config["mcrtsimulation"]["outputpath"], then scans for input/, peel/, and original/ subdirectories. Each valid subdirectory becomes a PhotonMCRTDataset entry in datasets.

Attributes

Attribute	Type	Description
`config`	`dict`	Full parsed YAML configuration
`datasets`	`dict[str, PhotonMCRTDataset]`	Photon datasets keyed by type (`"input"`, `"peel"`, `"original"`)

PhotonMCRTDataset

Represents a single photon dataset (one of input, peel, or original). Inherits from scida's DatasetWithUnitMixin for lazy Dask-backed HDF5 loading with unit support.

ds = mcrt.datasets["peel"]
print(ds.photon_type)   # "peel"
print(ds.linename)      # "Lya"
print(ds.lambda0)       # 1215.67
print(ds.cosmological)  # True

Constructor

PhotonMCRTDataset(
    outputpath: str,
    overwrite: bool = False,
    config: str | dict | None = None,
    default_los: int = 0,
    photon_type: str | None = None,
    **kwargs
)

Parameter	Description
`outputpath`	Path to the photon type subdirectory (e.g., `output/peel`)
`overwrite`	Overwrite cached dataset files
`config`	YAML config path or dict. Enables unit conversion and cosmological features.
`default_los`	Index into `peeling.lines_of_sight` for the default line of sight
`photon_type`	One of `"input"`, `"peel"`, `"original"`. Auto-detected from directory name if omitted.

Note

When used through PhotonMCRT, the constructor is called automatically with the correct config and photon_type.

Key Attributes

Attribute	Type	Description
`photon_type`	`str`	`"input"`, `"peel"`, or `"original"`
`cosmological`	`bool`	Whether the simulation is cosmological
`redshift`	`float`	Redshift (0.0 for non-cosmological)
`hubble_flow`	`float`	Hubble flow in km/s/Mpc
`linename`	`str`	Emission line name (e.g., `"Lya"`, `"MgII"`)
`lambda0`	`float`	Rest-frame wavelength in Angstrom
`los`	`ndarray` or `None`	Line-of-sight direction vector (shape `(3,)`)
`losdirs`	`ndarray` or `None`	LOS rotation matrix (shape `(3, 3)`)
`length_unit`	`pint.Quantity`	Physical length per code unit (e.g., `"35000.0 ckpc/h"`)
`domain`	`ndarray`	Domain bounds, shape `(3, 2)`: `[[x_min, x_max], ...]`
`center`	`ndarray`	Domain center in box coordinates
`center_los`	`ndarray`	Domain center in LOS coordinates
`data`	dict-like	Dask-backed field accessor (see Available Fields)

`return_projection`

Create a 2D surface brightness map by binning photon positions along the line of sight.

projection, extent = ds.return_projection(nbins=512, zoomfactor=2.0)

# With units
projection, extent, units_str = ds.return_projection(return_units=True)
# units_str: "erg s⁻¹ cm⁻² arcsec⁻²"

Parameter	Type	Default	Description
`nbins`	`int` or `None`	auto	Pixel resolution. Defaults to `ngrid` from config or 512.
`img_size`	`float` or `None`	auto	Image size in code units. Defaults to domain extent.
`img_center`	`ndarray` or `None`	auto	Image center in LOS coordinates. Defaults to domain center.
`zoomfactor`	`float`	`1.0`	Zoom factor (divides `img_size`)
`return_mccounts`	`bool`	`False`	Return raw photon counts instead of physical units
`return_units`	`bool`	`False`	Append units string to return tuple

Returns: (projection, extent) or (projection, extent, units_str).

projection: 2D array of shape (nbins, nbins) in surface brightness units
extent: [x_min, x_max, y_min, y_max] in physical kpc, centered on the domain

`return_datacube`

Create a 3D volumetric histogram of photon positions.

cube, ex, ey, ez = ds.return_datacube(nbins=128)

# In redshift space
cube, ex, ey, ez = ds.return_datacube(
    use_redshift_space=True,
    unitstr="luminosity_density"
)

Parameter	Type	Default	Description
`nbins`	`int` or `list[int]`	`128`	Bins per dimension
`ranges`	`ndarray` or `None`	auto	Min/max per dimension in code units, shape `(3, 2)`
`unitstr`	`str`	`"luminosity_density"`	Unit type: `"luminosity_density"`, `"luminosity"`, or `"none"`
`use_los`	`bool`	`True`	Bin in LOS coordinates
`use_redshift_space`	`bool`	`False`	Shift z-axis by spectral offset (requires `use_los=True`)
`return_units`	`bool`	`False`	Append units string

Returns: (datacube, edges_x, edges_y, edges_z) with edges in physical kpc.

`return_spectrum`

Extract a spectrum within a circular aperture.

hist, bin_centers = ds.return_spectrum(
    center_los=ds.center_los,
    aperture_radius_arcsec=2.0,
    bins=np.linspace(-10, 10, 100)
)

Parameter	Type	Default	Description
`center_los`	`ndarray`	required	Center position in LOS coordinates (shape `(3,)`)
`aperture_radius_arcsec`	`float`	required	Aperture radius in arcseconds
`bins`	`ndarray` or `None`	`linspace(-10, 10, 100)`	Bin edges for the wavelength axis (in Angstrom offset)

Returns: (hist, bin_centers) — histogram counts and bin center positions.

Other Methods

Method	Returns	Description
`is_cosmological()`	`bool`	Whether the simulation has cosmological settings
`get_cosmology_and_redshift()`	`(Cosmology, float)`	Astropy Planck15 cosmology and redshift
`get_boxsize_and_hubblefactor()`	`(float, float)`	Box size in ckpc/h and Hubble parameter \(h\)
`pos_to_los_pos(pos)`	`ndarray`	Rotate box coordinates to LOS frame
`los_pos_to_pos(pos_los)`	`ndarray`	Rotate LOS coordinates back to box frame
`arcsec_per_unitlength()`	`float`	Angular scale at the simulation redshift
`dlambda_to_vel(dlambda)`	`float`	Convert wavelength shift to velocity in km/s

Available Fields

Fields are accessed via ds.data["fieldname"] and return Dask arrays.

Stored Fields

Field	Shape	Units	Description
`position`	`(N, 3)`	code length	Photon position
`direction`	`(N, 3)`	dimensionless	Photon direction
`weight`	`(N,)`	\(10^{42}\) erg	Photon luminosity weight
`weight_peel`	`(N,)`	\(10^{42}\) erg	Peeling weight (peel datasets only)
`dlambda`	`(N,)`	Angstrom	Wavelength shift from line center
`tau_seen`	`(N,)`	dimensionless	Optical depth accumulated
`tau_target`	`(N,)`	dimensionless	Target optical depth
`column_density_seen`	`(N,)`	cm\(^{-2}\)	Column density traversed
`nscatterings`	`(N,)`	dimensionless	Number of scattering events
`nsteps`	`(N,)`	dimensionless	Number of simulation steps
`state`	`(N,)`	dimensionless	Photon state flag
`current_cell`	`(N,)`	dimensionless	Current cell index
`origin`	`(N, 3)`	code length	Emission position
`id`	`(N,)`	dimensionless	Photon ID
`lsp_dist`	`(N,)`	code length	Distance to last scattering point
`dlambda_peel_shifted`	`(N,)`	Angstrom	Peeling-shifted wavelength offset (peel datasets only)
`next_cell`	`(N,)`	dimensionless	Next cell index

Derived Fields

Computed on-the-fly from stored fields:

Field	Description
`position_los`	Position rotated to LOS frame
`position_rs_los`	Position in redshift space (LOS frame)
`lsp_los`	Last scattering position in LOS frame (backtracked by `lsp_dist` for peel photons)

Example: Projection with Matplotlib

import matplotlib.pyplot as plt
from thor.mcrt import PhotonMCRT

mcrt = PhotonMCRT("config.yaml")
ds = mcrt.datasets["peel"]

projection, extent, units = ds.return_projection(
    nbins=512, zoomfactor=1.0, return_units=True
)

fig, ax = plt.subplots()
im = ax.imshow(
    projection.T,
    origin="lower",
    extent=extent,
    cmap="inferno",
)
ax.set_xlabel("x [kpc]")
ax.set_ylabel("y [kpc]")
plt.colorbar(im, label=units)
plt.show()

Unit Conversion

Physical units are computed using utilities from thor.common.units:

from thor.common.units import get_units
import astropy.units as u

# Luminosity (10^42 erg/s per photon weight)
units = get_units("luminosity", cosmology, z)

# Surface brightness (area must be an astropy Quantity)
pixel_area = (0.01 * u.Mpc) ** 2
units = get_units("surface_brightness", cosmology, z, area=pixel_area)

# Other unit types: "flux", "luminosity_density"

Box size and redshift are auto-detected from the config file using thor.common.boxsize:

First checks the dataset section for boxsize (in cm), converts to physical kpc
Falls back to HDF5 header attributes from the simulation snapshot
Zoom factor from gadget.zoom_box or gadget.zoom_target is applied when consider_zoom=True