For positive integers n and c with c >= 2n (= the Coxeter number of the hyperoctahedral group of signed permutations of degree n), we present a finite-dimensional discrete orthogonality relation for Macdonald's three-parameter hyperoctahedral Hall-Littlewood polynomials of degree at most c in each of the n variables. These polynomials are labeled by partitions lambda that fit inside a rectangular box of shape c(n), i.e. the partitions in question are of length <= n, and have parts of size <= c. We employ coordinate patches around the vertices of the alcove of boxed partitions Lambda((n, c)) = {lambda subset of c(n)} to establish the self- adjointness of a one-parameter family of commuting discrete difference operators T-x((n, c)) acting on functions f : Lambda((n, c)) -> C. By construction, the basis of hyperoctahedral Hall-Littlewood polynomials constitutes a joint eigenbasis for T-X((n, c)) with simple spectrum, which gives rise to our discrete orthogonality relation. From the point of view of quantum integrable particle dynamics, the present geometric construction establishes the orthogonality of the Bethe Ansatz eigenfunctions for a recently studied q-boson system on a finite lattice with integrable open-end boundary conditions. The Bethe Ansatz equations enter in the geometric picture as compatibility conditions between coordinate patches stemming from distinct vertices. Two applications of the orthogonality relations are highlighted: (i) a cubature rule for the integration of symmetric functions with respect to the Haar measure on the compact symplectic group Sp(n), and (ii) a Verlinde formula for the structure constants of a deformation of the Wess-Zumino-Witten fusion ring associated with the affine Lie algebra of type (C) over cap (n) (= C-n((1))). (C) 2021 Elsevier Inc. All rights reserved.