How do you fit measured emission line data with Gaussian function in Python? (Atomic Spectroscopy)

Question

For a physics lab project, I am measuring various emission lines from various elements. High intensity peaks occur at certain wavelengths. My goal is to fit a Gaussian function in python in order to find at which wavelength the intensity is peaking.

I have already tried using the norm function from the scipy.stats library. Below is the code and the graph that is produced.

import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt

mean, std = norm.fit(he3888_1[:,0])
plt.plot(he3888_1[:,0], he3888_1[:,1], color='r')
x = np.linspace(min(he3888_1[:,0]), max(he3888_1[:,0]), 100)
y = norm.pdf(x, mean, std)
plt.plot(x, y)
plt.xlabel("Wavelength (Angstroms)")
plt.ylabel("Intensity")
plt.show()

Could this be because the intensity is low for a relatively long period prior to it?

Answer 1

Lmfit seems like a good option for your case. The code below simulates a Gaussian peak with a linear background added and shows how you can extract the parameters with lmfit. The latter has a number of other built-in models (Lorentzian, Voight, etc.) that can be easily combined with each other.

import numpy as np
from lmfit.models import Model, LinearModel
from lmfit.models import GaussianModel, LorentzianModel
import matplotlib.pyplot as plt

def generate_gaussian(amp, mu, sigma_sq, slope=0, const=0):
    x = np.linspace(mu-10*sigma_sq, mu+10*sigma_sq, num=200)
    y_gauss = (amp/np.sqrt(2*np.pi*sigma_sq))*np.exp(-0.5*(x-mu)**2/sigma_sq)
    y_linear = slope*x + const
    y = y_gauss + y_linear
    return x, y

# Gaussiand peak generation
amplitude = 6
center = 3884
variance = 4
slope = 0
intercept = 0.05
x, y = generate_gaussian(amplitude, center, variance, slope, intercept)


#Create a lmfit model: Gaussian peak + linear background
gaussian = GaussianModel()
background = LinearModel()
model = gaussian + background

#Find what model parameters you need to specify
print('parameter names: {}'.format(model.param_names))
print('independent variables: {}'.format(model.independent_vars))

#Model fit
result = model.fit(y, x=x, amplitude=3, center=3880,
                   sigma=3, slope=0, intercept=0.1)
y_fit = result.best_fit #the simulated intensity

result.best_values #the extracted peak parameters

# Comparison of the fitted spectrum with the original one
plt.plot(x, y, label='model spectrum')
plt.plot(x, y_fit, label='fitted spectrum')
plt.xlabel('wavelength, (angstroms')
plt.ylabel('intensity')
plt.legend()

Output:

parameter names: ['amplitude', 'center', 'sigma', 'slope', 'intercept']

independent variables: ['x']

result.best_values
Out[139]: 
{'slope': 2.261379140543626e-13,
 'intercept': 0.04999999912168238,
 'amplitude': 6.000000000000174,
 'center': 3883.9999999999977,
 'sigma': 2.0000000000013993}