Metadata-Version: 2.4
Name: nt_summary_stats
Version: 1.2
Summary: Fast computation of traditional summary statistics for neutrino telescopes
Author-email: "Felix J. Yu" <felixyu@g.harvard.edu>
Project-URL: Homepage, https://github.com/felixyu7/nt_summary_stats
Project-URL: Repository, https://github.com/felixyu7/nt_summary_stats
Project-URL: Issues, https://github.com/felixyu7/nt_summary_stats/issues
Keywords: neutrino,telescope,summary,statistics,icecube,physics
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Science/Research
Classifier: Topic :: Scientific/Engineering :: Physics
Classifier: License :: OSI Approved :: MIT License
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.8
Classifier: Programming Language :: Python :: 3.9
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: C++
Requires-Python: >=3.8
Description-Content-Type: text/markdown
Requires-Dist: numpy>=1.20.0
Provides-Extra: dev
Requires-Dist: pytest>=6.0; extra == "dev"
Requires-Dist: pytest-cov; extra == "dev"
Requires-Dist: black; extra == "dev"
Requires-Dist: flake8; extra == "dev"
Requires-Dist: pybind11>=2.11; extra == "dev"

# 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](https://arxiv.org/abs/2101.11589).

## Install

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

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

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

```python
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](https://arxiv.org/abs/2101.11589). 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.

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

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

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