diff --git a/PQAnalysis/analysis/rdf/rdf.py b/PQAnalysis/analysis/rdf/rdf.py
index ff8e37d1..dd6c3e1e 100644
--- a/PQAnalysis/analysis/rdf/rdf.py
+++ b/PQAnalysis/analysis/rdf/rdf.py
@@ -587,8 +587,8 @@ def _setup_legacy_bins(self, delta_r: PositiveReal):
         """Sets bins with the scalar semantics of ``thh_tools/RDF``."""
         self.r_min = 0.0
         self.delta_r = float(np.float32(delta_r))
-        self.n_bins = max(
-            int(max(cell.box_lengths) / 2.0 / self.delta_r)
+        self.n_bins = min(
+            int(min(cell.box_lengths) / 2.0 / self.delta_r)
             for cell in self.cells
         )
         self.r_max = self.delta_r * self.n_bins
diff --git a/tests/analysis/rdf/test_parity.py b/tests/analysis/rdf/test_parity.py
index 227777b2..b52802ed 100644
--- a/tests/analysis/rdf/test_parity.py
+++ b/tests/analysis/rdf/test_parity.py
@@ -9,6 +9,7 @@
 
 import hashlib
 import importlib
+import math
 from pathlib import Path
 
 import numpy as np
@@ -129,6 +130,61 @@ def test_float64_coordinate_controls_boundary_bin(
     assert np.array_equal(analysis.bins, np.array([1.0, 0.0, 0.0, 0.0]))
 
 
+@pytest.mark.parametrize(("frame_kernel", "batch_kernel"), KERNELS)
+def test_noncubic_box_bins_span_only_the_kernel_cutoff(
+    tmp_path,
+    monkeypatch,
+    frame_kernel,
+    batch_kernel,
+):
+    monkeypatch.setattr(
+        rdf_module, "legacy_rdf_frame_histogram", frame_kernel
+    )
+    monkeypatch.setattr(
+        rdf_module, "legacy_rdf_batch_histogram", batch_kernel
+    )
+
+    trajectory = tmp_path / "noncubic.xyz"
+    trajectory.write_text(
+        "2 20.0 20.0 60.0 90.0 90.0 90.0\n\n"
+        "X 0.0 0.0 0.0\n"
+        "Y 3.0 0.0 0.0\n",
+        encoding="utf-8",
+    )
+
+    analysis = RDF(
+        TrajectoryReader(str(trajectory)),
+        "X",
+        "Y",
+        delta_r=0.1,
+    )
+    result = analysis.run()
+
+    # The kernel drops every pair beyond min(20, 20, 60) / 2 = 10.0, so
+    # the histogram must cover int(10.0 / float32(0.1)) = 99 bins. Sizing
+    # it from the largest edge (299 bins) would normalize over a tail the
+    # kernel never counts.
+    delta_r = float(np.float32(0.1))
+
+    assert analysis._legacy_rdf  # pylint: disable=protected-access
+    assert analysis.n_bins == 99
+    assert analysis.r_max == 99 * delta_r
+
+    # The single pair sits at distance 3.0 -> bin floor(3.0 / delta_r) = 29.
+    expected_bins = np.zeros(99)
+    expected_bins[29] = 1.0
+    assert np.array_equal(analysis.bins, expected_bins)
+
+    # Hand-computed: one pair, one reference, one frame, number density
+    # 1 / (20 * 20 * 60), ideal-gas shell between 29 and 30 bin edges.
+    density = 1.0 / (20.0 * 20.0 * 60.0)
+    shell = 4.0 / 3.0 * math.pi * (
+        (30.0 * delta_r) ** 3 - (29.0 * delta_r) ** 3
+    )
+    assert result[1][29] == pytest.approx(1.0 / (density * shell), rel=1e-12)
+    assert result[2][-1] == 1.0
+
+
 @pytest.mark.parametrize("example_dir", ["msd"], indirect=False)
 def test_api_output_round_trips_exact_values(test_with_data_dir):
     Path("rdf.in").write_text(