pearl-scienta-preprocess.ipf
Go to the documentation of this file.
#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
function prompt_int_linbg_reduction(param)
string &param
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
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
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
threadsafe function /wave int_linbg_reduction(source, param)
wave source
string &param
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 &param
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
threadsafe function /wave int_quadbg_reduction(source, param)
wave source
string &param
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
function prompt_redim_linbg_reduction(param)
string &param
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
threadsafe function /wave redim_linbg_reduction(source, param)
wave source
string &param
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
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
function prompt_gauss4_reduction(param)
string &param
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
threadsafe function /wave gauss4_reduction(source, param)
wave source
string &param
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
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
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 &param
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
threadsafe function /wave Shockley_anglefit(source, param)
wave source
string &param
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
MultiGaussLinBG
threadsafe variable MultiGaussLinBG(wave w, variable x)
multiple gaussian peaks on a linear background fit function.
Definition: pearl-fitfuncs.ipf:44
MultiGaussLinBG_AO
threadsafe variable MultiGaussLinBG_AO(wave pw, wave yw, wave xw)
multiple gaussian peaks on a linear background fit function (all at once).
Definition: pearl-fitfuncs.ipf:79
prompt_gauss4_reduction
variable prompt_gauss4_reduction(string *param)
prompt for the gauss4_reduction parameters
Definition: pearl-scienta-preprocess.ipf:580
redim_linbg_reduction
threadsafe wave redim_linbg_reduction(wave source, string *param)
linear background reduction function for incorrectly dimensioned scienta image
Definition: pearl-scienta-preprocess.ipf:526
test_gauss4_reduction
variable test_gauss4_reduction(wave image)
apply the gauss4_reduction function to a single image
Definition: pearl-scienta-preprocess.ipf:546
ad_profile_x_w
threadsafe wave ad_profile_x_w(wave dataset, variable q1, variable q2, wave destwave, variable noavg=defaultValue)
1D cut through 2D dataset along X dimension, existing destination wave.
Definition: pearl-area-profiles.ipf:504
adh5_setup_profile
threadsafe variable adh5_setup_profile(wave image, wave profile, variable dim)
set up a one-dimensional wave for a line profile based on a 2D original wave.
Definition: pearl-area-import.ipf:1133
prompt_int_quadbg_reduction
variable prompt_int_quadbg_reduction(string *param)
Definition: pearl-scienta-preprocess.ipf:305
int_linbg_reduction
threadsafe wave int_linbg_reduction(wave source, string *param)
linear-background subtracted integration reduction function.
Definition: pearl-scienta-preprocess.ipf:214
gauss4_reduction
threadsafe wave gauss4_reduction(wave source, string *param)
fit horizontal cuts of an image with up to four gaussian peaks on a linear background ...
Definition: pearl-scienta-preprocess.ipf:672
prompt_redim_linbg_reduction
variable prompt_redim_linbg_reduction(string *param)
parameter dialog for the redim_linbg_reduction() function
Definition: pearl-scienta-preprocess.ipf:463
ad_profile_y_w
threadsafe wave ad_profile_y_w(wave dataset, variable p1, variable p2, wave destwave, variable noavg=defaultValue)
1D cut through 2D dataset along X dimension, existing destination wave.
Definition: pearl-area-profiles.ipf:567
int_quadbg_reduction
threadsafe wave int_quadbg_reduction(wave source, string *param)
integrate peak area minus a quadratic background
Definition: pearl-scienta-preprocess.ipf:367
prompt_int_linbg_reduction
variable prompt_int_linbg_reduction(string *param)
prompt the user for integrate on linear background reduction parameters.
Definition: pearl-scienta-preprocess.ipf:35
csr_int_linbg_reduction
string csr_int_linbg_reduction(string win)
set reduction parameters from cursors in a graph.
Definition: pearl-scienta-preprocess.ipf:89
capture_int_linbg_cursors
string capture_int_linbg_cursors()
this function is for testing only, until we implement a proper interface
Definition: pearl-scienta-preprocess.ipf:67