The anisotropic ion-ion correlation function g(R₁, R₂), where R₁ and R₂ are the coordinates of ions 1 and 2, is calculated for double layer systems composed of counterions between two charged, planar walls. The primitive model of the electrolyte is used. The Widom formula for the excess chemical potential is applied in the neighbourhood of a confined particle in the double layer to obtain g(R₁, R₂) from Monte Carlo (MC) simulations. The pair correlations are also calculated in the anisotropic RHNC approximation, where the reference hypernetted chain (RHNC) closure is applied at the pair (ion-ion) level and where the anisotropic bridge function is taken from an inhomogeneous reference system composed of a hard sphere fluid with the same density profile as the double-layer system. The results of the MC simulations and RHNC calculations are generally in very good agreement with each other. This demonstrates the fact that the anisotropic RHNC approximation is a very accurate method for calculating the various double-layer properties. It also follows that the Widom technique is a very versatile tool in simulation studies of double-layer systems. To obtain 15 points for the pair correlation function by simulation takes about 60 times longer than to obtain 200 000 points by anisotropic RHNC calculations (the latter takes a few minutes on an IBM RS-6000/590 workstation).