A Monte Carlo code was written to determine the saturation thickness for multiply scattered gamma rays from aluminium targets. Interactions of incident gamma rays with the energies of 123, 279, 360, 511, 662, 1115, and 1250 keV were simulated. This work aims to design a convenient code which can be used in investigations on gamma backscattering. Obtained results for saturation thickness values have been compared with experimental ones and the Monte Carlo N-particle (MCNP) code results, and showed good agreement. Also, based on the similar behavior of number of multiple scattered photons between these three methods, the expected spectrum of singly or multiply scattered photons which is not possible to observe with experiment has been presented.