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.
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/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/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/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/docs/source/user-guide/applications.md b/docs/source/user-guide/applications.md
index 76c277c7de..b8d38cacba 100644
--- a/docs/source/user-guide/applications.md
+++ b/docs/source/user-guide/applications.md
@@ -49,6 +49,7 @@ 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>
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 6bd86626a6..d253530718 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 = "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.py b/src/simtools/applications/compare_productions.py
new file mode 100644
index 0000000000..df6517b627
--- /dev/null
+++ b/src/simtools/applications/compare_productions.py
@@ -0,0 +1,75 @@
+#!/usr/bin/python3
+
+r"""Compare simulation productions at different comparison levels.
+
+The application accepts repeated production descriptors and dispatches to
+comparison-level specific implementations.
+
+Command line arguments
+----------------------
+production (repeated, required)
+    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.
+"""
+
+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 as needed for multiple productions."
+        ),
+    )
+    parser.add_argument(
+        "--comparison_level",
+        choices=["events", "signals", "compute"],
+        default="events",
+        help="Comparison level to execute.",
+    )
+
+
+def main():
+    """See CLI description."""
+    app_context = build_application(
+        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":
+        metrics_per_production = collect_production_metrics(production_descriptors)
+        plot_event_level_production_comparison.plot(
+            metrics_per_production,
+            output_path=output_directory,
+        )
+    else:
+        raise NotImplementedError(f"Comparison level '{comparison_level}' is not implemented yet.")
+
+
+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..866ec6f5a8
--- /dev/null
+++ b/src/simtools/sim_events/production_comparison.py
@@ -0,0 +1,374 @@
+"""Utilities for event-level comparison across multiple simulation productions."""
+
+import logging
+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
+
+
+@dataclass
+class ProductionDescriptor:
+    """Descriptor for one production input provided."""
+
+    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
+    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):
+        """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 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.
+
+    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__)
+
+    metric_accumulators = _initialize_metric_accumulators()
+    for event_data_file in production_descriptor.event_data_files:
+        _collect_metrics_from_event_data_file(event_data_file, telescope_list, metric_accumulators)
+
+    logger.info(
+        f"Collected production {production_descriptor.label}: "
+        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=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.
+    """
+    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.
+
+    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 = {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:
+        subset_counts["mixed_type"] = len(telescopes)
+    return subset_counts
+
+
+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:
+        return np.array([], dtype=dtype)
+    return np.concatenate(arrays)
diff --git a/src/simtools/statistics.py b/src/simtools/statistics.py
new file mode 100644
index 0000000000..1e78cc80a2
--- /dev/null
+++ b/src/simtools/statistics.py
@@ -0,0 +1,69 @@
+"""Generic statistical comparison helpers."""
+
+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 compare_samples_with_statistics(baseline_samples, candidate_samples, bin_edges):
+    """Compute KS 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
+
+    Returns
+    -------
+    dict
+        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,
+        "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)
+    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()
+    return result
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..69044a625c
--- /dev/null
+++ b/src/simtools/visualization/plot_event_level_production_comparison.py
@@ -0,0 +1,895 @@
+"""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.statistics import (
+    compare_samples_with_statistics,
+)
+from simtools.utils import names
+
+_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.
+    """
+    comparison_statistics = _initialize_comparison_statistics(metrics_per_production)
+
+    for plot_name, plot_function in _TOP_LEVEL_STAT_PLOTS:
+        _record_plot_statistics(
+            comparison_statistics,
+            plot_name,
+            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(
+                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):
+            if cumulative is None:
+                continue
+            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,
+                ),
+            )
+
+
+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
+            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:
+            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)
+
+
+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_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
+    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 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(
+            ax,
+            bin_edges,
+            fractions,
+            errors,
+            color=_artist_color(stairs_artist),
+        )
+
+    ax.set_xlabel("Triggered Telescopes per Event")
+    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(bin_edges[:-1])
+    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):
+    """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 None
+
+    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]
+
+    statistics = _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",
+    )
+    if statistics is not None:
+        statistics["plot_name"] = "trigger_combination"
+    return statistics
+
+
+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 None
+
+    statistics = _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",
+    )
+    if statistics is not None:
+        statistics["plot_name"] = "single_telescope_trigger_distribution"
+    return statistics
+
+
+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 None
+
+    statistics = _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",
+    )
+    if statistics is not None:
+        statistics["plot_name"] = "mixed_trigger_combinations"
+    return statistics
+
+
+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=(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(
+            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(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()
+
+    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):
+    """Return True only for mixed signatures 1+1 and 1+2 (max multiplicity 3)."""
+    type_counts, _ = _type_counts_from_combination(combination)
+    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> + ...'."""
+    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:
+        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
+    return type_counts, telescope_names
+
+
+def _plot_triggered_vs_quantity(
+    metrics_per_production,
+    output_path,
+    quantity_name,
+    x_label,
+    x_scale,
+    bins,
+    suffix="",
+):
+    """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:
+        simulated, triggered = _get_quantity_arrays(metrics, quantity_name)
+        if simulated.size == 0:
+            continue
+        plotted = True
+
+        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,
+            )
+        _plot_series(ax, bin_edges, efficiency, metrics.label, quantity_name, linewidth=2)
+
+    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")
+    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()
+
+    _save_figure(fig, output_path, f"triggered_fraction_vs_{quantity_name}{suffix}.png")
+
+
+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 _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))
+    bin_edges = _get_global_quantity_bin_edges(
+        metrics_per_production,
+        quantity_name=quantity_name,
+        x_scale=x_scale,
+        bins=bins,
+    )
+    plotted = False
+    counts_by_label = {}
+    for metrics in metrics_per_production:
+        simulated, triggered = _get_quantity_arrays(metrics, quantity_name)
+        if simulated.size == 0:
+            continue
+        plotted = True
+
+        sim_counts, _ = np.histogram(simulated, bins=bin_edges)
+        trig_counts, _ = np.histogram(triggered, bins=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,
+            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 None
+    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)
+
+    statistics = _build_quantity_distribution_statistics(
+        metrics_per_production,
+        quantity_name,
+        counts_by_label,
+        baseline_bin_edges=bin_edges,
+        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 _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(
+        {
+            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 None
+
+    statistics = _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",
+    )
+    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 = {
+            "candidate_label": metrics.label,
+            "baseline_counts": np.asarray(baseline_counts, dtype=int).tolist(),
+            "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,
+    baseline_bin_edges,
+    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]
+
+    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 = compare_samples_with_statistics(
+            baseline_data["simulated_samples"],
+            candidate_data["simulated_samples"],
+            baseline_bin_edges,
+        )
+        trig_stats = compare_samples_with_statistics(
+            baseline_data["triggered_samples"],
+            candidate_data["triggered_samples"],
+            baseline_bin_edges,
+        )
+
+        comparison_record = {
+            "candidate_label": metrics.label,
+            "simulated": sim_stats,
+            "triggered": trig_stats,
+            "baseline_bin_edges": np.asarray(baseline_bin_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 s/t={_format_statistic_value(sim_ks)}/"
+                f"{_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):
+    """Return two cumulative slots, using None to disable one for specific quantities."""
+    if quantity_name == "energy":
+        return (False, None)
+    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, _metrics=None):
+    """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"}
+
+_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,
+    *[
+        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)
+        if cum is not None
+    ],
+]
diff --git a/tests/integration_tests/config/compare_productions_run.yml b/tests/integration_tests/config/compare_productions_run.yml
new file mode 100644
index 0000000000..5d5e805862
--- /dev/null
+++ b/tests/integration_tests/config/compare_productions_run.yml
@@ -0,0 +1,18 @@
+---
+applications:
+- application: simtools-compare-productions
+  configuration:
+    comparison_level: events
+    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_multiplicity.png
+  - output_file: trigger_combination.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.py b/tests/unit_tests/applications/test_compare_productions.py
new file mode 100644
index 0000000000..37cb9b01ad
--- /dev/null
+++ b/tests/unit_tests/applications/test_compare_productions.py
@@ -0,0 +1,130 @@
+from pathlib import Path
+from types import SimpleNamespace
+from unittest.mock import MagicMock, patch
+
+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):
+    """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"]],
+            "comparison_level": "events",
+        },
+        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.build_application",
+            return_value=app_context,
+        ),
+        patch(
+            "simtools.applications.compare_productions.parse_production_arguments",
+            return_value=parsed_productions,
+        ) as mock_parse,
+        patch(
+            "simtools.applications.compare_productions.collect_production_metrics",
+            return_value=collected_metrics,
+        ) as mock_collect,
+        patch(
+            "simtools.applications.compare_productions.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)
+    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.sim_events.production_comparison.resolve_file_patterns",
+        side_effect=lambda patterns: patterns,
+    )
+
+    descriptors = production_comparison.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 = production_comparison.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"):
+        production_comparison.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.sim_events.production_comparison.resolve_file_patterns",
+        side_effect=lambda patterns: patterns,
+    )
+
+    with pytest.raises(ValueError, match=error_match):
+        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.sim_events.production_comparison.resolve_file_patterns", return_value=[])
+
+    with pytest.raises(ValueError, match="does not resolve to any files"):
+        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.sim_events.production_comparison.resolve_file_patterns",
+        side_effect=lambda patterns: patterns,
+    )
+
+    descriptors = production_comparison.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
new file mode 100644
index 0000000000..928be224d4
--- /dev/null
+++ b/tests/unit_tests/sim_events/test_production_comparison.py
@@ -0,0 +1,111 @@
+from types import SimpleNamespace
+
+import numpy as np
+import pytest
+
+from simtools.sim_events.production_comparison import (
+    collect_production_metrics,
+)
+
+
+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]),
+        angular_distance=np.array([0.5, 1.0, 1.5]),
+    )
+    triggered_shower_data = SimpleNamespace(
+        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"]),
+        ]
+    )
+
+    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 == 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, 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].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(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, 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, [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.2])
+
+
+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 == pytest.approx(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 == 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))
diff --git a/tests/unit_tests/test_statistics.py b/tests/unit_tests/test_statistics.py
new file mode 100644
index 0000000000..a1437dc8e2
--- /dev/null
+++ b/tests/unit_tests/test_statistics.py
@@ -0,0 +1,31 @@
+import numpy as np
+
+from simtools.statistics import compare_samples_with_statistics
+
+
+def test_compare_samples_with_statistics_returns_statistics():
+    """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])
+
+    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["valid"]
+    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
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..3ebd9a6de9
--- /dev/null
+++ b/tests/unit_tests/visualization/test_plot_event_level_production_comparison.py
@@ -0,0 +1,295 @@
+import json
+from collections import Counter
+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
+
+
+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]),
+        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}),
+        simulated_event_count=4,
+        triggered_event_count=2,
+    )
+
+
+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)
+    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 = [
+        "comparison_statistics.json",
+        "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",
+    ]
+    _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."""
+    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",
+    ]
+    _assert_files_exist(output_path, per_type_files)
+
+
+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"
+    )
+
+
+@pytest.mark.parametrize(
+    ("quantity_name", "expected"),
+    [
+        ("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."""
+    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)
+
+    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_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)
+    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()
+
+
+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
+    )