Factorizing the cross section for single hadron production in $e^+e^-$ annihilations, differential in $z_h$, $P_T$ and thrust, is a highly non trivial task. We have devised a factorization scheme that allows us to recast the $e^+e^- to hX$ cross section in the convolution of a perturbatively calculable coefficient and a universal Transverse Momentum Dependent (TMD) Fragmentation Function (FF). The predictions obtained from our NLO-NLL perturbative computation, together with a simple ansatz to model the non-perturbative part of the TMD, are applied to the experimental measurements of the BELLE Collaboration for the phenomenological extraction of this process independent TMD FF.