# A fit function for QENS data.

andyfaff

Wed, 04/06/2011 - 04:38 pm

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

//!!!!!!!!!!!!!!

//ASSUMES:

//1) THAT THE POINT SPACING IN XW1 IS LINEAR!!!!!!!!!!!!!!!!!

//2) THE RESOLUTION WAVE MUST HAVE THE SAME NUMBER OF POINTS AS XW1 AND BE THE WAVE ROOT:RES

//3) DO NOT TRY TO INTERPRET THE FUNCTION AT ANY OTHER POINTS THAN THE INPUT XW1 WAVE (IT DOESN'T MAKE SENSE).

//!!!!!!!!!!!!!!

//explanatory

//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

endfor

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

End

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

//!!!!!!!!!!!!!!

//ASSUMES:

//1) THAT THE POINT SPACING IN XW1 IS LINEAR!!!!!!!!!!!!!!!!!

//2) THE RESOLUTION WAVE MUST HAVE THE SAME NUMBER OF POINTS AS XW1 AND BE THE WAVE ROOT:RES

//3) DO NOT TRY TO INTERPRET THE FUNCTION AT ANY OTHER POINTS THAN THE INPUT XW1 WAVE (IT DOESN'T MAKE SENSE).

//!!!!!!!!!!!!!!

//explanatory

//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

endfor

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

End

Forum

Support

Gallery

### Igor Pro 8

Learn More

### Igor XOP Toolkit

Learn More

### Igor NIDAQ Tools MX

Learn More

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!"//ASSUMES THAT THE POINT SPACING IN XW1 IS LINEAR!!!!!!!!!!!!!!!!!"

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.

support@wavemetrics.com

April 7, 2011 at 09:27 am - Permalink

If you know a better way to do the convolution please let me know.

April 7, 2011 at 08:57 pm - Permalink

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.

support@wavemetrics.com

April 8, 2011 at 09:53 am - Permalink

March 1, 2014 at 01:54 pm - Permalink

//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.

Suggestions?

Dean

March 1, 2014 at 04:41 pm - Permalink