OpticalMap ≡ NamedTuple of 3D histograms keyed by channel Symbol

Handy type alias for LEGEND optical maps, i.e. three-dimensional histograms for each SiPM channel. The field names are the channel symbols (e.g. :S030).

_λ0_model_bulk_ops()

Low-level implementation of λ0_model using bulk array programming, Reactant / CUDA compatible.

_λ0_model_loops()

Low-level implementation of λ0_model using for loops.

analytical_λ0_curve(log_p0_nominal, x0_random_coin, channel; eps_range=range(0.0,1.0,length=101))

Compute the per-channel no-light probability λ₀ as a function of a single channel's efficiency.

This analytical model computes the expectation per event λ₀jk = mean(p0jk * x0jk) with p0jk = exp(logp0jk * ε) and x0_jk the random coincidence indicators.

Arguments

• log_p0_nominal::Table: table of log p0 values (events × channels).
• x0_random_coin::Table: table of random-coincidence indicators (events × channels).
• channel::Symbol: Symbol identifying the channel (must be present in columnnames(logp0nominal)).
• eps_range: efficiency values to evaluate (vector or range).

Returns

• NamedTuple with fields :eps (vector of efficiencies) and :λ0 (vector of no-light probabilities evaluated at those efficiencies).

Notes

Can only be used if no multiplicity threshold is applied!

analytical_λ0_curve_all(log_p0_nominal, x0_random_coin; eps_range=range(0.0,1.0,length=101))

Compute λ₀(eps) curves for every channel present in log_p0_nominal.

Description

• For each channel (columns of log_p0_nominal) evaluate the no-light probability λ₀ across eps_range by calling compute_λ0_curve.
• The returned curves are suitable for locating the efficiency at which a measured λ₀ would be reproduced (e.g. by interpolation).

Arguments

• see above in description of compute_λ0_curve

Returns

• Dict{Symbol, NamedTuple} mapping each channel symbol to the NamedTuple returned by compute_λ0_curve (fields :eps and :λ0).

Notes

• Channel order is derived from columnnames(log_p0_nominal). Use that table to control which channels are computed and their ordering.

ar39_beta_energy_dist() -> MixtureModel{Uniform}

Energy distribution of the beta particle emitted in an Ar-39 nuclear decay.

Return a continuous probability distribution for the beta decay spectrum of Ar-39. The distribution is constructed from tabulated values from the IAEA BetaShape database, downloadable at this link.

detection_prob(h, coords...)

Return the bin content of histogram h at the given coordinates coords (with units). Coordinates must be inside the histogram bounds, otherwise an error is thrown.

detection_prob_bcast(h, xss, yss, zss)

Apply [detection_prob_event] to each triple of vectors (xs, ys, zs) from (xss, yss, zss).

This is useful when xss, yss, and zss are VectorOfVectors, e.g. the coordinates of all hits for many events. Returns a vector of vectors of bin contents.

estimate_efficiencies_from_curves(curves, λ0_d)

Given a mapping curves from some key (e.g. fiber symbol) to a NamedTuple with fields :eps and :λ0 (as produced by λ0_vs_efficiency_all or similar), and a dictionary λ0_d mapping channel symbols to observed no-light probabilities, compute the efficiency at which the simulated curve crosses the observed λ₀ for each channel.

Arguments

• curves::AbstractDict{K,NamedTuple}: maps a key (fiber or channel) to a NamedTuple (eps, λ0) where eps is a vector of efficiencies and λ0 the corresponding no-light probabilities.
• λ0_d::NamedTuple: observed λ₀ per channel symbol.

Returns

• Dict{Symbol,Float64} mapping channel symbol -> estimated efficiency.

The routine finds a bracket where the curve crosses the horizontal line λ0 = λ0_d[channel] and linearly interpolates between neighbouring points. If no sign change is found, it falls back to the nearest sampled point.

load_optical_map(filename, runsel) -> OpticalMap

Load a LEGEND-200 optical map from file.

Examples

load_optical_map("./optmap.lh5", (:p13, :r001))

log_p0_nominal(sim_data, optical_map; ...)

Logarithm of no-light-probability for events simulated by remage.

This can be used to compute inputs for λ0_model.

Arguments

• sim_data: output table (stepping data) for a scintillation detector (the same the optical map refers to) from remage. Must contain the fields xloc, yloc, zloc and edep.
• optmap: liquid argon optical map as loaded by load_optical_map.
• n_events: number of Ar-39 decays to simulate.
• light_yield: liquid argon scintillation yield.

Returns

A table of log(p0) for each channel (columns) and event (rows).

Examples

# load an optical map
optmap = load_optical_map("map.lh5", (:p13, :r001))

# load some remage simulation data
sim_data = lh5open("th228.lh5") do h5
    return h5["stp/lar"][:]
end

# call the function
log_p0_nominal(sim_data, optmap)

log_p0_nominal_ar39(optical_map, n_events; ...)

Logarithm of no-light-probability for simulated Ar-39 events.

Uses the energy distribution of the beta particle emitted in the decay of an Ar-39 nucleus and the detection probability in liquid argon (through the optical map) to estimate the probability of seeing no light for a set of random events. This can be used to compute inputs for λ0_model.

Arguments

• optmap: liquid argon optical map as loaded by load_optical_map.
• n_events: number of Ar-39 decays to simulate.
• light_yield: liquid argon scintillation yield.
• rand_voxel_kwargs...: optional keyword arguments forwarded to rand_voxel.

Returns

A table of log(p0) for each channel (columns) and event (rows).

make_efficiencies_prior(channels)

Construct a hierarchical prior for channels efficiencies.

• a global hyperprior scale ~ Beta(7, 5) sets the typical efficiency scale.
• each channel efficiency has a Beta prior with mean scale * shift and concentration concentration, enforcing positivity and preference for lower values.

Arguments

• channels: A list of channel names.

make_λ0_likelihood(x0, log_p0_nominal, x0_random_coin; multiplicity_thr=0, n_rands=10, smear_factor=0, device=CPUDevice()) -> DensityFunction

Construct the likelihood of no-light fractions per channel.

We model the fraction of events with no detected light (λ0) as follows:

• For each channel, the expected no-light probability λ0_model comes from the simulation (log_p0_nominal) combined with random coincidences (x0_random_coin) and scaled by per-channel efficiencies (the parameters of the model).

• The observed no-light fraction in data is λ0_data = N0 / N_data, where N0 is the number of no-light events passing a multiplicity threshold.

• Since N_data is large, the binomial distribution N0 ~ Binomial(N_data, λ0_model) can be approximated by a normal distribution: λ0_data ~ Normal(μ = λ0_model, σ² = λ0_model (1 - λ0_model) / N_data).

The likelihood is the sum of log-probabilities across all channels.

Arguments

• x0: observed no-light indicators from data events.
• log_p0_nominal: logarithm of the probability to see no light for each (event, channel), typically from simulations.
• x0_random_coin: observed no-light indicators from random coincidence events.
• multiplicity_thr: discard events with multiplicity below this threshold (optional, defaults to 0).
• n_rands: average forward model results over this amount of random numbers.
• smear_factor: the width of the likelihood gaussian terms is increased by a factor smear_factor * mean.
• device: on which device to run the computation of the forward model. (default CPUDevice())

Examples

Get some data:

using LegendOpticalFits

runsel = (:p13, :r001)
nev_sim = 10_000
nev_data = 1_000
multiplicity_thr = 6

optmap = load_optical_map("./optmap-p13.lh5", runsel)
log_p0 = log_p0_nominal_ar39(optmap, nev_sim)

x0 = x0_data("l200-p13-r001-ath-tier_evt.lh5", runsel, max_events=nev_data)
x0_rc = x0_data("l200-p13-r001-ant-tier_evt.lh5", runsel, max_events=nev_sim)

Build the likelihood (on the CPU by default):

logl = make_λ0_likelihood(x0, lp0, x0_rc, multiplicity_thr=multiplicity_thr)

CUDA via Reactant/XLA:

using Reactant
Reactant.set_default_backend("cuda")

logl = make_λ0_likelihood(x0, lp0, x0_rc, multiplicity_thr=multiplicity_thr, device=ReactantDevice())

Returns

• A DensityFunction object representing the log-likelihood. It can be called with a parameter set.

rand_voxel(optmap::OpticalMap; xrange = nothing, yrange = nothing, zrange = nothing) -> (ix, iy, iz)

Sample a random valid voxel (bin indices) from an OpticalMap.

The function draws random voxel indices (ix, iy, iz) within the histogram domain of the optical map. The histogram of the first channel is used to determine the geometry (all channels share the same dimensions).

Arguments

• optmap: optical map (see load_optical_map.
• xrange, yrange, zrange: optional (min,max) in axis units. If nothing (default), the full axis range is used.

Returns

Tuple (ix, iy, iz) of voxel indices.

x0_data(filename, runsel; ...) -> Table

Load x0 values from LEGEND-200 SiPM data.

For each event and channel, this function records whether no photon above the high photoelectron threshold was detected within the [-1, 5] μs coincidence window. The output is a Table where each column corresponds to a SiPM channel and each row to an event.

Arguments

• filename: path to the pygama LH5 file containing event-tier data.
• runsel: (period, run) identifier, used to select the proper channel map.
• max_events: (optional) maximum number of events to read, default 10_000.

Returns

a Table of booleans with dimensions (max_events, n_channels). each entry is true if no qualifying photon was observed, false otherwise.

Examples

x0_data("l200-p13-r003-anp-20241217T094846Z-tier_evt.lh5", (:p13, :r003))

x0_mock(efficiencies_true, log_p0_nominal, x0_random_coin)

Generate mock x0 data (bernoulli draws folded with random-coincidence mask) from efficiencies_true and log_p0_nominal from model.

λ0_data(x0; multiplicity_thr=0) -> NamedTuple

Compute per–channel no-light fractions from a boolean Table of events.

Arguments

• x0: a Table where each column corresponds to a channel and each row to an event; entries are Bool indicating whether no qualifying photon was observed in the channel.
• multiplicity_thr: minimum number of channels with light per event required for the event to be considered. Defaults to 0.

Returns

A NamedTuple with one field per channel containing the fraction of selected events in which that channel had no light, and the total number of events passing the multiplicity threshold.

λ0_model(efficiencies, log_p0_nominal, x0_random_coin[, multiplicity_thr])

Expected fraction of events in which a SiPM channel sees no light.

Computes the expected fraction by looping over the rows (events) in log_p0_nominal. It samples the expected light/no-light observable by combining the probability to see no light p0 with the random coincidence data (at the same row index). If a multiplicity threshold is specified, events with lower multiplicity are discarded.

Arguments

All input data is keyed by detector name (a symbol)

• efficiencies: scaling factors for each SiPM channel.
• log_p0_nominal: logarithm of the probability to see no light for each (event, channel), typically from simulations.
• x0_random_coin: presence of light from random coincidences for each (event, channel). This is typically coming from a measurement.
• multiplicity_thr: discard events with multiplicity below this threshold (optional, defaults to 0).
• n_rands: average forward model results over this amount of random numbers.

Returns

• Vector of expectation values for each channel, ordered as the input data structures.