Gaussian Fit to Data in Python using Lmfit

Question

I'm trying to write a code that performs a Gaussian fit to a gamma ray calibration spectrum, i.e. multiple peaks. Here is my current code:

from numpy import loadtxt
import numpy as np


from lmfit.models import GaussianModel
import matplotlib.pyplot as plt



#Centre of each of the peaks we want to fit Gaussian to
peakCentroids = np.array([251, 398, 803, 908, 996, 1133, 1178, 2581, 3194, 3698, 4671])


#Import total data set
data = loadtxt('Phase_1_02.dat')
x = data[:, 0] #channel number/gamma energy
y = data[:, 1] #Count number

mod = GaussianModel()

pars = mod.guess(y, x=x)
out = mod.fit(y, pars, x=x)

print(out.fit_report(min_correl=0.25))



for currentPeakNumber in range(len(peakCentroids)):
     fig = plt.figure(figsize = (8, 8))
    
plt.plot(x, y, 'b')
plt.plot(x, out.init_fit, 'k--', label='initial fit')
plt.plot(x, out.best_fit, 'r-', label='best fit')
plt.legend(loc='best')
plt.show()

It's outputting the spectra for my data and printing the relevant parameters (e.g. center, sigma, fwhm etc.), but I'm having a bit of brain freeze in terms of fitting the Gaussian peak to each of the centroids specified! Currently the output spectra is only fitting to the first peak at a value of about 248. Is there anyone much better at coding in Python than me that can shed some light on how to go about this and if it's possible using Lmfit please? Thanks in advance!! :)

take a look here it gives you the possibility to fit with Gaussian (or the Lorentzian, or the mix of two), while sett ing the initial parameters visually. — shcrela

M Newville M Newville · Accepted Answer · 2021-02-27T02:39:32

If I understand the question correctly, you are looking to model the data you have with a series of Gaussian line shapes, centered at the many (10 or more) values you have.

If that is the case, the model should be constructed from the 10 or more Gaussians, but your model only builds one Gaussian. You'll want to build a model with something like

import numpy as np
from lmfit.models import GaussianModel

peakCentroids = np.array([251.0, 398, 803, 908, 996, 1133, 1178, 2581, 3194,
                          3698, 4671.0])

mod = None 

for i, cen in enumerate(peakCentroids):
     thispeak = GaussianModel(prefix='p%d_' %(i+1)) 
     if mod is None:  
          mod = thispeak
     else:
          mod = mod + thispeak

pars = mod.make_params()

for i, cen in enumerate(peakCentroids):
    pars['p%d_center' % (i+1)].value = cen
    pars['p%d_amplitude' % (i+1)].value = 1.0 # is that a sensible initial value?
    pars['p%d_sigma' % (i+1)].value = 1.0 # is that a sensible initial value?

out = mod.fit(y, pars, x=x)
print(out.fit_report(min_correl=0.25))

or at least that might be a good place to start...

Gaussian Fit to Data in Python using Lmfit

1 Answers