LightHit

LightHit computes the expected optical-module response to light propagating in a homogeneous scattering medium. It solves the time-dependent radiative transfer equation (RTE), including unscattered, once-scattered and multiply-scattered light, and returns expected photoelectrons per module and per time bin.

The medium, detector and source are ordinary Python objects. There is no global initialisation state, so several configurations can safely coexist in one program.

LightHit is currently an alpha release. Validate numerical settings for your medium and geometry before using its results in an analysis.

Installation

python -m pip install lighthit

Numba acceleration is strongly recommended for tracks, showers and large detector arrays:

python -m pip install 'lighthit[accelerate]'

For the installed synthetic demonstration, viewer and wheel self-tests:

python -m pip install 'lighthit[demo,test]'
lighthit-selftest

Python 3.11 or newer is required.

To run the detector-specific scripts in this repository, install the checkout into that environment first: python -m pip install -e '.[demo,test]'. A script under examples/ does not automatically import the adjacent src/ tree; an older PyPI installation in the active environment will otherwise be used.

Quick start

The example below defines a spectral medium, two optical modules and a monochromatic isotropic flash.

import numpy as np
import lighthit as lh

# Optical properties tabulated versus wavelength.
medium = lh.SpectralMedium(
    wavelength_nm=[400.0, 450.0, 500.0],
    absorption_per_m=[0.030, 0.020, 0.040],
    scattering_per_m=[0.030, 0.022, 0.017],
    phase_index=[1.344, 1.339, 1.336],
    group_index=[1.386, 1.374, 1.367],
    g=0.9,
)

# Response functions receive NumPy arrays and return arrays of the same shape.
def angular_acceptance(head_on_cosine):
    return np.ones_like(np.asarray(head_on_cosine, dtype=float))

def spectral_efficiency(wavelength_nm):
    return np.full_like(np.asarray(wavelength_nm, dtype=float), 0.20)

detector = lh.DetectorArray(
    positions_m=[[20.0, 0.0, 0.0], [35.0, 0.0, 0.0]],
    orientations=[[-1.0, 0.0, 0.0], [-1.0, 0.0, 0.0]],
    effective_area_m2=0.05,
    angular_acceptance=angular_acceptance,
    spectral_efficiency=spectral_efficiency,
)

source = lh.IsotropicFlash.monochromatic(
    position_m=[0.0, 0.0, 0.0],
    photons=1.0e8,
    wavelength_nm=450.0,
)

# build() returns a reusable kernel and prepares only the tables this source
# needs. The same kernel can then transport more compatible events.
kernel = lh.build(medium, detector, source)
response = kernel.transport(source)

print(response.charge_pe)   # integrated expected photoelectrons, one per OM
print(response.bins_pe)     # expected photoelectrons, shape (OM, time bin)

The first build may take noticeably longer than later event calculations. Set KernelConfig(cache_directory=...) to reuse transport tables across processes.

Sources

The high-level dispatcher accepts:

flash = lh.IsotropicFlash.monochromatic(
    position_m=[0, 0, 0], photons=1e8, wavelength_nm=450
)

track = lh.CherenkovTrack(
    start_m=[0, 0, -20],
    direction=[0.2, 0.1, 0.97],
    length_m=40,
    beta=1.0,
)

shower = lh.G4Shower.from_hdf5("event.h5", event=0)

synthetic_shower = lh.SyntheticShower.gaussian(
    lh.SourcePose(
        position_m=[25.0, -10.0, -720.0],
        direction=[0.4, 0.2, -0.89442719],
        time_ns=0.0,
    ),
    charged_track_length_m=300.0,
)

For tracks and showers, kernel.transport(source) automatically selects the source engine and uses the detector's angular response. A compatible prepared wavelength-folded cache is used when present; set KernelConfig(spectral_folded_cache=False) for the original reference path.

Applications can select any calculated signal without using the CLI:

response = kernel.transport(track)
prompt = response.select("0+1")  # unscattered plus once-scattered
print(prompt.charge_pe, prompt.bins_pe, prompt.rate_pe_per_ns)
# Other selectors: 0, 1, ">=2", "all".

This selects stored components; it does not skip the multiply-scattered solve or make the transport itself faster.

For a shower, transport_hybrid instead computes orders 0 and 1 with the direct physical-time prompt solver and order >=2 with the directional RTE. It requires the Numba accelerator and a prepared directional cache. With the standard 321-frequency cache, only 65 selected nodes are applied by default:

hybrid = kernel.transport_hybrid(shower, frequency_nodes=65)
print(hybrid.charge_pe, hybrid.bins_pe)

HybridTransportResponse.spectrum_ge2_pe contains only the >=2 spectrum on the computed frequency grid. Prompt orders have exact time bins but no stored nonzero-frequency spectrum. The active OMs receive the full prompt calculation; the integrated order-1 charge for screened, inactive OMs is a coarse estimate and is identified by metadata["prompt_order1_level"]. For a screening control, pass prompt_config=lh.PromptConfig(screen=False); this evaluates direct order 1 for every OM and costs more time.

To reduce the cache itself, build it on 65 frequencies and specify a separate 321-node readout grid. The readout grid is used only after the expensive solve:

config = lh.KernelConfig(
    omega_per_ns=np.linspace(0.0, 1.2, 65),
    readout_omega_per_ns=np.linspace(0.0, 1.2, 321),
    cache_directory="hybrid-cache",
    cache_policy="require",
)
kernel = lh.TransportKernel(medium, detector, config)
hybrid = kernel.transport_hybrid(shower, frequency_nodes=65)

Shower basis for order >=2

The order >=2 source of a shower can instead be represented on an axial-slice basis: cubic slices along the shower axis (0.25 m), a five-point transverse moment stencil and the exact angular harmonics of the kernel configuration. For compact showers, orders 0+1, the RTE cache, the active-OM rule and the readout are unchanged; only the order >=2 source differs:

basis = lh.ShowerBasis()          # cell_m=0.25, cubic slices, stencil
hybrid = kernel.transport_hybrid(shower, frequency_nodes=65, ge2_basis=basis)

The prepared order >=2 radial splines are kept on the kernel and reused by every later event. To start a job without rebuilding them, save them once and memory-map them afterwards:

kernel.ge2_kernel_store().save("ge2-pack")      # after the first event
kernel.ge2_kernel_store("ge2-pack")             # in a later process

On one BGVD cluster and a 100 GeV electron a warm event took 14 s instead of 628 s on the same 4-core machine (1.9 instead of 5.6 GiB), with the same active OMs, order >=2 charge within 2e-4 and total time-bin L1 within 0.2 % of the default hybrid; the basis source is closer to a converged source reference than the default one. See docs/research/shower-basis.md for the measured errors, including 1 TeV electrons and muons with showers.

For a long G4 muon source, the basis now uses a separate 3–300 m source window for each OM and sums the zero-frequency charge before selecting active modules. An explicitly selected nearly straight track can use the axial m=0 channel while its secondaries keep the shower basis:

hybrid = kernel.transport_hybrid(
    muon_event_source, frequency_nodes=65,
    ge2_basis=lh.ShowerBasis(track_uid=primary_muon_uid))

The track option is an approximation and remains opt-in. See docs/research/long-muon-radial-windows.md for the NTSim checks, cutoff comparison, timing, memory and current limitations.

The same ShowerBasis also supports cone_model="spectral": it applies the Cherenkov threshold and cone angle separately at every wavelength node. The existing RTE cache and prepared kernel pack are reused. On the checked 100 GeV NTSim electron it took 12.2 s versus 116 s for the spectral hybrid, with 26 matching active OMs. See docs/research/spectral-shower-basis.md for the formula, per-OM accuracy and NTSim run commands. The explicitly selected spectral muon route can additionally fuse wavelength bands and use 33 of the prepared 65 frequencies; on two checked 100 GeV G4 muons this preserved the active OMs and kept order >=2 time-bin L1 differences below 3% on bright OMs. See docs/research/muon-fast-spectral.md.

Configuring the calculation

Numerical settings are collected in KernelConfig:

config = lh.KernelConfig(
    wavelength_nodes=9,
    threshold_pe=0.01,
    cache_directory="lighthit-cache",
)

kernel = lh.build(medium, detector, track, config=config)
response = kernel.transport(track)

For production jobs, prepare the cache explicitly and forbid event processing from starting an expensive build unexpectedly:

build_config = lh.KernelConfig(cache_directory="lighthit-cache")
kernel = lh.TransportKernel(medium, detector, build_config)
kernel.build(source=track, progress="console")
print(kernel.last_build_report.as_dict())

run_config = lh.KernelConfig(
    cache_directory="lighthit-cache",
    cache_policy="require",
)
kernel = lh.TransportKernel(medium, detector, run_config)
response = kernel.transport(track)

Repeated build calls reuse compatible tables and report that there is nothing to do. Pass force=True only for an intentional rebuild.

For many independent Cherenkov events on the same medium and OM spectral response, an optional second build folds the wavelength quadrature into two directional kernels (one for a beta-one straight track). This avoids preparing nine large radial splines during each event:

kernel.build_folded_directional()
fast_config = lh.KernelConfig(
    cache_directory="lighthit-cache", cache_policy="require",
    spectral_folded_cache=True,
)
fast_kernel = lh.TransportKernel(medium, detector, fast_config)
response = fast_kernel.transport(track)

Use the same numerical configuration in both kernels. Once a compatible derived cache exists, it is selected automatically; set spectral_folded_cache=False to force the original nine-wavelength reference. Without a derived cache the original path remains available. The shared radial interpolation is a numerical approximation, so compare representative events against the reference before analysis. On the BGVD 321-frequency profile the three prepared arrays occupy about 10 GiB on disk.

Detector coordinates can be inspected without detector-specific assumptions:

geometry = lh.describe_geometry(detector)
print(geometry.centroid_m)
print(geometry.bounds_min_m, geometry.bounds_max_m)

For a point laser, cover every source-to-module distance when preparing its cache. This matters for a bright source: modules that could be safely skipped at low photon count may contribute above threshold_pe at high photon count.

laser_range = lh.point_source_radial_range(detector, flash.position_m)
laser_config = lh.KernelConfig(
    radial_range_m=laser_range,
    k_panel_per_m=min(0.04, 2 * np.pi / laser_range[1]),
    cache_directory="laser-cache",
)
laser_kernel = lh.build(medium, detector, flash, config=laser_config)

The returned range guarantees geometric coverage. The spatial panel rule resolves oscillations at the farthest OM better than the default panel; check convergence in k_panel_per_m, radial_nodes and the frequency grid for the medium and geometry being analysed. For several laser positions, choose one shared radial range that covers the signals of interest; the cache contains a radial transport kernel, not a source-position-specific result.

Important inputs are explicit:

• SpectralMedium contains absorption and scattering coefficients in m⁻¹, phase and group refractive indices, and the Henyey–Greenstein parameter g.
• DetectorArray contains OM positions, orientations, effective areas, angular acceptance and wavelength-dependent detection efficiency.
• A detector orientation points from the OM toward a head-on source; angular_acceptance(+1) is the head-on response.
• spectral_efficiency should include all wavelength-dependent detection factors required by the application, such as quantum efficiency and optical transmission.

Reading the response

kernel.transport(...) returns a TransportResponse.

Absolute bin edges for every module are

absolute_edges_ns = (
    response.time_origin_ns[:, None]
    + response.relative_time_edges_ns[None, :]
)

The integrated charge is evaluated independently at zero frequency. A finite time window and a finite frequency grid mean that bins_pe.sum(axis=1) need not equal charge_pe exactly. Signed ringing is reported rather than silently clipped or renormalised.

Responses can be stored as compressed NumPy files:

response.save("response.npz")

Interactive viewer

Install the optional viewer dependency:

python -m pip install 'lighthit[viewer]'

Turn a response into a standalone interactive HTML file:

payload = lh.viewer_payload(
    response,
    source=source,
    event_id="flash-450nm",
    label="450 nm calibration flash",
)

viewer_path = lh.write_event_viewer(payload, "lighthit-viewer.html")
print(viewer_path)

Compatible event payloads can share one viewer:

combined = lh.merge_event_viewers(laser_payload, track_payload, shower_payload)
lh.write_event_viewer(combined, "multi-event-viewer.html")

Open lighthit-viewer.html in a web browser. It contains the detector geometry, integrated charge by scattering order, a selectable per-OM time histogram and a time animation. A portable lighthit-viewer.json companion is written beside the HTML file. The Display selector offers raw bins, a 3 ns Gaussian view, or both overlaid on the OM waveform. Each scattering order and 0+1 has its own raw and smoothed curve. The Gaussian view redistributes already binned integrals under a uniform-within-bin assumption; it is for display and does not replace an exact detector timing response or modify saved raw data.

For tracks and showers the animation also shows the source itself: on the default detector clock a gold head runs along the muon track, or the shower spindle grows along its axis, in step with the time slider, while each OM shows its bin at time_origin + relative time. The Per-OM origin clock keeps the older view in which every OM is shown at its own time zero. The source timing is an optional timing block of an axis source (lighthit.viewer.source_timing builds it from light elements); payloads without it stay valid and static.

Installed synthetic demonstration

The wheel includes one public-data-only Baikal-like example. Its geometry and OM response are illustrative and are not a Baikal calibration:

lighthit-demo geometry
lighthit-demo cache --profile quick --cache ./demo-cache --sources all
lighthit-demo run --profile quick --cache ./demo-cache \
  --output ./demo-output --sources all

The repository, but not the wheel, contains detector-specific laser, track, multi-track and G4-shower examples with explicit source positions and directions.

Scope

The current model assumes a homogeneous, unbounded medium with elastic Henyey–Greenstein scattering. It does not include boundaries, layered media, structural shadowing, polarisation, detector electronics, trigger or noise.

License

LightHit is distributed under the BSD 3-Clause License.