  Factorizing the cross section for single hadron production in $e^+e^-$
annihilations is a highly non trivial task when the transverse momentum of the
outgoing hadron with respect to the thrust axis is taken into account. We work
in a scheme that allows to factorize the $e^+e^- \to HX$ cross section as a
convolution of a calculable hard coefficient and a Transverse Momentum
Dependent (TMD) fragmentation function. The result, differential in $z_h$,
$P_T$ and thrust, will be given to all orders in perturbation theory and
explicitly computed to Next to Leading Order (NLO) and Next to Leading Log
(NLL) accuracy. The predictions obtained from our computation, applying the
simplest and most natural ansatz to model the non-perturbative part of the TMD,
are in exceptional agreement with the experimental measurements of the BELLE
Collaboration. The factorization scheme we propose relates the TMD parton
densities defined in 1-hadron and 2-hadron processes, restoring the possibility
to perform global phenomenological studies of TMD physics including
experimental data from semi-inclusive deep inelastic scattering, Drell-Yan
processes, $e^+e^- \to H_1 H_2 X$ and $e^+e^- \to HX$ annihilations.
