A new formalism for the factorization of the cross section for single hadron production in e + e − annihilations, differential in z h , PT and thrust, is applied to the phenomenological analysis of data recently measured by the BELLE Collaboration. Within this scheme the e + e − → hX cross section can be recast in the convolution of a perturbatively calculable coefficient and a universal transverse momentum dependent fragmentation function. While performing a next-to-leading order calculation of the perturbative part of the process to next-to-leading logarithmic accuracy, we examine and thoroughly discuss the suitability of a number of possible ansatz to model the non-perturbative part of this universal transverse momentum dependent fragmentation function, showing the extent to which present experimental data can actually constrain its shape and functional form in terms of z h , PT and thrust.