Method of polarized semi-inclusive deep inelastic scattering (SIDIS) data analysis in the next to leading order (NLO) QCD is developed. Within the method one first directly extracts in NLO few first truncated (available to measurement) Mellin moments of the quark helicity distributions. Second, using these moments as an input to the proposed modification of the Jacobi polynomial expansion method (MJEM), one eventually reconstructs the local quark helicity distributions themselves. All numerical tests demonstrate that MJEM allows us to reproduce with the high precision the input local distributions even inside the narrow Bjorken $x$ region accessible for experiment. It is of importance that only four first input moments are sufficient to achieve a good quality of reconstruction. The application of the method to the simulated SIDIS data on the pion production is considered. The obtained results encourage one that the proposed NLO method can be successfully applied to the SIDIS data analysis. The analysis of HERMES data on pion production is performed. To this end the pion difference asymmetries are constructed from the measured by HERMES standard semi-inclusive spin asymmetries. The LO results of the valence distribution reconstruction are in a good accordance with the respective leading order SMC and HERMES results, while the NLO results are in agreement with the existing NLO parametrizations on these quantities.