Gnuplot: fitting and plotting peak with shoulder, but individual peaks need to be plotted too

Question

I would like to plot the subpeaks of my fitted peak. I have a composite peak with a main peak and a shoulder. I can fit and plot it, but i need the individual subpeaks which build up the spectra. How can I plot them too? My code is the following:

set term png
set fit quiet
set fit logfile "peakfitdatalog.txt"
do for[i = 0:20] for [j = 0:20] {
outfile = "testmap_Y".(sprintf('%02d_X%02d',j,i)).".png"
set output outfile

set xrange [900:1000]

x01=960
w1=10.64
a1=50000
mu1=0.90

x02=967
w2=10.64
a2=500
mu2=0.90

y0=350
y1=0.1

voigtfv(x)=y0+y1*x+a1*((2./pi)*w1/(4.*(x-x01)**2.+w1**2.))+a2*((2./pi)*w2/(4.*(x-x02)**2.+w1**2.))

fit voigtfv(x) "testmap_Y".(sprintf('%02d_X%02d',j,i)).".txt" via a1, a2, w1, w2, y0, y1

plot "testmap_Y".(sprintf('%02d_X%02d',j,i)).".txt", voigtfv(x)
set print "fit_param.dat" append
print sprintf("%i,%i,%i,%i",a1,a2,w1,w2)
set print
}

I have tried to plot each peaks alone after the plotting of the main fit and plot procedure but is is not working, I dont see the subpeaks only the result of the fitting procedure.

Answer 1

A bit late..., but maybe this is still helpful to you. Although, I don't understand how exactly your function is related to the Voigt profile (maybe an approximation?) and why it has some linear term y1*x. As I understand the Voigt-profile is a convolution of a Lorentz and a Gaussian distribution.

Actually, gnuplot has some Voigt-function implemented: check help voigt. As much as I understand voigt(x,0) is a Gauss profile and voigt(x,1) is a Lorentz profile.

Some comments:

• your data uses comma not point as decimal sign. Therefore, I had to use, e.g. set locale "German" and set decimalsign locale. Maybe there are other ways.
• it is important to get reasonable starting values for fitting. You might get errors: Undefined value during function evaluation. You also might find a way to get automated estimations, e.g. peak positions x1,x2 are easy, but b1, b2 and p might need more thoughts.
• for fitting the voigt(x,p) profile, parameter p has been kept the same for the two peaks.

Script:

### fit two voigt profiles
reset session

FILE = "SO77274181.dat"
set locale "German"
set decimalsign locale

v(x,x0,A,b,p) = A*voigt(b*(x-x0),p)
f1(x) = v(x,x1,A1,b1,p)
f2(x) = v(x,x2,A2,b2,p)
f(x)  = f1(x) + f2(x)
# initial values for fitting
x1=960
A1=10000
x2=970
A2=3000
b1 = 1.5
b2 = 0.5
p  = 0.5
fit f(x) FILE via x1,A1,b1,x2,A2,b2,p

set xrange [900:1000]
set samples 500

plot FILE u 1:2 w lp pt 7 ps 0.5, \
     f(x)  lc "red" lw 2, \
     f1(x) lc "web-green", \
     f2(x) lc "blue"
### end of script

Result:

Final set of parameters            Asymptotic Standard Error
=======================            ==========================
x1              = 961.328          +/- 0.08639      (0.008986%)
A1              = 12112.2          +/- 302.9        (2.501%)
b1              = 0.15852          +/- 0.003813     (2.406%)
x2              = 971.214          +/- 0.1677       (0.01727%)
A2              = 4018.12          +/- 173.4        (4.316%)
b2              = 0.246915         +/- 0.01145      (4.637%)
p               = 0.258716         +/- 0.02553      (9.867%)

Answer 2

First of all, 'w1' appears in the denominator of the fourth term of the voigtfv function is probably a typo for 'w2'.

General purpose non-linear fitting like Gnuplot's fit command may not converge to the solution you want, depending on the initial value. You can read the notes on initial values by typing 'help fit starting_values' on the gnuplot command line.

As long as your sample data, it worked well when the initial value of 'a2' was set to the same magnitude as 'a1'.

set xrange [900:1000]
set offset 0, 0, graph 0.3, 0

x01=960
w1=10.64
a1=50000

x02=967
w2=10.64
a2=50000  ### changed

y0=350
y1=0.1

voigtfv(x)=y0+y1*x\
          +a1*((2./pi)*w1/(4.*(x-x01)**2.+w1**2.))\
          +a2*((2./pi)*w2/(4.*(x-x02)**2.+w2**2.))

fit voigtfv(x) "testmap_Y.txt" via a1, a2, w1, w2, y0, y1

plot "testmap_Y.txt", voigtfv(x), \
                      a1*((2./pi)*w1/(4.*(x-x01)**2.+w1**2.)),\
                      a2*((2./pi)*w2/(4.*(x-x02)**2.+w2**2.))