LightHit

Spectral, time-dependent radiative transport and optical-module response. The production API folds wavelength-dependent water, Cherenkov emission, measured detector geometry, OM angular acceptance, quantum efficiency and glass/gel transmission into expected photoelectrons. The medium is homogeneous and unbounded; scattering uses Henyey–Greenstein with explicit g.

Production API

No global initialization state is used. A TransportKernel owns one medium, one detector and reusable wavelength-dependent caches; all sources and engines return the same TransportResponse contract.

import lighthit as lh

bgvd = lh.load_bgvd_model(
    "/private/path/bgvd-model-master",
    dataset="2021",                 # all 2304 OMs by default
)

kernel = bgvd.kernel(
    lh.KernelConfig(
        cache_directory="/local/cache/lighthit-bgvd",
        wavelength_nodes=9,
        threshold_pe=0.01,
    ),
).build()

laser = lh.IsotropicFlash.monochromatic(
    position_m=[0, 0, 0], photons=1e12, wavelength_nm=532,
)
laser_response = kernel.transport(laser)

track = lh.CherenkovTrack(
    start_m=[0, 0, -50], direction=[0.2, 0.1, 0.97], length_m=120,
)
track_response = kernel.transport(track, method="track")

shower = lh.G4Shower.from_hdf5("event.h5", event=5)
shower_response = kernel.transport(shower, method="axial")
assert shower_response.metadata["detector_angular_model"] == "exact_m_blocks"

A method name says how the source is handled. How the module is handled is a property of the module: a detector whose angular acceptance has a measured Legendre bandwidth gets the exact directional kernel of the m-block solver, and the result reports which model was used in metadata["detector_angular_model"]. method="auto" chooses isotropic, track, or axial from the source type; method="directional" forces the exact module model even for an isotropic detector, and the superseded source-centroid model is kept only as the experimental axial_centroid / track_centroid for regression.

build() builds the tables that choice will actually need --- call kernel.cache_kinds(method) to see which --- so a directional detector no longer pays for the multipole tables it will not read. It defaults to the element-source route; pass build(method="all") when both a flash and an event will be transported. build(source=...) resolves both the table and the wavelengths from the source, which matters for a monochromatic flash: it reads one line, not the quadrature, and building the quadrature for it produced two tables that were never touched and left the real one to the first transport.

Research variants (axial_full, with generic as its descriptive alias) share the same dispatcher but require explicit permission:

control = kernel.transport(
    shower, method="axial_full", allow_experimental=True,
)
print(kernel.available_methods())

The BGVD adapter reads the separately supplied private checkout at runtime. It uses the full median_om_coordinates_2021.csv, wavelength-dependent water, efficiency(lambda) * transmission_gel_glass(lambda), and the private angular OM polynomial. It deliberately ignores the private scattering indicatrix and uses HG with g=0.9. Neither the checkout nor derived private files are part of the LightHit distribution.

Cherenkov spectra use two source fields,

S0(lambda) = lambda^-2
S2(lambda) = lambda^-2 * n_phase(lambda)^-2

with spatial/angular coefficients stored independently. The complete derivation and direct mapping from every formula to code is in docs/chapters/12-rte-to-photoelectrons.qmd.

Minimal runnable integrations are examples/production_laser.py, production_track.py, and production_shower.py.

Low-level monochromatic solver

PointGreenSolver remains the reviewed low-level Green function of an instantaneous monochromatic point flash and a point isotropic receiver. A directed single photon and an isotropic single photon are supported.

What the calculation returns

PointGreenSolver.solve computes the complex response spectrum per unit effective area, m⁻². The zero-frequency component gives the integrated signal ("charge"). result.readout gives bin integrals, m⁻²; dividing by the bin width gives the mean registration rate, m⁻²·ns⁻¹.

Detection efficiency is 1, and sensitivity is the same in every direction. Area is never silently replaced by a bare number: multiplying by a small effective area gives the expected photon count. There is no OM surface, no first-entry condition, and no shadowing in this model.

A directed delta flash observed exactly on its own forward ray gives a singular response; that request is rejected. The API treats $\cos\theta\ge1-10^{-12}$ as singular. Geometries that close would need a finite aperture or a finite angular source distribution instead. In the shipped demonstration the angle is 60°, so the direct light is zero and the whole signal is scattered. For an isotropic flash the direct light is a delta function in time with a finite integral; bin integrals across it are exact.

Install into a clean checkout

Python 3.11 or newer is required. Commands below run from ~/Projects/LightHit:

python3 --version
python3 -m venv .venv
source .venv/bin/activate
python -m pip install --upgrade pip
python -m pip install -e '.[accelerate,dev]'
python -m pytest -q

After publication, the equivalent user installation is:

python -m pip install 'lighthit[accelerate]'

Add the notebook or viewer extra only when those interfaces are needed.

None of these commands creates a Git commit or publishes anything.

A first calculation without Jupyter

python -m lighthit --config examples/point-green.toml \
  --output .build/point-green --plots

With the optional acceleration dependency installed, exact first scattering uses its compiled scalar kernel automatically; select the compiled angular tail explicitly with --angular-backend numba. Use --single-backend numpy to force the uncompiled reference path for a comparison.

Output:

• .build/point-green/spectrum.npz: frequencies and the 0, 1, ≥2 components;
• .build/point-green/profiles.npz: unsmeared bins and bins with the given readout;
• .build/point-green/report.json: parameters, charge, timings, diagnostics;
• .build/point-green/figures/: separate profile and spectrum plots.

Charge without building a time spectrum:

python -m lighthit --charge-only --repeat 3 --output .build/charge

A quick pass and a stricter numerical setting:

python -m lighthit --preset quick --output .build/quick --plots
python -m lighthit --preset refined --output .build/refined --plots
python scripts/compare_runs.py .build/point-green .build/refined

Numerical library threads can be pinned for comparable timings:

OPENBLAS_NUM_THREADS=1 OMP_NUM_THREADS=1 VECLIB_MAXIMUM_THREADS=1 \
  python -m lighthit --repeat 3 --output .build/benchmark

--repeat repeats the calculation itself, not loading a saved result. All timings are kept in the report. This first release has no on-disk transport cache. Within one solve call, the angular solution is shared across every observation point.

Notebook

python -m jupyter lab notebooks/01_point_green.ipynb

Jupyter menu: Run → Run All Cells. The notebook calls the same package as the CLI. It shows normalization, geometry, the phase function, the spectral equations, charge, time profiles, error diagnostics, and timings.

Non-interactive execution and HTML:

mkdir -p .build/notebook
python -m jupyter nbconvert --to notebook --execute \
  notebooks/01_point_green.ipynb --output 01_point_green.executed \
  --output-dir .build/notebook --ExecutePreprocessor.timeout=600
python -m jupyter nbconvert --to html \
  .build/notebook/01_point_green.executed.ipynb

The source notebook is kept in Git without outputs. After running it, especially against private water, strip outputs before git add:

python scripts/strip_notebook_outputs.py notebooks/01_point_green.ipynb

.gitignore does not remove outputs already committed to a tracked notebook.

Interactive event viewer

The general event viewer shows several G4 showers and point flashes in one standalone HTML file: 3-D detector geometry, integrated charge by scattering order, a selectable per-OM time histogram, a full-array time heatmap, and frame/cumulative animation. Shower axes are drawn as directed spindles and muons as directed lines with their stored length. Signed Fourier-inversion bins are displayed with a signed-log colour scale and are never clipped or renormalized; adjacent bins can be combined ×2/×4 for display.

With the local files under g4_data, run all four stored samples (event 5) and an isotropic laser through the shared fast Numba cache:

python -m pip install -e '.[accelerate,viewer]'
python scripts/run_event_viewer.py --output .build/event-viewer --threads 4

The default triangular array has three clusters. Electron showers are placed obliquely near individual clusters, the long muons cross the array, and the laser sits at the triangle centre. Signals are converted to expected photoelectrons with explicit effective area and efficiency. Before computing a full spectrum, a cheap axial proxy plus an exact zero-frequency guard skips OMs below --threshold-pe 0.01; set the threshold to zero to disable this.

Open .build/event-viewer/viewer.html. The same viewer can load another portable viewer.json through its file button. Input HDF5 files and generated event arrays remain local and are not package data.

The screening construction and measured speedups are recorded in docs/research/signal-screening.md. The separate cutoff/node/bin-width study of time ringing is in docs/research/time-ringing.md.

Python API

import numpy as np
from lighthit import PointGreenSolver, SolverSettings, synthetic_medium

medium = synthetic_medium()  # test values close to Baikal water, not a calibration
solver = PointGreenSolver(medium, SolverSettings())
omega = np.linspace(0.0, 1.2, 241)  # rad/ns; include 0 for the charge
result = solver.solve(
    omega,
    displacement_m=[17.32050807568877, 0.0, 10.0],  # detector minus source, m
    direction=[0.0, 0.0, 1.0],
)
print(result.charge_per_m2)
print(result.timings_s)
front = result.front_time_ns[0]
edges = np.arange(front - 30, front + 651, 2.0)
raw = result.readout(edges, sigma_ns=0.0)
measured = result.readout(edges, sigma_ns=3.0)

For an isotropic flash, pass direction=None. For several detectors, pass displacement_m of shape (D,3). Emission time is set by emission_time_ns; the number of emitted photons by photons. A point isotropic detector needs no normal vector.

Private BGVD model

The distribution contains no private optical tables, OM coefficients, geometry rows, event files, or derived private responses. load_bgvd_model reads BaikalWater.py, OpticalModule.py, and the selected geometry CSV from a separately supplied checkout, recording content hashes rather than copying the inputs.

import lighthit as lh

bgvd = lh.load_bgvd_model(
    "/full/local/path/to/bgvd-model",
    dataset="2021",
    g=0.9,
)

The path may point at a checkout or its bgvd_model directory. If the private package is installed, the path can be omitted. The older providers.load_bgvd_water remains available for a low-level monochromatic PointGreenSolver, but intentionally has no detector response.

Keep private caches and results outside the source tree or under ignored local directories. The production adapter is checked against a synthetic double and locally against the separately supplied package; none of that package is written into the wheel or source distribution.

A run at $g=0.9$ needs its own scan of settings: numbers validated for a synthetic medium are not automatically a validation of real water, and the (scattering degree, spatial degree, $k_{\max}$) triple that converges at one radius does not automatically converge at another. The convergence variables and their separate roles belong to docs/chapters/15-validation-and-error-budget.qmd.

What exactly is solved, and where the formulas are

The math is written up as a short Quarto book in docs/:

Render it with Quarto (cd docs && quarto preview, or quarto render for a static copy in docs/_book/); nothing in the book needs private data or Geant4.

Every collision order is included. Orders $\ge2$ use HG coefficients only up to $L$, while the exact first order uses the full HG function — an explicit, stated composite approximation of the full HG response, whose difference from the full HG answer vanishes as $L\to\infty$. Free angular transport is not truncated at $L$: its infinite tail is removed analytically. The spatial angular inversion has its own, independent degree $J$.

The $k$ range and quadrature, the frequency band, and the frequency step are all finite as well. quick, balanced, refined change several spatial/angular settings at once; they do not fix a mathematical error bound. Unsmeared bins can oscillate right at the light front. The code keeps negative values and mass before the front; nothing is clipped or renormalized. The low-level point-source readout integrates orders 0 and 1 directly in physical time. The production track/shower route bins the ballistic term directly but Fourier-inverts both scattered components, so both are subject to the finite frequency grid. docs/chapters/07-time-and-bins.qmd keeps these routes separate.

Model boundaries

The production model includes measured OM orientations, angular acceptance, spectral efficiency, Cherenkov tracks and showers, wavelength integration and reusable radial caches. It does not model a finite OM surface, structural shadowing, boundaries, layered or inhomogeneous media, polarisation, detector electronics, trigger/noise, or GPU execution.

Sources and attribution for the mathematical construction are recorded in PROVENANCE.md. LightHit is distributed under BSD-3-Clause; see LICENSE.