Fitting two voigt curves, one after the other using lmfit

Question

I have the following emission spectra of Neon collected on a Raman (background subtracted data):

    x=np.array([[1114.120887, 1114.682293, 1115.243641, 1115.80493 , 1116.366161, 1116.927334, 1117.488449, 1118.049505, 1118.610503, 1119.171443, 1119.732324, 1120.293147, 1120.853912, 1121.414619, 1121.975267, 1122.535857, 1123.096389, 1123.656863, 1124.217278, 1124.777635, 1125.337934, 1125.898175, 1126.458357, 1127.018482, 1127.578548, 1128.138556, 1128.698505, 1129.258397, 1129.81823 , 1130.378005, 1130.937722, 1131.497381, 1132.056981]])

    y=np.array([[-4.89046878e+00, -4.90985832e+00, -5.92924587e+00, -3.28194437e+00, -1.96801488e+00, -3.32070938e+00, -5.34008887e+00, -3.59466330e-01, -2.04552879e+00, -1.06490224e+00,  8.24910035e+00,  5.32297309e+01, 1.11543677e+02,  8.98576241e+01,  2.18504948e+02,  7.15152212e+02, 7.62799601e+02,  2.89446870e+02,  7.24275144e+01,  1.94081610e+01, 1.72212272e+00,  7.02773412e-01, -3.16573861e-01,  4.99745483e+00, 7.97811157e+00,  6.25396305e-01,  6.27274408e+00, -4.41328018e+00, -7.76592840e+00,  3.88142539e+00,  6.52872017e+00,  1.50939096e+00, -8.43249208e-01]])

I have fitted a single Voigt function using lmfit, specifically:

    model = VoigtModel()+ ConstantModel()
    params=model.make_params(center=1123.096389, amplitude=1000, sigma=0.27)
    result = model.fit(y.flatten(), params, x=x.flatten())

There is a second peak on the LH shoulder (sorry can't post image)- people using commercial peak fitting software fit the first voigt, then add the second, and then it adjusts the fits of both. How would I do this in python?

A related question - is there a way to optimize how many points to include in the peak fit. Right now, I am only feeding x and y data covering a set spectral range to do the peak fitting. But commercial software optimizes how much range to include in a given peak fit (I presume using residuals). How would I recreate this?

Thanks!

Answer 1

You can do it manually as so:

import numpy as np
import matplotlib.pyplot as plt
from lmfit.models import VoigtModel, ConstantModel

x=np.array([1114.120887, 1114.682293, 1115.243641, 1115.80493 , 1116.366161, 1116.927334, 1117.488449, 1118.049505, 1118.610503, 1119.171443, 1119.732324, 1120.293147, 1120.853912, 1121.414619, 1121.975267, 1122.535857, 1123.096389, 1123.656863, 1124.217278, 1124.777635, 1125.337934, 1125.898175, 1126.458357, 1127.018482, 1127.578548, 1128.138556, 1128.698505, 1129.258397, 1129.81823 , 1130.378005, 1130.937722, 1131.497381, 1132.056981])
y=np.array([-4.89046878e+00, -4.90985832e+00, -5.92924587e+00, -3.28194437e+00, -1.96801488e+00, -3.32070938e+00, -5.34008887e+00, -3.59466330e-01, -2.04552879e+00, -1.06490224e+00,  8.24910035e+00,  5.32297309e+01, 1.11543677e+02,  8.98576241e+01,  2.18504948e+02,  7.15152212e+02, 7.62799601e+02,  2.89446870e+02,  7.24275144e+01,  1.94081610e+01, 1.72212272e+00,  7.02773412e-01, -3.16573861e-01,  4.99745483e+00, 7.97811157e+00,  6.25396305e-01,  6.27274408e+00, -4.41328018e+00, -7.76592840e+00,  3.88142539e+00,  6.52872017e+00,  1.50939096e+00, -8.43249208e-01])

model = VoigtModel() + ConstantModel()
params=model.make_params(center=1123.0, amplitude=1000, sigma=0.27)
result1 = model.fit(y.flatten(), params, x=x.flatten())

rest = y-result1.best_fit

model = VoigtModel() + ConstantModel()
params=model.make_params(center=1120.5, amplitude=200, sigma=0.27)
result2 = model.fit(rest, params, x=x.flatten())

rest -= result2.best_fit

plt.plot(x, y, label='Original')
plt.plot(x, result1.best_fit, label='1123.0')
plt.plot(x, result2.best_fit, label='1120.5')
plt.plot(x, rest, label='residual')
plt.legend()
plt.show()

You have to make sure that the residual makes sense. In this case, is quite close to 0, so I'd argue that it is fine.

lmfit does optimize the fit, so it is not necessary to pinpoint the exact value of the peak position. Also, it is important to point out that because of the resolution of this data (and spectroscopy in general), the highest points are not necessarily the centre of the peak. Additionally, because of the same, some shoulders might not be shoulders, though in this case looks like it is.

For your related question - judging by the documentation of lmfit it uses all the range you input. Residuals seem like not a solution since you fall in the same problem (what range to consider). I believe that the commercial SW you mention uses Multivariate Curve Resolution (MCR). These deconvolution problems have been a hot topic for decades. If you are interested in this kind of solution, I suggest reading about Multivariate Curve Resolution (MCR).