We present a comprehensive analysis and extraction of the unpolarized transverse momentum dependent (TMD) parton distribution functions, which are fundamental constituents of the TMD factorization theorem. We provide a general review of the theory of TMD distributions, and present a new scheme of scale fixation. This scheme, called the zeta-prescription, allows to minimize the impact of perturbative logarithms in a large range of scales and does not generate undesired power corrections. Within zeta-prescription we consistently include the perturbatively calculable parts up to next-to-next-to-leading order (NNLO), and perform the global fit of the Drell-Yan and Z-boson production, which include the data of E288, Tevatron and LHC experiments. The non-perturbative part of the TMDs are explored checking a variety of models. We support the obtained results by a study of theoretical uncertainties, perturbative convergence, and a dedicated study of the range of applicability of the TMD factorization theorem. The considered non-perturbative models present significant differences in the fitting behavior, which allow us to clearly disfavor most of them. The numerical evaluations are provided by the arTeMiDe code, which is introduced in this work and that can be used for current/future TMD phenomenology.