The factorization theorems for transverse momentum distributions of dilepton/boson production, recently formulated by Collins and Echevarria-Idilbi-Scimemi in terms of well-defined transverse momentum dependent distributions (TMDs), allows for a systematic and quantitative analysis of non-perturbative QCD effects of the cross sections involving these quantities. In this paper we perform a global fit using all current available data for Drell-Yan and Z-boson production at hadron colliders within this framework. The perturbative calculable pieces of our estimates are included using a complete resummation at next-to-next-to-leading-logarithmic accuracy. Performing the matching of transverse momentum distributions onto the standard collinear parton distribution functions and recalling that the corresponding matching coefficient can be partially exponentiated, we find that this exponentiated part is spin-independent and resummable. We argue that the inclusion of higher order perturbative pieces is necessary when data from lower energy scales are analyzed. We consider non-perturbative corrections both to the intrinsic nucleon structure and to the evolution kernel and find that the non-perturbative part of the TMDs could be parametrized in terms of a minimal set of parameters (namely 2-3). When all corrections are included the global fit so performed gives a χ2/d.o.f. ≲ 1 and a very precise prediction for vector boson production at the Large Hadron Collider (LHC).