In positron emission tomography (PET) imaging, statistical iterative reconstruction (IR) techniques appear particularly promising since they can provide accurate system model.The system model matrix which describes the relationship between image space and projection space is important to the image quality. It contains some factors such as geometrical component and blurring component. The blurring component is usually described by point spread function (PSF). A PSF matrix derived from the single photon incidence response function is studied. And then an IR method based on the system matrix containing the PSF is developed. More specifically, the gamma photon incidence on a crystal array is simulated by Monte Carlo (MC) simulation, and then the single photon incidence response functions are calculated. Subsequently, the single photon incidence response functions is used to compute the coincidence blurring factor according to the physical process of PET coincidence detection. Through weighting the ordinary system matrix response by the coincidence blurring factors, the IR system matrix containing PSF is finally established. Using this system matrix, the image is reconstructed by ordered subset expectation maximization (OSEM) algorithm. The experimental results show that the proposed system matrix can obviously improve the image radial resolution, contrast and noise property. Furthermore, the simulated single gamma-ray incidence response function only depends on the crystal configuration, so the method could be extended to any PET scanners with the same detector crystal configuration.