Science verification analysis

[zsselcuk]

This PR introduces the first version of the Science Verification Analysis, comprising multiple tests to verify the health and performance of the RNO-G stations. The main version uses the monitoring.root files to extract data from selected runs or time periods. We're keeping the outdated analysis based on the dataProviderRNOG module for reference, but it will no longer be maintained or updated. The analysis produces both a concise summary of the individual channel health and a more detailed set of diagnostic plots & test results in case a channel needs further investigation. It also provides a summary of any unreadable root files encountered during the analysis, which can help detect data storing issues early enough.
The analysis consists of the following tests (all applied on the FORCE trigger data):

• Spectral analysis test to determine if surface channels detect the galactic signal,
• The stability of the SNR values is validated through a z-score analysis (currently log-based, might be changed), which flags the outliers using channel-specific reference SNR values and reports how significant the deviation from the reference is. The reference values can be calculated for each channel/station using the dedicated script present under /expected_values,
• The KDE modality tests are performed on the RMS values to determine if the RMS distributions are unimodal, reports on the skewness/tail characterization as well in the detailed reports,
• The stability of the RMS values over time is tested again through a z-score test, where the outliers are flagged using channel-specific reference RMS values. The relative median shift between each run is also calculated and used to report on the RMS stability in addition to the z-score test. The reference values can be calculated for each channel/station using the dedicated script present under /expected_values,
• Performs a glitching analysis and determines the glitching health using a binomial test,
• Performs a block offset analysis using reference block offset values and reports about it (which is not accounted for in the overall channel health decision)
• Plots the trigger rates

The overall channel health combines these tests. More tests will be implemented in the future and existing tests will be improved.

[fschlueter]

Is that something which can be either defined globally or in a config file? Also I was wondering if it makes sense to generalize this a bit by only defining a dict with keys which act as label for the channel group and values which list the channel id's, e.g.:

channel_groups = {
   "phased_array": [0, 1, 2, 3],
   "surface": [12, 13, ...],
    ....
}

[zsselcuk]

Yes, I think I should put all of the station-specific stuff into a config file, including the analysis. This way it would also be easier for station experts to change some parameters if needed.

[fschlueter]

Can you specify that in an extra file? Makes the scripts easier to read.

[fschlueter]

can you combine these json files?

[zsselcuk]

Sure! For each station there will be only a single file.

[fschlueter]

Love that you gave it a doc-string. Can you move it up to the top of the file?

[fschlueter]

you stored these files in the repo - was that intended? Are those files meant to be created or used by other users? Does it make sense to add more meta data to the json files themself? Examples: station_id, timeperiod used (maybe with a list of runs, ...).

[zsselcuk]

So this script is only meant to be used by the station experts to calculate the expected snr values per station. The script then automatically saves these files in the repo so that the science_verification_analysis.py can use them to perform the snr analysis.
I think adding more metadata would be very nice, so that it is reproducably and we can also track if there is a change in the expected values over time.

[fschlueter]

This is a pretty monolitic file. I think it would improve readability to split the code up in several files. I gave a few suggestion in the other file. Maybe we can separate utility function and function calculating physics observables?

[zsselcuk]

Yes, we should definitely do that! Otherwise, it will be too complicated in the future.

[fschlueter]

do you want to use logging?

[zsselcuk]

Yes, I think that's a good idea

[zsselcuk]

Things to do:

• Add trigger rate analysis
• Complete doc strings
• Implement z-score analysis for RMS
• Test the RMS
• Freq vs time analysis for spectra
• Correct units
• Combine json files and add info on the time interval etc.
• Push the expected snr values for json
• Function to check GPS/file naming

[fschlueter]

Can you give this PR a description with a short overview over the analysis and steps involved?

[zsselcuk]

Updated the description

[fschlueter]

I like that (having a default and than station specific overrides)! Just thinking a head the new stations will all be in the 14 design. Maybe it is possible to have a dict or dicts of overrides, e.g., station_2024_geometry, and than a list/dict of stations and their overrides (blank if they fit the default)

[fschlueter]

hö, is that still necessary, what is with the previous json file?

[fschlueter]

Would it be easier / more readable if we put all these configs (block_offsets, glitching, rms, specral_analysis) in one json?

[fschlueter]

should we use the rnog idoem of specifying the data dir with an env variable ($RNO_G_DATA) (or via the cmd line)

[fschlueter]

Do we need the __init__ file at the top level? That would be only necessary if we import functions/modules from the sva elsewhere - is that the plan?

[fschlueter]

Hey zeynep, so I am running the analysis right now and will report (piece by piece) my experiences. First thing noted, we need the rnog_data repo (should talk to brian to get that into cvmfs)

[fschlueter]

Another thing, I think it would be good to store all output per executing. For example, I am executing python science_verification_analysis_main.py -st 14 --time_range 2026-07-16 2026-07-20. It would be nice if all results (plots, log file, text files, ...) are stored under the following path:
/some/base/directory/YY-MM-DD_station-SS_runXXXXXX-runXXXXXXX_RRRRRR The YY, MM, DD are for year, month and day, the SS for the station number, the XXX for the first and last run id and the RRRRRRR stand for a random 6 charather string (this allows to run it several time on the same data and compare). The idea behind this is the following: We are building a history in a central location and have versioning (in time) of the results. We could run the analysis automatically once per day/week, store the data there (in a directory people could access from their browser!). So shifters just have to inspect those directories (from the browser for example) compile the report or alert us.

What do you think?

[fschlueter]

I think we also have to make a prioritisation of the plots and maybe hide all other plots in a sub folder. For example we could have a few "standard plots" + plots which are associated with a failed condition (those which show a problem found by your analysis) and than all the rest. Separated by subdirectories, or the last one put in a subdirectory

[fschlueter]

What are those heatmaps, are they helpful?

[fschlueter]

We also need a top-level README.md which describes all the different directories + main function in detail.

[fschlueter]

This is an AI-generated suggestion:

Science Verification Analysis (SVA)

Purpose

The Science Verification Analysis checks whether an RNO-G station is behaving as expected, using a set of automated diagnostic tests run over a chosen set of runs. Each test looks at a different failure mode (galactic-noise sensitivity, SNR/RMS stability, channel glitching, ADC block offsets, trigger rates) and reports a per-channel health verdict (OK, !!, or X). The output is a concise CSV summary for every channel plus a set of detailed plots and text/JSON reports for anyone who needs to dig into a flagged channel.

The analysis is meant to be run regularly (e.g. after a station comes back online, or periodically to catch slow degradation) and to be usable by anyone on the collaboration, not just the analysis' original author — channel maps, thresholds, and plotting styles all live in config files rather than in code, and reference ("expected") values are calculated once per station from a known-good time period and then reused.

There are two ways to read the underlying data:

• monitoring.root (current, recommended): reads pre-computed per-event summary quantities (RMS, max amplitude, glitching test statistic, block offsets, spectra) from the station's monitoring stream. Fast, and works for all data taken with the monitoring pipeline.
• dataProviderRNOG (outdated/legacy): reads full waveforms via NuRadioReco for older runs that predate monitoring.root. Kept only for archival purposes — see the [Legacy path](#sva_dataproviderrnog-legacy) section below for its caveats.

Quick start / usage example

Run the main analysis for station 14 over an explicit list of runs:

python science_verification_analysis_main.py --station_id 14 --runs 260080 260090 260100

Or over a run range, excluding a couple of bad runs:

python science_verification_analysis_main.py --station_id 14 --run_range 260080 260150 --exclude-runs 260101 260102

Or over a date range (run numbers are looked up from the RNO-G run table):

python science_verification_analysis_main.py --station_id 14 --time_range 2026-04-26 2026-05-02

Add --debug_plot to also produce extra diagnostic plots, and --save_location /path/to/plots to redirect where plots are written (defaults to a plots/ directory next to the script). --runs, --run_range, and --time_range are mutually exclusive — pick exactly one way to select runs.

This produces, under the script directory by default:

• channel_health_summary/validation_summary_station<ID>_<run_label>.csv — the top-level per-channel verdict table.
• detailed_results/ — per-test text/JSON files with the numbers behind each verdict (spectral results, SNR/RMS outlier details, glitching stats, block offset stats, failed-run report).
• plots/ — all the diagnostic and debug plots.
• logs/logging_science_verification_analysis_station<ID>_<run_label>.log — a full log of the run, including warnings about missing files, invalid timestamps, and borderline channels.

Analysis overview

For a given station and set of runs, science_verification_analysis_main.py runs the following steps in order:

1. Read data — read_multiple_runs() reads monitoring.root and headers.root for every requested run, validates them against each other (matching event numbers, station/run IDs, non-overlapping trigger types), and concatenates everything into one combined dataset. Runs with missing or inconsistent files are skipped and reported rather than crashing the whole analysis.
2. Spectral (galactic noise) test — normalizes surface-channel spectra to a reference band and checks, in several frequency bands, whether upward-facing channels show more galactic excess than downward-facing ones (as expected).
3. SNR stability (z-score) test — compares each channel's log-SNR distribution for the current runs against a previously-computed reference (mean/std/k-value per channel) and flags statistically significant outlier events.
4. RMS/Vrms modality & tail test — uses a KDE of the RMS distribution per channel to check it's unimodal (not bimodal/flat, which would suggest a mis-biased or noisy channel), and separately checks for excess skew/tails.
5. RMS stability test — same z-score approach as SNR, plus a relative-median-shift metric across runs, combined into an overall stability decision per channel.
6. Glitching test — a one-sided binomial test on how often each channel's glitching test statistic is triggered, against an expected background rate.
7. Block offset test — summarizes each channel's ADC block-offset statistics (mean/median/std/IQR/P99) and compares them to reference values; informational only (not counted in the overall channel-health verdict).
8. Trigger rate plots — plots trigger rates over time per trigger type, no pass/fail verdict.
9. Summary CSV — create_result_csv_file() combines the SNR, spectral, RMS, RMS-stability, and glitching verdicts into one OK/!!/X per channel and writes the summary table.

Steps 2–7 each have their own config file (thresholds, band definitions, KDE parameters, etc. — see below) so they can be tuned without touching analysis code, and their own reference/"expected value" script under expected_values/ to (re-)generate reference numbers for a new station or after a known-good re-calibration period.

Directory structure

science_verification_analysis_main.py

The entry point described above. Also home to setup_logging() and the output-directory constants (PLOTS_DIR, RESULTS_DIR, CSV_DIR, LOGS_DIR, REFERENCE_DIR, CONFIG_DIR), all resolved relative to the script's own location so it can be run from anywhere.

config_files_sva/

All tunable parameters, one JSON file per topic (plus one small Python helper). Editing these does not require touching analysis code:

• config_station.json — channel maps (all_channels, surface_channels, deep_channels, upward_channels, downward_channels, vpol_channels, hpol_channels, phased_array_channels, reference_channels, reference_channels_galaxy) under a default_config, with a station_specific_adjustments block keyed by station ID (currently station 14 has a different channel layout).
• config_helper.py — get_station_config(station_id, default_config, station_specific_adjustments) merges the default config with any station-specific overrides. This is the only "code" in this directory; everything else is data.
• config_spectral_analysis.json — frequency bands used for the galactic-noise test (galactic_excess, plus a few RFI bands), the normalization band, and the significance thresholds (alpha_spec, ci_threshold_spec, log_ratio_thresholds_spec).
• config_rms.json — parameters for the KDE modality test (bandwidth, grid_points, peak_prominence, height_threshold), the tail/skewness test (percentiles, extreme_k, minimum event count), and the human-readable reporting thresholds (strong_skew, extreme_skew, etc.).
• config_glitching.json — the binomial test parameters for the glitching test (alpha, expected background pvalue, confidence level, and the CI thresholds that separate "weak"/"moderate"/"strong" excessive glitching).
• config_block_offsets.json — acceptable median/IQR block-offset limits.
• config_plotting.py — set_plot_style() centralizes all matplotlib rcParams (fonts, tick sizes, line widths, dpi) plus a shared 24-color palette (COLORS), so every plot in the analysis looks consistent.

analysis_functions_sva/

The statistics/physics behind each test, split by topic and independent of I/O:

• spectral_analysis_sva.py — normalize_channels(), find_amplitude_ratio_in_band[_specific_bkg](), excess_info_from_ratio[_specific_bkg](), validate_excess_in_bands(). Computes upward/downward amplitude ratios per frequency band and classifies them as NO/WEAK/MODERATE/STRONG excess (method differs slightly for monitoring vs. dataProviderRNOG data).
• z_score_analysis_sva.py — shared statistics machinery used by both the SNR and RMS tests: log-parameter statistics, z-scores against a reference mean/std, k-value derivation (find_k_value), outlier flagging/detail extraction, saving/loading reference values to/from JSON (including the metadata block), and a rolling-window z-score variant.
• vrms_analysis_sva.py — calculate_vrms() (from raw waveforms, dataProviderRNOG path) and get_rms_per_trigger_monitoring() (from precomputed RMS, monitoring path) both split values by trigger type; kde_modality() and tail_fraction_and_trimmed_skew_two_sided() implement the modality/tail tests; report_vrms_characteristics() turns those into human-readable labels.
• vrms_stability_analysis_sva.py — get_rms_per_run(), relative_median_shift() (pairwise relative shift in median RMS between runs), and decision_metric(), which combines outlier fraction, largest z-score excess, and median-shift into the final OK/!!/X RMS-stability verdict. Also contains a block of commented-out, currently-unused helpers (Wasserstein-distance and linear-regression-based stability metrics) kept for possible future use.
• glitching_analysis_sva.py — binomtest_glitch_fraction(), a one-sided binomial test per channel classifying glitching as NO/WEAK/MODERATE/STRONG excessive.
• block_offsets_analysis_sva_monitoring.py — get_force_block_offsets_monitoring(), block_offset_statistics_monitoring() (mean/median/std/IQR/P99/P95), and a violin-plot helper, for the monitoring.root path.
• block_offsets_analysis_sva_dataproviderrnog.py — equivalent block-offset functions for the legacy dataProviderRNOG path (computes offsets before and after removal, since that path removes block offsets from the waveform).

monitoring_data_functions_sva/

• get_monitoring_data_uproot.py — everything related to reading monitoring.root/headers.root with uproot: low-level readers (get_event_info_from_monitoring_file, get_run_summary_from_monitoring_file, get_info_from_header_file), trigger-type assignment and consistency checks (assign_trigger_types, check_event_numbers_according_to_trigger_types), SNR calculation, and the main read_multiple_runs() entry point that loops over runs, validates and concatenates everything, and reports failed runs instead of raising.

helper_functions/

Small, reusable utilities that don't belong to a specific test:

• config_helper.py is imported from config_files_sva/ (see above) but conceptually lives here too.
• output_writer.py — every "write results to disk" function used by the main script: failed-run CSV, spectral results text file, SNR/RMS outlier-detail text files, RMS modality text files, glitching results text file, block-offset results text file, and create_result_csv_file(), which assembles the final per-channel summary CSV and computes the combined channel_health() verdict.
• read_rnog_runtable.py — read_rnog_runtable(), a thin wrapper around rnog_data.runtable used to turn a --time_range into a list of run numbers.

plotting_functions_sva/

One module per plot family, all consuming already-computed arrays/dicts (no analysis logic lives here):

• plotting_sva_spectrum.py — time-integrated surface/deep spectrum plots (normalized and unnormalized).
• plotting_sva_snr.py — SNR-vs-time plots per channel, including the outlier flags and z-score/k-value bands; also choose_day_interval(), a small helper that picks a sensible tick spacing based on the time range.
• plotting_sva_vrms.py — RMS/Vrms-vs-time plots (all triggers together and per trigger), the z-score single-trigger plot, rolling-mean plots, and the relative-median-shift heatmap (create_heatmap_plot).
• plotting_sva_glitch.py — glitching violin plots and the 99th-percentile-glitching-over-time plot.
• plotting_sva_debug.py — extra diagnostic plots (amplitude-ratio distributions, raw/​z-scored SNR distributions, RMS distributions with the KDE overlay) enabled with --debug_plot.
• plotting_sva_trigger_rate.py — trigger-rate-over-time and trigger-rate-heatmap plots.

expected_values/

Scripts to (re-)generate the reference ("expected") values that the SNR and RMS stability tests compare against, plus the values themselves:

• expected_snr_values.py / expected_rms_values.py — standalone CLIs (same run-selection arguments as the main script, plus --save-values) that read a known-stable period for a station and compute per-channel k-value/mean/std, writing them to expected_snr/expected_snr_values_station<ID>.json / expected_rms/expected_rms_station<ID>.json. Each JSON file carries a metadata block (station ID, run range, excluded runs, event count, trigger type, start/end time, and any comment — including an automatic note when a k-value hit the min/max cap) alongside the values block, so a reference file is self-documenting about how and when it was produced.
• expected_snr/, expected_rms/ — one combined JSON file per station (not split by parameter as in earlier versions).
• expected_block_offsets/ — one JSON per station with the reference block-offset statistics (median/IQR/P95/P99, in both ADC counts and mV) and the simulation settings used to derive them.
• outdated/ — reference-value scripts and results for the legacy dataProviderRNOG method (expected_snr_values_dataproviderrnog.py, expected_rms_values_dataproviderrnog.py). Per its own README: "These scripts won't be updated anymore."

sva_dataproviderrnog/ (legacy)

• read_rnog_data_nuradio.py — reads full waveforms via NuRadioReco's readRNOGDataMattak from combined.root files.
• science_verification_analysis_dataprovider.py — the original, pre-monitoring.root version of the full analysis, kept only for stations/periods without monitoring.root files. Its own module docstring spells out the caveats: no reference-value files exist for it by default (you'd need to run the expected_values/outdated scripts first and update the file paths), and its RMS analysis is done in ADC units, which is explicitly flagged as incorrect. Prefer science_verification_analysis_main.py whenever monitoring.root is available.

Output directories (created at runtime, not checked in)

plots/, detailed_results/, channel_health_summary/, logs/ are created next to whichever script is run (science_verification_analysis_main.py, or the expected_values scripts, which get their own plots_reference/, logs_reference/, results_reference/). None of these need to exist beforehand — every script creates them with os.makedirs(..., exist_ok=True).