NT Summary Stats

Fast computation (C++) of traditional summary statistics for neutrino telescopes for Python. The summary stats computed are based off of this paper.

Install

pip install nt_summary_stats

The wheel ships with the C++ backend enabled by default; source installs will build the nt_summary_stats._native extension automatically (requires a C++17 compiler and pybind11). Set NTSS_DISABLE_NATIVE=1 if you need to fall back to the pure NumPy implementation.

Usage

import numpy as np
from nt_summary_stats import compute_summary_stats

# Basic usage (9 standard statistics; mode="standard" is the default)
times = np.array([10.0, 15.0, 25.0, 100.0])      # shape: (N,), dtype: float
charges = np.array([1.0, 2.0, 1.5, 0.5])         # shape: (N,), dtype: float
stats = compute_summary_stats(times, charges)     # returns: np.ndarray, shape (9,)

print(stats[0])  # total_charge: 5.0
print(stats[3])  # first_pulse_time: 10.0
print(stats[7])  # charge_weighted_mean_time: 26.0

# Minimal statistics (4 statistics)
minimal_stats = compute_summary_stats(times, charges, mode="minimal")  # shape (4,)
print(minimal_stats)  # [total_charge, first_pulse_time, cw_mean_time, cw_std_time]

# Extended statistics (25 total statistics)
extended_stats = compute_summary_stats(times, charges, mode="extended")  # returns: np.ndarray, shape (25,)
print(extended_stats.shape)  # (25,)
print(extended_stats[9])     # charge_5_percent_time
print(extended_stats[21])    # n_pulses: 4.0

Process detector events that expose the required photon-level keys:

from nt_summary_stats import process_event

# Input: event dictionary with required photon fields
event_data = {
    'photons': {
        'sensor_pos_x': [0.0, 0.0, 100.0],  # list[float], length M
        'sensor_pos_y': [0.0, 0.0, 0.0],    # list[float], length M
        'sensor_pos_z': [0.0, 0.0, 50.0],   # list[float], length M
        'string_id': [1, 1, 2],             # list[int], length M
        'sensor_id': [1, 1, 1],             # list[int], length M
        't': [10.0, 15.0, 20.0]             # list[float], length M
    }
}

# Default: no grouping (uses all hits as-is)
sensor_positions, sensor_stats = process_event(event_data)

# Optional: group hits within time windows
sensor_positions, sensor_stats = process_event(event_data, grouping_window_ns=2.0)
# sensor_positions: np.ndarray, shape (N_sensors, 3), dtype: float64
# sensor_stats: np.ndarray, shape (N_sensors, 9), dtype: float64
# Arrays are aligned: sensor_positions[i] corresponds to sensor_stats[i]

Process individual sensor data:

from nt_summary_stats import process_sensor_data

# Input: sensor hit data
sensor_times = [10.0, 10.5, 15.0, 100.0]    # list[float] or np.ndarray(N,)
sensor_charges = [1.0, 0.5, 2.0, 1.0]       # list[float] or np.ndarray(N,), optional

# Default: no grouping (uses all hits as-is)
stats = process_sensor_data(sensor_times, sensor_charges)  # returns: np.ndarray, shape (9,)

# Optional: group hits within time windows
stats = process_sensor_data(sensor_times, sensor_charges, grouping_window_ns=2.0)

Summary Statistics

Computes summary statistics for neutrino telescope sensors as described in the IceCube paper. All functions return numpy arrays with statistics in the following order:

Select the statistic set with the mode argument: "standard" (default, 9), "minimal" (4), or "extended" (25).

Minimal Statistics (4 stats)

A lightweight subset of the standard statistics (indices 0, 3, 7, 8). Because these are permutation-invariant, the input is not sorted unless grouping is requested.

stats = compute_summary_stats(times, charges, mode="minimal")  # shape: (4,)

# Array indices:
stats[0]  # total_charge: Total charge collected
stats[1]  # first_pulse_time: Time of first pulse
stats[2]  # charge_weighted_mean_time: Charge-weighted mean time
stats[3]  # charge_weighted_std_time: Charge-weighted standard deviation

Standard Statistics (9 stats, default)

stats = compute_summary_stats(times, charges)  # shape: (9,)

# Array indices:
stats[0]  # total_charge: Total charge collected
stats[1]  # charge_100ns: Charge within 100ns of first pulse
stats[2]  # charge_500ns: Charge within 500ns of first pulse
stats[3]  # first_pulse_time: Time of first pulse
stats[4]  # last_pulse_time: Time of last pulse
stats[5]  # charge_20_percent_time: Time at which 20% of charge is collected
stats[6]  # charge_50_percent_time: Time at which 50% of charge is collected
stats[7]  # charge_weighted_mean_time: Charge-weighted mean time
stats[8]  # charge_weighted_std_time: Charge-weighted standard deviation

Extended Statistics (25 stats, optional)

Pass mode="extended" to compute 16 additional statistics:

stats = compute_summary_stats(times, charges, mode="extended")  # shape: (25,)

# Includes all 9 standard statistics above, plus:
stats[9]   # charge_5_percent_time: Time at which 5% of charge is collected
stats[10]  # charge_10_percent_time: Time at which 10% of charge is collected
stats[11]  # charge_25_percent_time: Time at which 25% of charge is collected
stats[12]  # charge_75_percent_time: Time at which 75% of charge is collected
stats[13]  # charge_90_percent_time: Time at which 90% of charge is collected
stats[14]  # charge_95_percent_time: Time at which 95% of charge is collected
stats[15]  # charge_10ns: Charge within 10ns of first pulse
stats[16]  # charge_20ns: Charge within 20ns of first pulse
stats[17]  # charge_50ns: Charge within 50ns of first pulse
stats[18]  # charge_200ns: Charge within 200ns of first pulse
stats[19]  # charge_1000ns: Charge within 1000ns of first pulse
stats[20]  # charge_2000ns: Charge within 2000ns of first pulse
stats[21]  # n_pulses: Number of input pulses (pre-grouping count)
stats[22]  # q_max_frac: Peak charge fraction (max pulse charge / total charge)
stats[23]  # n_string_neighbors: HLC-style neighbor count (0 outside process_event)
stats[24]  # t_skewness: Charge-weighted time skewness (0 for < 3 pulses)

API

compute_summary_stats(times, charges, mode="standard")

Args:

• times: np.ndarray or list, shape (N,) - pulse arrival times in ns
• charges: np.ndarray or list, shape (N,) - pulse charges
• mode: str - "standard" (default, 9), "minimal" (4), or "extended" (25)

Returns: np.ndarray, shape (4,), (9,), or (25,) - array with summary statistics in order shown above

process_event(event_data, grouping_window_ns=None, mode="standard")

Args:

• event_data: dict holding photon-level arrays (either provide a top-level photons dictionary or store the required fields directly) with keys sensor_pos_x, sensor_pos_y, sensor_pos_z, string_id, sensor_id, t, and optional charge
• grouping_window_ns: float or None - time window for grouping hits (default: None, no grouping)
• mode: str - "standard" (default, 9), "minimal" (4), or "extended" (25) per sensor

Returns: tuple[np.ndarray, np.ndarray]

• sensor_positions: np.ndarray, shape (N_sensors, 3) - sensor positions
• sensor_stats: np.ndarray, shape (N_sensors, 4 | 9 | 25) - statistics for each sensor (aligned with positions)

process_sensor_data(sensor_times, sensor_charges=None, grouping_window_ns=None, mode="standard")

Args:

• sensor_times: np.ndarray or list, shape (N,) - hit times for sensor
• sensor_charges: np.ndarray or list, shape (N,) - hit charges (optional, defaults to 1.0)
• grouping_window_ns: float or None - time window for grouping hits (default: None, no grouping)
• mode: str - "standard" (default, 9), "minimal" (4), or "extended" (25)

Returns: np.ndarray, shape (4,), (9,), or (25,) - array with summary statistics in order shown above

License

MIT