A fit function for QENS data.

A generalised fit function for QENS data, fitting N lorentzians and a delta function, convolving with the instrumental resolution function. It wouldn't be hard to swap out lorentzians for gauss functions, etc.

Function fitUsingNLor(cw, yw, xw1) : FitFunc
    Wave cw, yw, xw1
    //An all at once fit function to fit QENS data.
    //function fits a linear background
    // with a delta function
    // and N lorentzians.
    // delta and lorentzians are centred at the same position (cw[2])
    // xw1 contains the x values
    // res contains the instrumental resolution function at the same x values.
    // the returned y wave consists of the sum of the delta and N lorentzian functions
    // convolved with the resolution function

        //the resolution function MUST be contained in the following wave
    Wave res = root:res

        //1) THAT THE POINT SPACING IN XW1 IS LINEAR!!!!!!!!!!!!!!!!!
    //cw[0] = a
    //cw[1] = b
    //cw[2] = PEAK POS
    //cw[3] = area - delta
    //cw[4] = FWHM - lor 0
    //cw[5] = AREA - lor 0
    //cw[2*n + 4] = FWHM - lor n
    //cw[2*n + 5]  = AREA - lor n
    variable timer
    variable numconts, ii, offset, area_res, conv_offset, rot_remaining, shift, thearea
    numconts = (numpnts(cw) - 4)/2
    yw = 0
    make/free/d/n=3 lor_cw = 0
    make/free/d/n=(numpnts(yw)) tempyw
    MarkPerfTestTime 0
    for(ii = 0 ; ii < numconts ; ii+=1)
        offset = 2*ii + 4
        lor_cw[0] = cw[2]
        lor_cw[1] = cw[offset]
        lor_cw[2] = cw[offset + 1]
        //multithread tempyw = 2*lor_cw[2]/pi * lor_cw[1]/(4*(xw1-lor_cw[0])^2 + lor_cw[1]^2)
        MPFXLorenzianPeak(lor_cw, tempyw, xw1)
        multithread yw += tempyw
    MarkPerfTestTime 1
    //convolve and account for the scaling
    convolve/a res, yw
    MarkPerfTestTime 2
        fastop yw = (xw1[1] - xw1[0]) *yw
    theArea = areaXY(xw1, res)
    fastop yw =  (1/theArea) * yw
    MarkPerfTestTime 3
    //have to rotate because of the convolution
    //the rotate operation gets you most of the way there
    conv_offset = binarysearchinterp(xw1, 0) - trunc(numpnts(xw1) / 2)
    rotate -trunc(conv_offset), yw
    duplicate/free yw, tempyw
    rot_remaining = conv_offset - trunc(conv_offset)
    yw = tempyw[p + rot_remaining]
    MarkPerfTestTime 4
    //add in scaling of the delta function
    duplicate/free res, tempyw
    shift = cw[2]/(xw1[1] - xw1[0])
    area_res = areaxy(W_energy, res)
    yw += tempyw[p-shift] * cw[3] / area_res
    //add in the background
    multithread yw += cw[0] + cw[1] * xw1
    MarkPerfTestTime 5
This is great, Andy- we like seeing folks post special-purpose code like this that can help others use Igor effectively, code that we don't have the expertise to write ourselves.

I wrote the following comments before I realized it was you... I'm sure you've thought about these issues... People reading my comments should not take them to indicate that they should avoid using this code!


It might also cause problems using the /D auto-destination feature.

You can most likely relieve that requirement by (carefully) making an intermediate wave to gather computed results, then using a wave assignment to read out linearly-interpolated values into the yw wave. See the second example here:

DisplayHelpTopic "All-At-Once Fitting Functions"

I am concerned about your use of a second independent variable to pass in the resolution function. That will make Igor think this is a multivariate fit function and /D flag will be made as a matrix. Your solution is easy to use, of course, since you can do it all in the Curve Fit dialog. But you would be better off to write a driver function and pass the resolution function via a wave that's looked up use a WAVE statement inside the function. You could also do this via a structure-based all-at-once fit function:

DisplayHelpTopic "Structure Fit Functions"

Since Igor is not able to make global structures, you *have* to make a driver function for the fit.

John Weeks
WaveMetrics, Inc.
As you mention the assumption that the code makes is that the xw1 wave is linearly spaced. The theoretical function has to be convolved with the resolution function (an experimentally measured curve). I had to assume linear spacing because the data I have is linearly spaced and I wanted to use the convolve operation. Depending on the number of points in the input data the convolution shifts data in its destination wave. I have to rotate those around manually (I spent ages trying to work this out)! I know this works if the x points are linearly spaced, not sure if it works if they are irregularly spaced.

If you know a better way to do the convolution please let me know.
As I said at the beginning of my comment, I wrote that before I realized who I was writing to. I know you well enough to know you thought this out...

My solution to convolution of non-evenly-spaced data is to make an intermediate wave that *is* evenly spaced and do all the computations using that. Then the last line is an interpolation into the yw wave. That's the way the convolution all-at-once complicated example works. It can be really tricky to figure out a heuristic for making an adequate intermediate wave.

John Weeks
WaveMetrics, Inc.
Hi Andy, I have been trying to implement this function but am only getting NANs out when fitting or using it as a function. I likely am missing something simple. I put the function in a procedure file with an .ipf extension. I have moved my resolution function to a wave in root which is named res. I created a coef wave 6 points long to use with a test of one lorentzian. I populated the coef wave with what seems to be reasonable starting points. I setup a curve fit session with my data wave and an x energy wave where the resolution wave is the same energies as the data waves. The function a returned a NAN apparently given the error in the fit. I try to plot the function as I change variables in the curve fit dialog box but it just generates NANs. Any suggestions on what I might be doing wrong? I am using igor 6.22 where it compiles without errors on a mac osx 10.9.2. I wonder if an example pxp file might be attached to help with new users of the function. Thanks! Dean
I think I have made some progress. There is a reference to an undefined attribute W_energy. After looking up areaway(), it seems to me like that should be the X energy wave so it can calculate the area of the resolution peak and then normalize the output wave. I changed W_energy to xw1 and it seems to work now for plotting. See below for a modified MarkPerfectTime 4 section.

    MarkPerfTestTime 4
    //add in scaling of the delta function
    duplicate/free res, tempyw
    shift = cw[2]/(xw1[1] - xw1[0])
    //daw print for debug
    //print W_energy
//  area_res = areaxy(W_energy, res)
    area_res = areaxy(xw1, res)
    yw += tempyw[p-shift] * cw[3] / area_res

Does this make sense?

So, at least I can now get the function to return something that seems workable. But, when I go to fit my QENS data what is returned as the fit wave has odd x scaling and the doesn't fit the data. It is broadened and too tall by 1-2 orders of magnitude. I have attached a image of the fit window with the resolution wave, test data wave, fit wave, and the wave I create with the fit parameters.

I think there is still something wrong.






Igor Pro 9

Learn More

Igor XOP Toolkit

Learn More

Igor NIDAQ Tools MX

Learn More