For multi-pinhole small animal SPECT imaging, system matrix (SM) is critical for accurate reconstruction. It is typically achieved by analytical calculation, experimental measurement, or Monte-Carlo (MC) simulation. MC simulation method is relatively accurate and avoids long-time measurements. However, to achieve high resolution for small animal reconstruction, SM should be extremely huge, and would have severe noise produced when simulation time is limited. In this paper, SM was factorized into two parts: geometrical projection center matrix and corresponding point spread function (PSF) matrix. PSFs were derived from MC simulation. Three models including circular Gaussian, standard Gaussian, and skew Gaussian functions were used to fit a MC simulated PSF. A series of detailed discussions and comparisons for typical PSF cases with above three models fitting were presented to explain and validate why standard Gaussian model is good for PSF fitting and accurate reconstruction, especially in oblique incidence cases, and propose the cases that circular and skew Gaussian is suitable. Simulated 10-point-source and the hot rod phantom experiments in an experimental ultra-high resolution small animal SPECT system demonstrated this method could achieve accurate image reconstruction quality with high compression ratio. Standard Gaussian model was verified to effectively reduce statistical noise and improve the accuracy of PSF even if the simulation time was not long enough, to greatly compress the SM storage, and be able to increase the sampling density of PSF. Thus, accurate image reconstruction for multi-pinhole small animal SPECT can be achieved based on standard Gaussian model fitted SM. This method could also be applied to SM measured from the point source experiments.