"""Calculate NSB trigger rates from reduced event-data HDF5 files."""
import logging
import re
from collections import defaultdict
from pathlib import Path
import astropy.units as u
import numpy as np
from astropy.table import Table
from simtools.io.file_type import product_value
from simtools.utils import general
_logger = logging.getLogger(__name__)
_THRESHOLD_RE = re.compile(r"_(?:[ad]sum|discriminator)(?P<threshold>\d+(?:\.\d+)?)(?=\.)")
_RUN_NUMBER_RE = re.compile(r"run(?P<run_number>\d+)", re.IGNORECASE)
_DEFAULT_REDUCED_EVENT_DATA_PATTERN = (
f"*{product_value('sim_telarray_event_data', 'default_suffix')}"
)
[docs]
def parse_nsb_hdf5_file(file_path):
"""
Parse one NSB reduced event-data HDF5 file.
Parameters
----------
file_path : Path or str
Path to the HDF5 file.
Returns
-------
dict or None
Dictionary with keys: 'run', 'threshold', 'triggers', 'events', 'file_path'.
Returns None if critical information cannot be extracted.
"""
file_path = Path(file_path)
try:
events_table = Table.read(file_path, path="SHOWERS")
except OSError, ValueError, KeyError:
_logger.exception(f"Failed to read SHOWERS table from {file_path}")
return None
try:
triggers_table = Table.read(file_path, path="TRIGGERS")
triggers = len(triggers_table)
except OSError, ValueError, KeyError:
_logger.info(f"No TRIGGERS table found in {file_path}; using 0 triggers")
triggers = 0
run_number = extract_run_number(file_path)
threshold = extract_threshold(file_path)
events = len(events_table)
if run_number is None or threshold is None:
_logger.warning(
f"Skipping {file_path}: missing critical info "
f"(run={run_number}, threshold={threshold}, triggers={triggers})"
)
return None
return {
"run": run_number,
"threshold": threshold,
"triggers": triggers,
"events": events,
"file_path": str(file_path),
}
[docs]
def parse_nsb_hdf5_files(file_list):
"""
Parse multiple NSB reduced event-data HDF5 files.
Parameters
----------
file_list : list of Path
List of HDF5 file paths to parse.
Returns
-------
list of dict
List of parsed data dictionaries. Failed parses are excluded.
Raises
------
ValueError
If no files could be parsed successfully.
"""
data = []
for file_path in file_list:
_logger.debug(f"Parsing {file_path}")
result = parse_nsb_hdf5_file(file_path)
if result is not None:
data.append(result)
if not data:
raise ValueError("No HDF5 files could be parsed successfully")
_logger.info(f"Parsed {len(data)} out of {len(file_list)} HDF5 files")
return data
[docs]
def find_hdf5_files(
root_dir,
pattern=_DEFAULT_REDUCED_EVENT_DATA_PATTERN,
):
"""Recursively find reduced event-data HDF5 files matching the pattern."""
return sorted(general.resolve_file_patterns(Path(root_dir) / "**" / pattern))
[docs]
def group_by_threshold_and_run(data):
"""
Organize parsed data into a structure grouped by threshold and run.
Parameters
----------
data : list of dict
List of parsed HDF5 file data.
Returns
-------
dict
Nested dictionary: {threshold: {run: {'triggers': int, 'events': int}}}
"""
grouped = defaultdict(dict)
source_files = defaultdict(dict)
for entry in data:
threshold = entry["threshold"]
run = entry["run"]
if run in grouped[threshold]:
previous = source_files[threshold][run]
raise ValueError(
f"Duplicate NSB input for threshold {threshold}, run {run}: "
f"{previous} and {entry.get('file_path', '<unknown>')}"
)
grouped[threshold][run] = {
"triggers": entry["triggers"],
"events": entry["events"],
}
source_files[threshold][run] = entry.get("file_path", "<unknown>")
return dict(grouped)
[docs]
def calculate_statistics(grouped_data, time_window):
"""
Calculate statistics for each threshold.
Parameters
----------
grouped_data : dict
Data grouped by threshold and run.
time_window : float, optional
Time window per event in seconds.
Returns
-------
dict
Statistics for each threshold with keys:
- 'runs': dict of {run_number: triggers}
- 'total_triggers': sum of all triggers
- 'total_events': sum of all events
- 'time_s': total time in seconds
- 'rate_hz': trigger rate in Hz
- 'rate_khz': trigger rate in kHz
- 'error_hz': exposure-weighted standard error of run rates in Hz
- 'num_runs': number of runs
"""
time_window = _time_window_quantity(time_window)
if not np.isfinite(time_window.value) or time_window <= 0 * u.s:
raise ValueError("time_window must be finite and positive.")
statistics = {}
for threshold, runs_data in grouped_data.items():
run_statistics = _collect_valid_run_statistics(runs_data, threshold, time_window)
statistics[threshold] = _calculate_threshold_statistics(run_statistics, time_window)
return statistics
def _time_window_quantity(time_window):
"""Return the time window as a duration in seconds."""
if not isinstance(time_window, u.Quantity):
time_window = time_window * u.s
return time_window.to(u.s)
def _collect_valid_run_statistics(runs_data, threshold, time_window):
"""Collect exposure and rate values for runs with valid event counts."""
run_statistics = {
"runs": {},
"triggers": [],
"events": [],
"rates": [],
"exposures": [],
}
for run_num, run_info in sorted(runs_data.items()):
triggers = run_info["triggers"]
events = run_info["events"]
if events is None or not np.isfinite(events) or events <= 0:
reason = (
"missing event count" if events is None else "nonpositive or non-finite exposure"
)
_logger.warning(f"Skipping run {run_num} for threshold {threshold}: {reason}")
continue
exposure = events * time_window
run_statistics["runs"][run_num] = triggers
run_statistics["triggers"].append(triggers)
run_statistics["events"].append(events)
run_statistics["rates"].append(triggers / exposure)
run_statistics["exposures"].append(exposure)
return run_statistics
def _calculate_threshold_statistics(run_statistics, time_window):
"""Calculate aggregate rate statistics for one threshold."""
triggers = run_statistics["triggers"]
events = run_statistics["events"]
total_triggers = np.sum(triggers)
total_events = np.sum(events) if events else 0
time_s = total_events * time_window if total_events > 0 else 0 * u.s
rate = total_triggers / time_s if time_s > 0 else 0 * u.Hz
return {
"runs": run_statistics["runs"],
"total_triggers": int(total_triggers),
"total_events": int(total_events),
"time_s": time_s.to_value(u.s),
"rate_hz": rate.to_value(u.Hz),
"rate_khz": rate.to_value(u.kHz),
"error_hz": _calculate_rate_error(run_statistics, time_s, rate),
"num_runs": len(triggers),
}
def _calculate_rate_error(run_statistics, time_s, weighted_mean):
"""Calculate exposure-weighted standard error of the valid run rates."""
if len(run_statistics["rates"]) < 2 or time_s <= 0:
return 0
exposures = u.Quantity(run_statistics["exposures"]).to(u.s)
rates = u.Quantity(run_statistics["rates"]).to(u.Hz)
weight_sum = np.sum(exposures)
effective_runs = weight_sum**2 / np.sum(exposures**2)
variance_denominator = weight_sum - np.sum(exposures**2) / weight_sum
if variance_denominator <= 0:
return 0
weighted_variance = np.sum(exposures * (rates - weighted_mean) ** 2)
weighted_variance /= variance_denominator
return np.sqrt(weighted_variance / effective_runs).to_value(u.Hz)
[docs]
def generate_ecsv_output(statistics, output_file, time_window):
"""
Generate ECSV table from statistics.
Parameters
----------
statistics : dict
Statistics dictionary from calculate_statistics().
output_file : Path or str
Output file path for ECSV table.
time_window : float, optional
Time window per event in seconds for metadata.
"""
output_file = Path(output_file)
if not statistics:
raise ValueError("No statistics to write")
all_runs = set()
for threshold_stats in statistics.values():
all_runs.update(threshold_stats["runs"].keys())
all_runs = sorted(all_runs)
threshold_col = []
run_cols = {run: [] for run in all_runs}
total_triggers_col = []
total_events_col = []
time_col = []
rate_hz_col = []
rate_khz_col = []
error_hz_col = []
for threshold in sorted(statistics.keys()):
stats = statistics[threshold]
threshold_col.append(threshold)
for run in all_runs:
if run in stats["runs"]:
run_cols[run].append(stats["runs"][run])
else:
run_cols[run].append(np.nan)
total_triggers_col.append(stats["total_triggers"])
total_events_col.append(stats["total_events"])
time_col.append(stats["time_s"])
rate_hz_col.append(stats["rate_hz"])
rate_khz_col.append(stats["rate_khz"])
error_hz_col.append(stats["error_hz"])
table_data = {"threshold": threshold_col}
for run in all_runs:
table_data[f"run{run}"] = run_cols[run]
table_data["Total trig'd"] = total_triggers_col
table_data["Out of"] = total_events_col
table_data["Time (s)"] = time_col
table_data["Rate (Hz)"] = rate_hz_col
table_data["Rate (kHz)"] = rate_khz_col
table_data["Error (Hz)"] = error_hz_col
table = Table(table_data)
total_events_sum = sum(total_events_col)
table.meta["comments"] = [
f"Total events: {total_events_sum}",
f"Time window per event: {time_window} s",
]
for col in ["Time (s)", "Rate (Hz)", "Rate (kHz)", "Error (Hz)"]:
table[col].format = ".2f"
output_file.parent.mkdir(parents=True, exist_ok=True)
table.write(output_file, format="ascii.ecsv", overwrite=True)
_logger.info(f"ECSV table written to {output_file}")
_logger.info(f"Table contains {len(threshold_col)} thresholds and {len(all_runs)} runs")
[docs]
def derive_nsb_triggers(args):
"""
Derive NSB trigger rates from reduced event-data HDF5 files.
Parameters
----------
args : dict
Configuration parameters with keys:
- root_dir: Root directory to search for reduced event-data HDF5 files
- pattern: Glob pattern for HDF5 files, default: ``*.reduced_event_data.hdf5``
- output: Output ECSV file path (optional, if None, no file is written)
- time_window: Time window per event in seconds (required; telescope-dependent)
- verbose: Enable verbose logging (optional)
Returns
-------
dict
Statistics dictionary with threshold as keys, containing:
- rate_hz: trigger rate in Hz
- rate_khz: trigger rate in kHz
- error_hz: standard error in Hz
- total_triggers: total number of triggers
- num_runs: number of runs
Raises
------
FileNotFoundError
If root directory doesn't exist or no matching HDF5 files are found.
ValueError
If time_window is missing/invalid or if no files could be parsed successfully.
"""
pattern = args.get("pattern", _DEFAULT_REDUCED_EVENT_DATA_PATTERN)
time_window = args.get("time_window")
if time_window is None:
raise ValueError(
"Missing required argument 'time_window'. "
"Provide telescope-specific time window in seconds."
)
try:
time_window = float(time_window)
except (TypeError, ValueError) as exc:
raise ValueError("Argument 'time_window' must be a number.") from exc
if not np.isfinite(time_window) or time_window <= 0:
raise ValueError("Argument 'time_window' must be finite and > 0.")
_logger.info("NSB Trigger Rate Calculator")
_logger.info(f"Root directory: {args['root_dir']}")
_logger.info(f"Time window: {time_window} s")
if args.get("output"):
_logger.info(f"Output file: {args['output']}")
_logger.info("Searching for reduced event-data HDF5 files...")
hdf5_files = find_hdf5_files(args["root_dir"], pattern)
_logger.info(f"Found {len(hdf5_files)} HDF5 files")
_logger.info("Parsing HDF5 files...")
parsed_data = parse_nsb_hdf5_files(hdf5_files)
_logger.info(f"Successfully parsed {len(parsed_data)} HDF5 files")
_logger.info("Grouping data by threshold and run...")
grouped_data = group_by_threshold_and_run(parsed_data)
thresholds = sorted(grouped_data.keys())
_logger.info(f"Found {len(thresholds)} thresholds: {thresholds}")
_logger.info("Calculating statistics...")
statistics = calculate_statistics(grouped_data, time_window)
_logger.info("Summary:")
for threshold in sorted(statistics.keys()):
stats = statistics[threshold]
_logger.info(
f" Threshold {threshold}: "
f"{stats['num_runs']} runs, "
f"{stats['total_triggers']} triggers, "
f"{stats['rate_hz']:.2f} Hz"
)
if args.get("output"):
_logger.info("Writing ECSV output...")
generate_ecsv_output(statistics, args["output"], time_window)
_logger.info(f"Output written to: {args['output']}")
_logger.info("NSB trigger rate calculation completed")
return statistics
