Custom peak function in MPF2 - what happened to the fit wave?

Share icon

I am using Multi Peak Fitting 2 to treat some x-ray photoemission spectroscopy (XPS) data. The model lineshape for the peaks is the so-called Doniac-Sunjic function, an asymmetric Lorenzian-like peak given analytically by

DS(a, G, A, x0) = A*( cos( Pi*a/2 + (1-a)*arctan((x-x0)/G) ) / (G^2 + (x-x0)^2)^(1/2-a/2)
Here A is the amplitude, G is the width, x0 is the centroid, and a is the asymmetry parameter. With a=0 the function simplifies to a Lorentzian, and with a between 0 and 1 it is asymmetric towards lower values of x. Typical fit values of a should be 0.1-0.3.

I believe I have successfully programmed this as an MPF2 fit function, and I can get the peaks to show up in the subcomponent axis. The fit appears to complete correctly, giving me reasonable-looking peaks. However, the fitted curve is not present on the graph - the original data is there, and the residuals look as if no fit is present (direct copy of the data, so it's subtracting zero from the data), see the attached image.

If I switch the peak type to Gaussian, it works correctly. The fit function reappears, and the residuals are calculated properly as well...

Have I messed up the peak definition? Why is it drawing the individual peaks, but not the fit function? Is it the arctangent function screwing me up?
It would be actually good to have your code to look at to get an idea what might be wrong.
The DS peak shape should be generally no problem for MPF2 as even non-analytic peaks are possible. Have you actually thought about peak broadening of your experiment (?) via Gaussian convolution or is this not relevant? Also, you want to be extra careful of what to give back as fit results, since especially the area in this case is useless.
Anyway, it would make things easier if you could provide your code or an example experiment (attached to you entry).
Maybe it would also be a good idea to work together to create a generalized DS peak snippet for the community to use as it is a relatively common peak shape.

January 26, 2012 at 02:23 am - Permalink

Of course the broadening is a concern - I had not thought much about how to do the convolution, I first just wanted to see if the fitting of multiple DS peaks could be handled on a relatively modern machine without too much slowdown.

I embarked upon this because I have a bunch of synchrotron XPS to fit. In the past I actually had a nice XOP compiled by a group where I did a postdoc, and this XOP included a gaussian convolution. Unfortunately that ran under windows only, and I am loath to switch operating systems back and forth just for running this one particular fitting. I wanted to see how fast (or slow) it would run without an XOP (I don't have the code, nor the XOP toolkit so I can't replicate this).

Here's the definitions for MPF2 - see if you get the same disappearing fit wave that I do?

Function/S DS_PeakFuncInfo(InfoDesired)
    Variable InfoDesired

    String info=""
       
    Switch (InfoDesired)
        case PeakFuncInfo_ParamNames:
            info = "Location;width;asymmetry;Height;"
            break;
        case PeakFuncInfo_PeakFName:
            info = "MPFDS"
            break;
        case PeakFuncInfo_GaussConvFName:
            info = "GausstoDSGuess"
            break;
        case PeakFuncInfo_ParameterFunc:
            info = "DSParams"
            break;
        case PeakFuncInfo_DerivedParamNames:
            info = "Location;Height;FWHM;Asymmetry;"
            break;
        default:
            break;
    endSwitch
   
    return info
end

Function MPFDS(w,yw,xw)
  wave w
  wave yw, xw

  yw = w[3]*cos(Pi*w[2]/2+(1-w[2])*atan((x-w[0])/w[1]))/((x-w[0])^2+w[1]^2)^((1-w[2])/2)
end

Function GaussToDSGuess(w)
  wave w
  variable height=w[2]
  Redimension/N=4 w
  w[2]=0.05 //just put in a small asymmetry to start the fit with.
  w[3]=height
  return 0
end

Function DSParams(cw, sw, outWave)
    Wave cw, sw, outWave
    Variable width = cw[1]
    Variable widthSigma = sqrt(sw[1][1])
    Variable asym=cw[2]
    Variable asymsigma=sqrt(sw[2][2])
    Variable amp = cw[3]
    Variable ampSigma = sqrt(sw[3][3])
   
    // location
    outWave[0][0] = cw[0]
    outWave[0][1] = sqrt(sw[0][0])
       
    // amplitude
    outWave[1][0] = amp
    outWave[1][1] = ampSigma
   
    // area    ..... haven't what to return here, the integral of DS diverges at large values of E, so the area is a bit meaningless.
   
    outWave[2][0] = nan
    outWave[2][1] = nan
   
    // FWHM
    outWave[3][0] = width*2*sqrt(ln(2))
    outWave[3][1] = widthSigma*2*sqrt(ln(2))
end

January 26, 2012 at 06:48 am - Permalink

If you can't figure out the problem, send a copy of your experiment file with the data to support@wavemetrics.com and I'll take a look at it. Also, make sure you're using the very latest version of Igor (and therefore of Multipeak Fit). I recall some bugs that look kind of like what you show here.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

January 26, 2012 at 09:26 am - Permalink

So John found my error instantly, it turns out that I defined the function correctly:
Function MPFDS(w,yw,xw)

But when I returned the result I used 'x' instead of 'xw'!
yw = w[3]*cos(Pi*w[2]/2+(1-w[2])*atan((x-w[0])/w[1]))/((x-w[0])^2+w[1]^2)^((1-w[2])/2) //this is wrong
yw = w[3]*cos(Pi*w[2]/2+(1-w[2])*atan((xw-w[0])/w[1]))/((xw-w[0])^2+w[1]^2)^((1-w[2])/2) //this is right!!

Works a champ now, but no broadening yet.

January 26, 2012 at 12:54 pm - Permalink

I sent Josh an e-mail in response to his message to support.

For the record here, the function MPFDS() needs to have "x" replaced with "xw" so that X values come from the input xw wave. One of those easy-to-make, hard-to-see mistakes.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

January 26, 2012 at 12:55 pm - Permalink

Here is my first version with added Gaussian broadening to the DS peak. I have messed a bit with some naming and the order of parameters in your function to make thinking easier for me. Sorry for that.
I used basically the approach described in ‘All-At-Once Fitting Functions’ in the manual with some modifications. That is, generating to waves with the peak and the Gaussian, convolving them, and then let the fit function pick the results.
There is now an additional GaussWidth parameter which determines the Gaussian broadening (set to zero for no broadening). Inside the fit function itself the approach is now to generate a wave which is more extended and fine-grained than the incoming peak-wave (controlled by a resolution factor inside the function). The DS-peak is written into this wave which can then be used for convolution. Looks like this:
...
Variable Resolution = 10                 // the sample resolution (higher = slower)
Variable Xdelta     = abs(xw[0] - xw[NumPnts(xw)-1]) / (Resolution * NumPnts(xw))    // calculate an high-res xdelta for scaling
Variable yPoints    = Resolution * NumPnts(yw)   // make the wave X times bigger than the original y wave
Make/FREE/D/N=(2*yPoints+1) yWave            // create a wave two times as big to have space before and after the peak
SetScale/P x xw[0]-yPoints/2*Xdelta, Xdelta, yWave   // reserve 1/2 size of space before the x wave starts (1/2 after)

yWave = w[2]*cos(Pi*w[3]/2+(1-w[3])*atan((x-w[0])/w[1]))/((x-w[0])^2+w[1]^2)^((1-w[3])/2)

if (w[4] != 0)                       // convolve with gaussian if width is set
    Variable width  = w[4]/Xdelta            // width normalized by x scaling
    yPoints         = 10*width       // create enough points to fit in the full gaussian
    Make/FREE/D/N=(yPoints + 1) GaussWave        // create a gaussian for the convolution
    GaussWave = Gauss(x,yPoints/2,width)
    Convolve/A GaussWave yWave           // convolve with the prepared gaussian
endif
yw = yWave(xw[p])                    // write the calculated values into the resultwave
...

There, your first version of the equation with the x (instead of xw) comes into play again. ;) The wave has to be more extended to avoid funny effects of the convolution at the edge in the final product. The waves created here might be overkill for most cases, but I don’t have any performance problems even with 5 of these peaks. Suggestions for optimizations are always welcome.
With the convolution in effect the calculated values for height, FWHM, and so on are rather meaningless. I use a function to calculate all these values numerically. I don’t have a clue how to get fit errors for these values, so I leave them blank. Maybe someone has a good idea for this.
I have also modified the guess function a bit to work better with automatic guesses. Don’t know if there are cases which let the guess fail this way (enormous asymmetry for example). Finally, the fit starts to look funny if the DS width gets close to zero, so one needs to look out for this case.
I have attached the full peak procedure. Please report any issues you run into.
DS Peak.ipf (5.59 KB)

January 27, 2012 at 10:36 am - Permalink

chozo wrote:
Here is my first version with added Gaussian broadening to the DS peak.


This works quite well for me! A little slower than non-convolved lineshapes, but quite useable.

I'm beginning to realize that what would be most useful would be 'global' asymmetry and gaussian width values - for example, I don't the gausswidth or asymmetry to be different when decomposing spectral data from a single peak (or doublet). I wonder if there is some hackery to be done by applying a constraint wave post hoc to fix these to be identically valued for all DS peaks, even though they are nominally being varied?

I was also thinking we could add a "DS_doublet" lineshape that produces two DS lineshapes with a preset energy splitting and intensity ratio - this would make fitting of closely spaced (say, Si 2p) XPS peaks easier than having a separate DS for each SO-split component and trying to fix the known energy splitting values. I'm sure this can be trivially implemented with chozo's code, I just don't have the time to do it myself right now - perhaps next week, but if anyone is feeling enthusiastic about working on this...

February 2, 2012 at 12:21 pm - Permalink

Great that it works for you. I’m not sure if there are still errors in the calculation, so you may want to double-check results before you go over to a daily use. And maybe someone comes up with solution for the missing errors.

As mentioned, this thing is not so optimized. I opted for a relatively high resolution (ten times the original points) to represent even extremely thin peaks and lots of space (double the initial size) to fit in the convolution. I’m having a hard time spying on the doings of MPF2 with the debugger, but usually there is a lot of unused space in the waves while going through the various steps. I can’t estimate the edge cases very well this way. I sure there are ways to adjust these things to the actual requirements of the particular fit. Also there are still some quirks of MPF2 which increase the difficulty here; like that the individual peak waves of the fit are not representing the correct size of the peaks (i.e. may be much too wide or too narrow).

As for the doublet peak, global asymmetry and so on, this would be possible with linked parameters (like global fit has for multiple sets). Maybe with some magic this can be implemented into MPF2 in the future. I still have it at the end of my crazy feature suggestion list. For now you could look out manually for the correct parameters while pushing the peaks around in the fit. Another lazy hack could be to synchronize parameters via a global variable or the Kx ones at runtime. Don’t know what this would break.

February 3, 2012 at 07:08 am - Permalink

I think it would be simple enough to write a peak shape function that implemented a doublet. The part that wouldn't work would be the automatic peak finder; I think it will be difficult to cull the extra peaks represented by the doublets (unless the individual peaks of the doublets are strongly overlapping). You might have to manually delete extra peaks.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

February 3, 2012 at 09:02 am - Permalink