We study from a theoretical point of view the Backus and Gilbert method in the case of Newtonian potential. If a mass distribution $m$ on a open set $\Omega$ creates a Newtonian potential $U^m$, which is known on an infinity of points $(y_n)_{n\in\mathbb N}$ out of $\overline{\Omega}$, we characterize the solution $m_0$, obtained as a generalization of the Backus and Gilbert method, as the projection of $m$ (for the scalar product of $L_2(\Omega)$) on a subspace of harmonic functions; this subspace may be the subspace of all harmonic, square-integrable functions (for example, if $\Omega$ is a starlike domain). Then we study the reproducing kernel $B$ associated to this projection, which satisfies $$ m_0(x)=\int\limits_{\Omega}B(x,y)m(y)\,dy $$ for any $m\in L_2(\Omega)$.