From 9e070e3db449fb712272d882e52a628980c9c57c Mon Sep 17 00:00:00 2001
From: Gernot Maier <gernot.maier@desy.de>
Date: Mon, 11 May 2026 16:33:12 +0200
Subject: [PATCH 01/27] Add event-level production comparison CLI

---
 PLAN-LongIntegrationTests.md                  | 106 ++++++++
 PR.md                                         |  24 ++
 docs/source/user-guide/applications.md        |   1 +
 ...ols-compare-productions-on-event-level.rst |   6 +
 pyproject.toml                                |   1 +
 .../compare_productions_on_event_level.py     |  79 ++++++
 .../sim_events/production_comparison.py       | 204 +++++++++++++++
 .../plot_event_level_production_comparison.py | 246 ++++++++++++++++++
 ...compare_productions_on_event_level_run.yml |  17 ++
 ...test_compare_productions_on_event_level.py |  45 ++++
 .../sim_events/test_production_comparison.py  |  84 ++++++
 ..._plot_event_level_production_comparison.py |  45 ++++
 12 files changed, 858 insertions(+)
 create mode 100644 PLAN-LongIntegrationTests.md
 create mode 100644 PR.md
 create mode 100644 docs/source/user-guide/applications/simtools-compare-productions-on-event-level.rst
 create mode 100644 src/simtools/applications/compare_productions_on_event_level.py
 create mode 100644 src/simtools/sim_events/production_comparison.py
 create mode 100644 src/simtools/visualization/plot_event_level_production_comparison.py
 create mode 100644 tests/integration_tests/config/compare_productions_on_event_level_run.yml
 create mode 100644 tests/unit_tests/applications/test_compare_productions_on_event_level.py
 create mode 100644 tests/unit_tests/sim_events/test_production_comparison.py
 create mode 100644 tests/unit_tests/visualization/test_plot_event_level_production_comparison.py

diff --git a/PLAN-LongIntegrationTests.md b/PLAN-LongIntegrationTests.md
new file mode 100644
index 0000000000..dd99858e13
--- /dev/null
+++ b/PLAN-LongIntegrationTests.md
@@ -0,0 +1,106 @@
+# Simtools Long Integration Tests
+
+Long integration tests (LITs) allow to compare simulations results obtained with different simtools version (includes CORSIKA, sim-telarray, simtools), production models, or hadronic interaction models. Goal is to
+
+- detect unintended physics regressions
+- quantify differences between model versions (expected and unexpected differences)
+- provide release qualification metrics
+- provide long-term monitoring of execution statistics
+
+Two types of tools are required - this work is about the implementation of the first type.
+
+1. configure and run simulation productions for long integration tests
+2. tools to compare two/several simulation production results
+
+## Simulation productions for LITs
+
+Assumption (for now) is that these productions are medium sized and run on the local HTCondor cluster. The future extension to a Grid-style setup should be taken into account (separation of production definition from execution backend).
+
+Start with fixed energy simulations for easier comparison (and quicker testing). This can easily be extended later to energy ranges.
+
+Comparison grid:
+
+- simtools versions
+- production model versions
+- interaction models
+
+Production grid:
+
+- zenith angles (e.g., 20, 40, 60 deg)
+- azimuth angles (e.g., 0, 180 deg)
+- energy (e.g., [0.03, 0.1, 0.3, 1., 3., 10., 30., 100. TeV])
+- primary (e.g. [‘gamma’, ‘proton’])
+
+The grid might be adjustable (testing with 1 grid point, release with many).
+
+Additional input:
+
+- event statistics per energy bin (number of runs x events per run) *fixed during prototyping*
+- scatter area (E) *fixed during prototyping*
+- adjust energy grid as function of zenith *fixed during prototyping*
+
+Implementation:
+
+- tool similar to `simtools-simulate-prod-htcondor-generator` (or adaption of this tool)
+- write one condor file per simtools versions (as these are defined by Apptainer times)
+- modify the configuration paramaeter `apptainer_image` such that it can digest a dict with (label: apptainer image path) entries, e.g.:
+```
+apptainer_image:
+  7.0.0: <apptainer image path for simtools 7.0.0>
+  6.3.0: <apptainer image path for simtools 6.3.0>
+```
+- write a config file with one line per production grid point (similar to VTSSimpipe) to be ingested to HTCondor (x number of runs)
+
+This means that there will be at most len(simtools versions) job submission steps.
+
+## Technical details
+
+Suggestion is to modify `simtools-simulate-prod-htcondor-generator`.
+
+Copy and adapt the example in tests/integration_tests/config/simulate_prod_htcondor_generator_gamma_20_deg_north.yml:
+
+- value listed in the comparison and production grid should be allowed to be lists in this configuration
+- one condor output file per apptainer_image (this is the `container_image` field in the condor file)
+
+simtools-simulate-prod-htcondor-generator should be modified that it doesn't write a file with the values from the grid included (like the example in htcondor_submit/simulate_prod.submit.sh) but :
+
+- a parameter file with one line per grid point (job) with the value listed:
+ 
+ Example for params.txt:
+
+ ```text
+ 0.1 20 0 gamma epos 7.0.0 <grid output path>
+ 0.1 20 0 gamma epos 6.3.0 <grid output path>
+ ```
+
+ - note that this is the combination of the comparison and production grid values (energy, zenith, azimuth, primary, interaction model), so this file might be quite long (depending on the grid size); the simtools version is handled by the container image fields.
+
+Condor should submit then with `queue name,value,tag from params.txt` (instead of `queue <nruns>`).
+
+<grid output path> is the parameter given with `--pack_for_grid_register` and should contain the apptainer label (see modification of `apptainer_image` field in the config file) to be used for the job execution. This label should also be used in condor files, parameter file, and the job execution script.
+
+Remember to keep the following for job execution:
+
+```text
+# Process ID used to generate run number
+process_id="$1"
+# Load environment variables (for DB access)
+set -a; source "$2"
+```
+
+For the implementation:
+
+- try to minimize code changes to the existing tool
+- even when running with the current functionality (no grid), the should should use a params.txt file
+- Scalar config values are internally normalized to lists of length 1.
+- all runs should start with the same random seed and run number (this might change later; let's start simple)
+  
+
+Do not implement:
+
+- Grid backend execution
+- adaptive statistics
+- energy-dependent scatter area
+- comparison tooling
+
+Make a plan with detailed suggestions. Goal is to have a work plan in markdown format that can be used for implementation. This should include technical details and code snippets where appropriate.
diff --git a/PR.md b/PR.md
new file mode 100644
index 0000000000..cdfdc11409
--- /dev/null
+++ b/PR.md
@@ -0,0 +1,24 @@
+## Summary
+
+Add grid-style HTCondor submission support to `simtools-simulate-prod-htcondor-generator`.
+
+This change expands selected production and comparison parameters into a `params.txt` file, writes HTCondor submit files that queue from those params, and updates the submit script to consume per-job values from HTCondor arguments instead of hardcoded values.
+
+## What Changed
+
+- allow grid axes such as `primary`, `azimuth_angle`, `zenith_angle`, `model_version`, and interaction models to be provided as lists
+- support `apptainer_image` as either a single image path or a label-to-image mapping
+- generate one HTCondor params file row per expanded job and one condor file per apptainer label
+- pass job-specific values, including resolved `array_layout_name`, through the params file into the submit script
+- preserve energy units in the params file and forward them correctly to `simtools-simulate-prod`
+- add unit and integration coverage for parser behavior and HTCondor grid generation
+
+## Example Config
+
+See the grid example in:
+
+`tests/integration_tests/config/simulate_prod_htcondor_generator_grid_example.yml`
+
+## Validation
+
+- `pytest -q tests/unit_tests/configuration/test_configurator.py tests/unit_tests/applications/test_simulate_prod_htcondor_generator.py tests/unit_tests/job_execution/test_htcondor_script_generator.py`
diff --git a/docs/source/user-guide/applications.md b/docs/source/user-guide/applications.md
index 76c277c7de..d96cc885f0 100644
--- a/docs/source/user-guide/applications.md
+++ b/docs/source/user-guide/applications.md
@@ -52,6 +52,7 @@ Parameters with the same functionality are named consistently the same among all
 simtools-convert-all-model-parameters-from-simtel <applications/simtools-convert-all-model-parameters-from-simtel>
 simtools-convert-geo-coordinates-of-array-elements <applications/simtools-convert-geo-coordinates-of-array-elements>
 simtools-convert-model-parameter-from-simtel <applications/simtools-convert-model-parameter-from-simtel>
+simtools-compare-productions-on-event-level <applications/simtools-compare-productions-on-event-level>
 simtools-db-add-file-to-db <applications/simtools-db-add-file-to-db>
 simtools-db-add-simulation-model-from-repository-to-db <applications/simtools-db-add-simulation-model-from-repository-to-db>
 simtools-db-add-value-from-json-to-db <applications/simtools-db-add-value-from-json-to-db>
diff --git a/docs/source/user-guide/applications/simtools-compare-productions-on-event-level.rst b/docs/source/user-guide/applications/simtools-compare-productions-on-event-level.rst
new file mode 100644
index 0000000000..d23ca2d9b3
--- /dev/null
+++ b/docs/source/user-guide/applications/simtools-compare-productions-on-event-level.rst
@@ -0,0 +1,6 @@
+
+simtools-compare-productions-on-event-level
+===========================================
+
+.. automodule:: compare_productions_on_event_level
+   :members:
diff --git a/pyproject.toml b/pyproject.toml
index 6bd86626a6..fae01d68b8 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -73,6 +73,7 @@ optional-dependencies.tests = [
 urls."bug tracker" = "https://github.com/gammasim/simtools/issues"
 urls.documentation = "https://gammasim.github.io/simtools/"
 urls.repository = "https://github.com/gammasim/simtools"
+scripts.simtools-compare-productions-on-event-level = "simtools.applications.compare_productions_on_event_level:main"
 scripts.simtools-convert-all-model-parameters-from-simtel = "simtools.applications.convert_all_model_parameters_from_simtel:main"
 scripts.simtools-convert-geo-coordinates-of-array-elements = "simtools.applications.convert_geo_coordinates_of_array_elements:main"
 scripts.simtools-convert-model-parameter-from-simtel = "simtools.applications.convert_model_parameter_from_simtel:main"
diff --git a/src/simtools/applications/compare_productions_on_event_level.py b/src/simtools/applications/compare_productions_on_event_level.py
new file mode 100644
index 0000000000..84dc05161c
--- /dev/null
+++ b/src/simtools/applications/compare_productions_on_event_level.py
@@ -0,0 +1,79 @@
+#!/usr/bin/python3
+
+r"""Compare multiple simulation productions on event level.
+
+The application accepts repeated production descriptors and creates event-level
+comparison plots using reduced event data files.
+
+Command line arguments
+----------------------
+production (repeated, required)
+    Production descriptor in two fields:
+    1) label
+    2) comma-separated event data file patterns
+output_path (str, required)
+    Output directory for generated comparison plots.
+
+Examples
+--------
+Compare two productions:
+
+.. code-block:: console
+
+    simtools-compare-productions-on-event-level \
+        --production baseline "data/baseline/*.h5" \
+        --production candidate "data/candidate/*.h5,data/candidate_extra/*.h5" \
+        --output_path simtools-output/
+"""
+
+from simtools.application_control import build_application
+from simtools.sim_events.production_comparison import (
+    collect_production_metrics,
+    parse_production_arguments,
+)
+from simtools.visualization import plot_event_level_production_comparison
+
+
+def _add_arguments(parser):
+    """Register application-specific command line arguments."""
+    parser.initialize_application_arguments(["output_path"])
+    parser.add_argument(
+        "--production",
+        action="append",
+        nargs="+",
+        metavar=("LABEL", "EVENT_DATA_FILES"),
+        required=True,
+        help=(
+            "Production descriptor. "
+            "Use as: --production <label> <comma-separated file patterns>. "
+            "Repeat this argument at least twice."
+        ),
+    )
+    parser.add_argument(
+        "--telescope_ids",
+        nargs="+",
+        default=None,
+        help="Optional telescope IDs to filter triggers.",
+    )
+
+
+def main():
+    """See CLI description."""
+    app_context = build_application(
+        initialization_kwargs={"db_config": False, "output": True},
+    )
+
+    production_descriptors = parse_production_arguments(app_context.args["production"])
+    metrics_per_production = collect_production_metrics(
+        production_descriptors,
+        telescope_list=app_context.args["telescope_ids"],
+    )
+
+    output_directory = app_context.io_handler.get_output_directory()
+    plot_event_level_production_comparison.plot(
+        metrics_per_production, output_path=output_directory
+    )
+
+
+if __name__ == "__main__":
+    main()
diff --git a/src/simtools/sim_events/production_comparison.py b/src/simtools/sim_events/production_comparison.py
new file mode 100644
index 0000000000..d7387832da
--- /dev/null
+++ b/src/simtools/sim_events/production_comparison.py
@@ -0,0 +1,204 @@
+"""Utilities for event-level comparison across multiple simulation productions."""
+
+import logging
+from collections import Counter
+from dataclasses import dataclass
+
+import numpy as np
+
+from simtools.sim_events.reader import EventDataReader
+from simtools.utils.general import resolve_file_patterns
+
+
+@dataclass
+class ProductionDescriptor:
+    """Descriptor for one production input provided by the user."""
+
+    label: str
+    event_data_files: list[str]
+
+
+@dataclass
+class ProductionEventMetrics:
+    """Aggregated event-level metrics for one production."""
+
+    label: str
+    simulated_energies: np.ndarray
+    triggered_energies: np.ndarray
+    simulated_core_distances: np.ndarray
+    triggered_core_distances: np.ndarray
+    trigger_multiplicity: np.ndarray
+    trigger_combinations: Counter
+    telescope_participation: Counter
+    simulated_event_count: int
+    triggered_event_count: int
+
+    @property
+    def trigger_fraction(self):
+        """Return triggered/simulated fraction."""
+        if self.simulated_event_count <= 0:
+            return 0.0
+        return self.triggered_event_count / self.simulated_event_count
+
+
+def parse_production_arguments(production_arguments):
+    """Parse repeated production arguments into validated descriptors.
+
+    Parameters
+    ----------
+    production_arguments : list[list[str]]
+        Repeated ``--production`` arguments in the shape ``[label, patterns]``.
+
+    Returns
+    -------
+    list[ProductionDescriptor]
+        Validated and normalized production descriptors.
+
+    Raises
+    ------
+    ValueError
+        If configuration is malformed or does not contain at least two productions.
+    """
+    parsed_productions = _normalize_production_arguments(production_arguments)
+    if len(parsed_productions) < 2:
+        raise ValueError("At least two productions are required.")
+
+    labels = [item[0] for item in parsed_productions]
+    if len(set(labels)) != len(labels):
+        raise ValueError("Production labels must be unique.")
+
+    descriptors = []
+    for label, pattern_list in parsed_productions:
+        patterns = [pattern.strip() for pattern in pattern_list.split(",") if pattern.strip()]
+        if len(patterns) == 0:
+            raise ValueError(f"Production '{label}' has no event_data_file pattern.")
+
+        resolved_files = [str(path) for path in resolve_file_patterns(patterns)]
+        if len(resolved_files) == 0:
+            raise ValueError(f"Production '{label}' does not resolve to any files.")
+
+        descriptors.append(ProductionDescriptor(label=label, event_data_files=resolved_files))
+
+    return descriptors
+
+
+def _normalize_production_arguments(production_arguments):
+    """Normalize raw production arguments into ``[(label, files), ...]``."""
+    if not production_arguments:
+        return []
+
+    if all(isinstance(item, str) for item in production_arguments):
+        if len(production_arguments) % 2 != 0:
+            raise ValueError("Production arguments must be provided as label/file pairs.")
+        return [
+            (production_arguments[index], production_arguments[index + 1])
+            for index in range(0, len(production_arguments), 2)
+        ]
+
+    normalized = []
+    for item in production_arguments:
+        if isinstance(item, list | tuple):
+            if len(item) == 2 and all(isinstance(value, str) for value in item):
+                normalized.append((item[0], item[1]))
+                continue
+
+            if all(isinstance(value, str) for value in item):
+                normalized.extend(_normalize_production_arguments(list(item)))
+                continue
+
+        raise ValueError("Production arguments must be provided as label/file pairs.")
+
+    return normalized
+
+
+def collect_production_metrics(production_descriptors, telescope_list=None):
+    """Collect event-level comparison metrics for each production.
+
+    Parameters
+    ----------
+    production_descriptors : list[ProductionDescriptor]
+        Input descriptor for each production.
+    telescope_list : list[str], optional
+        Telescope IDs to filter triggered events.
+
+    Returns
+    -------
+    list[ProductionEventMetrics]
+        Aggregated metrics per production.
+    """
+    return [
+        _collect_single_production_metrics(descriptor, telescope_list=telescope_list)
+        for descriptor in production_descriptors
+    ]
+
+
+def _collect_single_production_metrics(production_descriptor, telescope_list=None):
+    """Collect event-level metrics for one production descriptor."""
+    logger = logging.getLogger(__name__)
+
+    simulated_energies = []
+    triggered_energies = []
+    simulated_core_distances = []
+    triggered_core_distances = []
+    trigger_multiplicity = []
+    trigger_combinations = Counter()
+    telescope_participation = Counter()
+    simulated_event_count = 0
+    triggered_event_count = 0
+
+    for event_data_file in production_descriptor.event_data_files:
+        reader = EventDataReader(event_data_file, telescope_list=telescope_list)
+        for data_set in reader.data_sets:
+            _, shower_data, triggered_shower_data, triggered_data = reader.read_event_data(
+                event_data_file,
+                table_name_map=data_set,
+            )
+
+            if shower_data is None:
+                continue
+
+            simulated_energies.append(np.asarray(shower_data.simulated_energy))
+            simulated_core_distances.append(np.asarray(shower_data.core_distance_shower))
+            simulated_event_count += len(shower_data.simulated_energy)
+
+            if triggered_data is None or triggered_shower_data is None:
+                continue
+
+            triggered_energies.append(np.asarray(triggered_shower_data.simulated_energy))
+            triggered_core_distances.append(np.asarray(triggered_shower_data.core_distance_shower))
+
+            multiplicity = np.array([len(tel_list) for tel_list in triggered_data.telescope_list])
+            trigger_multiplicity.append(multiplicity)
+            triggered_event_count += len(multiplicity)
+
+            for tel_list in triggered_data.telescope_list:
+                telescopes = tuple(sorted(str(telescope) for telescope in tel_list))
+                combination_key = ",".join(telescopes)
+                trigger_combinations[combination_key] += 1
+                for telescope in set(telescopes):
+                    telescope_participation[telescope] += 1
+
+    logger.info(
+        f"Collected production {production_descriptor.label}: "
+        f"simulated={simulated_event_count}, triggered={triggered_event_count}"
+    )
+
+    return ProductionEventMetrics(
+        label=production_descriptor.label,
+        simulated_energies=_safe_concat(simulated_energies),
+        triggered_energies=_safe_concat(triggered_energies),
+        simulated_core_distances=_safe_concat(simulated_core_distances),
+        triggered_core_distances=_safe_concat(triggered_core_distances),
+        trigger_multiplicity=_safe_concat(trigger_multiplicity, dtype=int),
+        trigger_combinations=trigger_combinations,
+        telescope_participation=telescope_participation,
+        simulated_event_count=simulated_event_count,
+        triggered_event_count=triggered_event_count,
+    )
+
+
+def _safe_concat(arrays, dtype=float):
+    """Concatenate list of arrays, returning an empty array for empty input."""
+    if len(arrays) == 0:
+        return np.array([], dtype=dtype)
+    return np.concatenate(arrays)
diff --git a/src/simtools/visualization/plot_event_level_production_comparison.py b/src/simtools/visualization/plot_event_level_production_comparison.py
new file mode 100644
index 0000000000..ec04708c7b
--- /dev/null
+++ b/src/simtools/visualization/plot_event_level_production_comparison.py
@@ -0,0 +1,246 @@
+"""Plot event-level comparisons across multiple simulation productions."""
+
+import logging
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+_logger = logging.getLogger(__name__)
+
+
+def plot(metrics_per_production, output_path, bins=40):
+    """Create all event-level production comparison plots.
+
+    Parameters
+    ----------
+    metrics_per_production : list[ProductionEventMetrics]
+        Aggregated metrics per production.
+    output_path : pathlib.Path
+        Output directory for generated figures.
+    bins : int, optional
+        Number of bins for 1D histograms.
+    """
+    _plot_trigger_fraction(metrics_per_production, output_path)
+    _plot_trigger_multiplicity(metrics_per_production, output_path)
+    _plot_trigger_combinations(metrics_per_production, output_path)
+    _plot_triggered_vs_quantity(
+        metrics_per_production,
+        output_path,
+        quantity_name="energy",
+        x_label="Primary Energy (TeV)",
+        x_scale="log",
+        bins=bins,
+    )
+    _plot_triggered_vs_quantity(
+        metrics_per_production,
+        output_path,
+        quantity_name="core_distance",
+        x_label="Core Distance (m)",
+        x_scale="linear",
+        bins=bins,
+    )
+    _plot_telescope_participation(metrics_per_production, output_path)
+
+
+def _save_figure(fig, output_path, filename):
+    """Save figure and close it."""
+    output_file = output_path / filename
+    _logger.info(f"Saving comparison plot: {output_file}")
+    fig.savefig(output_file, dpi=300, bbox_inches="tight")
+    plt.close(fig)
+
+
+def _plot_trigger_fraction(metrics_per_production, output_path):
+    """Plot triggered/simulated event fractions by production."""
+    labels = [metrics.label for metrics in metrics_per_production]
+    values = [metrics.trigger_fraction for metrics in metrics_per_production]
+
+    fig, ax = plt.subplots(figsize=(9, 6))
+    bars = ax.bar(labels, values, color="tab:blue", alpha=0.85)
+    ax.set_ylabel("Trigger Fraction")
+    ax.set_title("Trigger Fraction by Production")
+    ax.set_ylim(0.0, max(1.0, max(values, default=0.0) * 1.15))
+    ax.grid(axis="y", alpha=0.25)
+    for bar, value in zip(bars, values):
+        ax.text(
+            bar.get_x() + bar.get_width() / 2,
+            bar.get_height(),
+            f"{value:.3f}",
+            ha="center",
+            va="bottom",
+            fontsize=9,
+        )
+
+    _save_figure(fig, output_path, "trigger_fraction_comparison.png")
+
+
+def _plot_trigger_multiplicity(metrics_per_production, output_path):
+    """Plot triggered telescope multiplicity distributions."""
+    fig, ax = plt.subplots(figsize=(9, 6))
+
+    global_max = 0
+    for metrics in metrics_per_production:
+        if metrics.trigger_multiplicity.size == 0:
+            continue
+        global_max = max(global_max, int(np.max(metrics.trigger_multiplicity)))
+
+    if global_max == 0:
+        _logger.warning("Skipping trigger multiplicity plot, no triggered events available.")
+        plt.close(fig)
+        return
+
+    bins = np.arange(1, global_max + 2)
+    for metrics in metrics_per_production:
+        if metrics.trigger_multiplicity.size == 0:
+            continue
+        counts, _ = np.histogram(metrics.trigger_multiplicity, bins=bins)
+        fractions = counts / counts.sum() if counts.sum() > 0 else counts
+        ax.step(bins[:-1], fractions, where="post", linewidth=2, label=metrics.label)
+
+    ax.set_xlabel("Triggered Telescopes per Event")
+    ax.set_ylabel("Fraction of Triggered Events")
+    ax.set_title("Trigger Multiplicity Distribution")
+    ax.set_xticks(bins[:-1])
+    ax.grid(alpha=0.25)
+    ax.legend()
+
+    _save_figure(fig, output_path, "trigger_multiplicity_comparison.png")
+
+
+def _plot_trigger_combinations(metrics_per_production, output_path, top_n=12):
+    """Plot trigger combination distributions for each production."""
+    combinations = set()
+    for metrics in metrics_per_production:
+        combinations.update(metrics.trigger_combinations.keys())
+
+    if len(combinations) == 0:
+        _logger.warning("Skipping trigger combination plot, no combinations available.")
+        return
+
+    totals = {}
+    for combination in combinations:
+        totals[combination] = sum(
+            metrics.trigger_combinations.get(combination, 0) for metrics in metrics_per_production
+        )
+
+    selected = [name for name, _ in sorted(totals.items(), key=lambda item: item[1], reverse=True)]
+    selected = selected[:top_n]
+
+    x_values = np.arange(len(selected))
+    width = 0.8 / max(1, len(metrics_per_production))
+    fig, ax = plt.subplots(figsize=(max(10, len(selected) * 1.1), 6))
+
+    for index, metrics in enumerate(metrics_per_production):
+        event_norm = metrics.triggered_event_count if metrics.triggered_event_count > 0 else 1
+        y_values = [metrics.trigger_combinations.get(name, 0) / event_norm for name in selected]
+        offset = (index - (len(metrics_per_production) - 1) / 2.0) * width
+        ax.bar(x_values + offset, y_values, width=width, label=metrics.label)
+
+    ax.set_xticks(x_values)
+    ax.set_xticklabels(selected, rotation=45, ha="right")
+    ax.set_ylabel("Fraction of Triggered Events")
+    ax.set_title("Top Trigger Combinations")
+    ax.grid(axis="y", alpha=0.25)
+    ax.legend()
+
+    _save_figure(fig, output_path, "trigger_combination_comparison.png")
+
+
+def _plot_triggered_vs_quantity(
+    metrics_per_production,
+    output_path,
+    quantity_name,
+    x_label,
+    x_scale,
+    bins,
+):
+    """Plot simulated vs triggered distributions for one quantity."""
+    fig, ax = plt.subplots(figsize=(9, 6))
+
+    for metrics in metrics_per_production:
+        if quantity_name == "energy":
+            simulated = metrics.simulated_energies
+            triggered = metrics.triggered_energies
+        else:
+            simulated = metrics.simulated_core_distances
+            triggered = metrics.triggered_core_distances
+
+        if simulated.size == 0:
+            continue
+
+        bin_edges = _get_bin_edges(simulated, x_scale=x_scale, bins=bins)
+        sim_counts, _ = np.histogram(simulated, bins=bin_edges)
+        trig_counts, _ = np.histogram(triggered, bins=bin_edges)
+        with np.errstate(divide="ignore", invalid="ignore"):
+            efficiency = np.divide(
+                trig_counts,
+                sim_counts,
+                out=np.zeros_like(trig_counts, dtype=float),
+                where=sim_counts > 0,
+            )
+
+        x_values = 0.5 * (bin_edges[:-1] + bin_edges[1:])
+        ax.plot(x_values, efficiency, linewidth=2, label=metrics.label)
+
+    ax.set_xlabel(x_label)
+    ax.set_ylabel("Triggered / Simulated")
+    ax.set_title(f"Triggered Event Fraction vs {x_label}")
+    ax.set_xscale(x_scale)
+    ax.set_ylim(0.0, 1.05)
+    ax.grid(alpha=0.25)
+    ax.legend()
+
+    filename = f"triggered_fraction_vs_{quantity_name}.png"
+    _save_figure(fig, output_path, filename)
+
+
+def _get_bin_edges(values, x_scale, bins):
+    """Return robust bin edges for the given values."""
+    if len(values) == 0:
+        return np.linspace(0.0, 1.0, bins + 1)
+
+    min_value = float(np.min(values))
+    max_value = float(np.max(values))
+    if np.isclose(min_value, max_value):
+        max_value = min_value + 1.0
+
+    if x_scale == "log":
+        min_value = max(min_value, np.finfo(float).tiny)
+        return np.logspace(np.log10(min_value), np.log10(max_value), bins + 1)
+    return np.linspace(min_value, max_value, bins + 1)
+
+
+def _plot_telescope_participation(metrics_per_production, output_path):
+    """Plot telescope participation fractions by production."""
+    telescopes = sorted(
+        {
+            telescope
+            for metrics in metrics_per_production
+            for telescope in metrics.telescope_participation.keys()
+        }
+    )
+    if len(telescopes) == 0:
+        _logger.warning("Skipping telescope participation plot, no triggered events available.")
+        return
+
+    x_values = np.arange(len(telescopes))
+    width = 0.8 / max(1, len(metrics_per_production))
+    fig, ax = plt.subplots(figsize=(max(10, len(telescopes) * 0.4), 6))
+
+    for index, metrics in enumerate(metrics_per_production):
+        event_norm = metrics.triggered_event_count if metrics.triggered_event_count > 0 else 1
+        fractions = [
+            metrics.telescope_participation.get(telescope, 0) / event_norm
+            for telescope in telescopes
+        ]
+        offset = (index - (len(metrics_per_production) - 1) / 2.0) * width
+        ax.bar(x_values + offset, fractions, width=width, label=metrics.label)
+
+    ax.set_xticks(x_values)
+    ax.set_xticklabels(telescopes, rotation=90)
+    ax.set_ylabel("Participation Fraction")
+    ax.set_title("Telescope Participation in Triggered Events")
+    ax.grid(axis="y", alpha=0.25)
+    ax.legend()
+
+    _save_figure(fig, output_path, "telescope_participation_fraction.png")
diff --git a/tests/integration_tests/config/compare_productions_on_event_level_run.yml b/tests/integration_tests/config/compare_productions_on_event_level_run.yml
new file mode 100644
index 0000000000..9aa1c408f7
--- /dev/null
+++ b/tests/integration_tests/config/compare_productions_on_event_level_run.yml
@@ -0,0 +1,17 @@
+---
+applications:
+- application: simtools-compare-productions-on-event-level
+  configuration:
+    production:
+    - baseline
+    - tests/resources/proton_za20deg_azm000deg_North_alpha_6.0.0_reduced_event_data.hdf5
+    - candidate
+    - tests/resources/proton_za20deg_azm000deg_North_alpha_6.0.0_reduced_event_data.hdf5
+    output_path: simtools-output
+  integration_tests:
+  - output_file: trigger_fraction_comparison.png
+  - output_file: trigger_multiplicity_comparison.png
+  - output_file: telescope_participation_fraction.png
+  test_name: run
+schema_name: application_workflow.metaschema
+schema_version: 0.4.0
diff --git a/tests/unit_tests/applications/test_compare_productions_on_event_level.py b/tests/unit_tests/applications/test_compare_productions_on_event_level.py
new file mode 100644
index 0000000000..9e6d3725e4
--- /dev/null
+++ b/tests/unit_tests/applications/test_compare_productions_on_event_level.py
@@ -0,0 +1,45 @@
+from pathlib import Path
+from types import SimpleNamespace
+from unittest.mock import MagicMock, patch
+
+import simtools.applications.compare_productions_on_event_level as app
+
+
+def test_main_collects_metrics_and_plots(tmp_test_directory):
+    """Test application orchestration from CLI args to plotting output."""
+    output_dir = Path(tmp_test_directory) / "plots"
+    app_context = SimpleNamespace(
+        args={
+            "production": [["baseline", "base_*.h5"], ["candidate", "cand_*.h5"]],
+            "telescope_ids": ["LSTN-01"],
+        },
+        io_handler=MagicMock(),
+    )
+    app_context.io_handler.get_output_directory.return_value = output_dir
+
+    parsed_productions = [MagicMock(), MagicMock()]
+    collected_metrics = [MagicMock(), MagicMock()]
+
+    with (
+        patch(
+            "simtools.applications.compare_productions_on_event_level.build_application",
+            return_value=app_context,
+        ),
+        patch(
+            "simtools.applications.compare_productions_on_event_level.parse_production_arguments",
+            return_value=parsed_productions,
+        ) as mock_parse,
+        patch(
+            "simtools.applications.compare_productions_on_event_level.collect_production_metrics",
+            return_value=collected_metrics,
+        ) as mock_collect,
+        patch(
+            "simtools.applications.compare_productions_on_event_level."
+            "plot_event_level_production_comparison.plot"
+        ) as mock_plot,
+    ):
+        app.main()
+
+    mock_parse.assert_called_once_with(app_context.args["production"])
+    mock_collect.assert_called_once_with(parsed_productions, telescope_list=["LSTN-01"])
+    mock_plot.assert_called_once_with(collected_metrics, output_path=output_dir)
diff --git a/tests/unit_tests/sim_events/test_production_comparison.py b/tests/unit_tests/sim_events/test_production_comparison.py
new file mode 100644
index 0000000000..12a55be408
--- /dev/null
+++ b/tests/unit_tests/sim_events/test_production_comparison.py
@@ -0,0 +1,84 @@
+from types import SimpleNamespace
+
+import numpy as np
+import pytest
+
+from simtools.sim_events.production_comparison import (
+    collect_production_metrics,
+    parse_production_arguments,
+)
+
+
+def test_parse_production_arguments_requires_two_productions():
+    """Test parser rejects less than two productions."""
+    with pytest.raises(ValueError, match="At least two productions"):
+        parse_production_arguments([["baseline", "base.h5"]])
+
+
+def test_parse_production_arguments_resolves_flattened_pairs(mocker):
+    """Test parser supports flattened label/file list from configuration files."""
+    mocker.patch(
+        "simtools.sim_events.production_comparison.resolve_file_patterns",
+        side_effect=lambda patterns: patterns,
+    )
+
+    descriptors = parse_production_arguments(["baseline", "base_*.h5", "candidate", "cand_*.h5"])
+
+    assert [descriptor.label for descriptor in descriptors] == ["baseline", "candidate"]
+    assert descriptors[0].event_data_files == ["base_*.h5"]
+    assert descriptors[1].event_data_files == ["cand_*.h5"]
+
+
+def test_parse_production_arguments_rejects_duplicate_labels(mocker):
+    """Test parser rejects duplicated production labels."""
+    mocker.patch(
+        "simtools.sim_events.production_comparison.resolve_file_patterns",
+        side_effect=lambda patterns: patterns,
+    )
+
+    with pytest.raises(ValueError, match="labels must be unique"):
+        parse_production_arguments([["same", "a.h5"], ["same", "b.h5"]])
+
+
+def test_collect_production_metrics_aggregates_event_data(mocker):
+    """Test production metrics aggregation across all event-level quantities."""
+    descriptor = SimpleNamespace(label="baseline", event_data_files=["baseline_file.h5"])
+
+    shower_data = SimpleNamespace(
+        simulated_energy=np.array([1.0, 2.0, 3.0]),
+        core_distance_shower=np.array([10.0, 20.0, 30.0]),
+    )
+    triggered_shower_data = SimpleNamespace(
+        simulated_energy=np.array([1.0, 3.0]),
+        core_distance_shower=np.array([10.0, 30.0]),
+    )
+    triggered_data = SimpleNamespace(
+        telescope_list=[
+            np.array(["LSTN-01", "MSTN-01"]),
+            np.array(["MSTN-01"]),
+        ]
+    )
+
+    mock_reader = mocker.patch("simtools.sim_events.production_comparison.EventDataReader")
+    mock_reader.return_value.data_sets = [{"SHOWERS": "SHOWERS", "TRIGGERS": "TRIGGERS"}]
+    mock_reader.return_value.read_event_data.return_value = (
+        SimpleNamespace(),
+        shower_data,
+        triggered_shower_data,
+        triggered_data,
+    )
+
+    metrics = collect_production_metrics([descriptor])
+
+    assert len(metrics) == 1
+    assert metrics[0].label == "baseline"
+    assert metrics[0].simulated_event_count == 3
+    assert metrics[0].triggered_event_count == 2
+    np.testing.assert_array_equal(metrics[0].simulated_energies, np.array([1.0, 2.0, 3.0]))
+    np.testing.assert_array_equal(metrics[0].triggered_energies, np.array([1.0, 3.0]))
+    np.testing.assert_array_equal(metrics[0].trigger_multiplicity, np.array([2, 1]))
+    assert metrics[0].trigger_combinations["LSTN-01,MSTN-01"] == 1
+    assert metrics[0].trigger_combinations["MSTN-01"] == 1
+    assert metrics[0].telescope_participation["MSTN-01"] == 2
+    assert metrics[0].telescope_participation["LSTN-01"] == 1
+    assert metrics[0].trigger_fraction == pytest.approx(2.0 / 3.0)
diff --git a/tests/unit_tests/visualization/test_plot_event_level_production_comparison.py b/tests/unit_tests/visualization/test_plot_event_level_production_comparison.py
new file mode 100644
index 0000000000..edaba6452c
--- /dev/null
+++ b/tests/unit_tests/visualization/test_plot_event_level_production_comparison.py
@@ -0,0 +1,45 @@
+from collections import Counter
+from pathlib import Path
+
+import numpy as np
+
+from simtools.sim_events.production_comparison import ProductionEventMetrics
+from simtools.visualization import plot_event_level_production_comparison
+
+
+def _build_metrics(label, simulated_scale, triggered_scale):
+    """Build a compact metrics fixture for plotting tests."""
+    return ProductionEventMetrics(
+        label=label,
+        simulated_energies=np.array([0.1, 0.3, 1.0, 3.0]) * simulated_scale,
+        triggered_energies=np.array([0.3, 1.0]) * triggered_scale,
+        simulated_core_distances=np.array([10.0, 20.0, 30.0, 40.0]),
+        triggered_core_distances=np.array([20.0, 40.0]),
+        trigger_multiplicity=np.array([2, 1]),
+        trigger_combinations=Counter({"LSTN-01,MSTN-01": 1, "MSTN-01": 1}),
+        telescope_participation=Counter({"LSTN-01": 1, "MSTN-01": 2}),
+        simulated_event_count=4,
+        triggered_event_count=2,
+    )
+
+
+def test_plot_writes_event_level_comparison_figures(tmp_test_directory):
+    """Test full comparison plotting writes all expected output files."""
+    output_path = Path(tmp_test_directory)
+    metrics = [
+        _build_metrics("baseline", simulated_scale=1.0, triggered_scale=1.0),
+        _build_metrics("candidate", simulated_scale=1.2, triggered_scale=1.1),
+    ]
+
+    plot_event_level_production_comparison.plot(metrics, output_path=output_path, bins=8)
+
+    expected_files = [
+        "trigger_fraction_comparison.png",
+        "trigger_multiplicity_comparison.png",
+        "trigger_combination_comparison.png",
+        "triggered_fraction_vs_energy.png",
+        "triggered_fraction_vs_core_distance.png",
+        "telescope_participation_fraction.png",
+    ]
+    for file_name in expected_files:
+        assert (output_path / file_name).exists()

From 6caa7c86d04fcbef07545eb87fd0a95b967fc907 Mon Sep 17 00:00:00 2001
From: Gernot Maier <gernot.maier@desy.de>
Date: Mon, 11 May 2026 16:34:15 +0200
Subject: [PATCH 02/27] Remove unrelated planning files from PR

---
 PLAN-LongIntegrationTests.md | 106 -----------------------------------
 PR.md                        |  24 --------
 2 files changed, 130 deletions(-)
 delete mode 100644 PLAN-LongIntegrationTests.md
 delete mode 100644 PR.md

diff --git a/PLAN-LongIntegrationTests.md b/PLAN-LongIntegrationTests.md
deleted file mode 100644
index dd99858e13..0000000000
--- a/PLAN-LongIntegrationTests.md
+++ /dev/null
@@ -1,106 +0,0 @@
-# Simtools Long Integration Tests
-
-Long integration tests (LITs) allow to compare simulations results obtained with different simtools version (includes CORSIKA, sim-telarray, simtools), production models, or hadronic interaction models. Goal is to
-
-- detect unintended physics regressions
-- quantify differences between model versions (expected and unexpected differences)
-- provide release qualification metrics
-- provide long-term monitoring of execution statistics
-
-Two types of tools are required - this work is about the implementation of the first type.
-
-1. configure and run simulation productions for long integration tests
-2. tools to compare two/several simulation production results
-
-## Simulation productions for LITs
-
-Assumption (for now) is that these productions are medium sized and run on the local HTCondor cluster. The future extension to a Grid-style setup should be taken into account (separation of production definition from execution backend).
-
-Start with fixed energy simulations for easier comparison (and quicker testing). This can easily be extended later to energy ranges.
-
-Comparison grid:
-
-- simtools versions
-- production model versions
-- interaction models
-
-Production grid:
-
-- zenith angles (e.g., 20, 40, 60 deg)
-- azimuth angles (e.g., 0, 180 deg)
-- energy (e.g., [0.03, 0.1, 0.3, 1., 3., 10., 30., 100. TeV])
-- primary (e.g. [‘gamma’, ‘proton’])
-
-The grid might be adjustable (testing with 1 grid point, release with many).
-
-Additional input:
-
-- event statistics per energy bin (number of runs x events per run) *fixed during prototyping*
-- scatter area (E) *fixed during prototyping*
-- adjust energy grid as function of zenith *fixed during prototyping*
-
-Implementation:
-
-- tool similar to `simtools-simulate-prod-htcondor-generator` (or adaption of this tool)
-- write one condor file per simtools versions (as these are defined by Apptainer times)
-- modify the configuration paramaeter `apptainer_image` such that it can digest a dict with (label: apptainer image path) entries, e.g.:
-```
-apptainer_image:
-  7.0.0: <apptainer image path for simtools 7.0.0>
-  6.3.0: <apptainer image path for simtools 6.3.0>
-```
-- write a config file with one line per production grid point (similar to VTSSimpipe) to be ingested to HTCondor (x number of runs)
-
-This means that there will be at most len(simtools versions) job submission steps.
-
-## Technical details
-
-Suggestion is to modify `simtools-simulate-prod-htcondor-generator`.
-
-Copy and adapt the example in tests/integration_tests/config/simulate_prod_htcondor_generator_gamma_20_deg_north.yml:
-
-- value listed in the comparison and production grid should be allowed to be lists in this configuration
-- one condor output file per apptainer_image (this is the `container_image` field in the condor file)
-
-simtools-simulate-prod-htcondor-generator should be modified that it doesn't write a file with the values from the grid included (like the example in htcondor_submit/simulate_prod.submit.sh) but :
-
-- a parameter file with one line per grid point (job) with the value listed:
- 
- Example for params.txt:
-
- ```text
- 0.1 20 0 gamma epos 7.0.0 <grid output path>
- 0.1 20 0 gamma epos 6.3.0 <grid output path>
- ```
-
- - note that this is the combination of the comparison and production grid values (energy, zenith, azimuth, primary, interaction model), so this file might be quite long (depending on the grid size); the simtools version is handled by the container image fields.
-
-Condor should submit then with `queue name,value,tag from params.txt` (instead of `queue <nruns>`).
-
-<grid output path> is the parameter given with `--pack_for_grid_register` and should contain the apptainer label (see modification of `apptainer_image` field in the config file) to be used for the job execution. This label should also be used in condor files, parameter file, and the job execution script.
-
-Remember to keep the following for job execution:
-
-```text
-# Process ID used to generate run number
-process_id="$1"
-# Load environment variables (for DB access)
-set -a; source "$2"
-```
-
-For the implementation:
-
-- try to minimize code changes to the existing tool
-- even when running with the current functionality (no grid), the should should use a params.txt file
-- Scalar config values are internally normalized to lists of length 1.
-- all runs should start with the same random seed and run number (this might change later; let's start simple)
-  
-
-Do not implement:
-
-- Grid backend execution
-- adaptive statistics
-- energy-dependent scatter area
-- comparison tooling
-
-Make a plan with detailed suggestions. Goal is to have a work plan in markdown format that can be used for implementation. This should include technical details and code snippets where appropriate.
diff --git a/PR.md b/PR.md
deleted file mode 100644
index cdfdc11409..0000000000
--- a/PR.md
+++ /dev/null
@@ -1,24 +0,0 @@
-## Summary
-
-Add grid-style HTCondor submission support to `simtools-simulate-prod-htcondor-generator`.
-
-This change expands selected production and comparison parameters into a `params.txt` file, writes HTCondor submit files that queue from those params, and updates the submit script to consume per-job values from HTCondor arguments instead of hardcoded values.
-
-## What Changed
-
-- allow grid axes such as `primary`, `azimuth_angle`, `zenith_angle`, `model_version`, and interaction models to be provided as lists
-- support `apptainer_image` as either a single image path or a label-to-image mapping
-- generate one HTCondor params file row per expanded job and one condor file per apptainer label
-- pass job-specific values, including resolved `array_layout_name`, through the params file into the submit script
-- preserve energy units in the params file and forward them correctly to `simtools-simulate-prod`
-- add unit and integration coverage for parser behavior and HTCondor grid generation
-
-## Example Config
-
-See the grid example in:
-
-`tests/integration_tests/config/simulate_prod_htcondor_generator_grid_example.yml`
-
-## Validation
-
-- `pytest -q tests/unit_tests/configuration/test_configurator.py tests/unit_tests/applications/test_simulate_prod_htcondor_generator.py tests/unit_tests/job_execution/test_htcondor_script_generator.py`

From 500d4a55ece5021147252e09e26cf97a937ba114 Mon Sep 17 00:00:00 2001
From: Gernot Maier <gernot.maier@desy.de>
Date: Mon, 11 May 2026 16:35:38 +0200
Subject: [PATCH 03/27] Add changelog for PR 2181

---
 docs/changes/2181.feature.md | 1 +
 1 file changed, 1 insertion(+)
 create mode 100644 docs/changes/2181.feature.md

diff --git a/docs/changes/2181.feature.md b/docs/changes/2181.feature.md
new file mode 100644
index 0000000000..ca5cae6df3
--- /dev/null
+++ b/docs/changes/2181.feature.md
@@ -0,0 +1 @@
+Add an event-level production comparison application with multi-production trigger-distribution plots.

From 8990ae1ca970ae55f026c6c0176723983fa2f68b Mon Sep 17 00:00:00 2001
From: Gernot Maier <gernot.maier@desy.de>
Date: Wed, 13 May 2026 10:27:22 +0200
Subject: [PATCH 04/27] improved plots

---
 .../compare_productions_on_event_level.py     |   2 +-
 .../sim_events/production_comparison.py       | 266 +++++++--
 .../plot_event_level_production_comparison.py | 515 ++++++++++++++++--
 .../sim_events/test_production_comparison.py  |  35 +-
 ..._plot_event_level_production_comparison.py |  62 ++-
 5 files changed, 766 insertions(+), 114 deletions(-)

diff --git a/src/simtools/applications/compare_productions_on_event_level.py b/src/simtools/applications/compare_productions_on_event_level.py
index 84dc05161c..933413e074 100644
--- a/src/simtools/applications/compare_productions_on_event_level.py
+++ b/src/simtools/applications/compare_productions_on_event_level.py
@@ -46,7 +46,7 @@ def _add_arguments(parser):
         help=(
             "Production descriptor. "
             "Use as: --production <label> <comma-separated file patterns>. "
-            "Repeat this argument at least twice."
+            "Repeat this argument as needed for multiple productions."
         ),
     )
     parser.add_argument(
diff --git a/src/simtools/sim_events/production_comparison.py b/src/simtools/sim_events/production_comparison.py
index d7387832da..167c52abe1 100644
--- a/src/simtools/sim_events/production_comparison.py
+++ b/src/simtools/sim_events/production_comparison.py
@@ -1,12 +1,13 @@
 """Utilities for event-level comparison across multiple simulation productions."""
 
 import logging
-from collections import Counter
-from dataclasses import dataclass
+from collections import Counter, defaultdict
+from dataclasses import dataclass, field
 
 import numpy as np
 
 from simtools.sim_events.reader import EventDataReader
+from simtools.utils import names
 from simtools.utils.general import resolve_file_patterns
 
 
@@ -32,6 +33,9 @@ class ProductionEventMetrics:
     telescope_participation: Counter
     simulated_event_count: int
     triggered_event_count: int
+    simulated_angular_distances: np.ndarray = field(default_factory=lambda: np.array([]))
+    triggered_angular_distances: np.ndarray = field(default_factory=lambda: np.array([]))
+    per_type: dict = field(default_factory=dict)
 
     @property
     def trigger_fraction(self):
@@ -57,11 +61,11 @@ def parse_production_arguments(production_arguments):
     Raises
     ------
     ValueError
-        If configuration is malformed or does not contain at least two productions.
+        If configuration is malformed or does not contain any production.
     """
     parsed_productions = _normalize_production_arguments(production_arguments)
-    if len(parsed_productions) < 2:
-        raise ValueError("At least two productions are required.")
+    if len(parsed_productions) == 0:
+        raise ValueError("At least one production is required.")
 
     labels = [item[0] for item in parsed_productions]
     if len(set(labels)) != len(labels):
@@ -136,67 +140,219 @@ def _collect_single_production_metrics(production_descriptor, telescope_list=Non
     """Collect event-level metrics for one production descriptor."""
     logger = logging.getLogger(__name__)
 
-    simulated_energies = []
-    triggered_energies = []
-    simulated_core_distances = []
-    triggered_core_distances = []
-    trigger_multiplicity = []
-    trigger_combinations = Counter()
-    telescope_participation = Counter()
-    simulated_event_count = 0
-    triggered_event_count = 0
-
+    metric_accumulators = _initialize_metric_accumulators()
     for event_data_file in production_descriptor.event_data_files:
-        reader = EventDataReader(event_data_file, telescope_list=telescope_list)
-        for data_set in reader.data_sets:
-            _, shower_data, triggered_shower_data, triggered_data = reader.read_event_data(
-                event_data_file,
-                table_name_map=data_set,
-            )
-
-            if shower_data is None:
-                continue
-
-            simulated_energies.append(np.asarray(shower_data.simulated_energy))
-            simulated_core_distances.append(np.asarray(shower_data.core_distance_shower))
-            simulated_event_count += len(shower_data.simulated_energy)
-
-            if triggered_data is None or triggered_shower_data is None:
-                continue
-
-            triggered_energies.append(np.asarray(triggered_shower_data.simulated_energy))
-            triggered_core_distances.append(np.asarray(triggered_shower_data.core_distance_shower))
-
-            multiplicity = np.array([len(tel_list) for tel_list in triggered_data.telescope_list])
-            trigger_multiplicity.append(multiplicity)
-            triggered_event_count += len(multiplicity)
-
-            for tel_list in triggered_data.telescope_list:
-                telescopes = tuple(sorted(str(telescope) for telescope in tel_list))
-                combination_key = ",".join(telescopes)
-                trigger_combinations[combination_key] += 1
-                for telescope in set(telescopes):
-                    telescope_participation[telescope] += 1
+        _collect_metrics_from_event_data_file(event_data_file, telescope_list, metric_accumulators)
 
     logger.info(
         f"Collected production {production_descriptor.label}: "
-        f"simulated={simulated_event_count}, triggered={triggered_event_count}"
+        f"simulated={metric_accumulators['simulated_event_count']}, "
+        f"triggered={metric_accumulators['triggered_event_count']}"
+    )
+
+    all_simulated_energies = _safe_concat(metric_accumulators["simulated_energies"])
+    all_simulated_core_distances = _safe_concat(metric_accumulators["simulated_core_distances"])
+    all_simulated_angular_distances = _safe_concat(
+        metric_accumulators["simulated_angular_distances"]
+    )
+    per_type = _build_per_type_metrics(
+        production_descriptor.label,
+        all_simulated_energies,
+        all_simulated_core_distances,
+        all_simulated_angular_distances,
+        metric_accumulators,
     )
 
     return ProductionEventMetrics(
         label=production_descriptor.label,
-        simulated_energies=_safe_concat(simulated_energies),
-        triggered_energies=_safe_concat(triggered_energies),
-        simulated_core_distances=_safe_concat(simulated_core_distances),
-        triggered_core_distances=_safe_concat(triggered_core_distances),
-        trigger_multiplicity=_safe_concat(trigger_multiplicity, dtype=int),
-        trigger_combinations=trigger_combinations,
-        telescope_participation=telescope_participation,
-        simulated_event_count=simulated_event_count,
-        triggered_event_count=triggered_event_count,
+        simulated_energies=all_simulated_energies,
+        triggered_energies=_safe_concat(metric_accumulators["triggered_energies"]),
+        simulated_core_distances=all_simulated_core_distances,
+        triggered_core_distances=_safe_concat(metric_accumulators["triggered_core_distances"]),
+        simulated_angular_distances=all_simulated_angular_distances,
+        triggered_angular_distances=_safe_concat(
+            metric_accumulators["triggered_angular_distances"]
+        ),
+        trigger_multiplicity=_safe_concat(metric_accumulators["trigger_multiplicity"], dtype=int),
+        trigger_combinations=metric_accumulators["trigger_combinations"],
+        telescope_participation=metric_accumulators["telescope_participation"],
+        simulated_event_count=metric_accumulators["simulated_event_count"],
+        triggered_event_count=metric_accumulators["triggered_event_count"],
+        per_type=per_type,
     )
 
 
+def _initialize_metric_accumulators():
+    """Initialize mutable accumulators used while scanning event files."""
+    return {
+        "simulated_energies": [],
+        "triggered_energies": [],
+        "simulated_core_distances": [],
+        "triggered_core_distances": [],
+        "simulated_angular_distances": [],
+        "triggered_angular_distances": [],
+        "trigger_multiplicity": [],
+        "trigger_combinations": Counter(),
+        "telescope_participation": Counter(),
+        "simulated_event_count": 0,
+        "triggered_event_count": 0,
+        "subset_accumulators": {
+            "energies": defaultdict(list),
+            "core_distances": defaultdict(list),
+            "angular_distances": defaultdict(list),
+            "multiplicities": defaultdict(list),
+            "counts": defaultdict(int),
+        },
+    }
+
+
+def _collect_metrics_from_event_data_file(event_data_file, telescope_list, metric_accumulators):
+    """Collect metrics from all data sets in one event data file."""
+    reader = EventDataReader(event_data_file, telescope_list=telescope_list)
+    for data_set in reader.data_sets:
+        _, shower_data, triggered_shower_data, triggered_data = reader.read_event_data(
+            event_data_file,
+            table_name_map=data_set,
+        )
+        if shower_data is None:
+            continue
+
+        _accumulate_simulated_events(shower_data, metric_accumulators)
+        if triggered_data is None or triggered_shower_data is None:
+            continue
+        _accumulate_triggered_events(triggered_shower_data, triggered_data, metric_accumulators)
+
+
+def _accumulate_simulated_events(shower_data, metric_accumulators):
+    """Accumulate simulated event quantities."""
+    metric_accumulators["simulated_energies"].append(np.asarray(shower_data.simulated_energy))
+    metric_accumulators["simulated_core_distances"].append(
+        np.asarray(shower_data.core_distance_shower)
+    )
+    metric_accumulators["simulated_angular_distances"].append(
+        np.asarray(shower_data.angular_distance)
+    )
+    metric_accumulators["simulated_event_count"] += len(shower_data.simulated_energy)
+
+
+def _accumulate_triggered_events(triggered_shower_data, triggered_data, metric_accumulators):
+    """Accumulate triggered event quantities and subset counters."""
+    trig_energies_arr = np.asarray(triggered_shower_data.simulated_energy)
+    trig_core_dist_arr = np.asarray(triggered_shower_data.core_distance_shower)
+    trig_angular_dist_arr = np.asarray(triggered_shower_data.angular_distance)
+    metric_accumulators["triggered_energies"].append(trig_energies_arr)
+    metric_accumulators["triggered_core_distances"].append(trig_core_dist_arr)
+    metric_accumulators["triggered_angular_distances"].append(trig_angular_dist_arr)
+
+    multiplicity = np.array([len(tel_list) for tel_list in triggered_data.telescope_list])
+    metric_accumulators["trigger_multiplicity"].append(multiplicity)
+    metric_accumulators["triggered_event_count"] += len(multiplicity)
+
+    for event_idx, tel_list in enumerate(triggered_data.telescope_list):
+        telescopes = tuple(sorted(str(telescope) for telescope in tel_list))
+        combination_key = ",".join(telescopes)
+        metric_accumulators["trigger_combinations"][combination_key] += 1
+        for telescope in set(telescopes):
+            metric_accumulators["telescope_participation"][telescope] += 1
+        _accumulate_per_subset(
+            telescopes,
+            float(trig_energies_arr[event_idx]),
+            float(trig_core_dist_arr[event_idx]),
+            float(trig_angular_dist_arr[event_idx]),
+            metric_accumulators["subset_accumulators"],
+        )
+
+
+def _build_per_type_metrics(
+    label,
+    all_simulated_energies,
+    all_simulated_core_distances,
+    all_simulated_angular_distances,
+    metric_accumulators,
+):
+    """Build per-subset metrics from collected accumulators."""
+    subset_accumulators = metric_accumulators["subset_accumulators"]
+    return {
+        key: ProductionEventMetrics(
+            label=label,
+            simulated_energies=all_simulated_energies,
+            triggered_energies=np.array(subset_accumulators["energies"][key]),
+            simulated_core_distances=all_simulated_core_distances,
+            triggered_core_distances=np.array(subset_accumulators["core_distances"][key]),
+            simulated_angular_distances=all_simulated_angular_distances,
+            triggered_angular_distances=np.array(subset_accumulators["angular_distances"][key]),
+            trigger_multiplicity=np.array(subset_accumulators["multiplicities"][key], dtype=int),
+            trigger_combinations=Counter(),
+            telescope_participation=Counter(),
+            simulated_event_count=metric_accumulators["simulated_event_count"],
+            triggered_event_count=subset_accumulators["counts"][key],
+        )
+        for key in subset_accumulators["energies"]
+    }
+
+
+def _accumulate_per_subset(telescopes, energy, core_dist, angular_dist, accumulators):
+    """Accumulate event data per telescope-type subset and trigger topology.
+
+    Parameters
+    ----------
+    telescopes : tuple[str]
+        Sorted telescope identifiers that triggered for this event.
+    energy : float
+        Primary energy of the simulated shower.
+    core_dist : float
+        Core distance of the simulated shower.
+    angular_dist : float
+        Angular distance of the shower from the pointing direction.
+    accumulators : dict
+        Mapping of quantity lists keyed by subset name.
+    """
+    type_counts = Counter()
+    valid_types = set()
+    for telescope in telescopes:
+        try:
+            tel_type = names.get_array_element_type_from_name(telescope)
+            type_counts[tel_type] += 1
+            valid_types.add(tel_type)
+        except ValueError:
+            pass
+    for tel_type, count in type_counts.items():
+        _accumulate_event_for_key(tel_type, count, energy, core_dist, angular_dist, accumulators)
+    if len(telescopes) == 1:
+        _accumulate_event_for_key(
+            "single_telescope", 1, energy, core_dist, angular_dist, accumulators
+        )
+    elif len(valid_types) > 1:
+        _accumulate_event_for_key(
+            "mixed_type", len(telescopes), energy, core_dist, angular_dist, accumulators
+        )
+
+
+def _accumulate_event_for_key(key, count, energy, core_dist, angular_dist, accumulators):
+    """Accumulate one event's data into a named subset key.
+
+    Parameters
+    ----------
+    key : str
+        Subset key (telescope type, ``"single_telescope"``, or ``"mixed_type"``).
+    count : int
+        Number of telescopes (of this type, or total) that triggered.
+    energy : float
+        Primary energy of the shower.
+    core_dist : float
+        Core distance of the shower.
+    angular_dist : float
+        Angular distance of the shower.
+    accumulators : dict
+        Mapping of quantity lists to append to.
+    """
+    accumulators["energies"][key].append(energy)
+    accumulators["core_distances"][key].append(core_dist)
+    accumulators["angular_distances"][key].append(angular_dist)
+    accumulators["multiplicities"][key].append(count)
+    accumulators["counts"][key] += 1
+
+
 def _safe_concat(arrays, dtype=float):
     """Concatenate list of arrays, returning an empty array for empty input."""
     if len(arrays) == 0:
diff --git a/src/simtools/visualization/plot_event_level_production_comparison.py b/src/simtools/visualization/plot_event_level_production_comparison.py
index ec04708c7b..fcaa6b7194 100644
--- a/src/simtools/visualization/plot_event_level_production_comparison.py
+++ b/src/simtools/visualization/plot_event_level_production_comparison.py
@@ -1,10 +1,13 @@
 """Plot event-level comparisons across multiple simulation productions."""
 
+import functools
 import logging
 
 import matplotlib.pyplot as plt
 import numpy as np
 
+from simtools.utils import names
+
 _logger = logging.getLogger(__name__)
 
 
@@ -20,27 +23,50 @@ def plot(metrics_per_production, output_path, bins=40):
     bins : int, optional
         Number of bins for 1D histograms.
     """
-    _plot_trigger_fraction(metrics_per_production, output_path)
     _plot_trigger_multiplicity(metrics_per_production, output_path)
     _plot_trigger_combinations(metrics_per_production, output_path)
-    _plot_triggered_vs_quantity(
-        metrics_per_production,
-        output_path,
-        quantity_name="energy",
-        x_label="Primary Energy (TeV)",
-        x_scale="log",
-        bins=bins,
-    )
-    _plot_triggered_vs_quantity(
-        metrics_per_production,
-        output_path,
-        quantity_name="core_distance",
-        x_label="Core Distance (m)",
-        x_scale="linear",
-        bins=bins,
-    )
+    _plot_single_telescope_trigger_frequencies(metrics_per_production, output_path)
+    _plot_mixed_trigger_combinations(metrics_per_production, output_path)
     _plot_telescope_participation(metrics_per_production, output_path)
 
+    for quantity_name, x_label, x_scale in _QUANTITY_CONFIGS:
+        if quantity_name in _TRIGGERED_FRACTION_QUANTITIES:
+            _plot_triggered_vs_quantity(
+                metrics_per_production,
+                output_path,
+                quantity_name=quantity_name,
+                x_label=x_label,
+                x_scale=x_scale,
+                bins=bins,
+            )
+        for cumulative in _distribution_cumulative_variants(quantity_name):
+            _plot_quantity_distribution(
+                metrics_per_production,
+                output_path,
+                quantity_name=quantity_name,
+                x_label=x_label,
+                x_scale=x_scale,
+                bins=bins,
+                cumulative=cumulative,
+            )
+
+    all_types = sorted(
+        {
+            tel_type
+            for metrics in metrics_per_production
+            for tel_type in metrics.per_type
+            if tel_type not in _SPECIAL_TRIGGER_SUBSETS
+        }
+    )
+    for tel_type in all_types:
+        type_metrics = [
+            metrics.per_type[tel_type]
+            for metrics in metrics_per_production
+            if tel_type in metrics.per_type
+        ]
+        for plot_fn in _PER_TYPE_PLOT_FNS:
+            plot_fn(type_metrics, output_path, suffix=f"_{tel_type}", bins=bins)
+
 
 def _save_figure(fig, output_path, filename):
     """Save figure and close it."""
@@ -50,7 +76,7 @@ def _save_figure(fig, output_path, filename):
     plt.close(fig)
 
 
-def _plot_trigger_fraction(metrics_per_production, output_path):
+def _plot_trigger_fraction(metrics_per_production, output_path, suffix=""):
     """Plot triggered/simulated event fractions by production."""
     labels = [metrics.label for metrics in metrics_per_production]
     values = [metrics.trigger_fraction for metrics in metrics_per_production]
@@ -58,23 +84,24 @@ def _plot_trigger_fraction(metrics_per_production, output_path):
     fig, ax = plt.subplots(figsize=(9, 6))
     bars = ax.bar(labels, values, color="tab:blue", alpha=0.85)
     ax.set_ylabel("Trigger Fraction")
-    ax.set_title("Trigger Fraction by Production")
+    type_label = f" ({suffix.lstrip('_').replace('_', ' ')})" if suffix else ""
+    ax.set_title(f"Trigger Fraction {type_label}")
     ax.set_ylim(0.0, max(1.0, max(values, default=0.0) * 1.15))
     ax.grid(axis="y", alpha=0.25)
-    for bar, value in zip(bars, values):
+    for bar_patch, value in zip(bars, values):
         ax.text(
-            bar.get_x() + bar.get_width() / 2,
-            bar.get_height(),
+            bar_patch.get_x() + bar_patch.get_width() / 2,
+            bar_patch.get_height(),
             f"{value:.3f}",
             ha="center",
             va="bottom",
             fontsize=9,
         )
 
-    _save_figure(fig, output_path, "trigger_fraction_comparison.png")
+    _save_figure(fig, output_path, f"trigger_fraction{suffix}.png")
 
 
-def _plot_trigger_multiplicity(metrics_per_production, output_path):
+def _plot_trigger_multiplicity(metrics_per_production, output_path, suffix="", bins=40):
     """Plot triggered telescope multiplicity distributions."""
     fig, ax = plt.subplots(figsize=(9, 6))
 
@@ -94,17 +121,25 @@ def _plot_trigger_multiplicity(metrics_per_production, output_path):
         if metrics.trigger_multiplicity.size == 0:
             continue
         counts, _ = np.histogram(metrics.trigger_multiplicity, bins=bins)
-        fractions = counts / counts.sum() if counts.sum() > 0 else counts
-        ax.step(bins[:-1], fractions, where="post", linewidth=2, label=metrics.label)
+        fractions, errors = _fraction_with_poisson_errors(counts)
+        stairs_artist = ax.stairs(fractions, bins, linewidth=1.5, label=metrics.label)
+        _plot_histogram_error_bars(
+            ax,
+            bins,
+            fractions,
+            errors,
+            color=_artist_color(stairs_artist),
+        )
 
     ax.set_xlabel("Triggered Telescopes per Event")
     ax.set_ylabel("Fraction of Triggered Events")
-    ax.set_title("Trigger Multiplicity Distribution")
+    type_label = f" ({suffix.lstrip('_').replace('_', ' ')})" if suffix else ""
+    ax.set_title(f"Trigger Multiplicity {type_label}")
     ax.set_xticks(bins[:-1])
     ax.grid(alpha=0.25)
     ax.legend()
 
-    _save_figure(fig, output_path, "trigger_multiplicity_comparison.png")
+    _save_figure(fig, output_path, f"trigger_multiplicity{suffix}.png")
 
 
 def _plot_trigger_combinations(metrics_per_production, output_path, top_n=12):
@@ -131,10 +166,19 @@ def _plot_trigger_combinations(metrics_per_production, output_path, top_n=12):
     fig, ax = plt.subplots(figsize=(max(10, len(selected) * 1.1), 6))
 
     for index, metrics in enumerate(metrics_per_production):
+        raw_counts = np.array([metrics.trigger_combinations.get(name, 0) for name in selected])
         event_norm = metrics.triggered_event_count if metrics.triggered_event_count > 0 else 1
-        y_values = [metrics.trigger_combinations.get(name, 0) / event_norm for name in selected]
+        y_values = raw_counts / event_norm
+        y_errors = np.sqrt(raw_counts) / event_norm
         offset = (index - (len(metrics_per_production) - 1) / 2.0) * width
-        ax.bar(x_values + offset, y_values, width=width, label=metrics.label)
+        ax.bar(
+            x_values + offset,
+            y_values,
+            width=width,
+            label=metrics.label,
+            yerr=y_errors,
+            error_kw={"elinewidth": 0.7, "capsize": 1, "capthick": 0.7},
+        )
 
     ax.set_xticks(x_values)
     ax.set_xticklabels(selected, rotation=45, ha="right")
@@ -143,7 +187,138 @@ def _plot_trigger_combinations(metrics_per_production, output_path, top_n=12):
     ax.grid(axis="y", alpha=0.25)
     ax.legend()
 
-    _save_figure(fig, output_path, "trigger_combination_comparison.png")
+    _save_figure(fig, output_path, "trigger_combination.png")
+
+
+def _plot_single_telescope_trigger_frequencies(metrics_per_production, output_path):
+    """Plot single-telescope trigger frequency distributions by telescope name."""
+    telescope_names = sorted(
+        {
+            combination
+            for metrics in metrics_per_production
+            for combination in metrics.trigger_combinations
+            if "," not in combination
+        }
+    )
+
+    if len(telescope_names) == 0:
+        _logger.warning("Skipping single-telescope trigger frequency plot, no data available.")
+        return
+
+    x_values = np.arange(len(telescope_names))
+    width = 0.8 / max(1, len(metrics_per_production))
+    fig, ax = plt.subplots(figsize=(max(10, len(telescope_names) * 0.45), 6))
+
+    for index, metrics in enumerate(metrics_per_production):
+        counts = np.array(
+            [
+                metrics.trigger_combinations.get(telescope_name, 0)
+                for telescope_name in telescope_names
+            ],
+            dtype=float,
+        )
+        fractions, errors = _fraction_with_poisson_errors(counts)
+        offset = (index - (len(metrics_per_production) - 1) / 2.0) * width
+        ax.bar(
+            x_values + offset,
+            fractions,
+            width=width,
+            label=metrics.label,
+            yerr=errors,
+            error_kw={"elinewidth": 0.7, "capsize": 1, "capthick": 0.7},
+        )
+
+    ax.set_xticks(x_values)
+    ax.set_xticklabels(telescope_names, rotation=90)
+    ax.set_ylabel("Fraction of Single-Telescope Triggers")
+    ax.set_title("Single-Telescope Trigger Distribution")
+    ax.grid(axis="y", alpha=0.25)
+    ax.legend()
+
+    _save_figure(fig, output_path, "single_telescope_trigger_distribution.png")
+
+
+def _plot_mixed_trigger_combinations(metrics_per_production, output_path):
+    """Plot mixed-type trigger combinations with multiplicity signatures and telescope names."""
+    mixed_labels = sorted(
+        {
+            _format_mixed_combination_label(combination)
+            for metrics in metrics_per_production
+            for combination in metrics.trigger_combinations
+            if _is_mixed_type_combination(combination)
+        }
+    )
+    if len(mixed_labels) == 0:
+        _logger.warning("Skipping mixed trigger combination plot, no mixed-type combinations.")
+        return
+
+    x_values = np.arange(len(mixed_labels))
+    width = 0.8 / max(1, len(metrics_per_production))
+    fig, ax = plt.subplots(figsize=(max(12, len(mixed_labels) * 0.8), 6))
+
+    for index, metrics in enumerate(metrics_per_production):
+        counts = np.array(
+            [
+                sum(
+                    count
+                    for combination, count in metrics.trigger_combinations.items()
+                    if _is_mixed_type_combination(combination)
+                    and _format_mixed_combination_label(combination) == label
+                )
+                for label in mixed_labels
+            ],
+            dtype=float,
+        )
+        fractions, errors = _fraction_with_poisson_errors(counts)
+        offset = (index - (len(metrics_per_production) - 1) / 2.0) * width
+        ax.bar(
+            x_values + offset,
+            fractions,
+            width=width,
+            label=metrics.label,
+            yerr=errors,
+            error_kw={"elinewidth": 0.7, "capsize": 1, "capthick": 0.7},
+        )
+
+    ax.set_xticks(x_values)
+    ax.set_xticklabels(mixed_labels, rotation=45, ha="right")
+    ax.set_ylabel("Fraction of Mixed-Type Triggers")
+    ax.set_title("Mixed-Type Trigger Combinations")
+    ax.grid(axis="y", alpha=0.25)
+    ax.legend()
+
+    _save_figure(fig, output_path, "mixed_trigger_combinations.png")
+
+
+def _is_mixed_type_combination(combination):
+    """Return True only for mixed signatures 1+1 and 1+2 (max multiplicity 3)."""
+    telescope_names = [name for name in combination.split(",") if name]
+    type_counts = {}
+    for telescope_name in telescope_names:
+        try:
+            tel_type = names.get_array_element_type_from_name(telescope_name)
+        except ValueError:
+            continue
+        type_counts[tel_type] = type_counts.get(tel_type, 0) + 1
+    if len(type_counts) < 2:
+        return False
+    signature = tuple(sorted(type_counts.values()))
+    return signature in {(1, 1), (1, 2)}
+
+
+def _format_mixed_combination_label(combination):
+    """Format mixed trigger combination as '<signature> | <tel1> + <tel2> + ...'."""
+    telescope_names = [name for name in combination.split(",") if name]
+    type_counts = {}
+    for telescope_name in telescope_names:
+        try:
+            tel_type = names.get_array_element_type_from_name(telescope_name)
+        except ValueError:
+            continue
+        type_counts[tel_type] = type_counts.get(tel_type, 0) + 1
+    signature = "+".join(str(count) for _, count in sorted(type_counts.items()))
+    telescopes_label = " + ".join(telescope_names)
+    return f"{signature} | {telescopes_label}"
 
 
 def _plot_triggered_vs_quantity(
@@ -153,20 +328,17 @@ def _plot_triggered_vs_quantity(
     x_label,
     x_scale,
     bins,
+    suffix="",
 ):
     """Plot simulated vs triggered distributions for one quantity."""
     fig, ax = plt.subplots(figsize=(9, 6))
 
+    plotted = False
     for metrics in metrics_per_production:
-        if quantity_name == "energy":
-            simulated = metrics.simulated_energies
-            triggered = metrics.triggered_energies
-        else:
-            simulated = metrics.simulated_core_distances
-            triggered = metrics.triggered_core_distances
-
+        simulated, triggered = _get_quantity_arrays(metrics, quantity_name)
         if simulated.size == 0:
             continue
+        plotted = True
 
         bin_edges = _get_bin_edges(simulated, x_scale=x_scale, bins=bins)
         sim_counts, _ = np.histogram(simulated, bins=bin_edges)
@@ -178,20 +350,22 @@ def _plot_triggered_vs_quantity(
                 out=np.zeros_like(trig_counts, dtype=float),
                 where=sim_counts > 0,
             )
+        _plot_series(ax, bin_edges, efficiency, metrics.label, quantity_name, linewidth=2)
 
-        x_values = 0.5 * (bin_edges[:-1] + bin_edges[1:])
-        ax.plot(x_values, efficiency, linewidth=2, label=metrics.label)
-
+    if not plotted:
+        _logger.warning(f"Skipping triggered fraction plot for {quantity_name}, no data available.")
+        plt.close(fig)
+        return
     ax.set_xlabel(x_label)
     ax.set_ylabel("Triggered / Simulated")
-    ax.set_title(f"Triggered Event Fraction vs {x_label}")
+    type_label = f" ({suffix.lstrip('_').replace('_', ' ')})" if suffix else ""
+    ax.set_title(f"Triggered Event Fraction vs {x_label}{type_label}")
     ax.set_xscale(x_scale)
     ax.set_ylim(0.0, 1.05)
     ax.grid(alpha=0.25)
     ax.legend()
 
-    filename = f"triggered_fraction_vs_{quantity_name}.png"
-    _save_figure(fig, output_path, filename)
+    _save_figure(fig, output_path, f"triggered_fraction_vs_{quantity_name}{suffix}.png")
 
 
 def _get_bin_edges(values, x_scale, bins):
@@ -210,6 +384,106 @@ def _get_bin_edges(values, x_scale, bins):
     return np.linspace(min_value, max_value, bins + 1)
 
 
+def _get_quantity_arrays(metrics, quantity_name):
+    """Return (simulated, triggered) arrays for a named quantity.
+
+    Parameters
+    ----------
+    metrics : ProductionEventMetrics
+        Metrics for one production.
+    quantity_name : str
+        One of ``"energy"``, ``"core_distance"``, or ``"angular_distance"``.
+
+    Returns
+    -------
+    tuple[np.ndarray, np.ndarray]
+        Simulated and triggered arrays for the requested quantity.
+    """
+    if quantity_name == "energy":
+        return metrics.simulated_energies, metrics.triggered_energies
+    if quantity_name == "core_distance":
+        return metrics.simulated_core_distances, metrics.triggered_core_distances
+    return metrics.simulated_angular_distances, metrics.triggered_angular_distances
+
+
+def _plot_quantity_distribution(
+    metrics_per_production,
+    output_path,
+    quantity_name,
+    x_label,
+    x_scale,
+    bins=40,
+    suffix="",
+    cumulative=False,
+):
+    """Plot simulated and triggered distributions for one quantity.
+
+    Parameters
+    ----------
+    metrics_per_production : list[ProductionEventMetrics]
+        Aggregated metrics per production.
+    output_path : pathlib.Path
+        Output directory for generated figures.
+    quantity_name : str
+        One of ``"energy"``, ``"core_distance"``, or ``"angular_distance"``.
+    x_label : str
+        Axis label for the quantity.
+    x_scale : str
+        Axis scale (``"log"`` or ``"linear"``).
+    bins : int, optional
+        Number of histogram bins.
+    suffix : str, optional
+        Filename and title suffix for per-type variants.
+    cumulative : bool, optional
+        Whether to plot cumulative distributions.
+    """
+    fig, ax = plt.subplots(figsize=(9, 6))
+    plotted = False
+    for metrics in metrics_per_production:
+        simulated, triggered = _get_quantity_arrays(metrics, quantity_name)
+        if simulated.size == 0:
+            continue
+        plotted = True
+
+        bin_edges = _get_bin_edges(simulated, x_scale=x_scale, bins=bins)
+        sim_counts, _ = np.histogram(simulated, bins=bin_edges)
+        trig_counts, _ = np.histogram(triggered, bins=bin_edges)
+        _plot_distribution_series(
+            ax,
+            bin_edges,
+            sim_counts,
+            label=f"{metrics.label} (simulated)",
+            quantity_name=quantity_name,
+            cumulative=cumulative,
+            linewidth=1.5,
+            linestyle="--",
+        )
+        _plot_distribution_series(
+            ax,
+            bin_edges,
+            trig_counts,
+            label=f"{metrics.label} (triggered)",
+            quantity_name=quantity_name,
+            cumulative=cumulative,
+            linewidth=2,
+        )
+
+    if not plotted:
+        _logger.warning(f"Skipping distribution plot for {quantity_name}, no data available.")
+        plt.close(fig)
+        return
+    cum_label = "Cumulative " if cumulative else ""
+    type_label = f" ({suffix.lstrip('_').replace('_', ' ')})" if suffix else ""
+    ax.set_xlabel(x_label)
+    ax.set_ylabel(f"{cum_label}Fraction of Events")
+    ax.set_title(f"{cum_label}Distribution: {x_label}{type_label}")
+    ax.set_xscale(x_scale)
+    ax.grid(alpha=0.25)
+    ax.legend(fontsize=8)
+    cum_str = "_cumulative" if cumulative else ""
+    _save_figure(fig, output_path, f"distribution_{quantity_name}{cum_str}{suffix}.png")
+
+
 def _plot_telescope_participation(metrics_per_production, output_path):
     """Plot telescope participation fractions by production."""
     telescopes = sorted(
@@ -229,12 +503,21 @@ def _plot_telescope_participation(metrics_per_production, output_path):
 
     for index, metrics in enumerate(metrics_per_production):
         event_norm = metrics.triggered_event_count if metrics.triggered_event_count > 0 else 1
-        fractions = [
-            metrics.telescope_participation.get(telescope, 0) / event_norm
-            for telescope in telescopes
-        ]
+        counts = np.array(
+            [metrics.telescope_participation.get(telescope, 0) for telescope in telescopes],
+            dtype=float,
+        )
+        fractions = counts / event_norm
+        errors = np.sqrt(counts) / event_norm
         offset = (index - (len(metrics_per_production) - 1) / 2.0) * width
-        ax.bar(x_values + offset, fractions, width=width, label=metrics.label)
+        ax.bar(
+            x_values + offset,
+            fractions,
+            width=width,
+            label=metrics.label,
+            yerr=errors,
+            error_kw={"elinewidth": 0.7, "capsize": 1, "capthick": 0.7},
+        )
 
     ax.set_xticks(x_values)
     ax.set_xticklabels(telescopes, rotation=90)
@@ -244,3 +527,137 @@ def _plot_telescope_participation(metrics_per_production, output_path):
     ax.legend()
 
     _save_figure(fig, output_path, "telescope_participation_fraction.png")
+
+
+def _distribution_cumulative_variants(quantity_name):
+    """Return enabled cumulative variants for one quantity."""
+    if quantity_name == "energy":
+        return (False,)
+    return (False, True)
+
+
+def _plot_series(
+    ax,
+    bin_edges,
+    values,
+    label,
+    quantity_name,
+    linewidth=2,
+    linestyle="-",
+    force_histogram=False,
+):
+    """Plot one series either as histogram stairs or as x/y line values."""
+    if force_histogram or quantity_name in _HISTOGRAM_STYLE_QUANTITIES:
+        return ax.stairs(values, bin_edges, linewidth=linewidth, linestyle=linestyle, label=label)
+    x_values = 0.5 * (bin_edges[:-1] + bin_edges[1:])
+    (line,) = ax.plot(x_values, values, linewidth=linewidth, linestyle=linestyle, label=label)
+    return line
+
+
+def _normalized_histogram_values(counts, cumulative=False):
+    """Return normalized bin values and Poisson errors for histogram counts."""
+    counts = np.asarray(counts, dtype=float)
+    total = counts.sum()
+    if total <= 0:
+        zeros = np.zeros_like(counts, dtype=float)
+        return zeros, zeros
+    if cumulative:
+        values = np.cumsum(counts) / total
+        errors = np.zeros_like(values, dtype=float)
+        return values, errors
+    values = counts / total
+    errors = np.sqrt(counts) / total
+    return values, errors
+
+
+def _plot_distribution_series(
+    ax,
+    bin_edges,
+    counts,
+    label,
+    quantity_name,
+    cumulative,
+    linewidth,
+    linestyle="-",
+):
+    """Plot one normalized distribution series and optional Poisson error bars."""
+    values, errors = _normalized_histogram_values(counts, cumulative=cumulative)
+    artist = _plot_series(
+        ax,
+        bin_edges,
+        values,
+        label,
+        quantity_name,
+        linewidth=linewidth,
+        linestyle=linestyle,
+        force_histogram=True,
+    )
+    if not cumulative:
+        _plot_histogram_error_bars(ax, bin_edges, values, errors, color=_artist_color(artist))
+
+
+def _fraction_with_poisson_errors(counts):
+    """Return normalized bin fractions and Poisson errors for counts."""
+    counts = np.asarray(counts, dtype=float)
+    total = counts.sum()
+    if total <= 0:
+        return np.zeros_like(counts, dtype=float), np.zeros_like(counts, dtype=float)
+    fractions = counts / total
+    errors = np.sqrt(counts) / total
+    return fractions, errors
+
+
+def _artist_color(artist):
+    """Return a representative color from a matplotlib artist."""
+    if artist is None:
+        return "black"
+    if hasattr(artist, "get_edgecolor"):
+        color = artist.get_edgecolor()
+        if isinstance(color, np.ndarray):
+            return color[0] if color.ndim > 1 else color
+        return color
+    if hasattr(artist, "get_color"):
+        return artist.get_color()
+    return "black"
+
+
+def _plot_histogram_error_bars(ax, bin_edges, values, errors, color):
+    """Overlay thin error bars on histogram-like series."""
+    centers = 0.5 * (bin_edges[:-1] + bin_edges[1:])
+    ax.errorbar(
+        centers,
+        values,
+        yerr=errors,
+        fmt="none",
+        ecolor=color,
+        elinewidth=0.7,
+        capsize=0,
+        alpha=0.9,
+    )
+
+
+_QUANTITY_CONFIGS = [
+    ("energy", "Primary Energy (TeV)", "log"),
+    ("core_distance", "Core Distance (m)", "linear"),
+    ("angular_distance", "Angular Distance (deg)", "linear"),
+]
+
+_HISTOGRAM_STYLE_QUANTITIES = {"energy"}
+_TRIGGERED_FRACTION_QUANTITIES = set()
+_SPECIAL_TRIGGER_SUBSETS = {"single_telescope", "mixed_type"}
+
+_PER_TYPE_PLOT_FNS = [
+    _plot_trigger_multiplicity,
+    *[
+        functools.partial(_plot_triggered_vs_quantity, quantity_name=q, x_label=lbl, x_scale=sc)
+        for q, lbl, sc in _QUANTITY_CONFIGS
+        if q in _TRIGGERED_FRACTION_QUANTITIES
+    ],
+    *[
+        functools.partial(
+            _plot_quantity_distribution, quantity_name=q, x_label=lbl, x_scale=sc, cumulative=cum
+        )
+        for q, lbl, sc in _QUANTITY_CONFIGS
+        for cum in _distribution_cumulative_variants(q)
+    ],
+]
diff --git a/tests/unit_tests/sim_events/test_production_comparison.py b/tests/unit_tests/sim_events/test_production_comparison.py
index 12a55be408..0a08ceb29b 100644
--- a/tests/unit_tests/sim_events/test_production_comparison.py
+++ b/tests/unit_tests/sim_events/test_production_comparison.py
@@ -9,10 +9,18 @@
 )
 
 
-def test_parse_production_arguments_requires_two_productions():
-    """Test parser rejects less than two productions."""
-    with pytest.raises(ValueError, match="At least two productions"):
-        parse_production_arguments([["baseline", "base.h5"]])
+def test_parse_production_arguments_accepts_single_production(mocker):
+    """Test parser accepts a single production descriptor."""
+    mocker.patch(
+        "simtools.sim_events.production_comparison.resolve_file_patterns",
+        side_effect=lambda patterns: patterns,
+    )
+
+    descriptors = parse_production_arguments([["baseline", "base.h5"]])
+
+    assert len(descriptors) == 1
+    assert descriptors[0].label == "baseline"
+    assert descriptors[0].event_data_files == ["base.h5"]
 
 
 def test_parse_production_arguments_resolves_flattened_pairs(mocker):
@@ -47,10 +55,12 @@ def test_collect_production_metrics_aggregates_event_data(mocker):
     shower_data = SimpleNamespace(
         simulated_energy=np.array([1.0, 2.0, 3.0]),
         core_distance_shower=np.array([10.0, 20.0, 30.0]),
+        angular_distance=np.array([0.5, 1.0, 1.5]),
     )
     triggered_shower_data = SimpleNamespace(
         simulated_energy=np.array([1.0, 3.0]),
         core_distance_shower=np.array([10.0, 30.0]),
+        angular_distance=np.array([0.5, 1.5]),
     )
     triggered_data = SimpleNamespace(
         telescope_list=[
@@ -82,3 +92,20 @@ def test_collect_production_metrics_aggregates_event_data(mocker):
     assert metrics[0].telescope_participation["MSTN-01"] == 2
     assert metrics[0].telescope_participation["LSTN-01"] == 1
     assert metrics[0].trigger_fraction == pytest.approx(2.0 / 3.0)
+    np.testing.assert_array_equal(metrics[0].simulated_angular_distances, np.array([0.5, 1.0, 1.5]))
+    np.testing.assert_array_equal(metrics[0].triggered_angular_distances, np.array([0.5, 1.5]))
+    assert set(metrics[0].per_type.keys()) == {"LSTN", "MSTN", "single_telescope", "mixed_type"}
+    assert metrics[0].per_type["LSTN"].triggered_event_count == 1
+    assert metrics[0].per_type["MSTN"].triggered_event_count == 2
+    np.testing.assert_array_equal(metrics[0].per_type["LSTN"].triggered_energies, [1.0])
+    np.testing.assert_array_equal(metrics[0].per_type["MSTN"].triggered_energies, [1.0, 3.0])
+    np.testing.assert_array_equal(metrics[0].per_type["LSTN"].trigger_multiplicity, [1])
+    np.testing.assert_array_equal(metrics[0].per_type["MSTN"].trigger_multiplicity, [1, 1])
+    assert metrics[0].per_type["single_telescope"].triggered_event_count == 1
+    assert metrics[0].per_type["mixed_type"].triggered_event_count == 1
+    np.testing.assert_array_equal(metrics[0].per_type["single_telescope"].triggered_energies, [3.0])
+    np.testing.assert_array_equal(metrics[0].per_type["mixed_type"].triggered_energies, [1.0])
+    np.testing.assert_array_equal(metrics[0].per_type["LSTN"].triggered_angular_distances, [0.5])
+    np.testing.assert_array_equal(
+        metrics[0].per_type["MSTN"].triggered_angular_distances, [0.5, 1.5]
+    )
diff --git a/tests/unit_tests/visualization/test_plot_event_level_production_comparison.py b/tests/unit_tests/visualization/test_plot_event_level_production_comparison.py
index edaba6452c..192bcaaf67 100644
--- a/tests/unit_tests/visualization/test_plot_event_level_production_comparison.py
+++ b/tests/unit_tests/visualization/test_plot_event_level_production_comparison.py
@@ -15,6 +15,8 @@ def _build_metrics(label, simulated_scale, triggered_scale):
         triggered_energies=np.array([0.3, 1.0]) * triggered_scale,
         simulated_core_distances=np.array([10.0, 20.0, 30.0, 40.0]),
         triggered_core_distances=np.array([20.0, 40.0]),
+        simulated_angular_distances=np.array([0.1, 0.2, 0.3, 0.4]),
+        triggered_angular_distances=np.array([0.2, 0.3]),
         trigger_multiplicity=np.array([2, 1]),
         trigger_combinations=Counter({"LSTN-01,MSTN-01": 1, "MSTN-01": 1}),
         telescope_participation=Counter({"LSTN-01": 1, "MSTN-01": 2}),
@@ -34,12 +36,62 @@ def test_plot_writes_event_level_comparison_figures(tmp_test_directory):
     plot_event_level_production_comparison.plot(metrics, output_path=output_path, bins=8)
 
     expected_files = [
-        "trigger_fraction_comparison.png",
-        "trigger_multiplicity_comparison.png",
-        "trigger_combination_comparison.png",
-        "triggered_fraction_vs_energy.png",
-        "triggered_fraction_vs_core_distance.png",
+        "trigger_multiplicity.png",
+        "trigger_combination.png",
+        "distribution_energy.png",
+        "distribution_core_distance.png",
+        "distribution_core_distance_cumulative.png",
+        "distribution_angular_distance.png",
+        "distribution_angular_distance_cumulative.png",
         "telescope_participation_fraction.png",
     ]
     for file_name in expected_files:
         assert (output_path / file_name).exists()
+
+
+def test_plot_writes_per_type_comparison_figures(tmp_test_directory):
+    """Test per-array-element-type comparison plots are written for each type."""
+    output_path = Path(tmp_test_directory)
+    per_type_lstn = _build_metrics("baseline", simulated_scale=1.0, triggered_scale=1.0)
+    per_type_mstn = _build_metrics("baseline", simulated_scale=1.0, triggered_scale=0.8)
+    metrics = [
+        ProductionEventMetrics(
+            label="baseline",
+            simulated_energies=np.array([0.1, 0.3, 1.0, 3.0]),
+            triggered_energies=np.array([0.3, 1.0]),
+            simulated_core_distances=np.array([10.0, 20.0, 30.0, 40.0]),
+            triggered_core_distances=np.array([20.0, 40.0]),
+            trigger_multiplicity=np.array([2, 1]),
+            trigger_combinations=Counter({"LSTN-01,MSTN-01": 1, "MSTN-01": 1}),
+            telescope_participation=Counter({"LSTN-01": 1, "MSTN-01": 2}),
+            simulated_event_count=4,
+            triggered_event_count=2,
+            per_type={"LSTN": per_type_lstn, "MSTN": per_type_mstn},
+        )
+    ]
+
+    plot_event_level_production_comparison.plot(metrics, output_path=output_path, bins=8)
+
+    per_type_files = [
+        "trigger_multiplicity_LSTN.png",
+        "trigger_multiplicity_MSTN.png",
+        "distribution_energy_LSTN.png",
+        "distribution_angular_distance_LSTN.png",
+        "distribution_angular_distance_cumulative_LSTN.png",
+    ]
+    for file_name in per_type_files:
+        assert (output_path / file_name).exists()
+
+
+def test_mixed_trigger_filter_accepts_only_allowed_signatures():
+    """Test mixed trigger selector keeps only 1+1 and 1+2/2+1 patterns."""
+    assert plot_event_level_production_comparison._is_mixed_type_combination("LSTN-01,MSTN-01")
+    assert plot_event_level_production_comparison._is_mixed_type_combination(
+        "LSTN-01,LSTN-02,MSTN-01"
+    )
+    assert not plot_event_level_production_comparison._is_mixed_type_combination(
+        "LSTN-01,MSTN-01,SSTS-01"
+    )
+    assert not plot_event_level_production_comparison._is_mixed_type_combination(
+        "LSTN-01,LSTN-02,LSTN-03,MSTN-01"
+    )

From 0e8060a5ed5227fe32599024b83207154fa6164f Mon Sep 17 00:00:00 2001
From: Gernot Maier <gernot.maier@desy.de>
Date: Wed, 13 May 2026 10:31:40 +0200
Subject: [PATCH 05/27] compare productions

---
 docs/source/user-guide/applications.md        |  2 +-
 ...ols-compare-productions-on-event-level.rst |  6 --
 .../simtools-compare-productions.rst          |  6 ++
 pyproject.toml                                |  3 +-
 ..._event_level.py => compare_productions.py} | 68 +++++++++++++------
 ...ductions_on_event_level_epos_urqmd_run.yml | 23 +++++++
 ...compare_productions_on_event_level_run.yml |  9 ++-
 ...test_compare_productions_on_event_level.py | 12 ++--
 8 files changed, 91 insertions(+), 38 deletions(-)
 delete mode 100644 docs/source/user-guide/applications/simtools-compare-productions-on-event-level.rst
 create mode 100644 docs/source/user-guide/applications/simtools-compare-productions.rst
 rename src/simtools/applications/{compare_productions_on_event_level.py => compare_productions.py} (59%)
 create mode 100644 tests/integration_tests/config/compare_productions_on_event_level_epos_urqmd_run.yml

diff --git a/docs/source/user-guide/applications.md b/docs/source/user-guide/applications.md
index d96cc885f0..e875bf3973 100644
--- a/docs/source/user-guide/applications.md
+++ b/docs/source/user-guide/applications.md
@@ -52,7 +52,7 @@ Parameters with the same functionality are named consistently the same among all
 simtools-convert-all-model-parameters-from-simtel <applications/simtools-convert-all-model-parameters-from-simtel>
 simtools-convert-geo-coordinates-of-array-elements <applications/simtools-convert-geo-coordinates-of-array-elements>
 simtools-convert-model-parameter-from-simtel <applications/simtools-convert-model-parameter-from-simtel>
-simtools-compare-productions-on-event-level <applications/simtools-compare-productions-on-event-level>
+simtools-compare-productions <applications/simtools-compare-productions>
 simtools-db-add-file-to-db <applications/simtools-db-add-file-to-db>
 simtools-db-add-simulation-model-from-repository-to-db <applications/simtools-db-add-simulation-model-from-repository-to-db>
 simtools-db-add-value-from-json-to-db <applications/simtools-db-add-value-from-json-to-db>
diff --git a/docs/source/user-guide/applications/simtools-compare-productions-on-event-level.rst b/docs/source/user-guide/applications/simtools-compare-productions-on-event-level.rst
deleted file mode 100644
index d23ca2d9b3..0000000000
--- a/docs/source/user-guide/applications/simtools-compare-productions-on-event-level.rst
+++ /dev/null
@@ -1,6 +0,0 @@
-
-simtools-compare-productions-on-event-level
-===========================================
-
-.. automodule:: compare_productions_on_event_level
-   :members:
diff --git a/docs/source/user-guide/applications/simtools-compare-productions.rst b/docs/source/user-guide/applications/simtools-compare-productions.rst
new file mode 100644
index 0000000000..76acf76fbe
--- /dev/null
+++ b/docs/source/user-guide/applications/simtools-compare-productions.rst
@@ -0,0 +1,6 @@
+
+simtools-compare-productions
+============================
+
+.. automodule:: compare_productions
+   :members:
diff --git a/pyproject.toml b/pyproject.toml
index fae01d68b8..3846752a64 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -73,7 +73,8 @@ optional-dependencies.tests = [
 urls."bug tracker" = "https://github.com/gammasim/simtools/issues"
 urls.documentation = "https://gammasim.github.io/simtools/"
 urls.repository = "https://github.com/gammasim/simtools"
-scripts.simtools-compare-productions-on-event-level = "simtools.applications.compare_productions_on_event_level:main"
+scripts.simtools-compare-productions = "simtools.applications.compare_productions:main"
+scripts.simtools-compare-productions-on-event-level = "simtools.applications.compare_productions:main"
 scripts.simtools-convert-all-model-parameters-from-simtel = "simtools.applications.convert_all_model_parameters_from_simtel:main"
 scripts.simtools-convert-geo-coordinates-of-array-elements = "simtools.applications.convert_geo_coordinates_of_array_elements:main"
 scripts.simtools-convert-model-parameter-from-simtel = "simtools.applications.convert_model_parameter_from_simtel:main"
diff --git a/src/simtools/applications/compare_productions_on_event_level.py b/src/simtools/applications/compare_productions.py
similarity index 59%
rename from src/simtools/applications/compare_productions_on_event_level.py
rename to src/simtools/applications/compare_productions.py
index 933413e074..7506ef338a 100644
--- a/src/simtools/applications/compare_productions_on_event_level.py
+++ b/src/simtools/applications/compare_productions.py
@@ -1,9 +1,9 @@
 #!/usr/bin/python3
 
-r"""Compare multiple simulation productions on event level.
+r"""Compare simulation productions at different comparison levels.
 
-The application accepts repeated production descriptors and creates event-level
-comparison plots using reduced event data files.
+The application accepts repeated production descriptors and dispatches to
+comparison-level specific implementations.
 
 Command line arguments
 ----------------------
@@ -11,19 +11,13 @@
     Production descriptor in two fields:
     1) label
     2) comma-separated event data file patterns
+comparison_level (str, optional)
+    Comparison level selector. Supported values are:
+    - events
+    - signals
+    - compute
 output_path (str, required)
     Output directory for generated comparison plots.
-
-Examples
---------
-Compare two productions:
-
-.. code-block:: console
-
-    simtools-compare-productions-on-event-level \
-        --production baseline "data/baseline/*.h5" \
-        --production candidate "data/candidate/*.h5,data/candidate_extra/*.h5" \
-        --output_path simtools-output/
 """
 
 from simtools.application_control import build_application
@@ -49,6 +43,12 @@ def _add_arguments(parser):
             "Repeat this argument as needed for multiple productions."
         ),
     )
+    parser.add_argument(
+        "--comparison_level",
+        choices=["events", "signals", "compute"],
+        default="events",
+        help="Comparison level to execute.",
+    )
     parser.add_argument(
         "--telescope_ids",
         nargs="+",
@@ -57,12 +57,8 @@ def _add_arguments(parser):
     )
 
 
-def main():
-    """See CLI description."""
-    app_context = build_application(
-        initialization_kwargs={"db_config": False, "output": True},
-    )
-
+def _run_event_comparison(app_context):
+    """Run event-level production comparison."""
     production_descriptors = parse_production_arguments(app_context.args["production"])
     metrics_per_production = collect_production_metrics(
         production_descriptors,
@@ -71,9 +67,39 @@ def main():
 
     output_directory = app_context.io_handler.get_output_directory()
     plot_event_level_production_comparison.plot(
-        metrics_per_production, output_path=output_directory
+        metrics_per_production,
+        output_path=output_directory,
+    )
+
+
+def _run_placeholder_comparison(comparison_level):
+    """Warn for future comparison-level implementations.
+
+    Notes
+    -----
+    Future work tracked in:
+    - signals level: issue #2183
+    - compute level: issue #2184
+    """
+    raise NotImplementedError(
+        f"Comparison level '{comparison_level}' is not implemented yet "
+        "(planned via issues #2183 and #2184)."
+    )
+
+
+def main():
+    """See CLI description."""
+    app_context = build_application(
+        initialization_kwargs={"db_config": False, "output": True},
     )
 
+    comparison_level = app_context.args["comparison_level"]
+    if comparison_level == "events":
+        _run_event_comparison(app_context)
+        return
+
+    _run_placeholder_comparison(comparison_level)
+
 
 if __name__ == "__main__":
     main()
diff --git a/tests/integration_tests/config/compare_productions_on_event_level_epos_urqmd_run.yml b/tests/integration_tests/config/compare_productions_on_event_level_epos_urqmd_run.yml
new file mode 100644
index 0000000000..2215763f4e
--- /dev/null
+++ b/tests/integration_tests/config/compare_productions_on_event_level_epos_urqmd_run.yml
@@ -0,0 +1,23 @@
+---
+applications:
+- application: simtools-compare-productions
+  configuration:
+    comparison_level: events
+    production:
+    - epos
+    - <EPOS_EVENT_DATA_GLOB_PATTERN>
+    - urqmd
+    - <URQMD_EVENT_DATA_GLOB_PATTERN>
+    output_path: simtools-output
+  integration_tests:
+  - output_file: trigger_multiplicity.png
+  - output_file: trigger_combination.png
+  - output_file: single_telescope_trigger_distribution.png
+  - output_file: mixed_trigger_combinations.png
+  - output_file: distribution_energy.png
+  - output_file: distribution_core_distance.png
+  - output_file: distribution_angular_distance.png
+  - output_file: telescope_participation_fraction.png
+  test_name: run_epos_vs_urqmd
+schema_name: application_workflow.metaschema
+schema_version: 0.4.0
diff --git a/tests/integration_tests/config/compare_productions_on_event_level_run.yml b/tests/integration_tests/config/compare_productions_on_event_level_run.yml
index 9aa1c408f7..2d202ea5f9 100644
--- a/tests/integration_tests/config/compare_productions_on_event_level_run.yml
+++ b/tests/integration_tests/config/compare_productions_on_event_level_run.yml
@@ -1,7 +1,8 @@
 ---
 applications:
-- application: simtools-compare-productions-on-event-level
+- application: simtools-compare-productions
   configuration:
+    comparison_level: events
     production:
     - baseline
     - tests/resources/proton_za20deg_azm000deg_North_alpha_6.0.0_reduced_event_data.hdf5
@@ -9,8 +10,10 @@ applications:
     - tests/resources/proton_za20deg_azm000deg_North_alpha_6.0.0_reduced_event_data.hdf5
     output_path: simtools-output
   integration_tests:
-  - output_file: trigger_fraction_comparison.png
-  - output_file: trigger_multiplicity_comparison.png
+  - output_file: trigger_multiplicity.png
+  - output_file: trigger_combination.png
+  - output_file: single_telescope_trigger_distribution.png
+  - output_file: mixed_trigger_combinations.png
   - output_file: telescope_participation_fraction.png
   test_name: run
 schema_name: application_workflow.metaschema
diff --git a/tests/unit_tests/applications/test_compare_productions_on_event_level.py b/tests/unit_tests/applications/test_compare_productions_on_event_level.py
index 9e6d3725e4..af2236a546 100644
--- a/tests/unit_tests/applications/test_compare_productions_on_event_level.py
+++ b/tests/unit_tests/applications/test_compare_productions_on_event_level.py
@@ -2,7 +2,7 @@
 from types import SimpleNamespace
 from unittest.mock import MagicMock, patch
 
-import simtools.applications.compare_productions_on_event_level as app
+import simtools.applications.compare_productions as app
 
 
 def test_main_collects_metrics_and_plots(tmp_test_directory):
@@ -11,6 +11,7 @@ def test_main_collects_metrics_and_plots(tmp_test_directory):
     app_context = SimpleNamespace(
         args={
             "production": [["baseline", "base_*.h5"], ["candidate", "cand_*.h5"]],
+            "comparison_level": "events",
             "telescope_ids": ["LSTN-01"],
         },
         io_handler=MagicMock(),
@@ -22,20 +23,19 @@ def test_main_collects_metrics_and_plots(tmp_test_directory):
 
     with (
         patch(
-            "simtools.applications.compare_productions_on_event_level.build_application",
+            "simtools.applications.compare_productions.build_application",
             return_value=app_context,
         ),
         patch(
-            "simtools.applications.compare_productions_on_event_level.parse_production_arguments",
+            "simtools.applications.compare_productions.parse_production_arguments",
             return_value=parsed_productions,
         ) as mock_parse,
         patch(
-            "simtools.applications.compare_productions_on_event_level.collect_production_metrics",
+            "simtools.applications.compare_productions.collect_production_metrics",
             return_value=collected_metrics,
         ) as mock_collect,
         patch(
-            "simtools.applications.compare_productions_on_event_level."
-            "plot_event_level_production_comparison.plot"
+            "simtools.applications.compare_productions.plot_event_level_production_comparison.plot"
         ) as mock_plot,
     ):
         app.main()

From 05fe632a201900afb38848b3e56f8fbc7cc90817 Mon Sep 17 00:00:00 2001
From: Gernot Maier <gernot.maier@desy.de>
Date: Wed, 13 May 2026 10:43:23 +0200
Subject: [PATCH 06/27] unit tests

---
 .../sim_events/test_production_comparison.py  | 114 ++++++++++++++
 ..._plot_event_level_production_comparison.py | 147 +++++++++++++++++-
 2 files changed, 257 insertions(+), 4 deletions(-)

diff --git a/tests/unit_tests/sim_events/test_production_comparison.py b/tests/unit_tests/sim_events/test_production_comparison.py
index 0a08ceb29b..db172b0446 100644
--- a/tests/unit_tests/sim_events/test_production_comparison.py
+++ b/tests/unit_tests/sim_events/test_production_comparison.py
@@ -109,3 +109,117 @@ def test_collect_production_metrics_aggregates_event_data(mocker):
     np.testing.assert_array_equal(
         metrics[0].per_type["MSTN"].triggered_angular_distances, [0.5, 1.5]
     )
+
+
+@pytest.mark.parametrize(
+    ("arguments", "error_match"),
+    [
+        ([], "At least one production is required"),
+        (["baseline", "base.h5", "dangling"], "label/file pairs"),
+        ([["baseline", "  ,   "]], "has no event_data_file pattern"),
+        ([["baseline", "a.h5"], ["candidate", 1]], "label/file pairs"),
+    ],
+)
+def test_parse_production_arguments_error_paths(mocker, arguments, error_match):
+    """Test parser validation failures for malformed production arguments."""
+    mocker.patch(
+        "simtools.sim_events.production_comparison.resolve_file_patterns",
+        side_effect=lambda patterns: patterns,
+    )
+
+    with pytest.raises(ValueError, match=error_match):
+        parse_production_arguments(arguments)
+
+
+def test_parse_production_arguments_rejects_unresolved_patterns(mocker):
+    """Test parser rejects productions that resolve to no files."""
+    mocker.patch("simtools.sim_events.production_comparison.resolve_file_patterns", return_value=[])
+
+    with pytest.raises(ValueError, match="does not resolve to any files"):
+        parse_production_arguments([["baseline", "missing_*.h5"]])
+
+
+def test_parse_production_arguments_accepts_nested_flattened_strings(mocker):
+    """Test parser accepts nested flattened string groups."""
+    mocker.patch(
+        "simtools.sim_events.production_comparison.resolve_file_patterns",
+        side_effect=lambda patterns: patterns,
+    )
+
+    descriptors = parse_production_arguments([["baseline", "a.h5", "candidate", "b.h5"]])
+
+    assert [descriptor.label for descriptor in descriptors] == ["baseline", "candidate"]
+
+
+def test_collect_production_metrics_handles_empty_and_partial_datasets(mocker):
+    """Test collector skips missing shower and missing trigger data sets gracefully."""
+    descriptor = SimpleNamespace(label="baseline", event_data_files=["baseline_file.h5"])
+    shower_data = SimpleNamespace(
+        simulated_energy=np.array([2.0]),
+        core_distance_shower=np.array([20.0]),
+        angular_distance=np.array([1.0]),
+    )
+
+    mock_reader = mocker.patch("simtools.sim_events.production_comparison.EventDataReader")
+    mock_reader.return_value.data_sets = [{"A": "A"}, {"B": "B"}, {"C": "C"}]
+    mock_reader.return_value.read_event_data.side_effect = [
+        (SimpleNamespace(), None, None, None),
+        (SimpleNamespace(), shower_data, None, None),
+        (SimpleNamespace(), shower_data, shower_data, None),
+    ]
+
+    metrics = collect_production_metrics([descriptor])
+
+    assert metrics[0].simulated_event_count == 2
+    assert metrics[0].triggered_event_count == 0
+    assert metrics[0].trigger_fraction == 0.0
+    assert metrics[0].trigger_combinations == {}
+    np.testing.assert_array_equal(metrics[0].triggered_energies, np.array([]))
+
+
+def test_collect_production_metrics_empty_input_keeps_arrays_empty(mocker):
+    """Test empty shower input yields empty arrays and zero trigger fraction."""
+    descriptor = SimpleNamespace(label="baseline", event_data_files=["baseline_file.h5"])
+
+    mock_reader = mocker.patch("simtools.sim_events.production_comparison.EventDataReader")
+    mock_reader.return_value.data_sets = [{"SHOWERS": "SHOWERS"}]
+    mock_reader.return_value.read_event_data.return_value = (SimpleNamespace(), None, None, None)
+
+    metrics = collect_production_metrics([descriptor])
+
+    assert metrics[0].simulated_event_count == 0
+    assert metrics[0].triggered_event_count == 0
+    assert metrics[0].trigger_fraction == 0.0
+    np.testing.assert_array_equal(metrics[0].simulated_energies, np.array([]))
+    np.testing.assert_array_equal(metrics[0].trigger_multiplicity, np.array([], dtype=int))
+
+
+def test_collect_production_metrics_ignores_unknown_telescope_type(mocker):
+    """Test unknown telescope names are ignored for type bins but still counted as single triggers."""
+    descriptor = SimpleNamespace(label="baseline", event_data_files=["baseline_file.h5"])
+    shower_data = SimpleNamespace(
+        simulated_energy=np.array([1.0]),
+        core_distance_shower=np.array([10.0]),
+        angular_distance=np.array([0.5]),
+    )
+    triggered_shower_data = SimpleNamespace(
+        simulated_energy=np.array([1.0]),
+        core_distance_shower=np.array([10.0]),
+        angular_distance=np.array([0.5]),
+    )
+    triggered_data = SimpleNamespace(telescope_list=[np.array(["BAD-01"])])
+
+    mock_reader = mocker.patch("simtools.sim_events.production_comparison.EventDataReader")
+    mock_reader.return_value.data_sets = [{"SHOWERS": "SHOWERS", "TRIGGERS": "TRIGGERS"}]
+    mock_reader.return_value.read_event_data.return_value = (
+        SimpleNamespace(),
+        shower_data,
+        triggered_shower_data,
+        triggered_data,
+    )
+
+    metrics = collect_production_metrics([descriptor])
+
+    assert "single_telescope" in metrics[0].per_type
+    assert metrics[0].per_type["single_telescope"].triggered_event_count == 1
+    assert "BAD" not in metrics[0].per_type
diff --git a/tests/unit_tests/visualization/test_plot_event_level_production_comparison.py b/tests/unit_tests/visualization/test_plot_event_level_production_comparison.py
index 192bcaaf67..ef3b0058ff 100644
--- a/tests/unit_tests/visualization/test_plot_event_level_production_comparison.py
+++ b/tests/unit_tests/visualization/test_plot_event_level_production_comparison.py
@@ -2,6 +2,8 @@
 from pathlib import Path
 
 import numpy as np
+import pytest
+from matplotlib import pyplot as plt
 
 from simtools.sim_events.production_comparison import ProductionEventMetrics
 from simtools.visualization import plot_event_level_production_comparison
@@ -25,6 +27,12 @@ def _build_metrics(label, simulated_scale, triggered_scale):
     )
 
 
+def _assert_files_exist(output_path, file_names):
+    """Assert all expected output files exist in output path."""
+    for file_name in file_names:
+        assert (output_path / file_name).exists()
+
+
 def test_plot_writes_event_level_comparison_figures(tmp_test_directory):
     """Test full comparison plotting writes all expected output files."""
     output_path = Path(tmp_test_directory)
@@ -45,8 +53,7 @@ def test_plot_writes_event_level_comparison_figures(tmp_test_directory):
         "distribution_angular_distance_cumulative.png",
         "telescope_participation_fraction.png",
     ]
-    for file_name in expected_files:
-        assert (output_path / file_name).exists()
+    _assert_files_exist(output_path, expected_files)
 
 
 def test_plot_writes_per_type_comparison_figures(tmp_test_directory):
@@ -79,8 +86,7 @@ def test_plot_writes_per_type_comparison_figures(tmp_test_directory):
         "distribution_angular_distance_LSTN.png",
         "distribution_angular_distance_cumulative_LSTN.png",
     ]
-    for file_name in per_type_files:
-        assert (output_path / file_name).exists()
+    _assert_files_exist(output_path, per_type_files)
 
 
 def test_mixed_trigger_filter_accepts_only_allowed_signatures():
@@ -95,3 +101,136 @@ def test_mixed_trigger_filter_accepts_only_allowed_signatures():
     assert not plot_event_level_production_comparison._is_mixed_type_combination(
         "LSTN-01,LSTN-02,LSTN-03,MSTN-01"
     )
+
+
+@pytest.mark.parametrize(
+    ("quantity_name", "expected"),
+    [("energy", (False,)), ("core_distance", (False, True)), ("angular_distance", (False, True))],
+)
+def test_distribution_cumulative_variants(quantity_name, expected):
+    """Test cumulative variant selection by quantity."""
+    assert (
+        plot_event_level_production_comparison._distribution_cumulative_variants(quantity_name)
+        == expected
+    )
+
+
+@pytest.mark.parametrize("cumulative", [False, True])
+def test_normalized_histogram_values(cumulative):
+    """Test normalized histogram values and Poisson errors helper."""
+    values, errors = plot_event_level_production_comparison._normalized_histogram_values(
+        np.array([1.0, 3.0]), cumulative=cumulative
+    )
+    assert values.shape == (2,)
+    assert errors.shape == (2,)
+    if cumulative:
+        np.testing.assert_allclose(values, np.array([0.25, 1.0]))
+        np.testing.assert_array_equal(errors, np.array([0.0, 0.0]))
+    else:
+        np.testing.assert_allclose(values, np.array([0.25, 0.75]))
+        np.testing.assert_allclose(errors, np.sqrt(np.array([1.0, 3.0])) / 4.0)
+
+
+def test_normalized_histogram_values_zero_counts():
+    """Test zero-count handling in histogram normalization helper."""
+    values, errors = plot_event_level_production_comparison._normalized_histogram_values(
+        np.array([0.0, 0.0]), cumulative=False
+    )
+    np.testing.assert_array_equal(values, np.array([0.0, 0.0]))
+    np.testing.assert_array_equal(errors, np.array([0.0, 0.0]))
+
+
+def test_plot_series_and_artist_color():
+    """Test line and histogram rendering branches in _plot_series and artist color extraction."""
+    fig, ax = plt.subplots(figsize=(4, 3))
+    bin_edges = np.array([1.0, 2.0, 3.0])
+    values = np.array([0.4, 0.6])
+    line_artist = plot_event_level_production_comparison._plot_series(
+        ax, bin_edges, values, "line", "core_distance"
+    )
+    hist_artist = plot_event_level_production_comparison._plot_series(
+        ax, bin_edges, values, "hist", "core_distance", force_histogram=True
+    )
+    assert line_artist is not None
+    assert hist_artist is not None
+    assert plot_event_level_production_comparison._artist_color(None) == "black"
+    assert plot_event_level_production_comparison._artist_color(line_artist) is not None
+    assert plot_event_level_production_comparison._artist_color(hist_artist) is not None
+    plt.close(fig)
+
+
+def test_trigger_fraction_and_skip_paths(tmp_test_directory):
+    """Test trigger fraction output and skip branches with empty inputs."""
+    output_path = Path(tmp_test_directory)
+    metrics = [_build_metrics("baseline", simulated_scale=1.0, triggered_scale=1.0)]
+    plot_event_level_production_comparison._plot_trigger_fraction(metrics, output_path)
+    assert (output_path / "trigger_fraction.png").exists()
+
+    empty_metric = _build_metrics("empty", simulated_scale=1.0, triggered_scale=1.0)
+    empty_metric.trigger_multiplicity = np.array([], dtype=int)
+    empty_metric.trigger_combinations = Counter()
+    empty_metric.telescope_participation = Counter()
+    plot_event_level_production_comparison._plot_trigger_multiplicity([empty_metric], output_path)
+    plot_event_level_production_comparison._plot_trigger_combinations([empty_metric], output_path)
+    plot_event_level_production_comparison._plot_telescope_participation(
+        [empty_metric], output_path
+    )
+    assert not (output_path / "trigger_multiplicity_empty.png").exists()
+
+
+def test_plot_invokes_triggered_fraction_branch_when_enabled(tmp_test_directory, mocker):
+    """Test branch enabling triggered-fraction plotting in top-level plot function."""
+    output_path = Path(tmp_test_directory)
+    metrics = [_build_metrics("baseline", simulated_scale=1.0, triggered_scale=1.0)]
+    mock_trigger = mocker.patch(
+        "simtools.visualization.plot_event_level_production_comparison._plot_triggered_vs_quantity"
+    )
+    mocker.patch(
+        "simtools.visualization.plot_event_level_production_comparison._TRIGGERED_FRACTION_QUANTITIES",
+        {"core_distance"},
+    )
+
+    plot_event_level_production_comparison.plot(metrics, output_path=output_path, bins=8)
+
+    mock_trigger.assert_called_once()
+
+
+def test_triggered_vs_quantity_outputs_and_empty_skip(tmp_test_directory):
+    """Test triggered-vs-quantity output and empty-input skip path."""
+    output_path = Path(tmp_test_directory)
+    metrics = [_build_metrics("baseline", simulated_scale=1.0, triggered_scale=1.0)]
+    plot_event_level_production_comparison._plot_triggered_vs_quantity(
+        metrics,
+        output_path,
+        quantity_name="core_distance",
+        x_label="Core Distance (m)",
+        x_scale="linear",
+        bins=8,
+    )
+    assert (output_path / "triggered_fraction_vs_core_distance.png").exists()
+
+    empty_metric = _build_metrics("empty", simulated_scale=1.0, triggered_scale=1.0)
+    empty_metric.simulated_core_distances = np.array([], dtype=float)
+    plot_event_level_production_comparison._plot_triggered_vs_quantity(
+        [empty_metric],
+        output_path,
+        quantity_name="core_distance",
+        x_label="Core Distance (m)",
+        x_scale="linear",
+        bins=8,
+        suffix="_empty",
+    )
+    assert not (output_path / "triggered_fraction_vs_core_distance_empty.png").exists()
+
+
+def test_single_and_mixed_trigger_skip_paths(tmp_test_directory):
+    """Test skip branches when no single or mixed trigger combinations exist."""
+    output_path = Path(tmp_test_directory)
+    metric = _build_metrics("baseline", simulated_scale=1.0, triggered_scale=1.0)
+    metric.trigger_combinations = Counter({"LSTN-01,MSTN-01,SSTS-01": 2})
+    plot_event_level_production_comparison._plot_single_telescope_trigger_frequencies(
+        [metric], output_path
+    )
+    plot_event_level_production_comparison._plot_mixed_trigger_combinations([metric], output_path)
+    assert not (output_path / "single_telescope_trigger_distribution.png").exists()
+    assert not (output_path / "mixed_trigger_combinations.png").exists()

From 4229f0465bb27bc4e3b124ac147cfd661cbb01c9 Mon Sep 17 00:00:00 2001
From: Gernot Maier <gernot.maier@desy.de>
Date: Wed, 13 May 2026 10:56:15 +0200
Subject: [PATCH 07/27] complexity

---
 .../sim_events/production_comparison.py       |  64 ++---
 .../plot_event_level_production_comparison.py | 227 +++++++++---------
 2 files changed, 144 insertions(+), 147 deletions(-)

diff --git a/src/simtools/sim_events/production_comparison.py b/src/simtools/sim_events/production_comparison.py
index 167c52abe1..5d6f5d1f63 100644
--- a/src/simtools/sim_events/production_comparison.py
+++ b/src/simtools/sim_events/production_comparison.py
@@ -91,28 +91,35 @@ def _normalize_production_arguments(production_arguments):
     if not production_arguments:
         return []
 
+    normalized = []
     if all(isinstance(item, str) for item in production_arguments):
-        if len(production_arguments) % 2 != 0:
-            raise ValueError("Production arguments must be provided as label/file pairs.")
-        return [
-            (production_arguments[index], production_arguments[index + 1])
-            for index in range(0, len(production_arguments), 2)
-        ]
+        return _pairwise_label_file_arguments(production_arguments)
 
-    normalized = []
     for item in production_arguments:
-        if isinstance(item, list | tuple):
-            if len(item) == 2 and all(isinstance(value, str) for value in item):
-                normalized.append((item[0], item[1]))
-                continue
+        normalized.extend(_normalize_single_production_argument(item))
+
+    return normalized
 
-            if all(isinstance(value, str) for value in item):
-                normalized.extend(_normalize_production_arguments(list(item)))
-                continue
 
+def _pairwise_label_file_arguments(flat_arguments):
+    """Convert a flat list of strings into ``[(label, files), ...]`` pairs."""
+    if len(flat_arguments) % 2 != 0:
         raise ValueError("Production arguments must be provided as label/file pairs.")
+    return [
+        (flat_arguments[index], flat_arguments[index + 1])
+        for index in range(0, len(flat_arguments), 2)
+    ]
 
-    return normalized
+
+def _normalize_single_production_argument(argument):
+    """Normalize one nested production argument into label/file pairs."""
+    if not isinstance(argument, list | tuple):
+        raise ValueError("Production arguments must be provided as label/file pairs.")
+    if not all(isinstance(value, str) for value in argument):
+        raise ValueError("Production arguments must be provided as label/file pairs.")
+    if len(argument) == 2:
+        return [(argument[0], argument[1])]
+    return _pairwise_label_file_arguments(list(argument))
 
 
 def collect_production_metrics(production_descriptors, telescope_list=None):
@@ -307,25 +314,26 @@ def _accumulate_per_subset(telescopes, energy, core_dist, angular_dist, accumula
     accumulators : dict
         Mapping of quantity lists keyed by subset name.
     """
+    for key, count in _subset_counts_for_event(telescopes).items():
+        _accumulate_event_for_key(key, count, energy, core_dist, angular_dist, accumulators)
+
+
+def _subset_counts_for_event(telescopes):
+    """Return subset keys and multiplicities for one triggered event."""
     type_counts = Counter()
-    valid_types = set()
     for telescope in telescopes:
         try:
             tel_type = names.get_array_element_type_from_name(telescope)
-            type_counts[tel_type] += 1
-            valid_types.add(tel_type)
         except ValueError:
-            pass
-    for tel_type, count in type_counts.items():
-        _accumulate_event_for_key(tel_type, count, energy, core_dist, angular_dist, accumulators)
+            continue
+        type_counts[tel_type] += 1
+
+    subset_counts = dict(type_counts)
     if len(telescopes) == 1:
-        _accumulate_event_for_key(
-            "single_telescope", 1, energy, core_dist, angular_dist, accumulators
-        )
-    elif len(valid_types) > 1:
-        _accumulate_event_for_key(
-            "mixed_type", len(telescopes), energy, core_dist, angular_dist, accumulators
-        )
+        subset_counts["single_telescope"] = 1
+    elif len(type_counts) > 1:
+        subset_counts["mixed_type"] = len(telescopes)
+    return subset_counts
 
 
 def _accumulate_event_for_key(key, count, energy, core_dist, angular_dist, accumulators):
diff --git a/src/simtools/visualization/plot_event_level_production_comparison.py b/src/simtools/visualization/plot_event_level_production_comparison.py
index fcaa6b7194..466640ab31 100644
--- a/src/simtools/visualization/plot_event_level_production_comparison.py
+++ b/src/simtools/visualization/plot_event_level_production_comparison.py
@@ -161,33 +161,21 @@ def _plot_trigger_combinations(metrics_per_production, output_path, top_n=12):
     selected = [name for name, _ in sorted(totals.items(), key=lambda item: item[1], reverse=True)]
     selected = selected[:top_n]
 
-    x_values = np.arange(len(selected))
-    width = 0.8 / max(1, len(metrics_per_production))
-    fig, ax = plt.subplots(figsize=(max(10, len(selected) * 1.1), 6))
-
-    for index, metrics in enumerate(metrics_per_production):
-        raw_counts = np.array([metrics.trigger_combinations.get(name, 0) for name in selected])
-        event_norm = metrics.triggered_event_count if metrics.triggered_event_count > 0 else 1
-        y_values = raw_counts / event_norm
-        y_errors = np.sqrt(raw_counts) / event_norm
-        offset = (index - (len(metrics_per_production) - 1) / 2.0) * width
-        ax.bar(
-            x_values + offset,
-            y_values,
-            width=width,
-            label=metrics.label,
-            yerr=y_errors,
-            error_kw={"elinewidth": 0.7, "capsize": 1, "capthick": 0.7},
-        )
-
-    ax.set_xticks(x_values)
-    ax.set_xticklabels(selected, rotation=45, ha="right")
-    ax.set_ylabel("Fraction of Triggered Events")
-    ax.set_title("Top Trigger Combinations")
-    ax.grid(axis="y", alpha=0.25)
-    ax.legend()
-
-    _save_figure(fig, output_path, "trigger_combination.png")
+    _plot_grouped_fraction_bars(
+        metrics_per_production,
+        categories=selected,
+        counts_getter=lambda metrics, category_labels: [
+            metrics.trigger_combinations.get(name, 0) for name in category_labels
+        ],
+        normalization_fn=_fractions_per_triggered_events,
+        output_path=output_path,
+        filename="trigger_combination.png",
+        y_label="Fraction of Triggered Events",
+        title="Top Trigger Combinations",
+        figure_width=max(10, len(selected) * 1.1),
+        x_rotation=45,
+        x_ha="right",
+    )
 
 
 def _plot_single_telescope_trigger_frequencies(metrics_per_production, output_path):
@@ -205,37 +193,22 @@ def _plot_single_telescope_trigger_frequencies(metrics_per_production, output_pa
         _logger.warning("Skipping single-telescope trigger frequency plot, no data available.")
         return
 
-    x_values = np.arange(len(telescope_names))
-    width = 0.8 / max(1, len(metrics_per_production))
-    fig, ax = plt.subplots(figsize=(max(10, len(telescope_names) * 0.45), 6))
-
-    for index, metrics in enumerate(metrics_per_production):
-        counts = np.array(
-            [
-                metrics.trigger_combinations.get(telescope_name, 0)
-                for telescope_name in telescope_names
-            ],
-            dtype=float,
-        )
-        fractions, errors = _fraction_with_poisson_errors(counts)
-        offset = (index - (len(metrics_per_production) - 1) / 2.0) * width
-        ax.bar(
-            x_values + offset,
-            fractions,
-            width=width,
-            label=metrics.label,
-            yerr=errors,
-            error_kw={"elinewidth": 0.7, "capsize": 1, "capthick": 0.7},
-        )
-
-    ax.set_xticks(x_values)
-    ax.set_xticklabels(telescope_names, rotation=90)
-    ax.set_ylabel("Fraction of Single-Telescope Triggers")
-    ax.set_title("Single-Telescope Trigger Distribution")
-    ax.grid(axis="y", alpha=0.25)
-    ax.legend()
-
-    _save_figure(fig, output_path, "single_telescope_trigger_distribution.png")
+    _plot_grouped_fraction_bars(
+        metrics_per_production,
+        categories=telescope_names,
+        counts_getter=lambda metrics, category_labels: [
+            metrics.trigger_combinations.get(telescope_name, 0)
+            for telescope_name in category_labels
+        ],
+        normalization_fn=_fraction_with_poisson_errors,
+        output_path=output_path,
+        filename="single_telescope_trigger_distribution.png",
+        y_label="Fraction of Single-Telescope Triggers",
+        title="Single-Telescope Trigger Distribution",
+        figure_width=max(10, len(telescope_names) * 0.45),
+        x_rotation=90,
+        x_ha="center",
+    )
 
 
 def _plot_mixed_trigger_combinations(metrics_per_production, output_path):
@@ -252,24 +225,56 @@ def _plot_mixed_trigger_combinations(metrics_per_production, output_path):
         _logger.warning("Skipping mixed trigger combination plot, no mixed-type combinations.")
         return
 
-    x_values = np.arange(len(mixed_labels))
+    _plot_grouped_fraction_bars(
+        metrics_per_production,
+        categories=mixed_labels,
+        counts_getter=lambda metrics, category_labels: [
+            sum(
+                count
+                for combination, count in metrics.trigger_combinations.items()
+                if _is_mixed_type_combination(combination)
+                and _format_mixed_combination_label(combination) == label
+            )
+            for label in category_labels
+        ],
+        normalization_fn=_fraction_with_poisson_errors,
+        output_path=output_path,
+        filename="mixed_trigger_combinations.png",
+        y_label="Fraction of Mixed-Type Triggers",
+        title="Mixed-Type Trigger Combinations",
+        figure_width=max(12, len(mixed_labels) * 0.8),
+        x_rotation=45,
+        x_ha="right",
+    )
+
+
+def _fractions_per_triggered_events(counts, metrics):
+    """Normalize counts by triggered-event count and return Poisson errors."""
+    event_norm = metrics.triggered_event_count if metrics.triggered_event_count > 0 else 1
+    return counts / event_norm, np.sqrt(counts) / event_norm
+
+
+def _plot_grouped_fraction_bars(
+    metrics_per_production,
+    categories,
+    counts_getter,
+    normalization_fn,
+    output_path,
+    filename,
+    y_label,
+    title,
+    figure_width,
+    x_rotation,
+    x_ha,
+):
+    """Plot grouped bars with fractions and Poisson error bars for each production."""
+    x_values = np.arange(len(categories))
     width = 0.8 / max(1, len(metrics_per_production))
-    fig, ax = plt.subplots(figsize=(max(12, len(mixed_labels) * 0.8), 6))
+    fig, ax = plt.subplots(figsize=(figure_width, 6))
 
     for index, metrics in enumerate(metrics_per_production):
-        counts = np.array(
-            [
-                sum(
-                    count
-                    for combination, count in metrics.trigger_combinations.items()
-                    if _is_mixed_type_combination(combination)
-                    and _format_mixed_combination_label(combination) == label
-                )
-                for label in mixed_labels
-            ],
-            dtype=float,
-        )
-        fractions, errors = _fraction_with_poisson_errors(counts)
+        counts = np.asarray(counts_getter(metrics, categories), dtype=float)
+        fractions, errors = normalization_fn(counts, metrics)
         offset = (index - (len(metrics_per_production) - 1) / 2.0) * width
         ax.bar(
             x_values + offset,
@@ -281,25 +286,18 @@ def _plot_mixed_trigger_combinations(metrics_per_production, output_path):
         )
 
     ax.set_xticks(x_values)
-    ax.set_xticklabels(mixed_labels, rotation=45, ha="right")
-    ax.set_ylabel("Fraction of Mixed-Type Triggers")
-    ax.set_title("Mixed-Type Trigger Combinations")
+    ax.set_xticklabels(categories, rotation=x_rotation, ha=x_ha)
+    ax.set_ylabel(y_label)
+    ax.set_title(title)
     ax.grid(axis="y", alpha=0.25)
     ax.legend()
 
-    _save_figure(fig, output_path, "mixed_trigger_combinations.png")
+    _save_figure(fig, output_path, filename)
 
 
 def _is_mixed_type_combination(combination):
     """Return True only for mixed signatures 1+1 and 1+2 (max multiplicity 3)."""
-    telescope_names = [name for name in combination.split(",") if name]
-    type_counts = {}
-    for telescope_name in telescope_names:
-        try:
-            tel_type = names.get_array_element_type_from_name(telescope_name)
-        except ValueError:
-            continue
-        type_counts[tel_type] = type_counts.get(tel_type, 0) + 1
+    type_counts, _ = _type_counts_from_combination(combination)
     if len(type_counts) < 2:
         return False
     signature = tuple(sorted(type_counts.values()))
@@ -308,6 +306,14 @@ def _is_mixed_type_combination(combination):
 
 def _format_mixed_combination_label(combination):
     """Format mixed trigger combination as '<signature> | <tel1> + <tel2> + ...'."""
+    type_counts, telescope_names = _type_counts_from_combination(combination)
+    signature = "+".join(str(count) for _, count in sorted(type_counts.items()))
+    telescopes_label = " + ".join(telescope_names)
+    return f"{signature} | {telescopes_label}"
+
+
+def _type_counts_from_combination(combination):
+    """Return per-type multiplicities and original telescope names for a combination."""
     telescope_names = [name for name in combination.split(",") if name]
     type_counts = {}
     for telescope_name in telescope_names:
@@ -316,9 +322,7 @@ def _format_mixed_combination_label(combination):
         except ValueError:
             continue
         type_counts[tel_type] = type_counts.get(tel_type, 0) + 1
-    signature = "+".join(str(count) for _, count in sorted(type_counts.items()))
-    telescopes_label = " + ".join(telescope_names)
-    return f"{signature} | {telescopes_label}"
+    return type_counts, telescope_names
 
 
 def _plot_triggered_vs_quantity(
@@ -497,36 +501,21 @@ def _plot_telescope_participation(metrics_per_production, output_path):
         _logger.warning("Skipping telescope participation plot, no triggered events available.")
         return
 
-    x_values = np.arange(len(telescopes))
-    width = 0.8 / max(1, len(metrics_per_production))
-    fig, ax = plt.subplots(figsize=(max(10, len(telescopes) * 0.4), 6))
-
-    for index, metrics in enumerate(metrics_per_production):
-        event_norm = metrics.triggered_event_count if metrics.triggered_event_count > 0 else 1
-        counts = np.array(
-            [metrics.telescope_participation.get(telescope, 0) for telescope in telescopes],
-            dtype=float,
-        )
-        fractions = counts / event_norm
-        errors = np.sqrt(counts) / event_norm
-        offset = (index - (len(metrics_per_production) - 1) / 2.0) * width
-        ax.bar(
-            x_values + offset,
-            fractions,
-            width=width,
-            label=metrics.label,
-            yerr=errors,
-            error_kw={"elinewidth": 0.7, "capsize": 1, "capthick": 0.7},
-        )
-
-    ax.set_xticks(x_values)
-    ax.set_xticklabels(telescopes, rotation=90)
-    ax.set_ylabel("Participation Fraction")
-    ax.set_title("Telescope Participation in Triggered Events")
-    ax.grid(axis="y", alpha=0.25)
-    ax.legend()
-
-    _save_figure(fig, output_path, "telescope_participation_fraction.png")
+    _plot_grouped_fraction_bars(
+        metrics_per_production,
+        categories=telescopes,
+        counts_getter=lambda metrics, category_labels: [
+            metrics.telescope_participation.get(telescope, 0) for telescope in category_labels
+        ],
+        normalization_fn=_fractions_per_triggered_events,
+        output_path=output_path,
+        filename="telescope_participation_fraction.png",
+        y_label="Participation Fraction",
+        title="Telescope Participation in Triggered Events",
+        figure_width=max(10, len(telescopes) * 0.4),
+        x_rotation=90,
+        x_ha="center",
+    )
 
 
 def _distribution_cumulative_variants(quantity_name):
@@ -596,7 +585,7 @@ def _plot_distribution_series(
         _plot_histogram_error_bars(ax, bin_edges, values, errors, color=_artist_color(artist))
 
 
-def _fraction_with_poisson_errors(counts):
+def _fraction_with_poisson_errors(counts, _metrics=None):
     """Return normalized bin fractions and Poisson errors for counts."""
     counts = np.asarray(counts, dtype=float)
     total = counts.sum()

From 7ad1a01132a6dfb1c37c33c792c6a8d1357c27d2 Mon Sep 17 00:00:00 2001
From: Gernot Maier <gernot.maier@desy.de>
Date: Wed, 13 May 2026 10:57:21 +0200
Subject: [PATCH 08/27] missed application

---
 pyproject.toml | 1 -
 1 file changed, 1 deletion(-)

diff --git a/pyproject.toml b/pyproject.toml
index 3846752a64..d253530718 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -74,7 +74,6 @@ urls."bug tracker" = "https://github.com/gammasim/simtools/issues"
 urls.documentation = "https://gammasim.github.io/simtools/"
 urls.repository = "https://github.com/gammasim/simtools"
 scripts.simtools-compare-productions = "simtools.applications.compare_productions:main"
-scripts.simtools-compare-productions-on-event-level = "simtools.applications.compare_productions:main"
 scripts.simtools-convert-all-model-parameters-from-simtel = "simtools.applications.convert_all_model_parameters_from_simtel:main"
 scripts.simtools-convert-geo-coordinates-of-array-elements = "simtools.applications.convert_geo_coordinates_of_array_elements:main"
 scripts.simtools-convert-model-parameter-from-simtel = "simtools.applications.convert_model_parameter_from_simtel:main"

From 8d6a0d2a776f141cc5cc18f048901ef0f0a5e72d Mon Sep 17 00:00:00 2001
From: Gernot Maier <gernot.maier@desy.de>
Date: Wed, 13 May 2026 11:02:39 +0200
Subject: [PATCH 09/27] simplify

---
 .../applications/compare_productions.py       | 49 ++++---------------
 ...test_compare_productions_on_event_level.py |  3 +-
 2 files changed, 10 insertions(+), 42 deletions(-)

diff --git a/src/simtools/applications/compare_productions.py b/src/simtools/applications/compare_productions.py
index 7506ef338a..2ffb6eacb4 100644
--- a/src/simtools/applications/compare_productions.py
+++ b/src/simtools/applications/compare_productions.py
@@ -49,42 +49,6 @@ def _add_arguments(parser):
         default="events",
         help="Comparison level to execute.",
     )
-    parser.add_argument(
-        "--telescope_ids",
-        nargs="+",
-        default=None,
-        help="Optional telescope IDs to filter triggers.",
-    )
-
-
-def _run_event_comparison(app_context):
-    """Run event-level production comparison."""
-    production_descriptors = parse_production_arguments(app_context.args["production"])
-    metrics_per_production = collect_production_metrics(
-        production_descriptors,
-        telescope_list=app_context.args["telescope_ids"],
-    )
-
-    output_directory = app_context.io_handler.get_output_directory()
-    plot_event_level_production_comparison.plot(
-        metrics_per_production,
-        output_path=output_directory,
-    )
-
-
-def _run_placeholder_comparison(comparison_level):
-    """Warn for future comparison-level implementations.
-
-    Notes
-    -----
-    Future work tracked in:
-    - signals level: issue #2183
-    - compute level: issue #2184
-    """
-    raise NotImplementedError(
-        f"Comparison level '{comparison_level}' is not implemented yet "
-        "(planned via issues #2183 and #2184)."
-    )
 
 
 def main():
@@ -93,12 +57,17 @@ def main():
         initialization_kwargs={"db_config": False, "output": True},
     )
 
+    output_directory = app_context.io_handler.get_output_directory()
+    production_descriptors = parse_production_arguments(app_context.args["production"])
+
     comparison_level = app_context.args["comparison_level"]
     if comparison_level == "events":
-        _run_event_comparison(app_context)
-        return
-
-    _run_placeholder_comparison(comparison_level)
+        plot_event_level_production_comparison.plot(
+            collect_production_metrics(production_descriptors),
+            output_path=output_directory,
+        )
+    else:
+        raise NotImplementedError(f"Comparison level '{comparison_level}' is not implemented yet.")
 
 
 if __name__ == "__main__":
diff --git a/tests/unit_tests/applications/test_compare_productions_on_event_level.py b/tests/unit_tests/applications/test_compare_productions_on_event_level.py
index af2236a546..016dc4d8f7 100644
--- a/tests/unit_tests/applications/test_compare_productions_on_event_level.py
+++ b/tests/unit_tests/applications/test_compare_productions_on_event_level.py
@@ -12,7 +12,6 @@ def test_main_collects_metrics_and_plots(tmp_test_directory):
         args={
             "production": [["baseline", "base_*.h5"], ["candidate", "cand_*.h5"]],
             "comparison_level": "events",
-            "telescope_ids": ["LSTN-01"],
         },
         io_handler=MagicMock(),
     )
@@ -41,5 +40,5 @@ def test_main_collects_metrics_and_plots(tmp_test_directory):
         app.main()
 
     mock_parse.assert_called_once_with(app_context.args["production"])
-    mock_collect.assert_called_once_with(parsed_productions, telescope_list=["LSTN-01"])
+    mock_collect.assert_called_once_with(parsed_productions)
     mock_plot.assert_called_once_with(collected_metrics, output_path=output_dir)

From 176a0b99914b2a099473a7d3b84eb514c5c67c83 Mon Sep 17 00:00:00 2001
From: Gernot Maier <gernot.maier@desy.de>
Date: Wed, 13 May 2026 11:08:43 +0200
Subject: [PATCH 10/27] sonar

---
 src/simtools/sim_events/production_comparison.py      | 11 ++++++++---
 .../plot_event_level_production_comparison.py         |  2 +-
 .../sim_events/test_production_comparison.py          |  4 ++--
 3 files changed, 11 insertions(+), 6 deletions(-)

diff --git a/src/simtools/sim_events/production_comparison.py b/src/simtools/sim_events/production_comparison.py
index 5d6f5d1f63..4d72f95f31 100644
--- a/src/simtools/sim_events/production_comparison.py
+++ b/src/simtools/sim_events/production_comparison.py
@@ -104,7 +104,7 @@ def _normalize_production_arguments(production_arguments):
 def _pairwise_label_file_arguments(flat_arguments):
     """Convert a flat list of strings into ``[(label, files), ...]`` pairs."""
     if len(flat_arguments) % 2 != 0:
-        raise ValueError("Production arguments must be provided as label/file pairs.")
+        _raise_invalid_production_arguments()
     return [
         (flat_arguments[index], flat_arguments[index + 1])
         for index in range(0, len(flat_arguments), 2)
@@ -114,14 +114,19 @@ def _pairwise_label_file_arguments(flat_arguments):
 def _normalize_single_production_argument(argument):
     """Normalize one nested production argument into label/file pairs."""
     if not isinstance(argument, list | tuple):
-        raise ValueError("Production arguments must be provided as label/file pairs.")
+        _raise_invalid_production_arguments()
     if not all(isinstance(value, str) for value in argument):
-        raise ValueError("Production arguments must be provided as label/file pairs.")
+        _raise_invalid_production_arguments()
     if len(argument) == 2:
         return [(argument[0], argument[1])]
     return _pairwise_label_file_arguments(list(argument))
 
 
+def _raise_invalid_production_arguments():
+    """Raise a standardized parser error for malformed production arguments."""
+    raise ValueError("Production arguments must be provided as label/file pairs.")
+
+
 def collect_production_metrics(production_descriptors, telescope_list=None):
     """Collect event-level comparison metrics for each production.
 
diff --git a/src/simtools/visualization/plot_event_level_production_comparison.py b/src/simtools/visualization/plot_event_level_production_comparison.py
index 466640ab31..8b5f33037c 100644
--- a/src/simtools/visualization/plot_event_level_production_comparison.py
+++ b/src/simtools/visualization/plot_event_level_production_comparison.py
@@ -101,7 +101,7 @@ def _plot_trigger_fraction(metrics_per_production, output_path, suffix=""):
     _save_figure(fig, output_path, f"trigger_fraction{suffix}.png")
 
 
-def _plot_trigger_multiplicity(metrics_per_production, output_path, suffix="", bins=40):
+def _plot_trigger_multiplicity(metrics_per_production, output_path, suffix="", bins=None):
     """Plot triggered telescope multiplicity distributions."""
     fig, ax = plt.subplots(figsize=(9, 6))
 
diff --git a/tests/unit_tests/sim_events/test_production_comparison.py b/tests/unit_tests/sim_events/test_production_comparison.py
index db172b0446..129a66781e 100644
--- a/tests/unit_tests/sim_events/test_production_comparison.py
+++ b/tests/unit_tests/sim_events/test_production_comparison.py
@@ -172,7 +172,7 @@ def test_collect_production_metrics_handles_empty_and_partial_datasets(mocker):
 
     assert metrics[0].simulated_event_count == 2
     assert metrics[0].triggered_event_count == 0
-    assert metrics[0].trigger_fraction == 0.0
+    assert metrics[0].trigger_fraction == pytest.approx(0.0)
     assert metrics[0].trigger_combinations == {}
     np.testing.assert_array_equal(metrics[0].triggered_energies, np.array([]))
 
@@ -189,7 +189,7 @@ def test_collect_production_metrics_empty_input_keeps_arrays_empty(mocker):
 
     assert metrics[0].simulated_event_count == 0
     assert metrics[0].triggered_event_count == 0
-    assert metrics[0].trigger_fraction == 0.0
+    assert metrics[0].trigger_fraction == pytest.approx(0.0)
     np.testing.assert_array_equal(metrics[0].simulated_energies, np.array([]))
     np.testing.assert_array_equal(metrics[0].trigger_multiplicity, np.array([], dtype=int))
 

From b050ae71126d83a2341c20d732871a7d7c0eff79 Mon Sep 17 00:00:00 2001
From: Gernot Maier <gernot.maier@desy.de>
Date: Wed, 13 May 2026 11:14:39 +0200
Subject: [PATCH 11/27] simplification

---
 src/simtools/sim_events/production_comparison.py | 10 +++-------
 1 file changed, 3 insertions(+), 7 deletions(-)

diff --git a/src/simtools/sim_events/production_comparison.py b/src/simtools/sim_events/production_comparison.py
index 4d72f95f31..7678c31ced 100644
--- a/src/simtools/sim_events/production_comparison.py
+++ b/src/simtools/sim_events/production_comparison.py
@@ -13,7 +13,7 @@
 
 @dataclass
 class ProductionDescriptor:
-    """Descriptor for one production input provided by the user."""
+    """Descriptor for one production input provided."""
 
     label: str
     event_data_files: list[str]
@@ -64,7 +64,7 @@ def parse_production_arguments(production_arguments):
         If configuration is malformed or does not contain any production.
     """
     parsed_productions = _normalize_production_arguments(production_arguments)
-    if len(parsed_productions) == 0:
+    if not parsed_productions:
         raise ValueError("At least one production is required.")
 
     labels = [item[0] for item in parsed_productions]
@@ -80,7 +80,6 @@ def parse_production_arguments(production_arguments):
         resolved_files = [str(path) for path in resolve_file_patterns(patterns)]
         if len(resolved_files) == 0:
             raise ValueError(f"Production '{label}' does not resolve to any files.")
-
         descriptors.append(ProductionDescriptor(label=label, event_data_files=resolved_files))
 
     return descriptors
@@ -327,10 +326,7 @@ def _subset_counts_for_event(telescopes):
     """Return subset keys and multiplicities for one triggered event."""
     type_counts = Counter()
     for telescope in telescopes:
-        try:
-            tel_type = names.get_array_element_type_from_name(telescope)
-        except ValueError:
-            continue
+        tel_type = names.get_array_element_type_from_name(telescope)
         type_counts[tel_type] += 1
 
     subset_counts = dict(type_counts)

From e453194bde17349447a726c8964ef7be6ae15fa4 Mon Sep 17 00:00:00 2001
From: Gernot Maier <gernot.maier@desy.de>
Date: Wed, 13 May 2026 11:30:37 +0200
Subject: [PATCH 12/27] sonar and sphinx, my best friends

---
 docs/source/api-reference/sim_events.md       |  9 ++++++
 docs/source/api-reference/visualization.md    |  7 +++++
 docs/source/conf.py                           |  1 +
 .../plot_event_level_production_comparison.py |  7 +++--
 .../sim_events/test_production_comparison.py  | 31 -------------------
 ..._plot_event_level_production_comparison.py |  6 +++-
 6 files changed, 27 insertions(+), 34 deletions(-)

diff --git a/docs/source/api-reference/sim_events.md b/docs/source/api-reference/sim_events.md
index e046313b85..566c222578 100644
--- a/docs/source/api-reference/sim_events.md
+++ b/docs/source/api-reference/sim_events.md
@@ -49,3 +49,12 @@ Provide tools for reading, writing, and processing simulated shower and telescop
 .. automodule:: simtools.sim_events.output_validator
    :members:
 ```
+
+## Production comparison
+
+(sim_events_production_comparison-1)=
+
+```{eval-rst}
+.. automodule:: simtools.sim_events.production_comparison
+   :members:
+```
diff --git a/docs/source/api-reference/visualization.md b/docs/source/api-reference/visualization.md
index b2b36eea11..45caf71a98 100644
--- a/docs/source/api-reference/visualization.md
+++ b/docs/source/api-reference/visualization.md
@@ -114,3 +114,10 @@ the visualization module.
 .. automodule:: visualization.camera_plot_utils
    :members:
 ```
+
+## plot_event_level_production_comparison
+
+```{eval-rst}
+.. automodule:: visualization.plot_event_level_production_comparison
+   :members:
+```
diff --git a/docs/source/conf.py b/docs/source/conf.py
index 38f5e03393..00dd594748 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -85,6 +85,7 @@ def get_python_version_from_pyproject():
 nitpick_ignore = [
     ("py:class", "astropy.table.table.Table"),
     ("py:class", "astropy.units.quantity.Quantity"),
+    ("py:class", "collections.Counter"),
     ("py:class", "logging.Logger"),
     ("py:class", "numpy.float64"),
     ("py:class", "numpy.ndarray"),
diff --git a/src/simtools/visualization/plot_event_level_production_comparison.py b/src/simtools/visualization/plot_event_level_production_comparison.py
index 8b5f33037c..e5b2a9f1d4 100644
--- a/src/simtools/visualization/plot_event_level_production_comparison.py
+++ b/src/simtools/visualization/plot_event_level_production_comparison.py
@@ -40,6 +40,8 @@ def plot(metrics_per_production, output_path, bins=40):
                 bins=bins,
             )
         for cumulative in _distribution_cumulative_variants(quantity_name):
+            if cumulative is None:
+                continue
             _plot_quantity_distribution(
                 metrics_per_production,
                 output_path,
@@ -519,9 +521,9 @@ def _plot_telescope_participation(metrics_per_production, output_path):
 
 
 def _distribution_cumulative_variants(quantity_name):
-    """Return enabled cumulative variants for one quantity."""
+    """Return two cumulative slots, using None to disable one for specific quantities."""
     if quantity_name == "energy":
-        return (False,)
+        return (False, None)
     return (False, True)
 
 
@@ -648,5 +650,6 @@ def _plot_histogram_error_bars(ax, bin_edges, values, errors, color):
         )
         for q, lbl, sc in _QUANTITY_CONFIGS
         for cum in _distribution_cumulative_variants(q)
+        if cum is not None
     ],
 ]
diff --git a/tests/unit_tests/sim_events/test_production_comparison.py b/tests/unit_tests/sim_events/test_production_comparison.py
index 129a66781e..0605ee162d 100644
--- a/tests/unit_tests/sim_events/test_production_comparison.py
+++ b/tests/unit_tests/sim_events/test_production_comparison.py
@@ -192,34 +192,3 @@ def test_collect_production_metrics_empty_input_keeps_arrays_empty(mocker):
     assert metrics[0].trigger_fraction == pytest.approx(0.0)
     np.testing.assert_array_equal(metrics[0].simulated_energies, np.array([]))
     np.testing.assert_array_equal(metrics[0].trigger_multiplicity, np.array([], dtype=int))
-
-
-def test_collect_production_metrics_ignores_unknown_telescope_type(mocker):
-    """Test unknown telescope names are ignored for type bins but still counted as single triggers."""
-    descriptor = SimpleNamespace(label="baseline", event_data_files=["baseline_file.h5"])
-    shower_data = SimpleNamespace(
-        simulated_energy=np.array([1.0]),
-        core_distance_shower=np.array([10.0]),
-        angular_distance=np.array([0.5]),
-    )
-    triggered_shower_data = SimpleNamespace(
-        simulated_energy=np.array([1.0]),
-        core_distance_shower=np.array([10.0]),
-        angular_distance=np.array([0.5]),
-    )
-    triggered_data = SimpleNamespace(telescope_list=[np.array(["BAD-01"])])
-
-    mock_reader = mocker.patch("simtools.sim_events.production_comparison.EventDataReader")
-    mock_reader.return_value.data_sets = [{"SHOWERS": "SHOWERS", "TRIGGERS": "TRIGGERS"}]
-    mock_reader.return_value.read_event_data.return_value = (
-        SimpleNamespace(),
-        shower_data,
-        triggered_shower_data,
-        triggered_data,
-    )
-
-    metrics = collect_production_metrics([descriptor])
-
-    assert "single_telescope" in metrics[0].per_type
-    assert metrics[0].per_type["single_telescope"].triggered_event_count == 1
-    assert "BAD" not in metrics[0].per_type
diff --git a/tests/unit_tests/visualization/test_plot_event_level_production_comparison.py b/tests/unit_tests/visualization/test_plot_event_level_production_comparison.py
index ef3b0058ff..2fc9597f60 100644
--- a/tests/unit_tests/visualization/test_plot_event_level_production_comparison.py
+++ b/tests/unit_tests/visualization/test_plot_event_level_production_comparison.py
@@ -105,7 +105,11 @@ def test_mixed_trigger_filter_accepts_only_allowed_signatures():
 
 @pytest.mark.parametrize(
     ("quantity_name", "expected"),
-    [("energy", (False,)), ("core_distance", (False, True)), ("angular_distance", (False, True))],
+    [
+        ("energy", (False, None)),
+        ("core_distance", (False, True)),
+        ("angular_distance", (False, True)),
+    ],
 )
 def test_distribution_cumulative_variants(quantity_name, expected):
     """Test cumulative variant selection by quantity."""

From 5f2753182ee389d322b4e3704067ea060e766fc7 Mon Sep 17 00:00:00 2001
From: Gernot Maier <gernot.maier@desy.de>
Date: Wed, 13 May 2026 11:32:30 +0200
Subject: [PATCH 13/27] obsolete test

---
 ...ductions_on_event_level_epos_urqmd_run.yml | 23 -------------------
 1 file changed, 23 deletions(-)
 delete mode 100644 tests/integration_tests/config/compare_productions_on_event_level_epos_urqmd_run.yml

diff --git a/tests/integration_tests/config/compare_productions_on_event_level_epos_urqmd_run.yml b/tests/integration_tests/config/compare_productions_on_event_level_epos_urqmd_run.yml
deleted file mode 100644
index 2215763f4e..0000000000
--- a/tests/integration_tests/config/compare_productions_on_event_level_epos_urqmd_run.yml
+++ /dev/null
@@ -1,23 +0,0 @@
----
-applications:
-- application: simtools-compare-productions
-  configuration:
-    comparison_level: events
-    production:
-    - epos
-    - <EPOS_EVENT_DATA_GLOB_PATTERN>
-    - urqmd
-    - <URQMD_EVENT_DATA_GLOB_PATTERN>
-    output_path: simtools-output
-  integration_tests:
-  - output_file: trigger_multiplicity.png
-  - output_file: trigger_combination.png
-  - output_file: single_telescope_trigger_distribution.png
-  - output_file: mixed_trigger_combinations.png
-  - output_file: distribution_energy.png
-  - output_file: distribution_core_distance.png
-  - output_file: distribution_angular_distance.png
-  - output_file: telescope_participation_fraction.png
-  test_name: run_epos_vs_urqmd
-schema_name: application_workflow.metaschema
-schema_version: 0.4.0

From 44d6b3bac93f221cd439f0446ed0bcb2ea32fe87 Mon Sep 17 00:00:00 2001
From: Gernot Maier <gernot.maier@desy.de>
Date: Wed, 13 May 2026 11:38:49 +0200
Subject: [PATCH 14/27] sonar

---
 .../plot_event_level_production_comparison.py         | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/src/simtools/visualization/plot_event_level_production_comparison.py b/src/simtools/visualization/plot_event_level_production_comparison.py
index e5b2a9f1d4..c231b5a099 100644
--- a/src/simtools/visualization/plot_event_level_production_comparison.py
+++ b/src/simtools/visualization/plot_event_level_production_comparison.py
@@ -105,6 +105,7 @@ def _plot_trigger_fraction(metrics_per_production, output_path, suffix=""):
 
 def _plot_trigger_multiplicity(metrics_per_production, output_path, suffix="", bins=None):
     """Plot triggered telescope multiplicity distributions."""
+    del bins
     fig, ax = plt.subplots(figsize=(9, 6))
 
     global_max = 0
@@ -118,16 +119,16 @@ def _plot_trigger_multiplicity(metrics_per_production, output_path, suffix="", b
         plt.close(fig)
         return
 
-    bins = np.arange(1, global_max + 2)
+    bin_edges = np.arange(1, global_max + 2)
     for metrics in metrics_per_production:
         if metrics.trigger_multiplicity.size == 0:
             continue
-        counts, _ = np.histogram(metrics.trigger_multiplicity, bins=bins)
+        counts, _ = np.histogram(metrics.trigger_multiplicity, bins=bin_edges)
         fractions, errors = _fraction_with_poisson_errors(counts)
-        stairs_artist = ax.stairs(fractions, bins, linewidth=1.5, label=metrics.label)
+        stairs_artist = ax.stairs(fractions, bin_edges, linewidth=1.5, label=metrics.label)
         _plot_histogram_error_bars(
             ax,
-            bins,
+            bin_edges,
             fractions,
             errors,
             color=_artist_color(stairs_artist),
@@ -137,7 +138,7 @@ def _plot_trigger_multiplicity(metrics_per_production, output_path, suffix="", b
     ax.set_ylabel("Fraction of Triggered Events")
     type_label = f" ({suffix.lstrip('_').replace('_', ' ')})" if suffix else ""
     ax.set_title(f"Trigger Multiplicity {type_label}")
-    ax.set_xticks(bins[:-1])
+    ax.set_xticks(bin_edges[:-1])
     ax.grid(alpha=0.25)
     ax.legend()
 

From 831c6f77c76cb3d9be4712ffab4b21bc03de4c81 Mon Sep 17 00:00:00 2001
From: Gernot Maier <gernot.maier@desy.de>
Date: Wed, 13 May 2026 12:05:09 +0200
Subject: [PATCH 15/27] unit test

---
 .../config/compare_productions_on_event_level_run.yml           | 2 --
 1 file changed, 2 deletions(-)

diff --git a/tests/integration_tests/config/compare_productions_on_event_level_run.yml b/tests/integration_tests/config/compare_productions_on_event_level_run.yml
index 2d202ea5f9..5d5e805862 100644
--- a/tests/integration_tests/config/compare_productions_on_event_level_run.yml
+++ b/tests/integration_tests/config/compare_productions_on_event_level_run.yml
@@ -12,8 +12,6 @@ applications:
   integration_tests:
   - output_file: trigger_multiplicity.png
   - output_file: trigger_combination.png
-  - output_file: single_telescope_trigger_distribution.png
-  - output_file: mixed_trigger_combinations.png
   - output_file: telescope_participation_fraction.png
   test_name: run
 schema_name: application_workflow.metaschema

From de395364013ca0a9c97160cc161de06f3aec5af4 Mon Sep 17 00:00:00 2001
From: Gernot Maier <gernot.maier@desy.de>
Date: Wed, 13 May 2026 12:27:40 +0200
Subject: [PATCH 16/27] copilot

---
 docs/source/user-guide/applications.md        |  2 +-
 .../plot_event_level_production_comparison.py | 25 -------------------
 ...el_run.yml => compare_productions_run.yml} |  0
 ...t_level.py => test_compare_productions.py} |  0
 ..._plot_event_level_production_comparison.py |  3 ---
 5 files changed, 1 insertion(+), 29 deletions(-)
 rename tests/integration_tests/config/{compare_productions_on_event_level_run.yml => compare_productions_run.yml} (100%)
 rename tests/unit_tests/applications/{test_compare_productions_on_event_level.py => test_compare_productions.py} (100%)

diff --git a/docs/source/user-guide/applications.md b/docs/source/user-guide/applications.md
index e875bf3973..b8d38cacba 100644
--- a/docs/source/user-guide/applications.md
+++ b/docs/source/user-guide/applications.md
@@ -49,10 +49,10 @@ Parameters with the same functionality are named consistently the same among all
 :glob: true
 :maxdepth: 1
 
+simtools-compare-productions <applications/simtools-compare-productions>
 simtools-convert-all-model-parameters-from-simtel <applications/simtools-convert-all-model-parameters-from-simtel>
 simtools-convert-geo-coordinates-of-array-elements <applications/simtools-convert-geo-coordinates-of-array-elements>
 simtools-convert-model-parameter-from-simtel <applications/simtools-convert-model-parameter-from-simtel>
-simtools-compare-productions <applications/simtools-compare-productions>
 simtools-db-add-file-to-db <applications/simtools-db-add-file-to-db>
 simtools-db-add-simulation-model-from-repository-to-db <applications/simtools-db-add-simulation-model-from-repository-to-db>
 simtools-db-add-value-from-json-to-db <applications/simtools-db-add-value-from-json-to-db>
diff --git a/src/simtools/visualization/plot_event_level_production_comparison.py b/src/simtools/visualization/plot_event_level_production_comparison.py
index c231b5a099..2a511f2786 100644
--- a/src/simtools/visualization/plot_event_level_production_comparison.py
+++ b/src/simtools/visualization/plot_event_level_production_comparison.py
@@ -78,31 +78,6 @@ def _save_figure(fig, output_path, filename):
     plt.close(fig)
 
 
-def _plot_trigger_fraction(metrics_per_production, output_path, suffix=""):
-    """Plot triggered/simulated event fractions by production."""
-    labels = [metrics.label for metrics in metrics_per_production]
-    values = [metrics.trigger_fraction for metrics in metrics_per_production]
-
-    fig, ax = plt.subplots(figsize=(9, 6))
-    bars = ax.bar(labels, values, color="tab:blue", alpha=0.85)
-    ax.set_ylabel("Trigger Fraction")
-    type_label = f" ({suffix.lstrip('_').replace('_', ' ')})" if suffix else ""
-    ax.set_title(f"Trigger Fraction {type_label}")
-    ax.set_ylim(0.0, max(1.0, max(values, default=0.0) * 1.15))
-    ax.grid(axis="y", alpha=0.25)
-    for bar_patch, value in zip(bars, values):
-        ax.text(
-            bar_patch.get_x() + bar_patch.get_width() / 2,
-            bar_patch.get_height(),
-            f"{value:.3f}",
-            ha="center",
-            va="bottom",
-            fontsize=9,
-        )
-
-    _save_figure(fig, output_path, f"trigger_fraction{suffix}.png")
-
-
 def _plot_trigger_multiplicity(metrics_per_production, output_path, suffix="", bins=None):
     """Plot triggered telescope multiplicity distributions."""
     del bins
diff --git a/tests/integration_tests/config/compare_productions_on_event_level_run.yml b/tests/integration_tests/config/compare_productions_run.yml
similarity index 100%
rename from tests/integration_tests/config/compare_productions_on_event_level_run.yml
rename to tests/integration_tests/config/compare_productions_run.yml
diff --git a/tests/unit_tests/applications/test_compare_productions_on_event_level.py b/tests/unit_tests/applications/test_compare_productions.py
similarity index 100%
rename from tests/unit_tests/applications/test_compare_productions_on_event_level.py
rename to tests/unit_tests/applications/test_compare_productions.py
diff --git a/tests/unit_tests/visualization/test_plot_event_level_production_comparison.py b/tests/unit_tests/visualization/test_plot_event_level_production_comparison.py
index 2fc9597f60..ce31046a8b 100644
--- a/tests/unit_tests/visualization/test_plot_event_level_production_comparison.py
+++ b/tests/unit_tests/visualization/test_plot_event_level_production_comparison.py
@@ -166,9 +166,6 @@ def test_plot_series_and_artist_color():
 def test_trigger_fraction_and_skip_paths(tmp_test_directory):
     """Test trigger fraction output and skip branches with empty inputs."""
     output_path = Path(tmp_test_directory)
-    metrics = [_build_metrics("baseline", simulated_scale=1.0, triggered_scale=1.0)]
-    plot_event_level_production_comparison._plot_trigger_fraction(metrics, output_path)
-    assert (output_path / "trigger_fraction.png").exists()
 
     empty_metric = _build_metrics("empty", simulated_scale=1.0, triggered_scale=1.0)
     empty_metric.trigger_multiplicity = np.array([], dtype=int)

From e2b034203111290e6038c5b0e30286513d8e9a3f Mon Sep 17 00:00:00 2001
From: Gernot Maier <gernot.maier@desy.de>
Date: Mon, 18 May 2026 09:43:29 +0200
Subject: [PATCH 17/27] report comparision statistics [skip ci]

---
 .../applications/compare_productions.py       |   3 +-
 .../sim_events/production_comparison.py       | 111 +++++++
 .../plot_event_level_production_comparison.py | 294 ++++++++++++++++--
 .../sim_events/test_production_comparison.py  |  32 ++
 ..._plot_event_level_production_comparison.py |   9 +
 5 files changed, 424 insertions(+), 25 deletions(-)

diff --git a/src/simtools/applications/compare_productions.py b/src/simtools/applications/compare_productions.py
index 2ffb6eacb4..df6517b627 100644
--- a/src/simtools/applications/compare_productions.py
+++ b/src/simtools/applications/compare_productions.py
@@ -62,8 +62,9 @@ def main():
 
     comparison_level = app_context.args["comparison_level"]
     if comparison_level == "events":
+        metrics_per_production = collect_production_metrics(production_descriptors)
         plot_event_level_production_comparison.plot(
-            collect_production_metrics(production_descriptors),
+            metrics_per_production,
             output_path=output_directory,
         )
     else:
diff --git a/src/simtools/sim_events/production_comparison.py b/src/simtools/sim_events/production_comparison.py
index 7678c31ced..e02e733b0e 100644
--- a/src/simtools/sim_events/production_comparison.py
+++ b/src/simtools/sim_events/production_comparison.py
@@ -5,6 +5,7 @@
 from dataclasses import dataclass, field
 
 import numpy as np
+from scipy import stats
 
 from simtools.sim_events.reader import EventDataReader
 from simtools.utils import names
@@ -36,6 +37,7 @@ class ProductionEventMetrics:
     simulated_angular_distances: np.ndarray = field(default_factory=lambda: np.array([]))
     triggered_angular_distances: np.ndarray = field(default_factory=lambda: np.array([]))
     per_type: dict = field(default_factory=dict)
+    comparison_statistics: dict = field(default_factory=dict)
 
     @property
     def trigger_fraction(self):
@@ -367,3 +369,112 @@ def _safe_concat(arrays, dtype=float):
     if len(arrays) == 0:
         return np.array([], dtype=dtype)
     return np.concatenate(arrays)
+
+
+def compute_ks_chi2_for_samples(baseline_samples, candidate_samples, bin_edges):
+    """Compute KS and Chi2 statistics for two sample arrays.
+
+    Parameters
+    ----------
+    baseline_samples : np.ndarray
+        Baseline sample values.
+    candidate_samples : np.ndarray
+        Candidate sample values.
+    bin_edges : np.ndarray
+        Histogram bin edges used to derive Chi2 counts.
+
+    Returns
+    -------
+    dict
+        KS and Chi2 statistics including p-values and counts metadata.
+    """
+    baseline_samples = np.asarray(baseline_samples)
+    candidate_samples = np.asarray(candidate_samples)
+    result = {
+        "ks_statistic": None,
+        "ks_pvalue": None,
+        "chi2_statistic": None,
+        "chi2_pvalue": None,
+        "valid": False,
+        "reason": "insufficient_data",
+    }
+
+    if baseline_samples.size == 0 or candidate_samples.size == 0:
+        return result
+
+    ks_test = stats.ks_2samp(baseline_samples, candidate_samples)
+    result["ks_statistic"] = float(ks_test.statistic)
+    result["ks_pvalue"] = float(ks_test.pvalue)
+
+    baseline_counts, _ = np.histogram(baseline_samples, bins=bin_edges)
+    candidate_counts, _ = np.histogram(candidate_samples, bins=bin_edges)
+    chi2_result = compute_ks_chi2_for_aligned_histograms(baseline_counts, candidate_counts)
+    result["chi2_statistic"] = chi2_result["chi2_statistic"]
+    result["chi2_pvalue"] = chi2_result["chi2_pvalue"]
+    result["valid"] = chi2_result["valid"]
+    result["reason"] = chi2_result["reason"]
+    result["baseline_counts"] = baseline_counts.astype(int).tolist()
+    result["candidate_counts"] = candidate_counts.astype(int).tolist()
+    result["bin_edges"] = np.asarray(bin_edges, dtype=float).tolist()
+    return result
+
+
+def compute_ks_chi2_for_aligned_histograms(baseline_counts, candidate_counts):
+    """Compute KS and Chi2 statistics from aligned count arrays.
+
+    Parameters
+    ----------
+    baseline_counts : np.ndarray
+        Baseline histogram counts on aligned support.
+    candidate_counts : np.ndarray
+        Candidate histogram counts on aligned support.
+
+    Returns
+    -------
+    dict
+        KS and Chi2 statistics including p-values and validity metadata.
+    """
+    baseline_counts = np.asarray(baseline_counts, dtype=float)
+    candidate_counts = np.asarray(candidate_counts, dtype=float)
+
+    result = {
+        "ks_statistic": None,
+        "ks_pvalue": None,
+        "chi2_statistic": None,
+        "chi2_pvalue": None,
+        "valid": False,
+        "reason": "insufficient_data",
+    }
+
+    baseline_total = float(np.sum(baseline_counts))
+    candidate_total = float(np.sum(candidate_counts))
+    if baseline_total <= 0 or candidate_total <= 0:
+        return result
+
+    baseline_cdf = np.cumsum(baseline_counts) / baseline_total
+    candidate_cdf = np.cumsum(candidate_counts) / candidate_total
+    result["ks_statistic"] = float(np.max(np.abs(baseline_cdf - candidate_cdf)))
+
+    support_mask = baseline_counts > 0
+    if not np.any(support_mask):
+        result["reason"] = "zero_baseline_counts"
+        return result
+
+    observed = candidate_counts[support_mask]
+    observed_total = float(np.sum(observed))
+    baseline_supported_total = float(np.sum(baseline_counts[support_mask]))
+    if observed_total <= 0 or baseline_supported_total <= 0:
+        result["reason"] = "invalid_expected_counts"
+        return result
+
+    expected = baseline_counts[support_mask] * (observed_total / baseline_supported_total)
+    if np.sum(expected) <= 0:
+        result["reason"] = "invalid_expected_counts"
+        return result
+
+    chi2 = stats.chisquare(f_obs=observed, f_exp=expected)
+    result["chi2_statistic"] = float(chi2.statistic)
+    result["chi2_pvalue"] = float(chi2.pvalue)
+    result["valid"] = True
+    result["reason"] = "ok"
+    return result
diff --git a/src/simtools/visualization/plot_event_level_production_comparison.py b/src/simtools/visualization/plot_event_level_production_comparison.py
index 2a511f2786..ab4504603c 100644
--- a/src/simtools/visualization/plot_event_level_production_comparison.py
+++ b/src/simtools/visualization/plot_event_level_production_comparison.py
@@ -1,11 +1,16 @@
 """Plot event-level comparisons across multiple simulation productions."""
 
 import functools
+import json
 import logging
 
 import matplotlib.pyplot as plt
 import numpy as np
 
+from simtools.sim_events.production_comparison import (
+    compute_ks_chi2_for_aligned_histograms,
+    compute_ks_chi2_for_samples,
+)
 from simtools.utils import names
 
 _logger = logging.getLogger(__name__)
@@ -23,11 +28,33 @@ def plot(metrics_per_production, output_path, bins=40):
     bins : int, optional
         Number of bins for 1D histograms.
     """
-    _plot_trigger_multiplicity(metrics_per_production, output_path)
-    _plot_trigger_combinations(metrics_per_production, output_path)
-    _plot_single_telescope_trigger_frequencies(metrics_per_production, output_path)
-    _plot_mixed_trigger_combinations(metrics_per_production, output_path)
-    _plot_telescope_participation(metrics_per_production, output_path)
+    comparison_statistics = _initialize_comparison_statistics(metrics_per_production)
+
+    _record_plot_statistics(
+        comparison_statistics,
+        "trigger_multiplicity",
+        _plot_trigger_multiplicity(metrics_per_production, output_path),
+    )
+    _record_plot_statistics(
+        comparison_statistics,
+        "trigger_combination",
+        _plot_trigger_combinations(metrics_per_production, output_path),
+    )
+    _record_plot_statistics(
+        comparison_statistics,
+        "single_telescope_trigger_distribution",
+        _plot_single_telescope_trigger_frequencies(metrics_per_production, output_path),
+    )
+    _record_plot_statistics(
+        comparison_statistics,
+        "mixed_trigger_combinations",
+        _plot_mixed_trigger_combinations(metrics_per_production, output_path),
+    )
+    _record_plot_statistics(
+        comparison_statistics,
+        "telescope_participation_fraction",
+        _plot_telescope_participation(metrics_per_production, output_path),
+    )
 
     for quantity_name, x_label, x_scale in _QUANTITY_CONFIGS:
         if quantity_name in _TRIGGERED_FRACTION_QUANTITIES:
@@ -42,14 +69,19 @@ def plot(metrics_per_production, output_path, bins=40):
         for cumulative in _distribution_cumulative_variants(quantity_name):
             if cumulative is None:
                 continue
-            _plot_quantity_distribution(
-                metrics_per_production,
-                output_path,
-                quantity_name=quantity_name,
-                x_label=x_label,
-                x_scale=x_scale,
-                bins=bins,
-                cumulative=cumulative,
+            cumulative_suffix = "_cumulative" if cumulative else ""
+            _record_plot_statistics(
+                comparison_statistics,
+                f"distribution_{quantity_name}{cumulative_suffix}",
+                _plot_quantity_distribution(
+                    metrics_per_production,
+                    output_path,
+                    quantity_name=quantity_name,
+                    x_label=x_label,
+                    x_scale=x_scale,
+                    bins=bins,
+                    cumulative=cumulative,
+                ),
             )
 
     all_types = sorted(
@@ -67,7 +99,13 @@ def plot(metrics_per_production, output_path, bins=40):
             if tel_type in metrics.per_type
         ]
         for plot_fn in _PER_TYPE_PLOT_FNS:
-            plot_fn(type_metrics, output_path, suffix=f"_{tel_type}", bins=bins)
+            statistics = plot_fn(type_metrics, output_path, suffix=f"_{tel_type}", bins=bins)
+            if statistics is None:
+                continue
+            plot_name = _plot_name_from_statistics(statistics)
+            _record_plot_statistics(comparison_statistics, plot_name, statistics)
+
+    _write_comparison_statistics_json(comparison_statistics, output_path)
 
 
 def _save_figure(fig, output_path, filename):
@@ -92,13 +130,15 @@ def _plot_trigger_multiplicity(metrics_per_production, output_path, suffix="", b
     if global_max == 0:
         _logger.warning("Skipping trigger multiplicity plot, no triggered events available.")
         plt.close(fig)
-        return
+        return None
 
     bin_edges = np.arange(1, global_max + 2)
+    counts_per_label = {}
     for metrics in metrics_per_production:
         if metrics.trigger_multiplicity.size == 0:
             continue
         counts, _ = np.histogram(metrics.trigger_multiplicity, bins=bin_edges)
+        counts_per_label[metrics.label] = counts
         fractions, errors = _fraction_with_poisson_errors(counts)
         stairs_artist = ax.stairs(fractions, bin_edges, linewidth=1.5, label=metrics.label)
         _plot_histogram_error_bars(
@@ -117,7 +157,16 @@ def _plot_trigger_multiplicity(metrics_per_production, output_path, suffix="", b
     ax.grid(alpha=0.25)
     ax.legend()
 
+    statistics = _build_aligned_count_statistics(
+        metrics_per_production,
+        counts_per_label,
+        metadata={"bin_edges": bin_edges.astype(float).tolist()},
+    )
+    _annotate_ks_statistics(ax, statistics)
+
     _save_figure(fig, output_path, f"trigger_multiplicity{suffix}.png")
+    statistics["plot_name"] = f"trigger_multiplicity{suffix}"
+    return statistics
 
 
 def _plot_trigger_combinations(metrics_per_production, output_path, top_n=12):
@@ -128,7 +177,7 @@ def _plot_trigger_combinations(metrics_per_production, output_path, top_n=12):
 
     if len(combinations) == 0:
         _logger.warning("Skipping trigger combination plot, no combinations available.")
-        return
+        return None
 
     totals = {}
     for combination in combinations:
@@ -139,7 +188,7 @@ def _plot_trigger_combinations(metrics_per_production, output_path, top_n=12):
     selected = [name for name, _ in sorted(totals.items(), key=lambda item: item[1], reverse=True)]
     selected = selected[:top_n]
 
-    _plot_grouped_fraction_bars(
+    statistics = _plot_grouped_fraction_bars(
         metrics_per_production,
         categories=selected,
         counts_getter=lambda metrics, category_labels: [
@@ -154,6 +203,9 @@ def _plot_trigger_combinations(metrics_per_production, output_path, top_n=12):
         x_rotation=45,
         x_ha="right",
     )
+    if statistics is not None:
+        statistics["plot_name"] = "trigger_combination"
+    return statistics
 
 
 def _plot_single_telescope_trigger_frequencies(metrics_per_production, output_path):
@@ -169,9 +221,9 @@ def _plot_single_telescope_trigger_frequencies(metrics_per_production, output_pa
 
     if len(telescope_names) == 0:
         _logger.warning("Skipping single-telescope trigger frequency plot, no data available.")
-        return
+        return None
 
-    _plot_grouped_fraction_bars(
+    statistics = _plot_grouped_fraction_bars(
         metrics_per_production,
         categories=telescope_names,
         counts_getter=lambda metrics, category_labels: [
@@ -187,6 +239,9 @@ def _plot_single_telescope_trigger_frequencies(metrics_per_production, output_pa
         x_rotation=90,
         x_ha="center",
     )
+    if statistics is not None:
+        statistics["plot_name"] = "single_telescope_trigger_distribution"
+    return statistics
 
 
 def _plot_mixed_trigger_combinations(metrics_per_production, output_path):
@@ -201,9 +256,9 @@ def _plot_mixed_trigger_combinations(metrics_per_production, output_path):
     )
     if len(mixed_labels) == 0:
         _logger.warning("Skipping mixed trigger combination plot, no mixed-type combinations.")
-        return
+        return None
 
-    _plot_grouped_fraction_bars(
+    statistics = _plot_grouped_fraction_bars(
         metrics_per_production,
         categories=mixed_labels,
         counts_getter=lambda metrics, category_labels: [
@@ -224,6 +279,9 @@ def _plot_mixed_trigger_combinations(metrics_per_production, output_path):
         x_rotation=45,
         x_ha="right",
     )
+    if statistics is not None:
+        statistics["plot_name"] = "mixed_trigger_combinations"
+    return statistics
 
 
 def _fractions_per_triggered_events(counts, metrics):
@@ -249,9 +307,11 @@ def _plot_grouped_fraction_bars(
     x_values = np.arange(len(categories))
     width = 0.8 / max(1, len(metrics_per_production))
     fig, ax = plt.subplots(figsize=(figure_width, 6))
+    counts_per_label = {}
 
     for index, metrics in enumerate(metrics_per_production):
         counts = np.asarray(counts_getter(metrics, categories), dtype=float)
+        counts_per_label[metrics.label] = counts
         fractions, errors = normalization_fn(counts, metrics)
         offset = (index - (len(metrics_per_production) - 1) / 2.0) * width
         ax.bar(
@@ -270,7 +330,15 @@ def _plot_grouped_fraction_bars(
     ax.grid(axis="y", alpha=0.25)
     ax.legend()
 
+    statistics = _build_aligned_count_statistics(
+        metrics_per_production,
+        counts_per_label,
+        metadata={"categories": list(categories)},
+    )
+    _annotate_ks_statistics(ax, statistics)
+
     _save_figure(fig, output_path, filename)
+    return statistics
 
 
 def _is_mixed_type_combination(combination):
@@ -348,6 +416,7 @@ def _plot_triggered_vs_quantity(
     ax.legend()
 
     _save_figure(fig, output_path, f"triggered_fraction_vs_{quantity_name}{suffix}.png")
+    return
 
 
 def _get_bin_edges(values, x_scale, bins):
@@ -421,6 +490,8 @@ def _plot_quantity_distribution(
     """
     fig, ax = plt.subplots(figsize=(9, 6))
     plotted = False
+    bin_edges_by_label = {}
+    counts_by_label = {}
     for metrics in metrics_per_production:
         simulated, triggered = _get_quantity_arrays(metrics, quantity_name)
         if simulated.size == 0:
@@ -430,6 +501,13 @@ def _plot_quantity_distribution(
         bin_edges = _get_bin_edges(simulated, x_scale=x_scale, bins=bins)
         sim_counts, _ = np.histogram(simulated, bins=bin_edges)
         trig_counts, _ = np.histogram(triggered, bins=bin_edges)
+        bin_edges_by_label[metrics.label] = bin_edges
+        counts_by_label[metrics.label] = {
+            "simulated": sim_counts,
+            "triggered": trig_counts,
+            "simulated_samples": simulated,
+            "triggered_samples": triggered,
+        }
         _plot_distribution_series(
             ax,
             bin_edges,
@@ -453,7 +531,7 @@ def _plot_quantity_distribution(
     if not plotted:
         _logger.warning(f"Skipping distribution plot for {quantity_name}, no data available.")
         plt.close(fig)
-        return
+        return None
     cum_label = "Cumulative " if cumulative else ""
     type_label = f" ({suffix.lstrip('_').replace('_', ' ')})" if suffix else ""
     ax.set_xlabel(x_label)
@@ -462,8 +540,20 @@ def _plot_quantity_distribution(
     ax.set_xscale(x_scale)
     ax.grid(alpha=0.25)
     ax.legend(fontsize=8)
+
+    statistics = _build_quantity_distribution_statistics(
+        metrics_per_production,
+        quantity_name,
+        counts_by_label,
+        bin_edges_by_label,
+        cumulative=cumulative,
+    )
+    _annotate_ks_statistics(ax, statistics)
+
     cum_str = "_cumulative" if cumulative else ""
     _save_figure(fig, output_path, f"distribution_{quantity_name}{cum_str}{suffix}.png")
+    statistics["plot_name"] = f"distribution_{quantity_name}{cum_str}{suffix}"
+    return statistics
 
 
 def _plot_telescope_participation(metrics_per_production, output_path):
@@ -477,9 +567,9 @@ def _plot_telescope_participation(metrics_per_production, output_path):
     )
     if len(telescopes) == 0:
         _logger.warning("Skipping telescope participation plot, no triggered events available.")
-        return
+        return None
 
-    _plot_grouped_fraction_bars(
+    statistics = _plot_grouped_fraction_bars(
         metrics_per_production,
         categories=telescopes,
         counts_getter=lambda metrics, category_labels: [
@@ -494,6 +584,162 @@ def _plot_telescope_participation(metrics_per_production, output_path):
         x_rotation=90,
         x_ha="center",
     )
+    if statistics is not None:
+        statistics["plot_name"] = "telescope_participation_fraction"
+    return statistics
+
+
+def _build_aligned_count_statistics(metrics_per_production, counts_per_label, metadata):
+    """Build comparison statistics for plots represented by aligned count arrays."""
+    baseline_label = metrics_per_production[0].label if metrics_per_production else None
+    baseline_counts = counts_per_label.get(baseline_label)
+    stats_summary = {
+        "baseline_label": baseline_label,
+        "metric_type": "aligned_counts",
+        "metadata": metadata,
+        "comparisons": [],
+    }
+
+    if baseline_label is None or baseline_counts is None:
+        return stats_summary
+
+    for metrics in metrics_per_production[1:]:
+        candidate_counts = counts_per_label.get(metrics.label)
+        if candidate_counts is None:
+            continue
+        result = compute_ks_chi2_for_aligned_histograms(baseline_counts, candidate_counts)
+        result["candidate_label"] = metrics.label
+        result["baseline_counts"] = np.asarray(baseline_counts, dtype=int).tolist()
+        result["candidate_counts"] = np.asarray(candidate_counts, dtype=int).tolist()
+        stats_summary["comparisons"].append(result)
+    return stats_summary
+
+
+def _build_quantity_distribution_statistics(
+    metrics_per_production,
+    quantity_name,
+    counts_by_label,
+    bin_edges_by_label,
+    cumulative,
+):
+    """Build comparison statistics for quantity distribution plots."""
+    baseline_label = metrics_per_production[0].label if metrics_per_production else None
+    stats_summary = {
+        "baseline_label": baseline_label,
+        "metric_type": "quantity_distribution",
+        "quantity_name": quantity_name,
+        "cumulative": bool(cumulative),
+        "comparisons": [],
+    }
+    if baseline_label is None or baseline_label not in counts_by_label:
+        return stats_summary
+
+    baseline_data = counts_by_label[baseline_label]
+    baseline_edges = bin_edges_by_label[baseline_label]
+
+    for metrics in metrics_per_production[1:]:
+        if metrics.label not in counts_by_label:
+            continue
+        candidate_data = counts_by_label[metrics.label]
+
+        sim_stats = compute_ks_chi2_for_samples(
+            baseline_data["simulated_samples"],
+            candidate_data["simulated_samples"],
+            baseline_edges,
+        )
+        trig_stats = compute_ks_chi2_for_samples(
+            baseline_data["triggered_samples"],
+            candidate_data["triggered_samples"],
+            baseline_edges,
+        )
+
+        comparison_record = {
+            "candidate_label": metrics.label,
+            "simulated": sim_stats,
+            "triggered": trig_stats,
+            "baseline_bin_edges": np.asarray(baseline_edges, dtype=float).tolist(),
+        }
+        stats_summary["comparisons"].append(comparison_record)
+    return stats_summary
+
+
+def _annotate_ks_statistics(ax, statistics):
+    """Overlay KS values for baseline-vs-candidate comparisons on a plot."""
+    if statistics is None or not statistics.get("comparisons"):
+        return
+
+    metric_type = statistics.get("metric_type")
+    lines = [f"Baseline: {statistics.get('baseline_label', 'n/a')}"]
+    for comparison in statistics["comparisons"]:
+        label = comparison.get("candidate_label", "candidate")
+        if metric_type == "quantity_distribution":
+            sim_ks = comparison.get("simulated", {}).get("ks_statistic")
+            trig_ks = comparison.get("triggered", {}).get("ks_statistic")
+            lines.append(
+                f"{label}: KS(sim)={_format_statistic_value(sim_ks)}, "
+                f"KS(trg)={_format_statistic_value(trig_ks)}"
+            )
+            continue
+        lines.append(f"{label}: KS={_format_statistic_value(comparison.get('ks_statistic'))}")
+
+    ax.text(
+        0.02,
+        0.98,
+        "\n".join(lines),
+        transform=ax.transAxes,
+        va="top",
+        ha="left",
+        fontsize=7,
+        bbox={"boxstyle": "round", "facecolor": "white", "alpha": 0.75, "edgecolor": "0.8"},
+    )
+
+
+def _format_statistic_value(value):
+    """Format statistic values for plot text annotations."""
+    if value is None:
+        return "n/a"
+    return f"{float(value):.4f}"
+
+
+def _initialize_comparison_statistics(metrics_per_production):
+    """Initialize top-level comparison statistics structure."""
+    baseline = metrics_per_production[0] if metrics_per_production else None
+    return {
+        "baseline": {
+            "label": baseline.label if baseline else None,
+            "simulated_event_count": baseline.simulated_event_count if baseline else 0,
+            "triggered_event_count": baseline.triggered_event_count if baseline else 0,
+        },
+        "comparison_sets": [
+            {
+                "label": metrics.label,
+                "simulated_event_count": metrics.simulated_event_count,
+                "triggered_event_count": metrics.triggered_event_count,
+            }
+            for metrics in metrics_per_production[1:]
+        ],
+        "plot_statistics": {},
+    }
+
+
+def _record_plot_statistics(comparison_statistics, plot_name, statistics):
+    """Store plot-level statistics if available."""
+    if statistics is None:
+        return
+    comparison_statistics["plot_statistics"][plot_name] = statistics
+
+
+def _plot_name_from_statistics(statistics):
+    """Return plot name from statistics object or fallback placeholder."""
+    return statistics.get("plot_name", "unnamed_plot")
+
+
+def _write_comparison_statistics_json(comparison_statistics, output_path):
+    """Write comparison statistics JSON file in the output directory."""
+    output_file = output_path / "comparison_statistics.json"
+    _logger.info(f"Writing comparison statistics: {output_file}")
+    with output_file.open("w", encoding="utf-8") as file_handle:
+        json.dump(comparison_statistics, file_handle, indent=2, sort_keys=True)
 
 
 def _distribution_cumulative_variants(quantity_name):
diff --git a/tests/unit_tests/sim_events/test_production_comparison.py b/tests/unit_tests/sim_events/test_production_comparison.py
index 0605ee162d..98dcde575e 100644
--- a/tests/unit_tests/sim_events/test_production_comparison.py
+++ b/tests/unit_tests/sim_events/test_production_comparison.py
@@ -5,6 +5,8 @@
 
 from simtools.sim_events.production_comparison import (
     collect_production_metrics,
+    compute_ks_chi2_for_aligned_histograms,
+    compute_ks_chi2_for_samples,
     parse_production_arguments,
 )
 
@@ -192,3 +194,33 @@ def test_collect_production_metrics_empty_input_keeps_arrays_empty(mocker):
     assert metrics[0].trigger_fraction == pytest.approx(0.0)
     np.testing.assert_array_equal(metrics[0].simulated_energies, np.array([]))
     np.testing.assert_array_equal(metrics[0].trigger_multiplicity, np.array([], dtype=int))
+
+
+def test_compute_ks_chi2_for_samples_returns_statistics():
+    """Test KS/Chi2 statistics for sample-based comparisons."""
+    baseline = np.array([1.0, 2.0, 3.0, 4.0])
+    candidate = np.array([1.0, 2.0, 2.5, 4.0])
+    bin_edges = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
+
+    result = compute_ks_chi2_for_samples(baseline, candidate, bin_edges)
+
+    assert result["ks_statistic"] is not None
+    assert result["ks_pvalue"] is not None
+    assert result["chi2_statistic"] is not None
+    assert result["chi2_pvalue"] is not None
+    assert result["valid"]
+    assert result["reason"] == "ok"
+    assert result["bin_edges"] == [1.0, 2.0, 3.0, 4.0, 5.0]
+
+
+def test_compute_ks_chi2_for_aligned_histograms_handles_empty_counts():
+    """Test aligned histogram statistics return invalid for empty counts."""
+    baseline = np.array([0.0, 0.0, 0.0])
+    candidate = np.array([1.0, 0.0, 0.0])
+
+    result = compute_ks_chi2_for_aligned_histograms(baseline, candidate)
+
+    assert result["ks_statistic"] is None
+    assert result["chi2_statistic"] is None
+    assert result["valid"] is False
+    assert result["reason"] == "insufficient_data"
diff --git a/tests/unit_tests/visualization/test_plot_event_level_production_comparison.py b/tests/unit_tests/visualization/test_plot_event_level_production_comparison.py
index ce31046a8b..ae9cc44336 100644
--- a/tests/unit_tests/visualization/test_plot_event_level_production_comparison.py
+++ b/tests/unit_tests/visualization/test_plot_event_level_production_comparison.py
@@ -1,3 +1,4 @@
+import json
 from collections import Counter
 from pathlib import Path
 
@@ -44,6 +45,7 @@ def test_plot_writes_event_level_comparison_figures(tmp_test_directory):
     plot_event_level_production_comparison.plot(metrics, output_path=output_path, bins=8)
 
     expected_files = [
+        "comparison_statistics.json",
         "trigger_multiplicity.png",
         "trigger_combination.png",
         "distribution_energy.png",
@@ -55,6 +57,13 @@ def test_plot_writes_event_level_comparison_figures(tmp_test_directory):
     ]
     _assert_files_exist(output_path, expected_files)
 
+    with (output_path / "comparison_statistics.json").open(encoding="utf-8") as file_handle:
+        stats_payload = json.load(file_handle)
+    assert stats_payload["baseline"]["label"] == "baseline"
+    assert [item["label"] for item in stats_payload["comparison_sets"]] == ["candidate"]
+    assert "distribution_energy" in stats_payload["plot_statistics"]
+    assert "trigger_multiplicity" in stats_payload["plot_statistics"]
+
 
 def test_plot_writes_per_type_comparison_figures(tmp_test_directory):
     """Test per-array-element-type comparison plots are written for each type."""

From 3bfdc5f6bd096f815c7f4f8eebe9fb0e7483c0d6 Mon Sep 17 00:00:00 2001
From: Gernot Maier <gernot.maier@desy.de>
Date: Mon, 18 May 2026 10:02:09 +0200
Subject: [PATCH 18/27] chi2/n [skip ci]

---
 .../sim_events/production_comparison.py       | 11 ++++++--
 .../plot_event_level_production_comparison.py | 28 ++++++++++++++++---
 .../sim_events/test_production_comparison.py  | 16 ++++-------
 3 files changed, 38 insertions(+), 17 deletions(-)

diff --git a/src/simtools/sim_events/production_comparison.py b/src/simtools/sim_events/production_comparison.py
index e02e733b0e..d36837e256 100644
--- a/src/simtools/sim_events/production_comparison.py
+++ b/src/simtools/sim_events/production_comparison.py
@@ -325,13 +325,20 @@ def _accumulate_per_subset(telescopes, energy, core_dist, angular_dist, accumula
 
 
 def _subset_counts_for_event(telescopes):
-    """Return subset keys and multiplicities for one triggered event."""
+    """Return subset keys and multiplicities for one triggered event.
+
+    Per-type keys (e.g. ``LSTN``, ``MSTN``) are filled only for events with
+    a single participating telescope type. Mixed-type events are tracked in
+    ``mixed_type`` only.
+    """
     type_counts = Counter()
     for telescope in telescopes:
         tel_type = names.get_array_element_type_from_name(telescope)
         type_counts[tel_type] += 1
 
-    subset_counts = dict(type_counts)
+    subset_counts = {}
+    if len(type_counts) == 1:
+        subset_counts = dict(type_counts)
     if len(telescopes) == 1:
         subset_counts["single_telescope"] = 1
     elif len(type_counts) > 1:
diff --git a/src/simtools/visualization/plot_event_level_production_comparison.py b/src/simtools/visualization/plot_event_level_production_comparison.py
index ab4504603c..08e626e71a 100644
--- a/src/simtools/visualization/plot_event_level_production_comparison.py
+++ b/src/simtools/visualization/plot_event_level_production_comparison.py
@@ -664,7 +664,7 @@ def _build_quantity_distribution_statistics(
 
 
 def _annotate_ks_statistics(ax, statistics):
-    """Overlay KS values for baseline-vs-candidate comparisons on a plot."""
+    """Overlay KS and Chi2 values for baseline-vs-candidate comparisons on a plot."""
     if statistics is None or not statistics.get("comparisons"):
         return
 
@@ -675,12 +675,19 @@ def _annotate_ks_statistics(ax, statistics):
         if metric_type == "quantity_distribution":
             sim_ks = comparison.get("simulated", {}).get("ks_statistic")
             trig_ks = comparison.get("triggered", {}).get("ks_statistic")
+            sim_chi2_over_n = _format_chi2_over_n(comparison.get("simulated", {}))
+            trig_chi2_over_n = _format_chi2_over_n(comparison.get("triggered", {}))
             lines.append(
-                f"{label}: KS(sim)={_format_statistic_value(sim_ks)}, "
-                f"KS(trg)={_format_statistic_value(trig_ks)}"
+                f"{label}: KS s/t={_format_statistic_value(sim_ks)}/"
+                f"{_format_statistic_value(trig_ks)}, "
+                f"Chi2/N s/t={sim_chi2_over_n}/{trig_chi2_over_n}"
             )
             continue
-        lines.append(f"{label}: KS={_format_statistic_value(comparison.get('ks_statistic'))}")
+        chi2_over_n = _format_chi2_over_n(comparison)
+        lines.append(
+            f"{label}: KS={_format_statistic_value(comparison.get('ks_statistic'))}, "
+            f"Chi2/N={chi2_over_n}"
+        )
 
     ax.text(
         0.02,
@@ -701,6 +708,19 @@ def _format_statistic_value(value):
     return f"{float(value):.4f}"
 
 
+def _format_chi2_over_n(statistics):
+    """Format Chi2/N using bins that contribute to the Chi2 calculation."""
+    chi2_value = statistics.get("chi2_statistic")
+    baseline_counts = statistics.get("baseline_counts")
+    if chi2_value is None or baseline_counts is None:
+        return "n/a"
+
+    n_contributing_bins = int(np.sum(np.asarray(baseline_counts, dtype=float) > 0))
+    if n_contributing_bins <= 0:
+        return "n/a"
+    return f"{float(chi2_value) / n_contributing_bins:.4f}"
+
+
 def _initialize_comparison_statistics(metrics_per_production):
     """Initialize top-level comparison statistics structure."""
     baseline = metrics_per_production[0] if metrics_per_production else None
diff --git a/tests/unit_tests/sim_events/test_production_comparison.py b/tests/unit_tests/sim_events/test_production_comparison.py
index 98dcde575e..b9eaec9f4c 100644
--- a/tests/unit_tests/sim_events/test_production_comparison.py
+++ b/tests/unit_tests/sim_events/test_production_comparison.py
@@ -96,21 +96,15 @@ def test_collect_production_metrics_aggregates_event_data(mocker):
     assert metrics[0].trigger_fraction == pytest.approx(2.0 / 3.0)
     np.testing.assert_array_equal(metrics[0].simulated_angular_distances, np.array([0.5, 1.0, 1.5]))
     np.testing.assert_array_equal(metrics[0].triggered_angular_distances, np.array([0.5, 1.5]))
-    assert set(metrics[0].per_type.keys()) == {"LSTN", "MSTN", "single_telescope", "mixed_type"}
-    assert metrics[0].per_type["LSTN"].triggered_event_count == 1
-    assert metrics[0].per_type["MSTN"].triggered_event_count == 2
-    np.testing.assert_array_equal(metrics[0].per_type["LSTN"].triggered_energies, [1.0])
-    np.testing.assert_array_equal(metrics[0].per_type["MSTN"].triggered_energies, [1.0, 3.0])
-    np.testing.assert_array_equal(metrics[0].per_type["LSTN"].trigger_multiplicity, [1])
-    np.testing.assert_array_equal(metrics[0].per_type["MSTN"].trigger_multiplicity, [1, 1])
+    assert set(metrics[0].per_type.keys()) == {"MSTN", "single_telescope", "mixed_type"}
+    assert metrics[0].per_type["MSTN"].triggered_event_count == 1
+    np.testing.assert_array_equal(metrics[0].per_type["MSTN"].triggered_energies, [3.0])
+    np.testing.assert_array_equal(metrics[0].per_type["MSTN"].trigger_multiplicity, [1])
     assert metrics[0].per_type["single_telescope"].triggered_event_count == 1
     assert metrics[0].per_type["mixed_type"].triggered_event_count == 1
     np.testing.assert_array_equal(metrics[0].per_type["single_telescope"].triggered_energies, [3.0])
     np.testing.assert_array_equal(metrics[0].per_type["mixed_type"].triggered_energies, [1.0])
-    np.testing.assert_array_equal(metrics[0].per_type["LSTN"].triggered_angular_distances, [0.5])
-    np.testing.assert_array_equal(
-        metrics[0].per_type["MSTN"].triggered_angular_distances, [0.5, 1.5]
-    )
+    np.testing.assert_array_equal(metrics[0].per_type["MSTN"].triggered_angular_distances, [1.5])
 
 
 @pytest.mark.parametrize(

From 8841720080e7112fda6cdab7587e60ce02eec414 Mon Sep 17 00:00:00 2001
From: Gernot Maier <gernot.maier@desy.de>
Date: Mon, 18 May 2026 10:46:13 +0200
Subject: [PATCH 19/27] fixed binning

---
 .../sim_events/production_comparison.py       | 10 +++---
 .../plot_event_level_production_comparison.py | 34 ++++++++++++++-----
 .../sim_events/test_production_comparison.py  | 27 ++++++++-------
 ..._plot_event_level_production_comparison.py | 25 ++++++++++++++
 4 files changed, 69 insertions(+), 27 deletions(-)

diff --git a/src/simtools/sim_events/production_comparison.py b/src/simtools/sim_events/production_comparison.py
index d36837e256..8d5362048d 100644
--- a/src/simtools/sim_events/production_comparison.py
+++ b/src/simtools/sim_events/production_comparison.py
@@ -327,18 +327,16 @@ def _accumulate_per_subset(telescopes, energy, core_dist, angular_dist, accumula
 def _subset_counts_for_event(telescopes):
     """Return subset keys and multiplicities for one triggered event.
 
-    Per-type keys (e.g. ``LSTN``, ``MSTN``) are filled only for events with
-    a single participating telescope type. Mixed-type events are tracked in
-    ``mixed_type`` only.
+    Per-type keys (e.g. ``LSTN``, ``MSTN``) are filled when at least two
+    telescopes of that type triggered, independent of whether the event is
+    single-type or mixed-type.
     """
     type_counts = Counter()
     for telescope in telescopes:
         tel_type = names.get_array_element_type_from_name(telescope)
         type_counts[tel_type] += 1
 
-    subset_counts = {}
-    if len(type_counts) == 1:
-        subset_counts = dict(type_counts)
+    subset_counts = {tel_type: count for tel_type, count in type_counts.items() if count >= 2}
     if len(telescopes) == 1:
         subset_counts["single_telescope"] = 1
     elif len(type_counts) > 1:
diff --git a/src/simtools/visualization/plot_event_level_production_comparison.py b/src/simtools/visualization/plot_event_level_production_comparison.py
index 08e626e71a..7f343e4afb 100644
--- a/src/simtools/visualization/plot_event_level_production_comparison.py
+++ b/src/simtools/visualization/plot_event_level_production_comparison.py
@@ -489,8 +489,13 @@ def _plot_quantity_distribution(
         Whether to plot cumulative distributions.
     """
     fig, ax = plt.subplots(figsize=(9, 6))
+    bin_edges = _get_global_quantity_bin_edges(
+        metrics_per_production,
+        quantity_name=quantity_name,
+        x_scale=x_scale,
+        bins=bins,
+    )
     plotted = False
-    bin_edges_by_label = {}
     counts_by_label = {}
     for metrics in metrics_per_production:
         simulated, triggered = _get_quantity_arrays(metrics, quantity_name)
@@ -498,10 +503,8 @@ def _plot_quantity_distribution(
             continue
         plotted = True
 
-        bin_edges = _get_bin_edges(simulated, x_scale=x_scale, bins=bins)
         sim_counts, _ = np.histogram(simulated, bins=bin_edges)
         trig_counts, _ = np.histogram(triggered, bins=bin_edges)
-        bin_edges_by_label[metrics.label] = bin_edges
         counts_by_label[metrics.label] = {
             "simulated": sim_counts,
             "triggered": trig_counts,
@@ -545,7 +548,7 @@ def _plot_quantity_distribution(
         metrics_per_production,
         quantity_name,
         counts_by_label,
-        bin_edges_by_label,
+        baseline_bin_edges=bin_edges,
         cumulative=cumulative,
     )
     _annotate_ks_statistics(ax, statistics)
@@ -556,6 +559,20 @@ def _plot_quantity_distribution(
     return statistics
 
 
+def _get_global_quantity_bin_edges(metrics_per_production, quantity_name, x_scale, bins):
+    """Return one shared set of bin edges for a quantity across productions."""
+    all_simulated_values = []
+    for metrics in metrics_per_production:
+        simulated, _ = _get_quantity_arrays(metrics, quantity_name)
+        if simulated.size == 0:
+            continue
+        all_simulated_values.append(np.asarray(simulated))
+
+    if len(all_simulated_values) == 0:
+        return _get_bin_edges(np.array([]), x_scale=x_scale, bins=bins)
+    return _get_bin_edges(np.concatenate(all_simulated_values), x_scale=x_scale, bins=bins)
+
+
 def _plot_telescope_participation(metrics_per_production, output_path):
     """Plot telescope participation fractions by production."""
     telescopes = sorted(
@@ -619,7 +636,7 @@ def _build_quantity_distribution_statistics(
     metrics_per_production,
     quantity_name,
     counts_by_label,
-    bin_edges_by_label,
+    baseline_bin_edges,
     cumulative,
 ):
     """Build comparison statistics for quantity distribution plots."""
@@ -635,7 +652,6 @@ def _build_quantity_distribution_statistics(
         return stats_summary
 
     baseline_data = counts_by_label[baseline_label]
-    baseline_edges = bin_edges_by_label[baseline_label]
 
     for metrics in metrics_per_production[1:]:
         if metrics.label not in counts_by_label:
@@ -645,19 +661,19 @@ def _build_quantity_distribution_statistics(
         sim_stats = compute_ks_chi2_for_samples(
             baseline_data["simulated_samples"],
             candidate_data["simulated_samples"],
-            baseline_edges,
+            baseline_bin_edges,
         )
         trig_stats = compute_ks_chi2_for_samples(
             baseline_data["triggered_samples"],
             candidate_data["triggered_samples"],
-            baseline_edges,
+            baseline_bin_edges,
         )
 
         comparison_record = {
             "candidate_label": metrics.label,
             "simulated": sim_stats,
             "triggered": trig_stats,
-            "baseline_bin_edges": np.asarray(baseline_edges, dtype=float).tolist(),
+            "baseline_bin_edges": np.asarray(baseline_bin_edges, dtype=float).tolist(),
         }
         stats_summary["comparisons"].append(comparison_record)
     return stats_summary
diff --git a/tests/unit_tests/sim_events/test_production_comparison.py b/tests/unit_tests/sim_events/test_production_comparison.py
index b9eaec9f4c..e2b1ec16b0 100644
--- a/tests/unit_tests/sim_events/test_production_comparison.py
+++ b/tests/unit_tests/sim_events/test_production_comparison.py
@@ -60,14 +60,15 @@ def test_collect_production_metrics_aggregates_event_data(mocker):
         angular_distance=np.array([0.5, 1.0, 1.5]),
     )
     triggered_shower_data = SimpleNamespace(
-        simulated_energy=np.array([1.0, 3.0]),
-        core_distance_shower=np.array([10.0, 30.0]),
-        angular_distance=np.array([0.5, 1.5]),
+        simulated_energy=np.array([1.0, 3.0, 2.0]),
+        core_distance_shower=np.array([10.0, 30.0, 20.0]),
+        angular_distance=np.array([0.5, 1.5, 1.2]),
     )
     triggered_data = SimpleNamespace(
         telescope_list=[
             np.array(["LSTN-01", "MSTN-01"]),
             np.array(["MSTN-01"]),
+            np.array(["MSTN-01", "MSTN-02"]),
         ]
     )
 
@@ -85,26 +86,28 @@ def test_collect_production_metrics_aggregates_event_data(mocker):
     assert len(metrics) == 1
     assert metrics[0].label == "baseline"
     assert metrics[0].simulated_event_count == 3
-    assert metrics[0].triggered_event_count == 2
+    assert metrics[0].triggered_event_count == 3
     np.testing.assert_array_equal(metrics[0].simulated_energies, np.array([1.0, 2.0, 3.0]))
-    np.testing.assert_array_equal(metrics[0].triggered_energies, np.array([1.0, 3.0]))
-    np.testing.assert_array_equal(metrics[0].trigger_multiplicity, np.array([2, 1]))
+    np.testing.assert_array_equal(metrics[0].triggered_energies, np.array([1.0, 3.0, 2.0]))
+    np.testing.assert_array_equal(metrics[0].trigger_multiplicity, np.array([2, 1, 2]))
     assert metrics[0].trigger_combinations["LSTN-01,MSTN-01"] == 1
     assert metrics[0].trigger_combinations["MSTN-01"] == 1
-    assert metrics[0].telescope_participation["MSTN-01"] == 2
+    assert metrics[0].trigger_combinations["MSTN-01,MSTN-02"] == 1
+    assert metrics[0].telescope_participation["MSTN-01"] == 3
+    assert metrics[0].telescope_participation["MSTN-02"] == 1
     assert metrics[0].telescope_participation["LSTN-01"] == 1
-    assert metrics[0].trigger_fraction == pytest.approx(2.0 / 3.0)
+    assert metrics[0].trigger_fraction == pytest.approx(1.0)
     np.testing.assert_array_equal(metrics[0].simulated_angular_distances, np.array([0.5, 1.0, 1.5]))
-    np.testing.assert_array_equal(metrics[0].triggered_angular_distances, np.array([0.5, 1.5]))
+    np.testing.assert_array_equal(metrics[0].triggered_angular_distances, np.array([0.5, 1.5, 1.2]))
     assert set(metrics[0].per_type.keys()) == {"MSTN", "single_telescope", "mixed_type"}
     assert metrics[0].per_type["MSTN"].triggered_event_count == 1
-    np.testing.assert_array_equal(metrics[0].per_type["MSTN"].triggered_energies, [3.0])
-    np.testing.assert_array_equal(metrics[0].per_type["MSTN"].trigger_multiplicity, [1])
+    np.testing.assert_array_equal(metrics[0].per_type["MSTN"].triggered_energies, [2.0])
+    np.testing.assert_array_equal(metrics[0].per_type["MSTN"].trigger_multiplicity, [2])
     assert metrics[0].per_type["single_telescope"].triggered_event_count == 1
     assert metrics[0].per_type["mixed_type"].triggered_event_count == 1
     np.testing.assert_array_equal(metrics[0].per_type["single_telescope"].triggered_energies, [3.0])
     np.testing.assert_array_equal(metrics[0].per_type["mixed_type"].triggered_energies, [1.0])
-    np.testing.assert_array_equal(metrics[0].per_type["MSTN"].triggered_angular_distances, [1.5])
+    np.testing.assert_array_equal(metrics[0].per_type["MSTN"].triggered_angular_distances, [1.2])
 
 
 @pytest.mark.parametrize(
diff --git a/tests/unit_tests/visualization/test_plot_event_level_production_comparison.py b/tests/unit_tests/visualization/test_plot_event_level_production_comparison.py
index ae9cc44336..6931f315ca 100644
--- a/tests/unit_tests/visualization/test_plot_event_level_production_comparison.py
+++ b/tests/unit_tests/visualization/test_plot_event_level_production_comparison.py
@@ -244,3 +244,28 @@ def test_single_and_mixed_trigger_skip_paths(tmp_test_directory):
     plot_event_level_production_comparison._plot_mixed_trigger_combinations([metric], output_path)
     assert not (output_path / "single_telescope_trigger_distribution.png").exists()
     assert not (output_path / "mixed_trigger_combinations.png").exists()
+
+
+def test_quantity_distribution_uses_shared_bin_edges_for_statistics(tmp_test_directory):
+    """Test shared binning keeps candidate samples in-range for statistics."""
+    output_path = Path(tmp_test_directory)
+    baseline = _build_metrics("baseline", simulated_scale=1.0, triggered_scale=1.0)
+    candidate = _build_metrics("candidate", simulated_scale=10.0, triggered_scale=10.0)
+
+    stats = plot_event_level_production_comparison._plot_quantity_distribution(
+        [baseline, candidate],
+        output_path,
+        quantity_name="angular_distance",
+        x_label="Angular Distance (deg)",
+        x_scale="linear",
+        bins=8,
+        cumulative=False,
+    )
+
+    comparison = stats["comparisons"][0]
+    assert sum(comparison["simulated"]["candidate_counts"]) == len(
+        candidate.simulated_angular_distances
+    )
+    assert sum(comparison["triggered"]["candidate_counts"]) == len(
+        candidate.triggered_angular_distances
+    )

From dd275dff9f552b6d9e08887d9b31314a3955101a Mon Sep 17 00:00:00 2001
From: Gernot Maier <gernot.maier@desy.de>
Date: Mon, 18 May 2026 11:24:20 +0200
Subject: [PATCH 20/27] statistics module

---
 .../applications/compare_productions.py       |  84 +++++++-
 .../sim_events/production_comparison.py       | 193 ------------------
 src/simtools/statistics.py                    | 116 +++++++++++
 .../plot_event_level_production_comparison.py |  75 +++----
 .../applications/test_compare_productions.py  |  83 ++++++++
 .../sim_events/test_production_comparison.py  |  96 ++-------
 ..._plot_event_level_production_comparison.py |  24 +++
 7 files changed, 358 insertions(+), 313 deletions(-)
 create mode 100644 src/simtools/statistics.py

diff --git a/src/simtools/applications/compare_productions.py b/src/simtools/applications/compare_productions.py
index df6517b627..87292e5d39 100644
--- a/src/simtools/applications/compare_productions.py
+++ b/src/simtools/applications/compare_productions.py
@@ -22,9 +22,10 @@
 
 from simtools.application_control import build_application
 from simtools.sim_events.production_comparison import (
+    ProductionDescriptor,
     collect_production_metrics,
-    parse_production_arguments,
 )
+from simtools.utils.general import resolve_file_patterns
 from simtools.visualization import plot_event_level_production_comparison
 
 
@@ -51,6 +52,87 @@ def _add_arguments(parser):
     )
 
 
+def parse_production_arguments(production_arguments):
+    """Parse repeated production arguments into validated descriptors.
+
+    Parameters
+    ----------
+    production_arguments : list[list[str]]
+        Repeated ``--production`` arguments in the shape ``[label, patterns]``.
+
+    Returns
+    -------
+    list[ProductionDescriptor]
+        Validated and normalized production descriptors.
+
+    Raises
+    ------
+    ValueError
+        If configuration is malformed or does not contain any production.
+    """
+    parsed_productions = _normalize_production_arguments(production_arguments)
+    if not parsed_productions:
+        raise ValueError("At least one production is required.")
+
+    labels = [label for label, _ in parsed_productions]
+    if len(set(labels)) != len(labels):
+        raise ValueError("Production labels must be unique.")
+
+    descriptors = []
+    for label, pattern_list in parsed_productions:
+        patterns = [pattern.strip() for pattern in pattern_list.split(",") if pattern.strip()]
+        if len(patterns) == 0:
+            raise ValueError(f"Production '{label}' has no event_data_file pattern.")
+
+        resolved_files = [str(path) for path in resolve_file_patterns(patterns)]
+        if len(resolved_files) == 0:
+            raise ValueError(f"Production '{label}' does not resolve to any files.")
+        descriptors.append(ProductionDescriptor(label=label, event_data_files=resolved_files))
+
+    return descriptors
+
+
+def _normalize_production_arguments(production_arguments):
+    """Normalize raw production arguments into ``[(label, files), ...]``."""
+    if not production_arguments:
+        return []
+
+    normalized = []
+    if all(isinstance(item, str) for item in production_arguments):
+        return _pairwise_label_file_arguments(production_arguments)
+
+    for item in production_arguments:
+        normalized.extend(_normalize_single_production_argument(item))
+
+    return normalized
+
+
+def _pairwise_label_file_arguments(flat_arguments):
+    """Convert a flat list of strings into ``[(label, files), ...]`` pairs."""
+    if len(flat_arguments) % 2 != 0:
+        _raise_invalid_production_arguments()
+    return [
+        (flat_arguments[index], flat_arguments[index + 1])
+        for index in range(0, len(flat_arguments), 2)
+    ]
+
+
+def _normalize_single_production_argument(argument):
+    """Normalize one nested production argument into label/file pairs."""
+    if not isinstance(argument, list | tuple):
+        _raise_invalid_production_arguments()
+    if not all(isinstance(value, str) for value in argument):
+        _raise_invalid_production_arguments()
+    if len(argument) == 2:
+        return [(argument[0], argument[1])]
+    return _pairwise_label_file_arguments(list(argument))
+
+
+def _raise_invalid_production_arguments():
+    """Raise a standardized parser error for malformed production arguments."""
+    raise ValueError("Production arguments must be provided as label/file pairs.")
+
+
 def main():
     """See CLI description."""
     app_context = build_application(
diff --git a/src/simtools/sim_events/production_comparison.py b/src/simtools/sim_events/production_comparison.py
index 8d5362048d..630d6fd384 100644
--- a/src/simtools/sim_events/production_comparison.py
+++ b/src/simtools/sim_events/production_comparison.py
@@ -5,11 +5,9 @@
 from dataclasses import dataclass, field
 
 import numpy as np
-from scipy import stats
 
 from simtools.sim_events.reader import EventDataReader
 from simtools.utils import names
-from simtools.utils.general import resolve_file_patterns
 
 
 @dataclass
@@ -37,7 +35,6 @@ class ProductionEventMetrics:
     simulated_angular_distances: np.ndarray = field(default_factory=lambda: np.array([]))
     triggered_angular_distances: np.ndarray = field(default_factory=lambda: np.array([]))
     per_type: dict = field(default_factory=dict)
-    comparison_statistics: dict = field(default_factory=dict)
 
     @property
     def trigger_fraction(self):
@@ -47,87 +44,6 @@ def trigger_fraction(self):
         return self.triggered_event_count / self.simulated_event_count
 
 
-def parse_production_arguments(production_arguments):
-    """Parse repeated production arguments into validated descriptors.
-
-    Parameters
-    ----------
-    production_arguments : list[list[str]]
-        Repeated ``--production`` arguments in the shape ``[label, patterns]``.
-
-    Returns
-    -------
-    list[ProductionDescriptor]
-        Validated and normalized production descriptors.
-
-    Raises
-    ------
-    ValueError
-        If configuration is malformed or does not contain any production.
-    """
-    parsed_productions = _normalize_production_arguments(production_arguments)
-    if not parsed_productions:
-        raise ValueError("At least one production is required.")
-
-    labels = [item[0] for item in parsed_productions]
-    if len(set(labels)) != len(labels):
-        raise ValueError("Production labels must be unique.")
-
-    descriptors = []
-    for label, pattern_list in parsed_productions:
-        patterns = [pattern.strip() for pattern in pattern_list.split(",") if pattern.strip()]
-        if len(patterns) == 0:
-            raise ValueError(f"Production '{label}' has no event_data_file pattern.")
-
-        resolved_files = [str(path) for path in resolve_file_patterns(patterns)]
-        if len(resolved_files) == 0:
-            raise ValueError(f"Production '{label}' does not resolve to any files.")
-        descriptors.append(ProductionDescriptor(label=label, event_data_files=resolved_files))
-
-    return descriptors
-
-
-def _normalize_production_arguments(production_arguments):
-    """Normalize raw production arguments into ``[(label, files), ...]``."""
-    if not production_arguments:
-        return []
-
-    normalized = []
-    if all(isinstance(item, str) for item in production_arguments):
-        return _pairwise_label_file_arguments(production_arguments)
-
-    for item in production_arguments:
-        normalized.extend(_normalize_single_production_argument(item))
-
-    return normalized
-
-
-def _pairwise_label_file_arguments(flat_arguments):
-    """Convert a flat list of strings into ``[(label, files), ...]`` pairs."""
-    if len(flat_arguments) % 2 != 0:
-        _raise_invalid_production_arguments()
-    return [
-        (flat_arguments[index], flat_arguments[index + 1])
-        for index in range(0, len(flat_arguments), 2)
-    ]
-
-
-def _normalize_single_production_argument(argument):
-    """Normalize one nested production argument into label/file pairs."""
-    if not isinstance(argument, list | tuple):
-        _raise_invalid_production_arguments()
-    if not all(isinstance(value, str) for value in argument):
-        _raise_invalid_production_arguments()
-    if len(argument) == 2:
-        return [(argument[0], argument[1])]
-    return _pairwise_label_file_arguments(list(argument))
-
-
-def _raise_invalid_production_arguments():
-    """Raise a standardized parser error for malformed production arguments."""
-    raise ValueError("Production arguments must be provided as label/file pairs.")
-
-
 def collect_production_metrics(production_descriptors, telescope_list=None):
     """Collect event-level comparison metrics for each production.
 
@@ -374,112 +290,3 @@ def _safe_concat(arrays, dtype=float):
     if len(arrays) == 0:
         return np.array([], dtype=dtype)
     return np.concatenate(arrays)
-
-
-def compute_ks_chi2_for_samples(baseline_samples, candidate_samples, bin_edges):
-    """Compute KS and Chi2 statistics for two sample arrays.
-
-    Parameters
-    ----------
-    baseline_samples : np.ndarray
-        Baseline sample values.
-    candidate_samples : np.ndarray
-        Candidate sample values.
-    bin_edges : np.ndarray
-        Histogram bin edges used to derive Chi2 counts.
-
-    Returns
-    -------
-    dict
-        KS and Chi2 statistics including p-values and counts metadata.
-    """
-    baseline_samples = np.asarray(baseline_samples)
-    candidate_samples = np.asarray(candidate_samples)
-    result = {
-        "ks_statistic": None,
-        "ks_pvalue": None,
-        "chi2_statistic": None,
-        "chi2_pvalue": None,
-        "valid": False,
-        "reason": "insufficient_data",
-    }
-
-    if baseline_samples.size == 0 or candidate_samples.size == 0:
-        return result
-
-    ks_test = stats.ks_2samp(baseline_samples, candidate_samples)
-    result["ks_statistic"] = float(ks_test.statistic)
-    result["ks_pvalue"] = float(ks_test.pvalue)
-
-    baseline_counts, _ = np.histogram(baseline_samples, bins=bin_edges)
-    candidate_counts, _ = np.histogram(candidate_samples, bins=bin_edges)
-    chi2_result = compute_ks_chi2_for_aligned_histograms(baseline_counts, candidate_counts)
-    result["chi2_statistic"] = chi2_result["chi2_statistic"]
-    result["chi2_pvalue"] = chi2_result["chi2_pvalue"]
-    result["valid"] = chi2_result["valid"]
-    result["reason"] = chi2_result["reason"]
-    result["baseline_counts"] = baseline_counts.astype(int).tolist()
-    result["candidate_counts"] = candidate_counts.astype(int).tolist()
-    result["bin_edges"] = np.asarray(bin_edges, dtype=float).tolist()
-    return result
-
-
-def compute_ks_chi2_for_aligned_histograms(baseline_counts, candidate_counts):
-    """Compute KS and Chi2 statistics from aligned count arrays.
-
-    Parameters
-    ----------
-    baseline_counts : np.ndarray
-        Baseline histogram counts on aligned support.
-    candidate_counts : np.ndarray
-        Candidate histogram counts on aligned support.
-
-    Returns
-    -------
-    dict
-        KS and Chi2 statistics including p-values and validity metadata.
-    """
-    baseline_counts = np.asarray(baseline_counts, dtype=float)
-    candidate_counts = np.asarray(candidate_counts, dtype=float)
-
-    result = {
-        "ks_statistic": None,
-        "ks_pvalue": None,
-        "chi2_statistic": None,
-        "chi2_pvalue": None,
-        "valid": False,
-        "reason": "insufficient_data",
-    }
-
-    baseline_total = float(np.sum(baseline_counts))
-    candidate_total = float(np.sum(candidate_counts))
-    if baseline_total <= 0 or candidate_total <= 0:
-        return result
-
-    baseline_cdf = np.cumsum(baseline_counts) / baseline_total
-    candidate_cdf = np.cumsum(candidate_counts) / candidate_total
-    result["ks_statistic"] = float(np.max(np.abs(baseline_cdf - candidate_cdf)))
-
-    support_mask = baseline_counts > 0
-    if not np.any(support_mask):
-        result["reason"] = "zero_baseline_counts"
-        return result
-
-    observed = candidate_counts[support_mask]
-    observed_total = float(np.sum(observed))
-    baseline_supported_total = float(np.sum(baseline_counts[support_mask]))
-    if observed_total <= 0 or baseline_supported_total <= 0:
-        result["reason"] = "invalid_expected_counts"
-        return result
-
-    expected = baseline_counts[support_mask] * (observed_total / baseline_supported_total)
-    if np.sum(expected) <= 0:
-        result["reason"] = "invalid_expected_counts"
-        return result
-
-    chi2 = stats.chisquare(f_obs=observed, f_exp=expected)
-    result["chi2_statistic"] = float(chi2.statistic)
-    result["chi2_pvalue"] = float(chi2.pvalue)
-    result["valid"] = True
-    result["reason"] = "ok"
-    return result
diff --git a/src/simtools/statistics.py b/src/simtools/statistics.py
new file mode 100644
index 0000000000..ee7bc36af9
--- /dev/null
+++ b/src/simtools/statistics.py
@@ -0,0 +1,116 @@
+"""Generic statistical comparison helpers for simtools and other applications."""
+
+import numpy as np
+from scipy import stats
+
+
+def ks_test_samples(samples1, samples2):
+    """Compute the Kolmogorov-Smirnov (KS) statistic and p-value for two 1D sample arrays.
+
+    Parameters
+    ----------
+    samples1 : np.ndarray
+        First sample array.
+    samples2 : np.ndarray
+        Second sample array.
+
+    Returns
+    -------
+    ks_statistic : float
+        KS statistic value.
+    ks_pvalue : float
+        KS test p-value.
+    """
+    samples1 = np.asarray(samples1)
+    samples2 = np.asarray(samples2)
+    if samples1.size == 0 or samples2.size == 0:
+        return None, None
+    ks_test = stats.ks_2samp(samples1, samples2)
+    return float(ks_test.statistic), float(ks_test.pvalue)
+
+
+def chi2_test_histograms(counts1, counts2):
+    """Compute the Chi2 statistic and p-value for two aligned histogram count arrays.
+
+    Parameters
+    ----------
+    counts1 : np.ndarray
+        First histogram counts (expected/baseline).
+    counts2 : np.ndarray
+        Second histogram counts (observed/candidate).
+
+    Returns
+    -------
+    chi2_statistic : float
+        Chi2 statistic value.
+    chi2_pvalue : float
+        Chi2 test p-value.
+    valid : bool
+        Whether the test is valid (enough data, no zero division).
+    reason : str
+        Reason for invalidity if not valid.
+    """
+    counts1 = np.asarray(counts1, dtype=float)
+    counts2 = np.asarray(counts2, dtype=float)
+    total1 = float(np.sum(counts1))
+    total2 = float(np.sum(counts2))
+    if total1 <= 0 or total2 <= 0:
+        return None, None, False, "insufficient_data"
+    support_mask = counts1 > 0
+    if not np.any(support_mask):
+        return None, None, False, "zero_baseline_counts"
+    observed = counts2[support_mask]
+    observed_total = float(np.sum(observed))
+    baseline_supported_total = float(np.sum(counts1[support_mask]))
+    if observed_total <= 0 or baseline_supported_total <= 0:
+        return None, None, False, "invalid_expected_counts"
+    expected = counts1[support_mask] * (observed_total / baseline_supported_total)
+    if np.sum(expected) <= 0:
+        return None, None, False, "invalid_expected_counts"
+    chi2 = stats.chisquare(f_obs=observed, f_exp=expected)
+    return float(chi2.statistic), float(chi2.pvalue), True, "ok"
+
+
+def compare_samples_with_statistics(baseline_samples, candidate_samples, bin_edges):
+    """Compute KS and Chi2 statistics for two sample arrays using provided bin edges.
+
+    Parameters
+    ----------
+    baseline_samples : np.ndarray
+        Baseline sample values.
+    candidate_samples : np.ndarray
+        Candidate sample values.
+    bin_edges : np.ndarray
+        Histogram bin edges for Chi2 calculation.
+
+    Returns
+    -------
+    dict
+        KS and Chi2 statistics, p-values, and histogram metadata.
+    """
+    baseline_samples = np.asarray(baseline_samples)
+    candidate_samples = np.asarray(candidate_samples)
+    result = {
+        "ks_statistic": None,
+        "ks_pvalue": None,
+        "chi2_statistic": None,
+        "chi2_pvalue": None,
+        "valid": False,
+        "reason": "insufficient_data",
+    }
+    if baseline_samples.size == 0 or candidate_samples.size == 0:
+        return result
+    ks_stat, ks_pval = ks_test_samples(baseline_samples, candidate_samples)
+    result["ks_statistic"] = ks_stat
+    result["ks_pvalue"] = ks_pval
+    baseline_counts, _ = np.histogram(baseline_samples, bins=bin_edges)
+    candidate_counts, _ = np.histogram(candidate_samples, bins=bin_edges)
+    chi2_stat, chi2_pval, valid, reason = chi2_test_histograms(baseline_counts, candidate_counts)
+    result["chi2_statistic"] = chi2_stat
+    result["chi2_pvalue"] = chi2_pval
+    result["valid"] = valid
+    result["reason"] = reason
+    result["baseline_counts"] = baseline_counts.astype(int).tolist()
+    result["candidate_counts"] = candidate_counts.astype(int).tolist()
+    result["bin_edges"] = np.asarray(bin_edges, dtype=float).tolist()
+    return result
diff --git a/src/simtools/visualization/plot_event_level_production_comparison.py b/src/simtools/visualization/plot_event_level_production_comparison.py
index 7f343e4afb..25fe51e2b6 100644
--- a/src/simtools/visualization/plot_event_level_production_comparison.py
+++ b/src/simtools/visualization/plot_event_level_production_comparison.py
@@ -7,9 +7,9 @@
 import matplotlib.pyplot as plt
 import numpy as np
 
-from simtools.sim_events.production_comparison import (
-    compute_ks_chi2_for_aligned_histograms,
-    compute_ks_chi2_for_samples,
+from simtools.statistics import (
+    chi2_test_histograms,
+    compare_samples_with_statistics,
 )
 from simtools.utils import names
 
@@ -30,31 +30,12 @@ def plot(metrics_per_production, output_path, bins=40):
     """
     comparison_statistics = _initialize_comparison_statistics(metrics_per_production)
 
-    _record_plot_statistics(
-        comparison_statistics,
-        "trigger_multiplicity",
-        _plot_trigger_multiplicity(metrics_per_production, output_path),
-    )
-    _record_plot_statistics(
-        comparison_statistics,
-        "trigger_combination",
-        _plot_trigger_combinations(metrics_per_production, output_path),
-    )
-    _record_plot_statistics(
-        comparison_statistics,
-        "single_telescope_trigger_distribution",
-        _plot_single_telescope_trigger_frequencies(metrics_per_production, output_path),
-    )
-    _record_plot_statistics(
-        comparison_statistics,
-        "mixed_trigger_combinations",
-        _plot_mixed_trigger_combinations(metrics_per_production, output_path),
-    )
-    _record_plot_statistics(
-        comparison_statistics,
-        "telescope_participation_fraction",
-        _plot_telescope_participation(metrics_per_production, output_path),
-    )
+    for plot_name, plot_function in _TOP_LEVEL_STAT_PLOTS:
+        _record_plot_statistics(
+            comparison_statistics,
+            plot_name,
+            plot_function(metrics_per_production, output_path),
+        )
 
     for quantity_name, x_label, x_scale in _QUANTITY_CONFIGS:
         if quantity_name in _TRIGGERED_FRACTION_QUANTITIES:
@@ -99,7 +80,7 @@ def plot(metrics_per_production, output_path, bins=40):
             if tel_type in metrics.per_type
         ]
         for plot_fn in _PER_TYPE_PLOT_FNS:
-            statistics = plot_fn(type_metrics, output_path, suffix=f"_{tel_type}", bins=bins)
+            statistics = plot_fn(type_metrics, suffix=f"_{tel_type}", bins=bins)
             if statistics is None:
                 continue
             plot_name = _plot_name_from_statistics(statistics)
@@ -116,7 +97,7 @@ def _save_figure(fig, output_path, filename):
     plt.close(fig)
 
 
-def _plot_trigger_multiplicity(metrics_per_production, output_path, suffix="", bins=None):
+def _plot_trigger_multiplicity(metrics_per_production, suffix="", bins=None):
     """Plot triggered telescope multiplicity distributions."""
     del bins
     fig, ax = plt.subplots(figsize=(9, 6))
@@ -164,7 +145,7 @@ def _plot_trigger_multiplicity(metrics_per_production, output_path, suffix="", b
     )
     _annotate_ks_statistics(ax, statistics)
 
-    _save_figure(fig, output_path, f"trigger_multiplicity{suffix}.png")
+    # (No code here; removed misplaced block)
     statistics["plot_name"] = f"trigger_multiplicity{suffix}"
     return statistics
 
@@ -382,6 +363,12 @@ def _plot_triggered_vs_quantity(
 ):
     """Plot simulated vs triggered distributions for one quantity."""
     fig, ax = plt.subplots(figsize=(9, 6))
+    bin_edges = _get_global_quantity_bin_edges(
+        metrics_per_production,
+        quantity_name=quantity_name,
+        x_scale=x_scale,
+        bins=bins,
+    )
 
     plotted = False
     for metrics in metrics_per_production:
@@ -390,7 +377,6 @@ def _plot_triggered_vs_quantity(
             continue
         plotted = True
 
-        bin_edges = _get_bin_edges(simulated, x_scale=x_scale, bins=bins)
         sim_counts, _ = np.histogram(simulated, bins=bin_edges)
         trig_counts, _ = np.histogram(triggered, bins=bin_edges)
         with np.errstate(divide="ignore", invalid="ignore"):
@@ -624,7 +610,18 @@ def _build_aligned_count_statistics(metrics_per_production, counts_per_label, me
         candidate_counts = counts_per_label.get(metrics.label)
         if candidate_counts is None:
             continue
-        result = compute_ks_chi2_for_aligned_histograms(baseline_counts, candidate_counts)
+        chi2_stat, chi2_pval, valid, reason = chi2_test_histograms(
+            baseline_counts, candidate_counts
+        )
+        result = {
+            "chi2_statistic": chi2_stat,
+            "chi2_pvalue": chi2_pval,
+            "valid": valid,
+            "reason": reason,
+            "candidate_label": metrics.label,
+            "baseline_counts": np.asarray(baseline_counts, dtype=int).tolist(),
+            "candidate_counts": np.asarray(candidate_counts, dtype=int).tolist(),
+        }
         result["candidate_label"] = metrics.label
         result["baseline_counts"] = np.asarray(baseline_counts, dtype=int).tolist()
         result["candidate_counts"] = np.asarray(candidate_counts, dtype=int).tolist()
@@ -658,12 +655,12 @@ def _build_quantity_distribution_statistics(
             continue
         candidate_data = counts_by_label[metrics.label]
 
-        sim_stats = compute_ks_chi2_for_samples(
+        sim_stats = compare_samples_with_statistics(
             baseline_data["simulated_samples"],
             candidate_data["simulated_samples"],
             baseline_bin_edges,
         )
-        trig_stats = compute_ks_chi2_for_samples(
+        trig_stats = compare_samples_with_statistics(
             baseline_data["triggered_samples"],
             candidate_data["triggered_samples"],
             baseline_bin_edges,
@@ -895,6 +892,14 @@ def _plot_histogram_error_bars(ax, bin_edges, values, errors, color):
 _TRIGGERED_FRACTION_QUANTITIES = set()
 _SPECIAL_TRIGGER_SUBSETS = {"single_telescope", "mixed_type"}
 
+_TOP_LEVEL_STAT_PLOTS = [
+    ("trigger_multiplicity", _plot_trigger_multiplicity),
+    ("trigger_combination", _plot_trigger_combinations),
+    ("single_telescope_trigger_distribution", _plot_single_telescope_trigger_frequencies),
+    ("mixed_trigger_combinations", _plot_mixed_trigger_combinations),
+    ("telescope_participation_fraction", _plot_telescope_participation),
+]
+
 _PER_TYPE_PLOT_FNS = [
     _plot_trigger_multiplicity,
     *[
diff --git a/tests/unit_tests/applications/test_compare_productions.py b/tests/unit_tests/applications/test_compare_productions.py
index 016dc4d8f7..5e1041220b 100644
--- a/tests/unit_tests/applications/test_compare_productions.py
+++ b/tests/unit_tests/applications/test_compare_productions.py
@@ -2,6 +2,8 @@
 from types import SimpleNamespace
 from unittest.mock import MagicMock, patch
 
+import pytest
+
 import simtools.applications.compare_productions as app
 
 
@@ -42,3 +44,84 @@ def test_main_collects_metrics_and_plots(tmp_test_directory):
     mock_parse.assert_called_once_with(app_context.args["production"])
     mock_collect.assert_called_once_with(parsed_productions)
     mock_plot.assert_called_once_with(collected_metrics, output_path=output_dir)
+
+
+def test_parse_production_arguments_accepts_single_production(mocker):
+    """Test parser accepts a single production descriptor."""
+    mocker.patch(
+        "simtools.applications.compare_productions.resolve_file_patterns",
+        side_effect=lambda patterns: patterns,
+    )
+
+    descriptors = app.parse_production_arguments([["baseline", "base.h5"]])
+
+    assert len(descriptors) == 1
+    assert descriptors[0].label == "baseline"
+    assert descriptors[0].event_data_files == ["base.h5"]
+
+
+def test_parse_production_arguments_resolves_flattened_pairs(mocker):
+    """Test parser supports flattened label/file list from configuration files."""
+    mocker.patch(
+        "simtools.applications.compare_productions.resolve_file_patterns",
+        side_effect=lambda patterns: patterns,
+    )
+
+    descriptors = app.parse_production_arguments(
+        ["baseline", "base_*.h5", "candidate", "cand_*.h5"]
+    )
+
+    assert [descriptor.label for descriptor in descriptors] == ["baseline", "candidate"]
+    assert descriptors[0].event_data_files == ["base_*.h5"]
+    assert descriptors[1].event_data_files == ["cand_*.h5"]
+
+
+def test_parse_production_arguments_rejects_duplicate_labels(mocker):
+    """Test parser rejects duplicated production labels."""
+    mocker.patch(
+        "simtools.applications.compare_productions.resolve_file_patterns",
+        side_effect=lambda patterns: patterns,
+    )
+
+    with pytest.raises(ValueError, match="labels must be unique"):
+        app.parse_production_arguments([["same", "a.h5"], ["same", "b.h5"]])
+
+
+@pytest.mark.parametrize(
+    ("arguments", "error_match"),
+    [
+        ([], "At least one production is required"),
+        (["baseline", "base.h5", "dangling"], "label/file pairs"),
+        ([["baseline", "  ,   "]], "has no event_data_file pattern"),
+        ([["baseline", "a.h5"], ["candidate", 1]], "label/file pairs"),
+    ],
+)
+def test_parse_production_arguments_error_paths(mocker, arguments, error_match):
+    """Test parser validation failures for malformed production arguments."""
+    mocker.patch(
+        "simtools.applications.compare_productions.resolve_file_patterns",
+        side_effect=lambda patterns: patterns,
+    )
+
+    with pytest.raises(ValueError, match=error_match):
+        app.parse_production_arguments(arguments)
+
+
+def test_parse_production_arguments_rejects_unresolved_patterns(mocker):
+    """Test parser rejects productions that resolve to no files."""
+    mocker.patch("simtools.applications.compare_productions.resolve_file_patterns", return_value=[])
+
+    with pytest.raises(ValueError, match="does not resolve to any files"):
+        app.parse_production_arguments([["baseline", "missing_*.h5"]])
+
+
+def test_parse_production_arguments_accepts_nested_flattened_strings(mocker):
+    """Test parser accepts nested flattened string groups."""
+    mocker.patch(
+        "simtools.applications.compare_productions.resolve_file_patterns",
+        side_effect=lambda patterns: patterns,
+    )
+
+    descriptors = app.parse_production_arguments([["baseline", "a.h5", "candidate", "b.h5"]])
+
+    assert [descriptor.label for descriptor in descriptors] == ["baseline", "candidate"]
diff --git a/tests/unit_tests/sim_events/test_production_comparison.py b/tests/unit_tests/sim_events/test_production_comparison.py
index e2b1ec16b0..aeaa413d81 100644
--- a/tests/unit_tests/sim_events/test_production_comparison.py
+++ b/tests/unit_tests/sim_events/test_production_comparison.py
@@ -5,49 +5,11 @@
 
 from simtools.sim_events.production_comparison import (
     collect_production_metrics,
-    compute_ks_chi2_for_aligned_histograms,
-    compute_ks_chi2_for_samples,
-    parse_production_arguments,
 )
-
-
-def test_parse_production_arguments_accepts_single_production(mocker):
-    """Test parser accepts a single production descriptor."""
-    mocker.patch(
-        "simtools.sim_events.production_comparison.resolve_file_patterns",
-        side_effect=lambda patterns: patterns,
-    )
-
-    descriptors = parse_production_arguments([["baseline", "base.h5"]])
-
-    assert len(descriptors) == 1
-    assert descriptors[0].label == "baseline"
-    assert descriptors[0].event_data_files == ["base.h5"]
-
-
-def test_parse_production_arguments_resolves_flattened_pairs(mocker):
-    """Test parser supports flattened label/file list from configuration files."""
-    mocker.patch(
-        "simtools.sim_events.production_comparison.resolve_file_patterns",
-        side_effect=lambda patterns: patterns,
-    )
-
-    descriptors = parse_production_arguments(["baseline", "base_*.h5", "candidate", "cand_*.h5"])
-
-    assert [descriptor.label for descriptor in descriptors] == ["baseline", "candidate"]
-    assert descriptors[0].event_data_files == ["base_*.h5"]
-    assert descriptors[1].event_data_files == ["cand_*.h5"]
-
-
-def test_parse_production_arguments_rejects_duplicate_labels(mocker):
-    """Test parser rejects duplicated production labels."""
-    mocker.patch(
-        "simtools.sim_events.production_comparison.resolve_file_patterns",
-        side_effect=lambda patterns: patterns,
-    )
-
-    with pytest.raises(ValueError, match="labels must be unique"):
-        parse_production_arguments([["same", "a.h5"], ["same", "b.h5"]])
+from simtools.statistics import (
+    chi2_test_histograms,
+    compare_samples_with_statistics,
+)
 
 
 def test_collect_production_metrics_aggregates_event_data(mocker):
@@ -110,46 +72,6 @@ def test_collect_production_metrics_aggregates_event_data(mocker):
     np.testing.assert_array_equal(metrics[0].per_type["MSTN"].triggered_angular_distances, [1.2])
 
 
-@pytest.mark.parametrize(
-    ("arguments", "error_match"),
-    [
-        ([], "At least one production is required"),
-        (["baseline", "base.h5", "dangling"], "label/file pairs"),
-        ([["baseline", "  ,   "]], "has no event_data_file pattern"),
-        ([["baseline", "a.h5"], ["candidate", 1]], "label/file pairs"),
-    ],
-)
-def test_parse_production_arguments_error_paths(mocker, arguments, error_match):
-    """Test parser validation failures for malformed production arguments."""
-    mocker.patch(
-        "simtools.sim_events.production_comparison.resolve_file_patterns",
-        side_effect=lambda patterns: patterns,
-    )
-
-    with pytest.raises(ValueError, match=error_match):
-        parse_production_arguments(arguments)
-
-
-def test_parse_production_arguments_rejects_unresolved_patterns(mocker):
-    """Test parser rejects productions that resolve to no files."""
-    mocker.patch("simtools.sim_events.production_comparison.resolve_file_patterns", return_value=[])
-
-    with pytest.raises(ValueError, match="does not resolve to any files"):
-        parse_production_arguments([["baseline", "missing_*.h5"]])
-
-
-def test_parse_production_arguments_accepts_nested_flattened_strings(mocker):
-    """Test parser accepts nested flattened string groups."""
-    mocker.patch(
-        "simtools.sim_events.production_comparison.resolve_file_patterns",
-        side_effect=lambda patterns: patterns,
-    )
-
-    descriptors = parse_production_arguments([["baseline", "a.h5", "candidate", "b.h5"]])
-
-    assert [descriptor.label for descriptor in descriptors] == ["baseline", "candidate"]
-
-
 def test_collect_production_metrics_handles_empty_and_partial_datasets(mocker):
     """Test collector skips missing shower and missing trigger data sets gracefully."""
     descriptor = SimpleNamespace(label="baseline", event_data_files=["baseline_file.h5"])
@@ -199,7 +121,7 @@ def test_compute_ks_chi2_for_samples_returns_statistics():
     candidate = np.array([1.0, 2.0, 2.5, 4.0])
     bin_edges = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
 
-    result = compute_ks_chi2_for_samples(baseline, candidate, bin_edges)
+    result = compare_samples_with_statistics(baseline, candidate, bin_edges)
 
     assert result["ks_statistic"] is not None
     assert result["ks_pvalue"] is not None
@@ -215,7 +137,13 @@ def test_compute_ks_chi2_for_aligned_histograms_handles_empty_counts():
     baseline = np.array([0.0, 0.0, 0.0])
     candidate = np.array([1.0, 0.0, 0.0])
 
-    result = compute_ks_chi2_for_aligned_histograms(baseline, candidate)
+    chi2_stat, chi2_pval, valid, reason = chi2_test_histograms(baseline, candidate)
+    result = {
+        "chi2_statistic": chi2_stat,
+        "chi2_pvalue": chi2_pval,
+        "valid": valid,
+        "reason": reason,
+    }
 
     assert result["ks_statistic"] is None
     assert result["chi2_statistic"] is None
diff --git a/tests/unit_tests/visualization/test_plot_event_level_production_comparison.py b/tests/unit_tests/visualization/test_plot_event_level_production_comparison.py
index 6931f315ca..3ebd9a6de9 100644
--- a/tests/unit_tests/visualization/test_plot_event_level_production_comparison.py
+++ b/tests/unit_tests/visualization/test_plot_event_level_production_comparison.py
@@ -233,6 +233,30 @@ def test_triggered_vs_quantity_outputs_and_empty_skip(tmp_test_directory):
     assert not (output_path / "triggered_fraction_vs_core_distance_empty.png").exists()
 
 
+def test_triggered_vs_quantity_uses_global_bin_edges(tmp_test_directory, mocker):
+    """Test triggered-fraction plot derives one shared binning across productions."""
+    output_path = Path(tmp_test_directory)
+    metrics = [
+        _build_metrics("baseline", simulated_scale=1.0, triggered_scale=1.0),
+        _build_metrics("candidate", simulated_scale=2.0, triggered_scale=1.5),
+    ]
+    mock_global_bins = mocker.patch(
+        "simtools.visualization.plot_event_level_production_comparison._get_global_quantity_bin_edges",
+        return_value=np.linspace(10.0, 50.0, 9),
+    )
+
+    plot_event_level_production_comparison._plot_triggered_vs_quantity(
+        metrics,
+        output_path,
+        quantity_name="core_distance",
+        x_label="Core Distance (m)",
+        x_scale="linear",
+        bins=8,
+    )
+
+    mock_global_bins.assert_called_once()
+
+
 def test_single_and_mixed_trigger_skip_paths(tmp_test_directory):
     """Test skip branches when no single or mixed trigger combinations exist."""
     output_path = Path(tmp_test_directory)

From ca3582f664d732a024dc6c79cc191eaa8bb93445 Mon Sep 17 00:00:00 2001
From: Gernot Maier <gernot.maier@desy.de>
Date: Mon, 18 May 2026 11:34:00 +0200
Subject: [PATCH 21/27] unit tests

---
 .../plot_event_level_production_comparison.py |  4 +-
 .../sim_events/test_production_comparison.py  | 40 ---------------
 tests/unit_tests/test_statistics.py           | 49 +++++++++++++++++++
 3 files changed, 51 insertions(+), 42 deletions(-)
 create mode 100644 tests/unit_tests/test_statistics.py

diff --git a/src/simtools/visualization/plot_event_level_production_comparison.py b/src/simtools/visualization/plot_event_level_production_comparison.py
index 25fe51e2b6..a04dd281ad 100644
--- a/src/simtools/visualization/plot_event_level_production_comparison.py
+++ b/src/simtools/visualization/plot_event_level_production_comparison.py
@@ -97,7 +97,7 @@ def _save_figure(fig, output_path, filename):
     plt.close(fig)
 
 
-def _plot_trigger_multiplicity(metrics_per_production, suffix="", bins=None):
+def _plot_trigger_multiplicity(metrics_per_production, output_path, suffix="", bins=None):
     """Plot triggered telescope multiplicity distributions."""
     del bins
     fig, ax = plt.subplots(figsize=(9, 6))
@@ -145,7 +145,7 @@ def _plot_trigger_multiplicity(metrics_per_production, suffix="", bins=None):
     )
     _annotate_ks_statistics(ax, statistics)
 
-    # (No code here; removed misplaced block)
+    _save_figure(fig, output_path, f"trigger_multiplicity{suffix}.png")
     statistics["plot_name"] = f"trigger_multiplicity{suffix}"
     return statistics
 
diff --git a/tests/unit_tests/sim_events/test_production_comparison.py b/tests/unit_tests/sim_events/test_production_comparison.py
index aeaa413d81..928be224d4 100644
--- a/tests/unit_tests/sim_events/test_production_comparison.py
+++ b/tests/unit_tests/sim_events/test_production_comparison.py
@@ -6,10 +6,6 @@
 from simtools.sim_events.production_comparison import (
     collect_production_metrics,
 )
-from simtools.statistics import (
-    chi2_test_histograms,
-    compare_samples_with_statistics,
-)
 
 
 def test_collect_production_metrics_aggregates_event_data(mocker):
@@ -113,39 +109,3 @@ def test_collect_production_metrics_empty_input_keeps_arrays_empty(mocker):
     assert metrics[0].trigger_fraction == pytest.approx(0.0)
     np.testing.assert_array_equal(metrics[0].simulated_energies, np.array([]))
     np.testing.assert_array_equal(metrics[0].trigger_multiplicity, np.array([], dtype=int))
-
-
-def test_compute_ks_chi2_for_samples_returns_statistics():
-    """Test KS/Chi2 statistics for sample-based comparisons."""
-    baseline = np.array([1.0, 2.0, 3.0, 4.0])
-    candidate = np.array([1.0, 2.0, 2.5, 4.0])
-    bin_edges = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
-
-    result = compare_samples_with_statistics(baseline, candidate, bin_edges)
-
-    assert result["ks_statistic"] is not None
-    assert result["ks_pvalue"] is not None
-    assert result["chi2_statistic"] is not None
-    assert result["chi2_pvalue"] is not None
-    assert result["valid"]
-    assert result["reason"] == "ok"
-    assert result["bin_edges"] == [1.0, 2.0, 3.0, 4.0, 5.0]
-
-
-def test_compute_ks_chi2_for_aligned_histograms_handles_empty_counts():
-    """Test aligned histogram statistics return invalid for empty counts."""
-    baseline = np.array([0.0, 0.0, 0.0])
-    candidate = np.array([1.0, 0.0, 0.0])
-
-    chi2_stat, chi2_pval, valid, reason = chi2_test_histograms(baseline, candidate)
-    result = {
-        "chi2_statistic": chi2_stat,
-        "chi2_pvalue": chi2_pval,
-        "valid": valid,
-        "reason": reason,
-    }
-
-    assert result["ks_statistic"] is None
-    assert result["chi2_statistic"] is None
-    assert result["valid"] is False
-    assert result["reason"] == "insufficient_data"
diff --git a/tests/unit_tests/test_statistics.py b/tests/unit_tests/test_statistics.py
new file mode 100644
index 0000000000..1b0470f1bf
--- /dev/null
+++ b/tests/unit_tests/test_statistics.py
@@ -0,0 +1,49 @@
+import numpy as np
+
+from simtools.statistics import (
+    chi2_test_histograms,
+    compare_samples_with_statistics,
+)
+
+
+def test_compare_samples_with_statistics_returns_statistics():
+    """Test KS/Chi2 statistics for sample-based comparisons."""
+    baseline = np.array([1.0, 2.0, 3.0, 4.0])
+    candidate = np.array([1.0, 2.0, 2.5, 4.0])
+    bin_edges = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
+
+    result = compare_samples_with_statistics(baseline, candidate, bin_edges)
+
+    assert result["ks_statistic"] is not None
+    assert result["ks_pvalue"] is not None
+    assert result["chi2_statistic"] is not None
+    assert result["chi2_pvalue"] is not None
+    assert result["valid"]
+    assert result["reason"] == "ok"
+    assert result["bin_edges"] == [1.0, 2.0, 3.0, 4.0, 5.0]
+
+
+def test_chi2_test_histograms_handles_empty_counts():
+    """Test aligned histogram statistics return invalid for empty counts."""
+    baseline = np.array([0.0, 0.0, 0.0])
+    candidate = np.array([1.0, 0.0, 0.0])
+
+    chi2_stat, chi2_pval, valid, reason = chi2_test_histograms(baseline, candidate)
+
+    assert chi2_stat is None
+    assert chi2_pval is None
+    assert valid is False
+    assert reason == "insufficient_data"
+
+
+def test_chi2_test_histograms_handles_outside_support_counts():
+    """Test Chi2 invalid_expected_counts when candidate counts are outside baseline support."""
+    baseline = np.array([1.0, 0.0, 0.0])
+    candidate = np.array([0.0, 1.0, 0.0])
+
+    chi2_stat, chi2_pval, valid, reason = chi2_test_histograms(baseline, candidate)
+
+    assert chi2_stat is None
+    assert chi2_pval is None
+    assert valid is False
+    assert reason == "invalid_expected_counts"

From b1971aa1c8732c121783e27d5bde7c0dc2481cca Mon Sep 17 00:00:00 2001
From: Gernot Maier <gernot.maier@desy.de>
Date: Mon, 18 May 2026 11:37:07 +0200
Subject: [PATCH 22/27] unit tests

---
 .../plot_event_level_production_comparison.py       |  2 +-
 tests/unit_tests/test_statistics.py                 | 13 +++++++++++++
 2 files changed, 14 insertions(+), 1 deletion(-)

diff --git a/src/simtools/visualization/plot_event_level_production_comparison.py b/src/simtools/visualization/plot_event_level_production_comparison.py
index a04dd281ad..0775426a85 100644
--- a/src/simtools/visualization/plot_event_level_production_comparison.py
+++ b/src/simtools/visualization/plot_event_level_production_comparison.py
@@ -80,7 +80,7 @@ def plot(metrics_per_production, output_path, bins=40):
             if tel_type in metrics.per_type
         ]
         for plot_fn in _PER_TYPE_PLOT_FNS:
-            statistics = plot_fn(type_metrics, suffix=f"_{tel_type}", bins=bins)
+            statistics = plot_fn(type_metrics, output_path, suffix=f"_{tel_type}", bins=bins)
             if statistics is None:
                 continue
             plot_name = _plot_name_from_statistics(statistics)
diff --git a/tests/unit_tests/test_statistics.py b/tests/unit_tests/test_statistics.py
index 1b0470f1bf..b4336cd50f 100644
--- a/tests/unit_tests/test_statistics.py
+++ b/tests/unit_tests/test_statistics.py
@@ -22,6 +22,19 @@ def test_compare_samples_with_statistics_returns_statistics():
     assert result["reason"] == "ok"
     assert result["bin_edges"] == [1.0, 2.0, 3.0, 4.0, 5.0]
 
+    results_none = compare_samples_with_statistics([], [], bin_edges)
+    assert results_none["ks_statistic"] is None
+
+
+def test_ks_test_samples_handles_empty_samples():
+    """Test KS test returns None for empty samples."""
+    ks_stat, ks_pval = (
+        compare_samples_with_statistics([], [], [0.0, 1.0])["ks_statistic"],
+        compare_samples_with_statistics([], [], [0.0, 1.0])["ks_pvalue"],
+    )
+    assert ks_stat is None
+    assert ks_pval is None
+
 
 def test_chi2_test_histograms_handles_empty_counts():
     """Test aligned histogram statistics return invalid for empty counts."""

From 167c1b2d842f9622c4d62626329dc7162714f1db Mon Sep 17 00:00:00 2001
From: Gernot Maier <gernot.maier@desy.de>
Date: Mon, 18 May 2026 11:39:32 +0200
Subject: [PATCH 23/27] sonar

---
 .../plot_event_level_production_comparison.py       | 13 +++++++++++--
 1 file changed, 11 insertions(+), 2 deletions(-)

diff --git a/src/simtools/visualization/plot_event_level_production_comparison.py b/src/simtools/visualization/plot_event_level_production_comparison.py
index 0775426a85..20ee1db18d 100644
--- a/src/simtools/visualization/plot_event_level_production_comparison.py
+++ b/src/simtools/visualization/plot_event_level_production_comparison.py
@@ -37,6 +37,14 @@ def plot(metrics_per_production, output_path, bins=40):
             plot_function(metrics_per_production, output_path),
         )
 
+    _plot_quantity_comparisons(metrics_per_production, output_path, bins, comparison_statistics)
+    _plot_per_type_comparisons(metrics_per_production, output_path, bins, comparison_statistics)
+
+    _write_comparison_statistics_json(comparison_statistics, output_path)
+
+
+def _plot_quantity_comparisons(metrics_per_production, output_path, bins, comparison_statistics):
+    """Plot quantity-based comparisons and collect their statistics."""
     for quantity_name, x_label, x_scale in _QUANTITY_CONFIGS:
         if quantity_name in _TRIGGERED_FRACTION_QUANTITIES:
             _plot_triggered_vs_quantity(
@@ -65,6 +73,9 @@ def plot(metrics_per_production, output_path, bins=40):
                 ),
             )
 
+
+def _plot_per_type_comparisons(metrics_per_production, output_path, bins, comparison_statistics):
+    """Plot per-telescope-type comparisons and collect their statistics."""
     all_types = sorted(
         {
             tel_type
@@ -86,8 +97,6 @@ def plot(metrics_per_production, output_path, bins=40):
             plot_name = _plot_name_from_statistics(statistics)
             _record_plot_statistics(comparison_statistics, plot_name, statistics)
 
-    _write_comparison_statistics_json(comparison_statistics, output_path)
-
 
 def _save_figure(fig, output_path, filename):
     """Save figure and close it."""

From 8f8c29a8be56688cb0f3d45741751d3f905d0188 Mon Sep 17 00:00:00 2001
From: Gernot Maier <gernot.maier@desy.de>
Date: Mon, 18 May 2026 11:44:58 +0200
Subject: [PATCH 24/27] sonar

---
 .../visualization/plot_event_level_production_comparison.py      | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/simtools/visualization/plot_event_level_production_comparison.py b/src/simtools/visualization/plot_event_level_production_comparison.py
index 20ee1db18d..e8f254f8fa 100644
--- a/src/simtools/visualization/plot_event_level_production_comparison.py
+++ b/src/simtools/visualization/plot_event_level_production_comparison.py
@@ -411,7 +411,6 @@ def _plot_triggered_vs_quantity(
     ax.legend()
 
     _save_figure(fig, output_path, f"triggered_fraction_vs_{quantity_name}{suffix}.png")
-    return
 
 
 def _get_bin_edges(values, x_scale, bins):

From 765e76fa2da8d1553c810ea0aad9a897335d79a1 Mon Sep 17 00:00:00 2001
From: Gernot Maier <gernot.maier@desy.de>
Date: Mon, 18 May 2026 11:53:50 +0200
Subject: [PATCH 25/27] docs

---
 docs/source/api-reference/index.md      |  1 +
 docs/source/api-reference/statistics.md | 14 ++++++++++++++
 src/simtools/statistics.py              |  2 +-
 3 files changed, 16 insertions(+), 1 deletion(-)
 create mode 100644 docs/source/api-reference/statistics.md

diff --git a/docs/source/api-reference/index.md b/docs/source/api-reference/index.md
index 38691d832b..1526cfe8ec 100644
--- a/docs/source/api-reference/index.md
+++ b/docs/source/api-reference/index.md
@@ -25,6 +25,7 @@ reporting
 runners
 sim_events
 sim_telarray
+statistics
 simulators
 testing
 trigger
diff --git a/docs/source/api-reference/statistics.md b/docs/source/api-reference/statistics.md
new file mode 100644
index 0000000000..54e60c0b9b
--- /dev/null
+++ b/docs/source/api-reference/statistics.md
@@ -0,0 +1,14 @@
+(statistics)=
+
+# Statistics
+
+Statistical helper functions for comparing distributions and aligned histograms.
+
+## statistics
+
+(statisticsmodule)=
+
+```{eval-rst}
+.. automodule:: simtools.statistics
+   :members:
+```
diff --git a/src/simtools/statistics.py b/src/simtools/statistics.py
index ee7bc36af9..c3592fe262 100644
--- a/src/simtools/statistics.py
+++ b/src/simtools/statistics.py
@@ -1,4 +1,4 @@
-"""Generic statistical comparison helpers for simtools and other applications."""
+"""Generic statistical comparison helpers."""
 
 import numpy as np
 from scipy import stats

From f4d642258f6fe9e73b8c72af9d3a0e052c9991be Mon Sep 17 00:00:00 2001
From: Gernot Maier <gernot.maier@desy.de>
Date: Tue, 19 May 2026 17:04:33 +0200
Subject: [PATCH 26/27] remove chi2s; unit tests

---
 .../applications/compare_productions.py       | 84 +------------------
 .../sim_events/production_comparison.py       | 82 ++++++++++++++++++
 src/simtools/statistics.py                    | 55 +-----------
 .../plot_event_level_production_comparison.py | 37 +-------
 .../applications/test_compare_productions.py  | 27 +++---
 5 files changed, 105 insertions(+), 180 deletions(-)

diff --git a/src/simtools/applications/compare_productions.py b/src/simtools/applications/compare_productions.py
index 87292e5d39..df6517b627 100644
--- a/src/simtools/applications/compare_productions.py
+++ b/src/simtools/applications/compare_productions.py
@@ -22,10 +22,9 @@
 
 from simtools.application_control import build_application
 from simtools.sim_events.production_comparison import (
-    ProductionDescriptor,
     collect_production_metrics,
+    parse_production_arguments,
 )
-from simtools.utils.general import resolve_file_patterns
 from simtools.visualization import plot_event_level_production_comparison
 
 
@@ -52,87 +51,6 @@ def _add_arguments(parser):
     )
 
 
-def parse_production_arguments(production_arguments):
-    """Parse repeated production arguments into validated descriptors.
-
-    Parameters
-    ----------
-    production_arguments : list[list[str]]
-        Repeated ``--production`` arguments in the shape ``[label, patterns]``.
-
-    Returns
-    -------
-    list[ProductionDescriptor]
-        Validated and normalized production descriptors.
-
-    Raises
-    ------
-    ValueError
-        If configuration is malformed or does not contain any production.
-    """
-    parsed_productions = _normalize_production_arguments(production_arguments)
-    if not parsed_productions:
-        raise ValueError("At least one production is required.")
-
-    labels = [label for label, _ in parsed_productions]
-    if len(set(labels)) != len(labels):
-        raise ValueError("Production labels must be unique.")
-
-    descriptors = []
-    for label, pattern_list in parsed_productions:
-        patterns = [pattern.strip() for pattern in pattern_list.split(",") if pattern.strip()]
-        if len(patterns) == 0:
-            raise ValueError(f"Production '{label}' has no event_data_file pattern.")
-
-        resolved_files = [str(path) for path in resolve_file_patterns(patterns)]
-        if len(resolved_files) == 0:
-            raise ValueError(f"Production '{label}' does not resolve to any files.")
-        descriptors.append(ProductionDescriptor(label=label, event_data_files=resolved_files))
-
-    return descriptors
-
-
-def _normalize_production_arguments(production_arguments):
-    """Normalize raw production arguments into ``[(label, files), ...]``."""
-    if not production_arguments:
-        return []
-
-    normalized = []
-    if all(isinstance(item, str) for item in production_arguments):
-        return _pairwise_label_file_arguments(production_arguments)
-
-    for item in production_arguments:
-        normalized.extend(_normalize_single_production_argument(item))
-
-    return normalized
-
-
-def _pairwise_label_file_arguments(flat_arguments):
-    """Convert a flat list of strings into ``[(label, files), ...]`` pairs."""
-    if len(flat_arguments) % 2 != 0:
-        _raise_invalid_production_arguments()
-    return [
-        (flat_arguments[index], flat_arguments[index + 1])
-        for index in range(0, len(flat_arguments), 2)
-    ]
-
-
-def _normalize_single_production_argument(argument):
-    """Normalize one nested production argument into label/file pairs."""
-    if not isinstance(argument, list | tuple):
-        _raise_invalid_production_arguments()
-    if not all(isinstance(value, str) for value in argument):
-        _raise_invalid_production_arguments()
-    if len(argument) == 2:
-        return [(argument[0], argument[1])]
-    return _pairwise_label_file_arguments(list(argument))
-
-
-def _raise_invalid_production_arguments():
-    """Raise a standardized parser error for malformed production arguments."""
-    raise ValueError("Production arguments must be provided as label/file pairs.")
-
-
 def main():
     """See CLI description."""
     app_context = build_application(
diff --git a/src/simtools/sim_events/production_comparison.py b/src/simtools/sim_events/production_comparison.py
index 630d6fd384..866ec6f5a8 100644
--- a/src/simtools/sim_events/production_comparison.py
+++ b/src/simtools/sim_events/production_comparison.py
@@ -8,6 +8,7 @@
 
 from simtools.sim_events.reader import EventDataReader
 from simtools.utils import names
+from simtools.utils.general import resolve_file_patterns
 
 
 @dataclass
@@ -44,6 +45,87 @@ def trigger_fraction(self):
         return self.triggered_event_count / self.simulated_event_count
 
 
+def parse_production_arguments(production_arguments):
+    """Parse repeated production arguments into validated descriptors.
+
+    Parameters
+    ----------
+    production_arguments : list[list[str]]
+        Repeated ``--production`` arguments in the shape ``[label, patterns]``.
+
+    Returns
+    -------
+    list[ProductionDescriptor]
+        Validated and normalized production descriptors.
+
+    Raises
+    ------
+    ValueError
+        If configuration is malformed or does not contain any production.
+    """
+    parsed_productions = _normalize_production_arguments(production_arguments)
+    if not parsed_productions:
+        raise ValueError("At least one production is required.")
+
+    labels = [label for label, _ in parsed_productions]
+    if len(set(labels)) != len(labels):
+        raise ValueError("Production labels must be unique.")
+
+    descriptors = []
+    for label, pattern_list in parsed_productions:
+        patterns = [pattern.strip() for pattern in pattern_list.split(",") if pattern.strip()]
+        if len(patterns) == 0:
+            raise ValueError(f"Production '{label}' has no event_data_file pattern.")
+
+        resolved_files = [str(path) for path in resolve_file_patterns(patterns)]
+        if len(resolved_files) == 0:
+            raise ValueError(f"Production '{label}' does not resolve to any files.")
+        descriptors.append(ProductionDescriptor(label=label, event_data_files=resolved_files))
+
+    return descriptors
+
+
+def _normalize_production_arguments(production_arguments):
+    """Normalize raw production arguments into ``[(label, files), ...]``."""
+    if not production_arguments:
+        return []
+
+    normalized = []
+    if all(isinstance(item, str) for item in production_arguments):
+        return _pairwise_label_file_arguments(production_arguments)
+
+    for item in production_arguments:
+        normalized.extend(_normalize_single_production_argument(item))
+
+    return normalized
+
+
+def _pairwise_label_file_arguments(flat_arguments):
+    """Convert a flat list of strings into ``[(label, files), ...]`` pairs."""
+    if len(flat_arguments) % 2 != 0:
+        _raise_invalid_production_arguments()
+    return [
+        (flat_arguments[index], flat_arguments[index + 1])
+        for index in range(0, len(flat_arguments), 2)
+    ]
+
+
+def _normalize_single_production_argument(argument):
+    """Normalize one nested production argument into label/file pairs."""
+    if not isinstance(argument, list | tuple):
+        _raise_invalid_production_arguments()
+    if not all(isinstance(value, str) for value in argument):
+        _raise_invalid_production_arguments()
+    if len(argument) == 2:
+        return [(argument[0], argument[1])]
+    return _pairwise_label_file_arguments(list(argument))
+
+
+def _raise_invalid_production_arguments():
+    """Raise a standardized parser error for malformed production arguments."""
+    raise ValueError("Production arguments must be provided as label/file pairs.")
+
+
 def collect_production_metrics(production_descriptors, telescope_list=None):
     """Collect event-level comparison metrics for each production.
 
diff --git a/src/simtools/statistics.py b/src/simtools/statistics.py
index c3592fe262..99b78d8e0d 100644
--- a/src/simtools/statistics.py
+++ b/src/simtools/statistics.py
@@ -29,50 +29,8 @@ def ks_test_samples(samples1, samples2):
     return float(ks_test.statistic), float(ks_test.pvalue)
 
 
-def chi2_test_histograms(counts1, counts2):
-    """Compute the Chi2 statistic and p-value for two aligned histogram count arrays.
-
-    Parameters
-    ----------
-    counts1 : np.ndarray
-        First histogram counts (expected/baseline).
-    counts2 : np.ndarray
-        Second histogram counts (observed/candidate).
-
-    Returns
-    -------
-    chi2_statistic : float
-        Chi2 statistic value.
-    chi2_pvalue : float
-        Chi2 test p-value.
-    valid : bool
-        Whether the test is valid (enough data, no zero division).
-    reason : str
-        Reason for invalidity if not valid.
-    """
-    counts1 = np.asarray(counts1, dtype=float)
-    counts2 = np.asarray(counts2, dtype=float)
-    total1 = float(np.sum(counts1))
-    total2 = float(np.sum(counts2))
-    if total1 <= 0 or total2 <= 0:
-        return None, None, False, "insufficient_data"
-    support_mask = counts1 > 0
-    if not np.any(support_mask):
-        return None, None, False, "zero_baseline_counts"
-    observed = counts2[support_mask]
-    observed_total = float(np.sum(observed))
-    baseline_supported_total = float(np.sum(counts1[support_mask]))
-    if observed_total <= 0 or baseline_supported_total <= 0:
-        return None, None, False, "invalid_expected_counts"
-    expected = counts1[support_mask] * (observed_total / baseline_supported_total)
-    if np.sum(expected) <= 0:
-        return None, None, False, "invalid_expected_counts"
-    chi2 = stats.chisquare(f_obs=observed, f_exp=expected)
-    return float(chi2.statistic), float(chi2.pvalue), True, "ok"
-
-
 def compare_samples_with_statistics(baseline_samples, candidate_samples, bin_edges):
-    """Compute KS and Chi2 statistics for two sample arrays using provided bin edges.
+    """Compute KS statistics for two sample arrays using provided bin edges.
 
     Parameters
     ----------
@@ -86,15 +44,13 @@ def compare_samples_with_statistics(baseline_samples, candidate_samples, bin_edg
     Returns
     -------
     dict
-        KS and Chi2 statistics, p-values, and histogram metadata.
+        KS statistics, p-values, and histogram metadata.
     """
     baseline_samples = np.asarray(baseline_samples)
     candidate_samples = np.asarray(candidate_samples)
     result = {
         "ks_statistic": None,
         "ks_pvalue": None,
-        "chi2_statistic": None,
-        "chi2_pvalue": None,
         "valid": False,
         "reason": "insufficient_data",
     }
@@ -105,11 +61,8 @@ def compare_samples_with_statistics(baseline_samples, candidate_samples, bin_edg
     result["ks_pvalue"] = ks_pval
     baseline_counts, _ = np.histogram(baseline_samples, bins=bin_edges)
     candidate_counts, _ = np.histogram(candidate_samples, bins=bin_edges)
-    chi2_stat, chi2_pval, valid, reason = chi2_test_histograms(baseline_counts, candidate_counts)
-    result["chi2_statistic"] = chi2_stat
-    result["chi2_pvalue"] = chi2_pval
-    result["valid"] = valid
-    result["reason"] = reason
+    result["valid"] = True
+    result["reason"] = "ok"
     result["baseline_counts"] = baseline_counts.astype(int).tolist()
     result["candidate_counts"] = candidate_counts.astype(int).tolist()
     result["bin_edges"] = np.asarray(bin_edges, dtype=float).tolist()
diff --git a/src/simtools/visualization/plot_event_level_production_comparison.py b/src/simtools/visualization/plot_event_level_production_comparison.py
index e8f254f8fa..69044a625c 100644
--- a/src/simtools/visualization/plot_event_level_production_comparison.py
+++ b/src/simtools/visualization/plot_event_level_production_comparison.py
@@ -8,7 +8,6 @@
 import numpy as np
 
 from simtools.statistics import (
-    chi2_test_histograms,
     compare_samples_with_statistics,
 )
 from simtools.utils import names
@@ -618,21 +617,11 @@ def _build_aligned_count_statistics(metrics_per_production, counts_per_label, me
         candidate_counts = counts_per_label.get(metrics.label)
         if candidate_counts is None:
             continue
-        chi2_stat, chi2_pval, valid, reason = chi2_test_histograms(
-            baseline_counts, candidate_counts
-        )
         result = {
-            "chi2_statistic": chi2_stat,
-            "chi2_pvalue": chi2_pval,
-            "valid": valid,
-            "reason": reason,
             "candidate_label": metrics.label,
             "baseline_counts": np.asarray(baseline_counts, dtype=int).tolist(),
             "candidate_counts": np.asarray(candidate_counts, dtype=int).tolist(),
         }
-        result["candidate_label"] = metrics.label
-        result["baseline_counts"] = np.asarray(baseline_counts, dtype=int).tolist()
-        result["candidate_counts"] = np.asarray(candidate_counts, dtype=int).tolist()
         stats_summary["comparisons"].append(result)
     return stats_summary
 
@@ -685,7 +674,7 @@ def _build_quantity_distribution_statistics(
 
 
 def _annotate_ks_statistics(ax, statistics):
-    """Overlay KS and Chi2 values for baseline-vs-candidate comparisons on a plot."""
+    """Overlay KS values for baseline-vs-candidate comparisons on a plot."""
     if statistics is None or not statistics.get("comparisons"):
         return
 
@@ -696,19 +685,12 @@ def _annotate_ks_statistics(ax, statistics):
         if metric_type == "quantity_distribution":
             sim_ks = comparison.get("simulated", {}).get("ks_statistic")
             trig_ks = comparison.get("triggered", {}).get("ks_statistic")
-            sim_chi2_over_n = _format_chi2_over_n(comparison.get("simulated", {}))
-            trig_chi2_over_n = _format_chi2_over_n(comparison.get("triggered", {}))
             lines.append(
                 f"{label}: KS s/t={_format_statistic_value(sim_ks)}/"
-                f"{_format_statistic_value(trig_ks)}, "
-                f"Chi2/N s/t={sim_chi2_over_n}/{trig_chi2_over_n}"
+                f"{_format_statistic_value(trig_ks)}"
             )
             continue
-        chi2_over_n = _format_chi2_over_n(comparison)
-        lines.append(
-            f"{label}: KS={_format_statistic_value(comparison.get('ks_statistic'))}, "
-            f"Chi2/N={chi2_over_n}"
-        )
+        lines.append(f"{label}: KS={_format_statistic_value(comparison.get('ks_statistic'))}")
 
     ax.text(
         0.02,
@@ -729,19 +711,6 @@ def _format_statistic_value(value):
     return f"{float(value):.4f}"
 
 
-def _format_chi2_over_n(statistics):
-    """Format Chi2/N using bins that contribute to the Chi2 calculation."""
-    chi2_value = statistics.get("chi2_statistic")
-    baseline_counts = statistics.get("baseline_counts")
-    if chi2_value is None or baseline_counts is None:
-        return "n/a"
-
-    n_contributing_bins = int(np.sum(np.asarray(baseline_counts, dtype=float) > 0))
-    if n_contributing_bins <= 0:
-        return "n/a"
-    return f"{float(chi2_value) / n_contributing_bins:.4f}"
-
-
 def _initialize_comparison_statistics(metrics_per_production):
     """Initialize top-level comparison statistics structure."""
     baseline = metrics_per_production[0] if metrics_per_production else None
diff --git a/tests/unit_tests/applications/test_compare_productions.py b/tests/unit_tests/applications/test_compare_productions.py
index 5e1041220b..37cb9b01ad 100644
--- a/tests/unit_tests/applications/test_compare_productions.py
+++ b/tests/unit_tests/applications/test_compare_productions.py
@@ -5,6 +5,7 @@
 import pytest
 
 import simtools.applications.compare_productions as app
+import simtools.sim_events.production_comparison as production_comparison
 
 
 def test_main_collects_metrics_and_plots(tmp_test_directory):
@@ -49,11 +50,11 @@ def test_main_collects_metrics_and_plots(tmp_test_directory):
 def test_parse_production_arguments_accepts_single_production(mocker):
     """Test parser accepts a single production descriptor."""
     mocker.patch(
-        "simtools.applications.compare_productions.resolve_file_patterns",
+        "simtools.sim_events.production_comparison.resolve_file_patterns",
         side_effect=lambda patterns: patterns,
     )
 
-    descriptors = app.parse_production_arguments([["baseline", "base.h5"]])
+    descriptors = production_comparison.parse_production_arguments([["baseline", "base.h5"]])
 
     assert len(descriptors) == 1
     assert descriptors[0].label == "baseline"
@@ -63,11 +64,11 @@ def test_parse_production_arguments_accepts_single_production(mocker):
 def test_parse_production_arguments_resolves_flattened_pairs(mocker):
     """Test parser supports flattened label/file list from configuration files."""
     mocker.patch(
-        "simtools.applications.compare_productions.resolve_file_patterns",
+        "simtools.sim_events.production_comparison.resolve_file_patterns",
         side_effect=lambda patterns: patterns,
     )
 
-    descriptors = app.parse_production_arguments(
+    descriptors = production_comparison.parse_production_arguments(
         ["baseline", "base_*.h5", "candidate", "cand_*.h5"]
     )
 
@@ -79,12 +80,12 @@ def test_parse_production_arguments_resolves_flattened_pairs(mocker):
 def test_parse_production_arguments_rejects_duplicate_labels(mocker):
     """Test parser rejects duplicated production labels."""
     mocker.patch(
-        "simtools.applications.compare_productions.resolve_file_patterns",
+        "simtools.sim_events.production_comparison.resolve_file_patterns",
         side_effect=lambda patterns: patterns,
     )
 
     with pytest.raises(ValueError, match="labels must be unique"):
-        app.parse_production_arguments([["same", "a.h5"], ["same", "b.h5"]])
+        production_comparison.parse_production_arguments([["same", "a.h5"], ["same", "b.h5"]])
 
 
 @pytest.mark.parametrize(
@@ -99,29 +100,31 @@ def test_parse_production_arguments_rejects_duplicate_labels(mocker):
 def test_parse_production_arguments_error_paths(mocker, arguments, error_match):
     """Test parser validation failures for malformed production arguments."""
     mocker.patch(
-        "simtools.applications.compare_productions.resolve_file_patterns",
+        "simtools.sim_events.production_comparison.resolve_file_patterns",
         side_effect=lambda patterns: patterns,
     )
 
     with pytest.raises(ValueError, match=error_match):
-        app.parse_production_arguments(arguments)
+        production_comparison.parse_production_arguments(arguments)
 
 
 def test_parse_production_arguments_rejects_unresolved_patterns(mocker):
     """Test parser rejects productions that resolve to no files."""
-    mocker.patch("simtools.applications.compare_productions.resolve_file_patterns", return_value=[])
+    mocker.patch("simtools.sim_events.production_comparison.resolve_file_patterns", return_value=[])
 
     with pytest.raises(ValueError, match="does not resolve to any files"):
-        app.parse_production_arguments([["baseline", "missing_*.h5"]])
+        production_comparison.parse_production_arguments([["baseline", "missing_*.h5"]])
 
 
 def test_parse_production_arguments_accepts_nested_flattened_strings(mocker):
     """Test parser accepts nested flattened string groups."""
     mocker.patch(
-        "simtools.applications.compare_productions.resolve_file_patterns",
+        "simtools.sim_events.production_comparison.resolve_file_patterns",
         side_effect=lambda patterns: patterns,
     )
 
-    descriptors = app.parse_production_arguments([["baseline", "a.h5", "candidate", "b.h5"]])
+    descriptors = production_comparison.parse_production_arguments(
+        [["baseline", "a.h5", "candidate", "b.h5"]]
+    )
 
     assert [descriptor.label for descriptor in descriptors] == ["baseline", "candidate"]

From 0348eebed7a8e35c9f000e56fa2cd5495ad37296 Mon Sep 17 00:00:00 2001
From: Gernot Maier <gernot.maier@desy.de>
Date: Tue, 19 May 2026 17:31:34 +0200
Subject: [PATCH 27/27] unit tests

---
 src/simtools/statistics.py          |  2 +-
 tests/unit_tests/test_statistics.py | 35 ++---------------------------
 2 files changed, 3 insertions(+), 34 deletions(-)

diff --git a/src/simtools/statistics.py b/src/simtools/statistics.py
index 99b78d8e0d..1e78cc80a2 100644
--- a/src/simtools/statistics.py
+++ b/src/simtools/statistics.py
@@ -39,7 +39,7 @@ def compare_samples_with_statistics(baseline_samples, candidate_samples, bin_edg
     candidate_samples : np.ndarray
         Candidate sample values.
     bin_edges : np.ndarray
-        Histogram bin edges for Chi2 calculation.
+        Histogram bin edges
 
     Returns
     -------
diff --git a/tests/unit_tests/test_statistics.py b/tests/unit_tests/test_statistics.py
index b4336cd50f..a1437dc8e2 100644
--- a/tests/unit_tests/test_statistics.py
+++ b/tests/unit_tests/test_statistics.py
@@ -1,13 +1,10 @@
 import numpy as np
 
-from simtools.statistics import (
-    chi2_test_histograms,
-    compare_samples_with_statistics,
-)
+from simtools.statistics import compare_samples_with_statistics
 
 
 def test_compare_samples_with_statistics_returns_statistics():
-    """Test KS/Chi2 statistics for sample-based comparisons."""
+    """Test KS statistics for sample-based comparisons."""
     baseline = np.array([1.0, 2.0, 3.0, 4.0])
     candidate = np.array([1.0, 2.0, 2.5, 4.0])
     bin_edges = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
@@ -16,8 +13,6 @@ def test_compare_samples_with_statistics_returns_statistics():
 
     assert result["ks_statistic"] is not None
     assert result["ks_pvalue"] is not None
-    assert result["chi2_statistic"] is not None
-    assert result["chi2_pvalue"] is not None
     assert result["valid"]
     assert result["reason"] == "ok"
     assert result["bin_edges"] == [1.0, 2.0, 3.0, 4.0, 5.0]
@@ -34,29 +29,3 @@ def test_ks_test_samples_handles_empty_samples():
     )
     assert ks_stat is None
     assert ks_pval is None
-
-
-def test_chi2_test_histograms_handles_empty_counts():
-    """Test aligned histogram statistics return invalid for empty counts."""
-    baseline = np.array([0.0, 0.0, 0.0])
-    candidate = np.array([1.0, 0.0, 0.0])
-
-    chi2_stat, chi2_pval, valid, reason = chi2_test_histograms(baseline, candidate)
-
-    assert chi2_stat is None
-    assert chi2_pval is None
-    assert valid is False
-    assert reason == "insufficient_data"
-
-
-def test_chi2_test_histograms_handles_outside_support_counts():
-    """Test Chi2 invalid_expected_counts when candidate counts are outside baseline support."""
-    baseline = np.array([1.0, 0.0, 0.0])
-    candidate = np.array([0.0, 1.0, 0.0])
-
-    chi2_stat, chi2_pval, valid, reason = chi2_test_histograms(baseline, candidate)
-
-    assert chi2_stat is None
-    assert chi2_pval is None
-    assert valid is False
-    assert reason == "invalid_expected_counts"