The NA61/SHINE collaboration has presented experimental data on $${{\rho }^{0}}$$, $$\omega $$ and $${{K}^{{0 * }}}$$ meson productions in $${{\pi }^{ - }}{\text{C}}$$ interactions at $$Plab = 158$$ and 350 GeV/c. As shown by the collaboration, the existing Monte Carlo event generators badly describe the data. In this paper a possibility to describe the data is considered in the Quark-Gluon String (QGS) model of Geant4. QGS model was proposed by A.B. Kaidalov and K.A. Ter-Martirosyan. Within the model, A.B. Kaidalov found a behaviour of quark and diquark fragmentation functions for $$z \to 0$$ and $$z \to 1$$, and proposed interpolation formulae for the functions in the whole region of $$z$$. These functions must be a solution of the well-known system of the integral equations. A simplified Monte Carlo estimation of the functions, based on usage of the fragmentation functions at $$z \to 1$$ as the kernel functions of the system, does not reproduce Kaidalov’s results. An improvement of the Monte Carlo simulations is proposed in this paper. It can be implemented in Monte Carlo event generators such as Los Alamos QGSM, QGSJet-II and the Geant4 QGS model. It will improve a description of experimental data in the models, especially, the description of the latest NA61/SHINE Collaboration data. Description of the data is a problem in DPMJet, QGSJet, EPOS and Sibyll models.