The ground state correlations induced by a general pairing Hamiltonian in a finite system of like fermions are described in terms of four-body correlated structures (quartets). These are real superpositions of products of two pairs of particles in time-reversed states. Quartets are determined variationally through an iterative sequence of diagonalizations of the Hamiltonian in restricted model spaces and are, in principle, all distinct from one another. The ground state is represented as a product of quartets to which, depending on the number of particles (supposed to be even, in any case), an extra collective pair is added. The extra pair is also determined variationally. In case of pairing in a spherically symmetric mean field, both the quartets and the extra pair (if any) are characterized by a total angular momentum J=0. Realistic applications of the quartet formalism are carried out for the Sn isotopes with the valence neutrons in the 50-82 neutron shell. Exact ground state correlation energies, occupation numbers and pair transfer matrix elements are reproduced to a very high degree of precision. The formalism also lends itself to a straightforward and accurate description of the lowest seniority 0 and 2 excited states of the pairing Hamiltonian. A simplified representation of the ground state as a product of identical quartets is eventually discussed and found to improve considerably upon the more traditional particle-number projected-BCS approach.