We present the theoretical development of the 3D multipole probability tomography applied to the electric Self-Potential (SP) method of geophysical exploration. We assume that an SP dataset can be thought of as the response of an aggregation of poles, dipoles, quadrupoles and octopoles. These physical sources are used to reconstruct, without a priori assumptions, the most probable position and shape of the true SP buried sources, by determining the location of their centres and critical points of their boundaries, as corners, wedges and vertices. At first, a few synthetic cases with cubic bodies are examined in order to determine the resolution power of the new technique. Then, an experimental SP dataset collected in the Mt. Somma-Vesuvius volcanic district (Naples, Italy) is elaborated in order to define location and shape of the sources of two SP anomalies of opposite sign detected in the northwestern sector of the surveyed area. The modelled sources are interpreted as the polarization state induced by an intense hydrothermal convective flow mechanism within the volcanic apparatus, from the free surface down to about 3 km of depth b.s.l.