matthias muntwiler
fda49c3195
- new data reduction interface for more efficient multi-peak fitting. the new interface breaks compatibility with pre-2.0 data reduction functions. user-defined functions must be adapted to the new interface. - new angle scan processing panel for interactive data analysis.
1200 lines
37 KiB
Igor
1200 lines
37 KiB
Igor
#pragma rtGlobals=3 // Use modern global access method and strict wave access.
|
|
#pragma IgorVersion = 6.1
|
|
#pragma ModuleName = PearlScientaPreprocess
|
|
#include "pearl-fitfuncs"
|
|
|
|
// Copyright (c) 2013-18 Paul Scherrer Institut
|
|
//
|
|
// Licensed under the Apache License, Version 2.0 (the "License");
|
|
// you may not use this file except in compliance with the License.
|
|
// You may obtain a copy of the License at
|
|
// http://www.apache.org/licenses/LICENSE-2.0
|
|
|
|
/// @file
|
|
/// @brief preprocessing functions for Scienta detector images.
|
|
/// @ingroup ArpesPackage
|
|
///
|
|
/// this procedure contains functions for data reduction
|
|
/// and instrument-specific normalization.
|
|
///
|
|
/// @author matthias muntwiler, matthias.muntwiler@psi.ch
|
|
///
|
|
/// @copyright 2013-18 Paul Scherrer Institut @n
|
|
/// Licensed under the Apache License, Version 2.0 (the "License"); @n
|
|
/// you may not use this file except in compliance with the License. @n
|
|
/// You may obtain a copy of the License at
|
|
/// http://www.apache.org/licenses/LICENSE-2.0
|
|
|
|
/// @namespace PearlScientaPreprocess
|
|
/// @brief preprocessing functions for Scienta detector images.
|
|
///
|
|
/// PearlScientaPreprocess is declared in @ref pearl-scienta-preprocess.ipf.
|
|
|
|
/// prompt the user for integrate on linear background reduction parameters.
|
|
///
|
|
function prompt_int_linbg_reduction(param)
|
|
string ¶m
|
|
|
|
variable Lcrop = NumberByKey("Lcrop", param, "=", ";")
|
|
variable Lsize = NumberByKey("Lsize", param, "=", ";")
|
|
variable Hcrop = NumberByKey("Hcrop", param, "=", ";")
|
|
variable Hsize = NumberByKey("Hsize", param, "=", ";")
|
|
variable Cpos = NumberByKey("Cpos", param, "=", ";")
|
|
variable Csize = NumberByKey("Csize", param, "=", ";")
|
|
|
|
prompt Lcrop, "Lower cropping region"
|
|
prompt Hcrop, "Upper cropping region"
|
|
prompt Lsize, "Lower background region"
|
|
prompt Hsize, "Upper background region"
|
|
prompt Cpos, "Center position"
|
|
prompt Csize, "Center integration region"
|
|
|
|
doprompt "int_linbg_reduction Parameters", lcrop, hcrop, lsize, hsize, cpos, csize
|
|
if (v_flag == 0)
|
|
param = ReplaceNumberByKey("Lcrop", param, Lcrop, "=", ";")
|
|
param = ReplaceNumberByKey("Lsize", param, Lsize, "=", ";")
|
|
param = ReplaceNumberByKey("Hcrop", param, Hcrop, "=", ";")
|
|
param = ReplaceNumberByKey("Hsize", param, Hsize, "=", ";")
|
|
param = ReplaceNumberByKey("Cpos", param, Cpos, "=", ";")
|
|
param = ReplaceNumberByKey("Csize", param, Csize, "=", ";")
|
|
endif
|
|
|
|
return v_flag
|
|
end
|
|
|
|
/// this function is for testing only, until we implement a proper interface
|
|
///
|
|
function /s capture_int_linbg_cursors()
|
|
string param = csr_int_linbg_reduction("")
|
|
svar /z global_params = root:packages:pearl_explorer:s_reduction_params
|
|
if (svar_exists(global_params))
|
|
global_params = param
|
|
endif
|
|
return param
|
|
end
|
|
|
|
/// set reduction parameters from cursors in a graph.
|
|
///
|
|
/// PRELIMINARY - function arguments may change
|
|
///
|
|
/// sets reduction parameters from cursors in a graph.
|
|
/// an even number of cursors (2 or more) must be set on the image.
|
|
/// cursor names and order do not matter,
|
|
/// except that the alphabetically first cursor which is attached to an image selects the image.
|
|
/// the cursors mark the following positions, from innermost to outermost pair:
|
|
/// 1) low and high limits of peak region.
|
|
/// 2) peak-side boundary of lower and upper background region.
|
|
/// 3) lower and upper cropping region.
|
|
///
|
|
function /s csr_int_linbg_reduction(win)
|
|
string win
|
|
|
|
// read all cursor positions
|
|
variable ic
|
|
string sc
|
|
variable nc = 10
|
|
make /n=(nc) /free positions
|
|
variable np = 0
|
|
wave /z image = $""
|
|
string imagename = ""
|
|
string tracename = ""
|
|
string info
|
|
|
|
for (ic = 0; ic < nc; ic += 1)
|
|
sc = num2char(char2num("A") + ic)
|
|
wave /z wc = CsrWaveRef($sc, win)
|
|
info = CsrInfo($sc, win)
|
|
tracename = StringByKey("TNAME", info, ":", ";")
|
|
if (waveexists(wc) && (wavedims(wc) == 2))
|
|
if (!waveexists(image))
|
|
wave image = wc
|
|
imagename = tracename
|
|
endif
|
|
if (cmpstr(tracename, imagename) == 0)
|
|
positions[np] = pcsr($sc, win)
|
|
np += 1
|
|
endif
|
|
endif
|
|
endfor
|
|
|
|
np = floor(np / 2) * 2 // ignore odd cursor
|
|
redimension /n=(np) positions
|
|
sort positions, positions
|
|
// shift upper positions by one so that the rightmost pixel becomes 1.0
|
|
positions = p >= np / 2 ? positions + 1 : positions
|
|
positions = positions / dimsize(image, 0)
|
|
|
|
// map innermost cursor pair to peak center and size
|
|
variable ip2 = np / 2
|
|
variable ip1 = ip2 - 1
|
|
variable Cpos = (positions[ip1] + positions[ip2]) / 2
|
|
variable Csize = positions[ip2] - positions[ip1]
|
|
if (ip1 >= 0)
|
|
Cpos = (positions[ip1] + positions[ip2]) / 2
|
|
Csize = positions[ip2] - positions[ip1]
|
|
else
|
|
// default: a small region in the center
|
|
Cpos = 0.5
|
|
Csize = 0.2
|
|
endif
|
|
|
|
// background region
|
|
ip1 -= 1
|
|
ip2 += 1
|
|
variable Lsize
|
|
variable Hsize
|
|
if (ip1 >= 0)
|
|
Lsize = positions[ip1]
|
|
Hsize = 1 - positions[ip2]
|
|
else
|
|
// default: everything outside the peak region
|
|
Lsize = Cpos - Csize / 2
|
|
Hsize = 1 - (Cpos + Csize / 2)
|
|
endif
|
|
|
|
// crop region
|
|
ip1 -= 1
|
|
ip2 += 1
|
|
variable Lcrop
|
|
variable Hcrop
|
|
if (ip1 >= 0)
|
|
Lcrop = positions[ip1]
|
|
Hcrop = 1 - positions[ip2]
|
|
else
|
|
// default: dead corners of the EW4000 at PEARL
|
|
Lcrop = 0.11
|
|
Hcrop = 0.11
|
|
endif
|
|
Lsize = max(Lsize - Lcrop, 0)
|
|
Hsize = max(Hsize - Hcrop, 0)
|
|
|
|
string param = ""
|
|
param = ReplaceNumberByKey("Lcrop", param, Lcrop, "=", ";")
|
|
param = ReplaceNumberByKey("Lsize", param, Lsize, "=", ";")
|
|
param = ReplaceNumberByKey("Hcrop", param, Hcrop, "=", ";")
|
|
param = ReplaceNumberByKey("Hsize", param, Hsize, "=", ";")
|
|
param = ReplaceNumberByKey("Cpos", param, Cpos, "=", ";")
|
|
param = ReplaceNumberByKey("Csize", param, Csize, "=", ";")
|
|
|
|
return param
|
|
end
|
|
|
|
/// linear-background subtracted integration reduction function.
|
|
///
|
|
/// data reduction function for adh5_load_reduced_detector.
|
|
/// cf. @ref adh5_default_reduction for an explanation of reduction functions.
|
|
///
|
|
/// this function calculates the average pixel value of each angular slice
|
|
/// in one center and two background intervals.
|
|
/// a background value is calculated at the center position
|
|
/// by linear interpolation from the two background values.
|
|
/// returns the center minus linear background in dest1.
|
|
/// returns the Poisson one-sigma error in dest2.
|
|
///
|
|
/// typical values (peak centered on detector, FWHM ~ 20 % of image)
|
|
/// Lcrop=0.11;Hcrop=0.11;Lsize=0.2;Hsize=0.2;Cpos=0.5;Csize=0.2
|
|
///
|
|
/// @param source scienta detector image, energy axis along X, angle axis along Y.
|
|
/// two-dimensional intensity distribution (image).
|
|
/// the scales are carried over to the result waves.
|
|
///
|
|
/// @param param parameters in a key1=value1;key2=value2;... list.
|
|
/// all region parameters are relative to the image size (0...1).
|
|
/// @arg Lcrop = size of the lower cropping region
|
|
/// @arg Hcrop = size of the upper cropping region
|
|
/// @arg Lsize = size of the lower background integration region
|
|
/// @arg Hsize = size of the upper background integration region
|
|
/// @arg Cpos = center position of the of the peak integration region
|
|
/// @arg Csize = size of the peak integration region
|
|
///
|
|
/// @return free wave containing references of the two result waves.
|
|
/// the first wave is the integral minus linear background.
|
|
/// the second wave is the Poisson one-sigma error.
|
|
///
|
|
threadsafe function /wave int_linbg_reduction(source, param)
|
|
wave source
|
|
string ¶m
|
|
|
|
variable nx = dimsize(source, 0)
|
|
variable ny = dimsize(source, 1)
|
|
|
|
// read parameters
|
|
variable lcrop = NumberByKey("Lcrop", param, "=", ";")
|
|
variable lsize = NumberByKey("Lsize", param, "=", ";")
|
|
variable hcrop = NumberByKey("Hcrop", param, "=", ";")
|
|
variable hsize = NumberByKey("Hsize", param, "=", ";")
|
|
variable cpos = NumberByKey("Cpos", param, "=", ";")
|
|
variable csize = NumberByKey("Csize", param, "=", ";")
|
|
|
|
make /wave /free /n=2 result_waves
|
|
make /free /n=0 dest1, dest2
|
|
result_waves[0] = dest1
|
|
result_waves[1] = dest2
|
|
adh5_setup_profile(source, dest1, 1)
|
|
adh5_setup_profile(source, dest2, 1)
|
|
|
|
// validate parameters
|
|
// background parameters are optional, center parameter is required.
|
|
if (numtype(lcrop) != 0)
|
|
lcrop = 0
|
|
endif
|
|
if (numtype(lsize) != 0)
|
|
lsize = 0
|
|
endif
|
|
if (numtype(hcrop) != 0)
|
|
hcrop = 0
|
|
endif
|
|
if (numtype(hsize) != 0)
|
|
hsize = 0
|
|
endif
|
|
if (numtype(Cpos) != 0)
|
|
redimension /n=0 result_waves
|
|
return result_waves // Cpos parameter missing
|
|
endif
|
|
if (numtype(Csize) != 0)
|
|
redimension /n=0 result_waves
|
|
return result_waves // Csize parameter missing
|
|
endif
|
|
|
|
variable lpos = lcrop + lsize / 2
|
|
variable hpos = 1 - (hcrop + hsize / 2)
|
|
|
|
variable p0
|
|
variable p1
|
|
|
|
duplicate /free dest1, lbg, hbg
|
|
if (lsize > 0)
|
|
p0 = round(lcrop * nx)
|
|
p1 = round((lcrop + lsize) * nx)
|
|
ad_profile_y_w(source, p0, p1, lbg)
|
|
else
|
|
lbg = 0
|
|
endif
|
|
if (hsize > 0)
|
|
p0 = round((1 - hcrop - hsize) * nx)
|
|
p1 = round((1 - hcrop) * nx)
|
|
ad_profile_y_w(source, p0, p1, hbg)
|
|
else
|
|
hbg = 0
|
|
endif
|
|
if (csize > 0)
|
|
p0 = round((cpos - csize/2) * nx)
|
|
p1 = round((cpos + csize/2) * nx)
|
|
ad_profile_y_w(source, p0, p1, dest1)
|
|
else
|
|
dest1 = 0
|
|
endif
|
|
|
|
variable scale = (cpos - lpos) / (hpos - lpos)
|
|
dest2 = dest1
|
|
dest1 -= scale * (hbg - lbg) + lbg
|
|
dest2 = sqrt(dest2 + scale^2 * (hbg + lbg))
|
|
|
|
string s_note1
|
|
string s_note2
|
|
sprintf s_note1, "AxisLabelD=peak integral"
|
|
sprintf s_note2, "KineticEnergy=%.3f", cpos * nx * dimdelta(source, 0) + dimoffset(source, 0)
|
|
Note dest1, s_note1
|
|
Note dest1, s_note2
|
|
Note dest2, s_note1
|
|
Note dest2, s_note2
|
|
|
|
return result_waves
|
|
end
|
|
|
|
function prompt_int_quadbg_reduction(param)
|
|
string ¶m
|
|
|
|
variable Lcrop = NumberByKey("Lcrop", param, "=", ";")
|
|
variable Lsize = NumberByKey("Lsize", param, "=", ";")
|
|
variable Hcrop = NumberByKey("Hcrop", param, "=", ";")
|
|
variable Hsize = NumberByKey("Hsize", param, "=", ";")
|
|
variable Cpos = NumberByKey("Cpos", param, "=", ";")
|
|
variable Csize = NumberByKey("Csize", param, "=", ";")
|
|
|
|
prompt Lcrop, "Lower cropping region"
|
|
prompt Hcrop, "Upper cropping region"
|
|
prompt Lsize, "Lower background region"
|
|
prompt Hsize, "Upper background region"
|
|
prompt Cpos, "Center position"
|
|
prompt Csize, "Center integration region"
|
|
|
|
doprompt "int_quadbg_reduction Parameters", lcrop, hcrop, lsize, hsize, cpos, csize
|
|
if (v_flag == 0)
|
|
param = ReplaceNumberByKey("Lcrop", param, Lcrop, "=", ";")
|
|
param = ReplaceNumberByKey("Lsize", param, Lsize, "=", ";")
|
|
param = ReplaceNumberByKey("Hcrop", param, Hcrop, "=", ";")
|
|
param = ReplaceNumberByKey("Hsize", param, Hsize, "=", ";")
|
|
param = ReplaceNumberByKey("Cpos", param, Cpos, "=", ";")
|
|
param = ReplaceNumberByKey("Csize", param, Csize, "=", ";")
|
|
endif
|
|
|
|
return v_flag
|
|
end
|
|
|
|
/// integrate peak area minus a quadratic background
|
|
///
|
|
/// data reduction function for adh5_load_reduced_detector.
|
|
/// cf. @ref adh5_default_reduction for an explanation of reduction functions.
|
|
///
|
|
/// this function calculates the average pixel value of each angular slice
|
|
/// in one center and two background intervals.
|
|
/// a background value is calculated at the center position
|
|
/// by linear interpolation from the two background values.
|
|
/// returns the center minus linear background in dest1.
|
|
/// returns the Poisson one-sigma error in dest2.
|
|
///
|
|
/// typical values (peak centered on detector, FWHM ~ 20 % of image)
|
|
/// Lcrop=0.11;Hcrop=0.11;Lsize=0.2;Hsize=0.2;Cpos=0.5;Csize=0.2
|
|
///
|
|
/// @param source scienta detector image, energy axis along X, angle axis along Y.
|
|
/// two-dimensional intensity distribution (image).
|
|
/// the scales are carried over to the result waves.
|
|
///
|
|
/// @param param parameters in a key1=value1;key2=value2;... list.
|
|
/// all region parameters are relative to the image size (0...1).
|
|
/// @arg Lcrop = size of the lower cropping region
|
|
/// @arg Hcrop = size of the upper cropping region
|
|
/// @arg Lsize = size of the lower background integration region
|
|
/// @arg Hsize = size of the upper background integration region
|
|
/// @arg Cpos = center position of the of the peak integration region
|
|
/// @arg Csize = size of the peak integration region
|
|
///
|
|
/// @return free wave containing references of the two result waves.
|
|
/// the first wave is the integral minus linear background.
|
|
/// the second wave is the Poisson one-sigma error.
|
|
///
|
|
threadsafe function /wave int_quadbg_reduction(source, param)
|
|
wave source
|
|
string ¶m
|
|
|
|
variable nx = dimsize(source, 0)
|
|
variable ny = dimsize(source, 1)
|
|
|
|
// read parameters
|
|
variable lcrop = NumberByKey("Lcrop", param, "=", ";")
|
|
variable lsize = NumberByKey("Lsize", param, "=", ";")
|
|
variable hcrop = NumberByKey("Hcrop", param, "=", ";")
|
|
variable hsize = NumberByKey("Hsize", param, "=", ";")
|
|
variable cpos = NumberByKey("Cpos", param, "=", ";")
|
|
variable csize = NumberByKey("Csize", param, "=", ";")
|
|
|
|
make /wave /free /n=2 result_waves
|
|
make /free /n=0 dest1, dest2
|
|
result_waves[0] = dest1
|
|
result_waves[1] = dest2
|
|
adh5_setup_profile(source, dest1, 1)
|
|
adh5_setup_profile(source, dest2, 1)
|
|
|
|
// validate parameters
|
|
// background parameters are optional, center parameter is required.
|
|
if (numtype(lcrop) != 0)
|
|
lcrop = 0
|
|
endif
|
|
if (numtype(lsize) != 0)
|
|
lsize = 0
|
|
endif
|
|
if (numtype(hcrop) != 0)
|
|
hcrop = 0
|
|
endif
|
|
if (numtype(hsize) != 0)
|
|
hsize = 0
|
|
endif
|
|
if (numtype(Cpos) != 0)
|
|
redimension /n=0 result_waves
|
|
return result_waves // Cpos parameter missing
|
|
endif
|
|
if (numtype(Csize) != 0)
|
|
redimension /n=0 result_waves
|
|
return result_waves // Csize parameter missing
|
|
endif
|
|
|
|
// crop boundaries
|
|
variable pcl = round(lcrop * nx)
|
|
variable pch = round((1 - hcrop) * nx)
|
|
// fit boundaries
|
|
variable pfl = round((lcrop + lsize) * nx)
|
|
variable pfh = round((1 - hcrop - hsize) * nx)
|
|
// integration boundaries
|
|
variable pil = round((cpos - csize/2) * nx)
|
|
variable pih = round((cpos + csize/2) * nx)
|
|
|
|
// prepare intermediate data buffer
|
|
make /n=(nx) /d /free profile, mask, fit
|
|
setscale /p x dimoffset(source,0), dimdelta(source,0), waveunits(source,0), profile, mask, fit
|
|
mask = ((p >= pcl) && (p < pfl)) || ((p >= pfh) && (p < pch))
|
|
|
|
variable qq
|
|
variable sp, sf
|
|
variable xil = x2pnt(profile, pil)
|
|
variable xih = x2pnt(profile, pih)
|
|
|
|
make /n=3 /free /d w_coef
|
|
for (qq = 0; qq < ny; qq += 1)
|
|
profile = source[p][qq]
|
|
curvefit /Q /NTHR=1 /W=2 poly 3, kwCWave=w_coef, profile /M=mask
|
|
fit = poly(w_coef, x)
|
|
sp = sum(profile, xil, xih)
|
|
sf = sum(fit, xil, xih)
|
|
dest1[qq] = sp - sf
|
|
dest2[qq] = sqrt(sp)
|
|
endfor
|
|
|
|
string s_note1
|
|
string s_note2
|
|
sprintf s_note1, "AxisLabelD=peak integral"
|
|
sprintf s_note2, "KineticEnergy=%.3f", cpos * nx * dimdelta(source, 0) + dimoffset(source, 0)
|
|
Note dest1, s_note1
|
|
Note dest1, s_note2
|
|
Note dest2, s_note1
|
|
Note dest2, s_note2
|
|
|
|
return result_waves
|
|
end
|
|
|
|
/// parameter dialog for the redim_linbg_reduction() function
|
|
///
|
|
/// @param param parameter string in a key1=value1;key2=value2;... list.
|
|
/// the parameter string is passed by reference.
|
|
/// see redim_linbg_reduction() for a description of parameters.
|
|
///
|
|
/// @return zero if the user clicked OK, non-zero if the user clicked Cancel.
|
|
///
|
|
function prompt_redim_linbg_reduction(param)
|
|
string ¶m
|
|
|
|
variable Lcrop = NumberByKey("Lcrop", param, "=", ";")
|
|
variable Lsize = NumberByKey("Lsize", param, "=", ";")
|
|
variable Hcrop = NumberByKey("Hcrop", param, "=", ";")
|
|
variable Hsize = NumberByKey("Hsize", param, "=", ";")
|
|
variable Cpos = NumberByKey("Cpos", param, "=", ";")
|
|
variable Csize = NumberByKey("Csize", param, "=", ";")
|
|
|
|
prompt Lcrop, "Lower cropping region"
|
|
prompt Hcrop, "Upper cropping region"
|
|
prompt Lsize, "Lower background region"
|
|
prompt Hsize, "Upper background region"
|
|
prompt Cpos, "Center position"
|
|
prompt Csize, "Center integration region"
|
|
|
|
doprompt "redim_linbg_reduction Parameters", lcrop, hcrop, lsize, hsize, cpos, csize
|
|
if (v_flag == 0)
|
|
param = ReplaceNumberByKey("Lcrop", param, Lcrop, "=", ";")
|
|
param = ReplaceNumberByKey("Lsize", param, Lsize, "=", ";")
|
|
param = ReplaceNumberByKey("Hcrop", param, Hcrop, "=", ";")
|
|
param = ReplaceNumberByKey("Hsize", param, Hsize, "=", ";")
|
|
param = ReplaceNumberByKey("Cpos", param, Cpos, "=", ";")
|
|
param = ReplaceNumberByKey("Csize", param, Csize, "=", ";")
|
|
endif
|
|
|
|
return v_flag
|
|
end
|
|
|
|
/// linear background reduction function for incorrectly dimensioned scienta image
|
|
///
|
|
/// if the energy step size does not divide the energy range to an integer number,
|
|
/// the scienta image is exported with the wrong array size.
|
|
/// this can be fixed by redimensioning the array.
|
|
///
|
|
/// the current implementation works in the case where dimension 0 needs to be incremented.
|
|
/// the function may be generalized to dimension 1 and/or decrementing by additional parameters.
|
|
/// it is not known yet whether a generalization is needed or whether it can cover all cases.
|
|
///
|
|
/// background subtraction and peak integration is the same as by the int_linbg_reduction() function.
|
|
///
|
|
/// @param source source wave
|
|
/// Scienta detector image, energy axis along X, angle axis along Y
|
|
///
|
|
/// @param param parameter string in a key1=value1;key2=value2;... list.
|
|
/// the parameter string is passed by reference.
|
|
///
|
|
/// all region parameters are relative to the image size (0...1).
|
|
/// @arg Lcrop size of the lower cropping region
|
|
/// @arg Hcrop size of the upper cropping region
|
|
/// @arg Lsize size of the lower background integration region
|
|
/// @arg Hsize size of the upper background integration region
|
|
/// @arg Cpos center position of the of the peak integration region
|
|
/// @arg Csize size of the peak integration region
|
|
///
|
|
/// typical values (peak centered on detector, FWHM ~ 20 % of image)
|
|
/// Lcrop=0.11;Hcrop=0.11;Lsize=0.2;Hsize=0.2;Cpos=0.5;Csize=0.2
|
|
///
|
|
/// @return free wave containing references of the two result waves.
|
|
/// the first wave is the integral minus linear background.
|
|
/// the second wave is the Poisson one-sigma error.
|
|
///
|
|
threadsafe function /wave redim_linbg_reduction(source, param)
|
|
wave source
|
|
string ¶m
|
|
|
|
variable nx = dimsize(source, 0)
|
|
variable ny = dimsize(source, 1)
|
|
|
|
duplicate /free source, source_redim
|
|
redimension /n=(nx * ny) source_redim
|
|
nx += 1
|
|
redimension /n=(nx, ny) source_redim
|
|
|
|
return int_linbg_reduction(source_redim, param)
|
|
end
|
|
|
|
/// apply the gauss4_reduction function to a single image
|
|
///
|
|
/// useful for testing or manual processing.
|
|
/// to debug, (temporarily) remove the threadsafe attribute from the gauss2_reduction function.
|
|
///
|
|
function test_gauss4_reduction(image)
|
|
wave image
|
|
|
|
string param = ""
|
|
|
|
param = ReplaceNumberByKey("rngl", param, -inf, "=", ";")
|
|
param = ReplaceNumberByKey("rngh", param, inf, "=", ";")
|
|
param = ReplaceNumberByKey("npeaks", param, 4, "=", ";")
|
|
param = ReplaceNumberByKey("ybox", param, 1, "=", ";")
|
|
param = ReplaceNumberByKey("pos1", param, 11, "=", ";")
|
|
param = ReplaceNumberByKey("wid1", param, 0.1, "=", ";")
|
|
param = ReplaceNumberByKey("pos2", param, 12, "=", ";")
|
|
param = ReplaceNumberByKey("wid2", param, 0.2, "=", ";")
|
|
param = ReplaceNumberByKey("pos3", param, 13, "=", ";")
|
|
param = ReplaceNumberByKey("wid3", param, 0.3, "=", ";")
|
|
param = ReplaceNumberByKey("pos4", param, 14, "=", ";")
|
|
param = ReplaceNumberByKey("wid4", param, 0.4, "=", ";")
|
|
|
|
wave /wave results = gauss4_reduction(image, param)
|
|
|
|
variable npk = numpnts(results) / 2
|
|
variable ipk
|
|
string sw
|
|
for (ipk = 0; ipk < npk; ipk += 1)
|
|
sw = "test_int_" + num2str(ipk + 1)
|
|
duplicate /o results[ipk], $sw
|
|
sw = "test_sig_" + num2str(ipk + 1)
|
|
duplicate /o results[ipk + npk], $sw
|
|
endfor
|
|
end
|
|
|
|
/// prompt for the gauss4_reduction parameters
|
|
///
|
|
///
|
|
function prompt_gauss4_reduction(param)
|
|
string ¶m
|
|
|
|
variable rngl = NumberByKey("rngl", param, "=", ";")
|
|
variable rngh = NumberByKey("rngh", param, "=", ";")
|
|
variable pos1 = NumberByKey("pos1", param, "=", ";")
|
|
variable wid1 = NumberByKey("wid1", param, "=", ";")
|
|
variable pos2 = NumberByKey("pos2", param, "=", ";")
|
|
variable wid2 = NumberByKey("wid2", param, "=", ";")
|
|
variable pos3 = NumberByKey("pos3", param, "=", ";")
|
|
variable wid3 = NumberByKey("wid3", param, "=", ";")
|
|
variable pos4 = NumberByKey("pos4", param, "=", ";")
|
|
variable wid4 = NumberByKey("wid4", param, "=", ";")
|
|
variable npeaks = NumberByKey("npeaks", param, "=", ";")
|
|
variable ybox = NumberByKey("ybox", param, "=", ";")
|
|
|
|
prompt rngl, "range low"
|
|
prompt rngh, "range high"
|
|
prompt pos1, "position 1"
|
|
prompt wid1, "width 1"
|
|
prompt pos2, "position 2"
|
|
prompt wid2, "width 2"
|
|
prompt pos3, "position 3"
|
|
prompt wid3, "width 3"
|
|
prompt pos4, "position 4"
|
|
prompt wid4, "width 4"
|
|
prompt npeaks, "number of peaks"
|
|
prompt ybox, "ybox (1 or 3)"
|
|
|
|
doprompt "gauss4_reduction reduction parameters (1/2)", rngl, rngh, npeaks, ybox
|
|
if (v_flag == 0)
|
|
param = ReplaceNumberByKey("rngl", param, rngl, "=", ";")
|
|
param = ReplaceNumberByKey("rngh", param, rngh, "=", ";")
|
|
param = ReplaceNumberByKey("npeaks", param, npeaks, "=", ";")
|
|
param = ReplaceNumberByKey("ybox", param, ybox, "=", ";")
|
|
|
|
doprompt "gauss4_reduction reduction parameters (2/2)", pos1, wid1, pos2, wid2, pos3, wid3, pos4, wid4
|
|
if (v_flag == 0)
|
|
param = ReplaceNumberByKey("pos1", param, pos1, "=", ";")
|
|
param = ReplaceNumberByKey("wid1", param, wid1, "=", ";")
|
|
param = ReplaceNumberByKey("pos2", param, pos2, "=", ";")
|
|
param = ReplaceNumberByKey("wid2", param, wid2, "=", ";")
|
|
param = ReplaceNumberByKey("pos3", param, pos3, "=", ";")
|
|
param = ReplaceNumberByKey("wid3", param, wid3, "=", ";")
|
|
param = ReplaceNumberByKey("pos4", param, pos4, "=", ";")
|
|
param = ReplaceNumberByKey("wid4", param, wid4, "=", ";")
|
|
endif
|
|
endif
|
|
|
|
return v_flag
|
|
end
|
|
|
|
/// fit horizontal cuts of an image with up to four gaussian peaks on a linear background
|
|
///
|
|
/// the function fits each horizontal profile (EDC) with four gaussian peaks on a linear background.
|
|
/// the position and width of the peaks is kept fixed according to input parameters.
|
|
/// the peak amplitude is constrained to positive value.
|
|
///
|
|
/// the width parameter is defined as in Igor's gauss curve fit function
|
|
/// (standard deviation divided by the square root of two).
|
|
/// the return value in dest1 is the integrated peak of one of the peaks.
|
|
/// dest2 returns the corresponding error estimate.
|
|
///
|
|
/// @param source source wave.
|
|
/// two-dimensional distribution of counts.
|
|
/// for correct weighting and error estimation it is important
|
|
/// that the source wave contains actual counts (Poisson statistics).
|
|
///
|
|
/// @param param (in, out) semicolon-separated key=value list of processing parameters.
|
|
/// this is a pass-by-reference argument.
|
|
/// the following parameters are required.
|
|
/// position, width and limit parameters are on the x (energy) scale.
|
|
/// @arg rngl low limit of fit interval
|
|
/// @arg rngh high limit of fit interval
|
|
/// @arg pos1 position of peak 1
|
|
/// @arg wid1 width of peak 1
|
|
/// @arg pos2 position of peak 2
|
|
/// @arg wid2 width of peak 2
|
|
/// @arg pos3 position of peak 3
|
|
/// @arg wid3 width of peak 3
|
|
/// @arg pos4 position of peak 4
|
|
/// @arg wid4 width of peak 4
|
|
/// @arg npeaks number of peaks to fit: 1...4
|
|
/// the others are held at amplitude 0.
|
|
/// @arg ybox box size of averaging in y direction, must be 1 or 3.
|
|
/// other values lead to corrupt data.
|
|
///
|
|
/// @return free wave containing references of the result waves.
|
|
/// the number of waves is two times the number of peaks that are fit.
|
|
/// the first npeaks waves contain the peak integrals,
|
|
/// the second npeaks waves the corresponding error estimates.
|
|
///
|
|
threadsafe function /wave gauss4_reduction(source, param)
|
|
wave source
|
|
string ¶m
|
|
|
|
variable nx = dimsize(source, 0)
|
|
variable ny = dimsize(source, 1)
|
|
|
|
// read parameters
|
|
variable rngl = NumberByKey("rngl", param, "=", ";")
|
|
variable rngh = NumberByKey("rngh", param, "=", ";")
|
|
variable pos1 = NumberByKey("pos1", param, "=", ";")
|
|
variable wid1 = NumberByKey("wid1", param, "=", ";")
|
|
variable pos2 = NumberByKey("pos2", param, "=", ";")
|
|
variable wid2 = NumberByKey("wid2", param, "=", ";")
|
|
variable pos3 = NumberByKey("pos3", param, "=", ";")
|
|
variable wid3 = NumberByKey("wid3", param, "=", ";")
|
|
variable pos4 = NumberByKey("pos4", param, "=", ";")
|
|
variable wid4 = NumberByKey("wid4", param, "=", ";")
|
|
variable npeaks = NumberByKey("npeaks", param, "=", ";")
|
|
variable ybox = NumberByKey("ybox", param, "=", ";")
|
|
|
|
// prepare curve fit
|
|
variable ipk
|
|
make /free xprof
|
|
adh5_setup_profile(source, xprof, 0)
|
|
duplicate /free xprof, xprof_sig
|
|
variable pl = max(x2pnt(xprof, rngl), 0)
|
|
variable ph = min(x2pnt(xprof, rngh), numpnts(xprof) - 1)
|
|
|
|
make /free /n=(npeaks) peak_coef
|
|
peak_coef = p * 3 + 2
|
|
variable n_coef = npeaks * 3 + 2
|
|
make /free /d /n=14 w_coef, W_sigma
|
|
w_coef[0] = {0, 0, 1, pos1, wid1, 1, pos2, wid2, 1, pos3, wid3, 1, pos4, wid4}
|
|
redimension /n=(n_coef) w_coef, w_sigma
|
|
|
|
// text constraints cannot be used in threadsafe functions.
|
|
// the following matrix-vector forumlation is equivalent to:
|
|
// make /free /T /N=6 constraints
|
|
// constraints[0] = {"K2 >= 0", "K5 >= 0", "K8 >= 0", "K11 >= 0", "K1 <= 0", "K0 => 0"}
|
|
make /free /n=(npeaks + 2, numpnts(w_coef)) cmat
|
|
make /free /n=(npeaks + 2) cvec
|
|
cmat = 0
|
|
cmat[0][0] = -1
|
|
cmat[1][1] = 1
|
|
cvec = 0
|
|
|
|
string hold = "00"
|
|
for (ipk=0; ipk < npeaks; ipk += 1)
|
|
hold += "011"
|
|
cmat[2 + ipk][2 + ipk*3] = -1
|
|
endfor
|
|
|
|
// prepare output
|
|
make /free /n=(npeaks * 2) /wave result_waves
|
|
string s_note
|
|
for (ipk = 0; ipk < npeaks; ipk += 1)
|
|
make /free /n=0 pk_int
|
|
adh5_setup_profile(source, pk_int, 1)
|
|
pk_int = nan
|
|
sprintf s_note, "AxisLabelD=peak %u integral", ipk+1
|
|
Note pk_int, s_note
|
|
sprintf s_note, "KineticEnergy=%.3f", w_coef[3 + ipk * 3]
|
|
Note pk_int, s_note
|
|
result_waves[ipk] = pk_int
|
|
|
|
make /free /n=0 pk_sig
|
|
adh5_setup_profile(source, pk_sig, 1)
|
|
pk_sig = nan
|
|
sprintf s_note, "AxisLabelD=peak %u sigma", ipk+1
|
|
Note pk_sig, s_note
|
|
sprintf s_note, "KineticEnergy=%.3f", w_coef[3 + ipk * 3]
|
|
Note pk_sig, s_note
|
|
result_waves[ipk + npeaks] = pk_sig
|
|
|
|
waveclear pk_int, pk_sig
|
|
endfor
|
|
|
|
// loop over angle scale
|
|
variable p0 = 0
|
|
variable p1 = dimsize(source, 1) - 1
|
|
variable pp
|
|
variable wmin
|
|
variable wmax
|
|
if (ybox > 1)
|
|
p0 += ceil((ybox - 1) / 2)
|
|
p1 -= ceil((ybox - 1) / 2)
|
|
endif
|
|
variable V_FitNumIters
|
|
|
|
for (pp = p0; pp <= p1; pp += 1)
|
|
// box average
|
|
xprof = source[p][pp]
|
|
if (ybox > 1)
|
|
xprof += source[p][pp-1] + source[p][pp+1]
|
|
endif
|
|
xprof_sig = max(sqrt(xprof), 1)
|
|
xprof /= ybox
|
|
xprof_sig /= ybox
|
|
|
|
// generate guess
|
|
wmin = wavemin(xprof)
|
|
wmax = wavemax(xprof)
|
|
w_coef[0] = wmin
|
|
w_coef[1] = 0
|
|
|
|
for (ipk=0; ipk < npeaks; ipk += 1)
|
|
w_coef[2 + ipk*3] = wmax - wmin
|
|
endfor
|
|
FuncFit /H=hold /Q /NTHR=1 /N /W=2 MultiGaussLinBG_AO w_coef xprof[pl,ph] /C={cmat, cvec} /I=1 /W=xprof_sig[pl,ph]
|
|
wave w_sigma
|
|
|
|
// retrieve results, leave them at nan if the fit did not converge
|
|
if (V_FitNumIters < 40)
|
|
for (ipk = 0; ipk < npeaks; ipk += 1)
|
|
wave val = result_waves[ipk]
|
|
wave sig = result_waves[ipk + npeaks]
|
|
val[pp] = max(w_coef[peak_coef[ipk]], 0)
|
|
sig[pp] = max(w_sigma[peak_coef[ipk]], 0)
|
|
endfor
|
|
endif
|
|
endfor
|
|
|
|
// calculate integral
|
|
for (ipk = 0; ipk < npeaks; ipk += 1)
|
|
wave val = result_waves[ipk]
|
|
wave sig = result_waves[ipk + npeaks]
|
|
val *= w_coef[peak_coef[ipk] + 2] * sqrt(pi)
|
|
sig *= w_coef[peak_coef[ipk] + 2] * sqrt(pi)
|
|
endfor
|
|
|
|
return result_waves
|
|
end
|
|
|
|
|
|
/// find peak positions for the gauss-fit reduction function
|
|
///
|
|
/// the function fits a reference spectrum with up to four gaussian peaks on a linear background.
|
|
/// and returns the results in a parameter string for the @ref gauss4_reduction function.
|
|
///
|
|
/// @param spectrum reference spectrum (EDC),
|
|
/// e.g. a ScientaImage integrated over the scan axes:
|
|
/// @code ad_extract_rod(ScientaImage, nan, nan, -inf, inf, -inf, inf, "ScientaSpectrum").
|
|
///
|
|
/// @param peakpos initial guess of the peak position.
|
|
/// the length of this wave determines the number of peaks to fit.
|
|
///
|
|
/// @return semicolon-separated key=value list of processing parameters
|
|
/// for use with the reduction functions.
|
|
/// position, width and limit parameters are on the x (energy) scale.
|
|
/// @arg rngl low limit of fit interval
|
|
/// @arg rngh high limit of fit interval
|
|
/// @arg pos1 position of peak 1
|
|
/// @arg wid1 width of peak 1
|
|
/// @arg pos2 position of peak 2
|
|
/// @arg wid2 width of peak 2
|
|
/// @arg pos3 position of peak 3
|
|
/// @arg wid3 width of peak 3
|
|
/// @arg pos4 position of peak 4
|
|
/// @arg wid4 width of peak 4
|
|
/// @arg return select result to return, "amp1" ... "amp4"
|
|
/// @arg npeaks number of peaks to fit: 1...4
|
|
/// the others are held at amplitude 0.
|
|
/// @arg ybox box size of averaging in y direction, must be 1 or 3.
|
|
/// other values lead to corrupt data.
|
|
///
|
|
function /s find_gauss4_reduction_params(spectrum, peakpos)
|
|
wave spectrum
|
|
wave peakpos
|
|
string param = ""
|
|
|
|
variable wmin = wavemin(spectrum)
|
|
variable wmax = wavemax(spectrum)
|
|
|
|
// read parameters
|
|
variable rngl = dimoffset(spectrum, 0)
|
|
variable rngh = dimoffset(spectrum, 0) + dimdelta(spectrum, 0) * (dimsize(spectrum, 0) - 1)
|
|
make /n=4 /free positions, widths
|
|
variable npeaks = numpnts(peakpos)
|
|
variable ybox = 1
|
|
positions = 0
|
|
positions[0, npeaks-1] = peakpos[p]
|
|
widths = 0.2
|
|
|
|
variable n_coef = npeaks * 3 + 2
|
|
make /free /d /n=(n_coef) w_coef
|
|
w_coef = 0
|
|
w_coef[0] = wmin
|
|
w_coef[1] = 0
|
|
|
|
make /free /n=(2+npeaks, numpnts(w_coef)) cmat
|
|
make /free /n=(2+npeaks) cvec
|
|
cmat = 0
|
|
cmat[0][0] = -1
|
|
cmat[1][1] = 1
|
|
cvec = 0
|
|
|
|
variable ip
|
|
for (ip=0; ip < npeaks; ip += 1)
|
|
cmat[2 + ip][2 + ip*3] = -1
|
|
w_coef[2 + ip*3] = wmax - wmin
|
|
w_coef[3 + ip*3] = peakpos[ip]
|
|
w_coef[4 + ip*3] = widths[ip]
|
|
endfor
|
|
|
|
variable V_FitNumIters
|
|
FuncFit /Q /NTHR=1 /N MultiGaussLinBG w_coef spectrum /C={cmat, cvec}
|
|
|
|
for (ip=0; ip < npeaks; ip += 1)
|
|
positions[ip] = w_coef[3 + ip * 3]
|
|
widths[ip ] = abs(w_coef[4 + ip * 3])
|
|
endfor
|
|
for (ip=npeaks; ip < 4; ip += 1)
|
|
positions[ip] = 0
|
|
widths[ip] = 0.2
|
|
endfor
|
|
|
|
param = ReplaceNumberByKey("rngl", param, rngl, "=", ";")
|
|
param = ReplaceNumberByKey("rngh", param, rngh, "=", ";")
|
|
param = ReplaceNumberByKey("npeaks", param, npeaks, "=", ";")
|
|
param = ReplaceNumberByKey("ybox", param, ybox, "=", ";")
|
|
param = ReplaceNumberByKey("pos1", param, positions[0], "=", ";")
|
|
param = ReplaceNumberByKey("wid1", param, widths[0], "=", ";")
|
|
param = ReplaceNumberByKey("pos2", param, positions[1], "=", ";")
|
|
param = ReplaceNumberByKey("wid2", param, widths[1], "=", ";")
|
|
param = ReplaceNumberByKey("pos3", param, positions[2], "=", ";")
|
|
param = ReplaceNumberByKey("wid3", param, widths[2], "=", ";")
|
|
param = ReplaceNumberByKey("pos4", param, positions[3], "=", ";")
|
|
param = ReplaceNumberByKey("wid4", param, widths[3], "=", ";")
|
|
|
|
return param
|
|
end
|
|
|
|
/// apply the Shockley_anglefit function to a single image
|
|
///
|
|
/// useful for testing or manual processing
|
|
/// since the Shockley_anglefit function cannot be called from the command line directly.
|
|
///
|
|
/// @param branch -1 or +1: select negative (positive) angles for the fit interval, respectively
|
|
///
|
|
function test_shockley_anglefit(image, branch)
|
|
wave image
|
|
variable branch
|
|
|
|
string param = ""
|
|
param = ReplaceStringByKey("branch", param, num2str(branch), "=", ";")
|
|
|
|
string s_branch
|
|
if (branch >= 0)
|
|
s_branch = "p"
|
|
else
|
|
s_branch = "n"
|
|
endif
|
|
string pkpos_name = "saf_pkpos_" + s_branch
|
|
string pkwid_name = "saf_pkwid_" + s_branch
|
|
|
|
wave /wave results = shockley_anglefit(image, param)
|
|
duplicate results[0], $pkpos_name
|
|
duplicate results[1], $pkwid_name
|
|
end
|
|
|
|
function prompt_Shockley_anglefit(param)
|
|
string ¶m
|
|
|
|
variable branch = NumberByKey("branch", param, "=", ";")
|
|
|
|
prompt branch, "Branch (-1 or +1)"
|
|
|
|
doprompt "Shockley_anglefit_reduction Parameters", branch
|
|
if (v_flag == 0)
|
|
param = ReplaceNumberByKey("branch", param, branch, "=", ";")
|
|
endif
|
|
|
|
return v_flag
|
|
end
|
|
|
|
/// data reduction function for analysing the Cu(111) Shockley surface state.
|
|
///
|
|
/// do curve fitting of one branch of the surface state.
|
|
/// the result is peak position versus energy.
|
|
///
|
|
/// TODO: this function contains hard-coded parameters. please generalize as necessary.
|
|
///
|
|
/// @param source scienta detector image, energy axis along X, angle axis along Y.
|
|
/// two-dimensional intensity distribution (image).
|
|
/// the scales are carried over to the result waves.
|
|
/// for correct weighting and error estimation it is important
|
|
/// that the source wave contains actual counts (Poisson statistics).
|
|
///
|
|
/// @param param (in, out) semicolon-separated key=value list of processing parameters.
|
|
/// this is a pass-by-reference argument.
|
|
/// the following parameters are required.
|
|
/// @arg branch -1 or +1: select negative (positive) angles for the fit interval, respectively
|
|
///
|
|
/// @return free wave containing references of the result waves.
|
|
/// @arg result[0] peak position
|
|
/// @arg result[1] peak width (sigma)
|
|
///
|
|
threadsafe function /wave Shockley_anglefit(source, param)
|
|
wave source
|
|
string ¶m
|
|
|
|
variable nx = dimsize(source, 0)
|
|
variable ny = dimsize(source, 1)
|
|
|
|
// read parameters
|
|
variable branch = NumberByKey("branch", param, "=", ";")
|
|
|
|
// validate parameters
|
|
if (numtype(branch) != 0)
|
|
branch = +1
|
|
endif
|
|
|
|
// prepare output
|
|
make /wave /free /n=2 result_waves
|
|
make /free /n=0 dest1, dest2
|
|
result_waves[0] = dest1
|
|
result_waves[1] = dest2
|
|
adh5_setup_profile(source, dest1, 0)
|
|
adh5_setup_profile(source, dest2, 0)
|
|
dest1 = nan
|
|
dest2 = nan
|
|
|
|
// select angle range
|
|
// hard-coded for a particular measurement series
|
|
variable y0
|
|
variable y1
|
|
if (branch < 0)
|
|
y0 = -5
|
|
y1 = 0
|
|
else
|
|
y0 = 0
|
|
y1 = 5
|
|
endif
|
|
|
|
// select energy range
|
|
// start at the point of highest intensity and go up 0.45 eV
|
|
variable p0
|
|
variable p1
|
|
variable q0
|
|
variable q1
|
|
duplicate /free dest1, center
|
|
q0 = round((y0 - dimoffset(source, 1)) / dimdelta(source, 1))
|
|
q1 = round((y1 - dimoffset(source, 1)) / dimdelta(source, 1))
|
|
ad_profile_x_w(source, q0, q1, center)
|
|
wavestats /q/m=1 center
|
|
p0 = round((v_maxloc - dimoffset(source, 0)) / dimdelta(source, 0))
|
|
p1 = round((v_maxloc + 0.4 - dimoffset(source, 0)) / dimdelta(source, 0))
|
|
|
|
// prepare intermediate data buffer
|
|
make /n=(ny)/d/free profile
|
|
setscale /p x dimoffset(source,1), dimdelta(source,1), waveunits(source,1), profile
|
|
|
|
variable pp
|
|
for (pp = p0; pp <= p1; pp += 1)
|
|
profile = source[pp][p]
|
|
curvefit /Q /NTHR=1 /W=2 gauss profile(y0,y1)
|
|
wave w_coef
|
|
dest1[pp] = w_coef[2]
|
|
dest2[pp] = w_coef[3]
|
|
endfor
|
|
|
|
return result_waves
|
|
end
|
|
|
|
function scienta_norm(w, x): fitfunc
|
|
wave w
|
|
variable x
|
|
|
|
return w[0] * (x^2 - w[1]^2)
|
|
end
|
|
|
|
function /wave fit_scienta_ang_transm(data, params)
|
|
wave data // measured angular distribution (1D)
|
|
wave /z params
|
|
|
|
if (!waveexists(params))
|
|
make /n=12 /o params
|
|
endif
|
|
redimension /n=12/d params
|
|
|
|
variable h = wavemax(data) - wavemin(data)
|
|
params[0] = h / 2
|
|
params[1] = 0
|
|
params[2] = 7
|
|
params[3] = h / 4
|
|
params[4] = -23
|
|
params[5] = 4
|
|
params[6] = h / 4
|
|
params[7] = +23
|
|
params[8] = 4
|
|
params[9] = h / 2
|
|
params[10] = 0
|
|
params[11] = -0.001
|
|
FuncFit /NTHR=0 /q scienta_ang_transm params data
|
|
|
|
return params
|
|
end
|
|
|
|
threadsafe function scienta_ang_transm(w, x): fitfunc
|
|
// parameterized angular transmission function of the analyser
|
|
wave w // coefficients
|
|
// w[0] = amplitude gauss 1
|
|
// w[1] = position gauss 1
|
|
// w[2] = width gauss 1
|
|
// w[3] = amplitude gauss 2
|
|
// w[4] = position gauss 2
|
|
// w[5] = width gauss 2
|
|
// w[6] = amplitude gauss 3
|
|
// w[7] = position gauss 3
|
|
// w[8] = width gauss 3
|
|
// w[9] = constant background
|
|
// w[10] = linear background
|
|
// w[11] = quadratic background
|
|
variable x
|
|
|
|
make /free /n=4 /d w_int
|
|
w_int[0] = 0
|
|
w_int[1,] = w[p - 1]
|
|
variable pk1 = gauss1d(w_int, x)
|
|
w_int[1,] = w[p + 2]
|
|
variable pk2 = gauss1d(w_int, x)
|
|
w_int[1,] = w[p + 5]
|
|
variable pk3 = gauss1d(w_int, x)
|
|
w_int[0,2] = w[9 + p]
|
|
w_int[3] = 0
|
|
variable bg = poly(w_int, x)
|
|
|
|
return bg + pk1 + pk2 + pk3
|
|
end
|
|
|
|
function /wave fit_scienta_poly_bg(data, params, bgterms)
|
|
wave data // measured distribution (2D)
|
|
wave /z params // wave, will be redimensioned for the correct size
|
|
variable bgterms // number of terms in the polynomial background: 2, 3, or 4
|
|
|
|
if (!waveexists(params))
|
|
make /n=15 /o params
|
|
endif
|
|
redimension /n=15 /d params
|
|
|
|
variable wmax = wavemax(data)
|
|
variable wmin = wavemin(data)
|
|
params[0] = 0
|
|
params[1] = 7
|
|
params[2] = 1 / 2
|
|
params[3] = -23
|
|
params[4] = 4
|
|
params[5] = 1 / 2
|
|
params[6] = +23
|
|
params[7] = 4
|
|
params[8] = 1
|
|
params[9] = 0
|
|
params[10] = -0.001
|
|
params[11] = wmin
|
|
params[12] = (wmax - wmin) / dimdelta(data,1) / dimsize(data,1)
|
|
params[13] = 0
|
|
params[14] = 0
|
|
|
|
string h = "0000000000000"
|
|
if (bgterms < 3)
|
|
h = h + "1"
|
|
else
|
|
h = h + "0"
|
|
endif
|
|
if (bgterms < 4)
|
|
h = h + "1"
|
|
else
|
|
h = h + "0"
|
|
endif
|
|
FuncFitMD /NTHR=1 /q /h=h scienta_poly_bg params data
|
|
|
|
return params
|
|
end
|
|
|
|
function scienta_poly_bg(w, e, a): fitfunc
|
|
// polynomial background with
|
|
// parameterized angular transmission function of the analyser
|
|
wave w // coefficients
|
|
// angular transmission, varies with a
|
|
// amplitude of gauss 1 = 1 constant
|
|
// other peak amplitudes and linear terms are relative to gauss 1
|
|
// w[0] = position gauss 1
|
|
// w[1] = width gauss 1
|
|
// w[2] = amplitude gauss 2, relative to gauss 1
|
|
// w[3] = position gauss 2
|
|
// w[4] = width gauss 2
|
|
// w[5] = amplitude gauss 3, relative to gauss 1
|
|
// w[6] = position gauss 3
|
|
// w[7] = width gauss 3
|
|
// w[8] = constant term
|
|
// w[9] = linear term
|
|
// w[10] = quadratic term
|
|
// spectral background, varies with e
|
|
// w[11] = constant term
|
|
// w[12] = linear term
|
|
// w[13] = quadratic term
|
|
// w[14] = cubic term
|
|
variable a // detection angle
|
|
variable e // electron energy
|
|
|
|
make /free /n=4 /d w_int
|
|
variable p0 = 0
|
|
|
|
w_int[0] = 0
|
|
w_int[1] = 1
|
|
w_int[2,] = w[p0 + p - 2]
|
|
variable pk1 = gauss1d(w_int, a)
|
|
p0 += 2
|
|
|
|
w_int[1,] = w[p0 + p - 1]
|
|
variable pk2 = gauss1d(w_int, a)
|
|
p0 += 3
|
|
|
|
w_int[1,] = w[p0 + p - 1]
|
|
variable pk3 = gauss1d(w_int, a)
|
|
p0 += 3
|
|
|
|
w_int[0,2] = w[p0 + p]
|
|
w_int[3] = 0
|
|
variable base = poly(w_int, a)
|
|
p0 += 3
|
|
|
|
w_int[0,3] = w[p0 + p]
|
|
variable bg = poly(w_int, e)
|
|
|
|
return bg * (base + pk1 + pk2 + pk3)
|
|
end
|