Fitting the peaks of gamma spectra is a typical task in non-destructive analysis of spent fuel or neutron activation analysis. Since these applications are already "quite old", there is some standard software available, like Genie 2000. A paper Evaluation of Peak-Fitting Software for Gamma Spectrum Analysis from 2015 compares a number of these tools. However, most of them are commercial.
All methods are based on the concept of "base functions", typically Gaussian functions, that are used in a linear combination. The two main problems are then to estimate the number of peaks and to estimate their locations. Both are done using prior knowledge: what does the user think is present in the sample? Once this is done, a non-linear optimisation algorithm is used (Levenberg-Marquardt) to minimize the error between the fit and the measured data. But, as mentioned by EMP in a comment, a lot of this is in the eye of the beholder. Most commercial software have a GUI that allows the user to interact with the code to "tinkle until convergence" (where convergence means "it fits according to my human eye and according to what radioactive species I think I have in my material").
A paper describing a non-linear fitting procedure using the Levenberg-Marquardt method can be found here. This will give you some of the mathematical background. This paper "Peak fitting and identification software library for high resolution gamma-ray spectra" also gives some mathematical background. An older report describing this problem can be found here. Based on the latter and the available solvers in SciPy, you should be able to roll your own code.