0

I have been trying to fit a gaussian into my spectrum. (intensity v/s velocity spectrum)

spectrum

New fitted spectrum

I used the following code to fit the data to a gaussian profile. However, as seen in the result only one data point was included in the fit. Is there anything I can do so that i can include more points into the gaussian.

from numpy import exp, linspace, random
import matplotlib.pyplot as plt
def gaussian(x, amp, cen, wid):
    return amp * exp(-(x-cen)**2 / wid)

from scipy.optimize import curve_fit
x = velocity
y = data
print(x)
print(y)

init_vals = [0.00950554, 60000, 35]  # for [amp, cen, wid]
best_vals, covar = curve_fit(gaussian, x, y, p0=init_vals)
print(best_vals)
print(repr(covar))

ym = gaussian(x, best_vals[0], best_vals[1], best_vals[2])
fig = plt.figure()
ax = fig.add_subplot(211)
ax.plot(x,y)
ax = fig.add_subplot(212)
ax.plot(x, y, c='k', label='Function')
ax.plot(x, ym, c='r', label='Best fit')
plt.legend()
plt.show()

Data points:

x: [-5.99999993e+04 -4.99999993e+04 -3.99999993e+04 -2.99999993e+04
 -1.99999993e+04 -9.99999934e+03  6.65010004e-04  1.00000007e+04
  2.00000007e+04  3.00000007e+04  4.00000007e+04  5.00000007e+04
  6.00000007e+04  7.00000007e+04  8.00000007e+04  9.00000007e+04
  1.00000001e+05  1.10000001e+05  1.20000001e+05  1.30000001e+05
  1.40000001e+05]

y: [ 0.00056511 -0.00098584 -0.00325616 -0.00101042  0.00168894 -0.00097406
 -0.00134408  0.00128847 -0.00111633 -0.00151621  0.00299326  0.00916455
  0.00960554  0.00317363  0.00311124 -0.00080881  0.00215932  0.00596419
 -0.00192256 -0.00190138 -0.00013216]

These are the data points of the spectrum. Can anybody help me get a better gaussian fit into the data. I have been trying everything I can.

Thanks a lot.

Lidia
  • 27
  • 5

1 Answers1

1

The first step is always to plot the data, which you already did. The next is to guess initial values. Those for amp and cen look reasonably, if compared to the plot. But what about wid? It is 2 times the width of the distribution SQUARED. From the plot, the width itself must be of odrer of several thousand. If squared, it may reach 10^7, times 2 gives 2*10^7 as the initial value. Very far from your 35!

One possible solution:

amp = 0.0106283
cen = 55784
wid = 1.92911e+08

Plot:

Plot of itting function

zkoza
  • 2,644
  • 3
  • 16
  • 24
  • Thank you very much. Could you please tell me what changes I should make in the code to obtain the above plot? I tried but couldn't succeed in replicating the graph. Thank you again for your help. – Lidia Apr 16 '21 at 18:23
  • Change the initial value of `wid` to a value similar to the one I gave as a possible solution. For example, `1e8`. – zkoza Apr 16 '21 at 21:01
  • `init_vals = [0.00950554, 60000, 1e8]` works like a charm. – zkoza Apr 16 '21 at 21:02
  • Thanks again. I tried it and the obtained result is attached as 'New fitted spectrum'. Is there any way I can get a curve as smooth as the one you obtained? – Lidia Apr 16 '21 at 21:45
  • I don't know. I don't use python :-) – zkoza Apr 16 '21 at 21:53
  • Define another "x" vector, say, `xx`, but with far more elements. For example, at least 100 points with values 0...100000, step 1000. Use it as the argument to `ym = gaussian(xx, best_vals[0], best_vals[1], best_vals[2])`. Then use both `ym` and `xx` in a plot, `ax.plot(xx, ym, c='r', label='Best fit')` ` – zkoza Apr 16 '21 at 21:59
  • Okay. I'll keep trying again. I was instructed to measure the full width at half maximum line width when I fit the data to a gaussian. So, will that change anything in the obtained gaussian fit? Thanks again for the help. I have been working on this problem for weeks. – Lidia Apr 16 '21 at 22:02
  • No, it won't. One thing that you should have been given is the uncertainties of the values. It is clear they're pretty large. Without it it's hard to estimate uncertainties of the fit parameters. Also, the data suggest that they may be described by a superposition of 2 Gaussians, but that would need to have some theory justyfying such a hypothesis. – zkoza Apr 16 '21 at 22:41
  • If I want to find the RMS noise of the spectrum can I proceed by taking the standard deviation of the spectrum? – Lidia Apr 17 '21 at 16:26
  • Could you please address the question and let me know if i am in the right direction? https://stackoverflow.com/questions/67137643/finding-rms-noise-in-a-spectra – Lidia Apr 17 '21 at 16:27