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

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.

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%)