diff --git a/image_analysis/geom_refinement/BeamCenterFromBackground.cpp b/image_analysis/geom_refinement/BeamCenterFromBackground.cpp
index 9129589bc..b2c306b9e 100644
--- a/image_analysis/geom_refinement/BeamCenterFromBackground.cpp
+++ b/image_analysis/geom_refinement/BeamCenterFromBackground.cpp
@@ -127,6 +127,12 @@ FindBeamCenterFromBackground(const DiffractionExperiment &experiment, const Pixe
     std::vector<double> block_jy(static_cast<size_t>(BLOCKS) * n_cells);
     std::vector<int32_t> block_count(static_cast<size_t>(BLOCKS) * n_cells);
 
+    // Most pixels lie outside the band. Those clearly outside it in tan(2theta) = rho / lz - by a
+    // margin far above float rounding - skip the square root and both atan2 below; every pixel the exact
+    // test would keep still reaches it.
+    const double tan_lo = std::tan(static_cast<double>(tt_lo)) * (1.0 - 1e-3);
+    const double tan_hi = tt_hi < PI / 2 ? std::tan(static_cast<double>(tt_hi)) * (1.0 + 1e-3) : INFINITY;
+
     float step_x = 0.0f, step_y = 0.0f, sigma_x = 0.0f, sigma_y = 0.0f;
     float previous_x = 0.0f, previous_y = 0.0f;
     int reversals = 0;
@@ -153,7 +159,13 @@ FindBeamCenterFromBackground(const DiffractionExperiment &experiment, const Pixe
                     const float lx = rot[0] * u + rot[1] * v + rot[2] * distance;
                     const float ly = rot[3] * u + rot[4] * v + rot[5] * distance;
                     const float lz = rot[6] * u + rot[7] * v + rot[8] * distance;
-                    const float rho = std::sqrt(lx * lx + ly * ly);
+                    const float rho_sq = lx * lx + ly * ly;
+                    if (lz > 0.0f) {
+                        const double lz_sq = static_cast<double>(lz) * lz;
+                        if (rho_sq < tan_lo * tan_lo * lz_sq || rho_sq > tan_hi * tan_hi * lz_sq)
+                            continue;
+                    }
+                    const float rho = std::sqrt(rho_sq);
                     const float two_theta = std::atan2(rho, lz);
                     if (two_theta < tt_lo || two_theta >= tt_hi || rho == 0.0f)
                         continue;