# The Clear BSD License
#
# Copyright (c) 2026 Tobias Heibges
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted (subject to the limitations in the disclaimer
# below) provided that the following conditions are met:
#
# * Redistributions of source code must retain the above copyright notice,
# this list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright
# notice, this list of conditions and the following disclaimer in the
# documentation and/or other materials provided with the distribution.
#
# * Neither the name of the copyright holder nor the names of its
# contributors may be used to endorse or promote products derived from this
# software without specific prior written permission.
#
# NO EXPRESS OR IMPLIED LICENSES TO ANY PARTY'S PATENT RIGHTS ARE GRANTED BY
# THIS LICENSE. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
# CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
# CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
# BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
# IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
"""
Detector Analysis Utilities
Provides utility functions for analyzing detection events, including
angular and temporal distribution histograms.
Examples
--------
>>> from surface_roughness.detectors import (
... compute_angular_distribution,
... compute_time_distribution
... )
>>>
>>> angles, counts = compute_angular_distribution(
... detector.events,
... reference_direction=np.array([0, 0, 1])
... )
"""
import numpy as np
from numpy.typing import NDArray
from .base import DetectionEvent
[docs]
def compute_angular_distribution(
events: list[DetectionEvent],
reference_direction: NDArray[np.float32],
num_bins: int = 50,
) -> tuple[NDArray[np.float64], NDArray[np.int64]]:
"""
Compute angular distribution histogram.
Parameters
----------
events : list of DetectionEvent
Detection events to analyze.
reference_direction : ndarray, shape (3,)
Reference direction for angle calculation.
num_bins : int, optional
Number of histogram bins. Default is 50.
Returns
-------
bin_centers : ndarray, shape (num_bins,)
Bin centers in degrees.
counts : ndarray, shape (num_bins,)
Number of events in each bin.
Examples
--------
>>> angles, counts = compute_angular_distribution(
... detector.events,
... reference_direction=np.array([0, 0, 1]),
... num_bins=90 # 2-degree bins
... )
>>> import matplotlib.pyplot as plt
>>> plt.bar(angles, counts, width=2)
>>> plt.xlabel('Angle (degrees)')
>>> plt.ylabel('Count')
"""
if len(events) == 0:
return np.array([], dtype=np.float64), np.array([], dtype=np.int64)
ref = reference_direction / np.linalg.norm(reference_direction)
angles = []
for event in events:
dir_norm = event.direction / np.linalg.norm(event.direction)
cos_angle = np.dot(dir_norm, ref)
cos_angle = np.clip(cos_angle, -1.0, 1.0)
angles.append(np.degrees(np.arccos(cos_angle)))
angles = np.array(angles)
counts, bin_edges = np.histogram(angles, bins=num_bins, range=(0, 180))
bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2
return bin_centers, counts
[docs]
def compute_time_distribution(
events: list[DetectionEvent],
num_bins: int = 50,
time_range: tuple[float, float] | None = None,
) -> tuple[NDArray[np.float64], NDArray[np.int64]]:
"""
Compute arrival time distribution histogram.
Parameters
----------
events : list of DetectionEvent
Detection events to analyze.
num_bins : int, optional
Number of histogram bins. Default is 50.
time_range : tuple of float, optional
(min, max) time range for histogram. If None, uses data range.
Returns
-------
bin_centers : ndarray, shape (N,)
Bin centers in seconds.
counts : ndarray, shape (N,)
Number of events in each bin.
Notes
-----
Automatically handles cases where all times are very similar by
reducing the number of bins to avoid empty histograms.
Examples
--------
>>> times, counts = compute_time_distribution(
... detector.events,
... num_bins=100
... )
>>> print(f"Time spread: {times.max() - times.min():.3e} s")
"""
if len(events) == 0:
return np.array([], dtype=np.float64), np.array([], dtype=np.int64)
times = np.array([e.time for e in events])
if time_range is None:
time_range = (times.min(), times.max())
# Handle case where all times are very similar
time_span = time_range[1] - time_range[0]
if time_span < 1e-12: # Less than 1 picosecond spread
return np.array([times.mean()]), np.array([len(times)])
# Adjust bins if time span is very small
effective_bins = min(num_bins, max(1, int(time_span * 1e12)))
counts, bin_edges = np.histogram(times, bins=effective_bins, range=time_range)
bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2
return bin_centers, counts
[docs]
def compute_intensity_distribution(
events: list[DetectionEvent],
num_bins: int = 50,
) -> tuple[NDArray[np.float64], NDArray[np.int64]]:
"""
Compute intensity distribution histogram.
Parameters
----------
events : list of DetectionEvent
Detection events to analyze.
num_bins : int, optional
Number of histogram bins. Default is 50.
Returns
-------
bin_centers : ndarray, shape (num_bins,)
Bin centers (intensity values).
counts : ndarray, shape (num_bins,)
Number of events in each bin.
"""
if len(events) == 0:
return np.array([], dtype=np.float64), np.array([], dtype=np.int64)
intensities = np.array([e.intensity for e in events])
counts, bin_edges = np.histogram(intensities, bins=num_bins)
bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2
return bin_centers, counts
[docs]
def compute_wavelength_distribution(
events: list[DetectionEvent],
num_bins: int = 50,
) -> tuple[NDArray[np.float64], NDArray[np.int64]]:
"""
Compute wavelength distribution histogram.
Parameters
----------
events : list of DetectionEvent
Detection events to analyze.
num_bins : int, optional
Number of histogram bins. Default is 50.
Returns
-------
bin_centers : ndarray, shape (num_bins,)
Bin centers in meters.
counts : ndarray, shape (num_bins,)
Number of events in each bin.
"""
if len(events) == 0:
return np.array([], dtype=np.float64), np.array([], dtype=np.int64)
wavelengths = np.array([e.wavelength for e in events])
counts, bin_edges = np.histogram(wavelengths, bins=num_bins)
bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2
return bin_centers, counts
[docs]
def compute_statistics(events: list[DetectionEvent]) -> dict:
"""
Compute summary statistics for detection events.
Parameters
----------
events : list of DetectionEvent
Detection events to analyze.
Returns
-------
stats : dict
Dictionary containing:
- count: number of events
- total_intensity: sum of intensities
- mean_time: average arrival time
- std_time: arrival time standard deviation
- min_time: earliest arrival
- max_time: latest arrival
- mean_wavelength: average wavelength
- time_spread: max_time - min_time
Examples
--------
>>> stats = compute_statistics(detector.events)
>>> print(f"Detected {stats['count']} rays")
>>> print(f"Total intensity: {stats['total_intensity']:.3e}")
>>> print(f"Time spread: {stats['time_spread']:.3e} s")
"""
if len(events) == 0:
return {
"count": 0,
"total_intensity": 0.0,
"mean_time": 0.0,
"std_time": 0.0,
"min_time": 0.0,
"max_time": 0.0,
"mean_wavelength": 0.0,
"time_spread": 0.0,
}
times = np.array([e.time for e in events])
intensities = np.array([e.intensity for e in events])
wavelengths = np.array([e.wavelength for e in events])
return {
"count": len(events),
"total_intensity": float(np.sum(intensities)),
"mean_time": float(np.mean(times)),
"std_time": float(np.std(times)),
"min_time": float(np.min(times)),
"max_time": float(np.max(times)),
"mean_wavelength": float(np.mean(wavelengths)),
"time_spread": float(np.max(times) - np.min(times)),
}