Metadata-Version: 2.4
Name: lighthit
Version: 0.2.0a1
Summary: Spectral radiative-transfer response for neutrino-telescope optical modules
Author: Dmitry Naumov
Project-URL: Repository, https://github.com/NeutrinoHit/LightHit
Keywords: radiative transfer,neutrino telescope,Cherenkov,Baikal-GVD
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Science/Research
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Topic :: Scientific/Engineering :: Physics
Requires-Python: >=3.11
Description-Content-Type: text/markdown
Requires-Dist: numpy>=1.26
Requires-Dist: scipy>=1.11
Provides-Extra: accelerate
Requires-Dist: numba>=0.61; extra == "accelerate"
Provides-Extra: notebook
Requires-Dist: jupyterlab>=4; extra == "notebook"
Requires-Dist: ipykernel>=6; extra == "notebook"
Requires-Dist: nbformat>=5; extra == "notebook"
Requires-Dist: nbclient>=0.10; extra == "notebook"
Requires-Dist: nbconvert>=7; extra == "notebook"
Requires-Dist: matplotlib>=3.8; extra == "notebook"
Requires-Dist: plotly>=5; extra == "notebook"
Provides-Extra: viewer
Requires-Dist: plotly>=5; extra == "viewer"
Provides-Extra: dev
Requires-Dist: pytest>=8; extra == "dev"
Requires-Dist: build>=1; extra == "dev"

# 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.

```python
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")
```

`method="auto"` chooses `isotropic`, `track`, or `axial` from the source type.
Research variants (`axial_full`, with `generic` as its descriptive alias) share
the same dispatcher but require explicit permission:

```python
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,

```text
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/production-rte-to-pe.qmd`](docs/production-rte-to-pe.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`:

```bash
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:

```bash
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

```bash
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:

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

A quick pass and a stricter numerical setting:

```bash
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:

```bash
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

```bash
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:

```bash
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`:

```bash
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:

```bash
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

```python
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.

```python
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 the
synthetic $g=0.7$ examples are not a validation of $g=0.9$ real water, and
the (scattering degree, spatial degree, $k_{\max}$) triple that converges at
one radius does not automatically converge at another.
`docs/chapters/04-bgvd-water.qmd` works through this in detail: the
constraint that actually governs convergence is the spatial multipole
degree tracking $k_{\max}\cdot r$, not $g$ as such.

## What exactly is solved, and where the formulas are

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

| Chapter | Content |
|---|---|
| `docs/production-rte-to-pe.qmd` | Production chain from spectral RTE and two source fields to expected photoelectrons |
| `docs/chapters/01-rte.qmd` | The radiative transfer equation, from a photon balance |
| `docs/chapters/02-point-source.qmd` | Point source, point detector: the adjoint system, the exact free tail, spatial inversion |
| `docs/chapters/03-histogram-binning.qmd` | Why orders 0 and 1 are binned in physical time instead of Fourier-inverted |
| `docs/chapters/04-bgvd-water.qmd` | The solver run against measured Baikal water at 450 nm |
| `docs/appendices/notation.qmd` | Symbol table |
| `docs/VALIDATION.md` | Numbers from the checks that were actually run |

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.

| File | Contents |
|---|---|
| `src/lighthit/medium.py` | Parameters of one spectral node, m and ns |
| `src/lighthit/angular.py` | Tridiagonal adjoint problem, exact free tail |
| `src/lighthit/single.py` | Coordinate-space first order with the full HG function |
| `src/lighthit/single_fast.py` | Optional compiled scalar kernel for the same exact first order |
| `src/lighthit/green.py` | Radial inversion and the combined 0+1+≥2 spectrum |
| `src/lighthit/readout.py` | Bins and instrument smearing, independent of transport |
| `src/lighthit/viewer.py` | Standalone multi-event 3-D viewer, OM charges, time histograms and animation |
| `src/lighthit/providers.py` | Local private-water provider |
| `tests/` | Independent matrix, analytic, and geometric checks |

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. Orders 0 and 1 are integrated directly in time, so their
fronts do not depend on the chosen frequency cutoff; the inverse Fourier
transform is applied only to the $\ge2$ part (`docs/chapters/03-histogram-binning.qmd`
explains why).

## What is not in the first module

A finite OM sphere, real acceptance and a PDE for it, cones and track
segments, showers, spectral convolution, a general spatial cache,
SVD/NUFFT, GPU execution, boundary surfaces, and an inhomogeneous medium.
Measurements on the test medium do not confirm accuracy on private optics.

Sources for the mathematical construction and the license status are in
[PROVENANCE.md](PROVENANCE.md). Copyright holders and a license must be
agreed before any publication; this delivery assigns neither.
