new user distribution of PEARL procedures
This commit is contained in:
705
pearl/pearl-scienta-preprocess.ipf
Normal file
705
pearl/pearl-scienta-preprocess.ipf
Normal file
@ -0,0 +1,705 @@
|
||||
#pragma rtGlobals=3 // Use modern global access method and strict wave access.
|
||||
#pragma IgorVersion = 6.1
|
||||
#pragma ModuleName = PearlScientaPreprocess
|
||||
#pragma version = 1.02
|
||||
|
||||
// $Id$
|
||||
// author: matthias.muntwiler@psi.ch
|
||||
// Copyright (c) 2013-14 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-15 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.
|
||||
|
||||
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
|
||||
|
||||
function /s capture_int_linbg_cursors()
|
||||
// this function is for testing only, until we implement a proper interface
|
||||
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)
|
||||
// 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.
|
||||
|
||||
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
|
||||
|
||||
function test_int_linbg(image)
|
||||
// useful for testing or manual processing
|
||||
// since the int_linbg_reduction function cannot be called from the command line directly.
|
||||
wave image
|
||||
|
||||
string param = ""
|
||||
prompt_int_linbg_reduction(param)
|
||||
|
||||
string data1_name = "test_data1"
|
||||
string data2_name = "test_data2"
|
||||
duplicate /o image, $data1_name, $data2_name
|
||||
wave w_data1 = $data1_name
|
||||
wave w_data2 = $data2_name
|
||||
|
||||
int_linbg_reduction(image, w_data1, w_data2, param)
|
||||
end
|
||||
|
||||
threadsafe function int_linbg_reduction(source, dest1, dest2, param)
|
||||
// data reduction function for adh5_load_reduced_detector
|
||||
// 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.
|
||||
wave source // source wave
|
||||
// Scienta detector image, energy axis along X, angle axis along Y
|
||||
wave dest1, dest2 // destination waves
|
||||
// each wave is a one-dimensional intensity distribution
|
||||
// the function may redimension these waves to one of the image dimensions
|
||||
// (it must be clear to the user which dimension this is).
|
||||
// the meaning of dest1 and dest2 is up to the particular function,
|
||||
// e.g. dest1 could hold the mean value and dest2 the one-sigma error,
|
||||
// or dest1 could hold the X-profile, and dest2 the Y-profile.
|
||||
string ¶m // parameters in a key1=value1;key2=value2;... list
|
||||
// all region parameters are relative to the image size (0...1)
|
||||
// Lcrop = size of the lower cropping region
|
||||
// Hcrop = size of the upper cropping region
|
||||
// Lsize = size of the lower background integration region
|
||||
// Hsize = size of the upper background integration region
|
||||
// Cpos = center position of the of the peak integration region
|
||||
// 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
|
||||
|
||||
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, "=", ";")
|
||||
|
||||
// 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)
|
||||
return 1 // Cpos parameter missing
|
||||
endif
|
||||
if (numtype(Csize) != 0)
|
||||
return 2 // Csize parameter missing
|
||||
endif
|
||||
|
||||
variable lpos = lcrop + lsize / 2
|
||||
variable hpos = 1 - (hcrop + hsize / 2)
|
||||
|
||||
variable p0
|
||||
variable p1
|
||||
|
||||
adh5_setup_profile(source, dest1, 1)
|
||||
adh5_setup_profile(source, dest2, 1)
|
||||
|
||||
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))
|
||||
return 0 // return zero if successful, non-zero if an error occurs
|
||||
end
|
||||
|
||||
function test_shockley_anglefit(image, branch)
|
||||
// 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.
|
||||
wave image
|
||||
variable branch // +1 or -1
|
||||
|
||||
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
|
||||
duplicate /o image, $pkpos_name, $pkwid_name
|
||||
wave w_pkpos = $pkpos_name
|
||||
wave w_pkwid = $pkwid_name
|
||||
|
||||
shockley_anglefit(image, w_pkpos, w_pkwid, param)
|
||||
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
|
||||
|
||||
threadsafe function Shockley_anglefit(source, dest1, dest2, param)
|
||||
// data reduction function for adh5_load_reduced_detector
|
||||
// specialized 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.
|
||||
wave source // source wave
|
||||
// Scienta detector image, energy axis along X, angle axis along Y
|
||||
// the apex of the surface state must be at angle 0
|
||||
wave dest1, dest2 // destination waves
|
||||
// dest1: peak position
|
||||
// dest2: peak width (sigma)
|
||||
string ¶m // parameters in a key1=value1;key2=value2;... list
|
||||
// branch=-1 or +1: select negative (positive) angles for the fit interval, respectively
|
||||
|
||||
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
|
||||
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 0 // return zero if successful, non-zero if an error occurs
|
||||
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
|
||||
|
||||
function test_int_quadbg(image)
|
||||
// useful for testing or manual processing
|
||||
// since the int_quadbg_reduction function cannot be called from the command line directly.
|
||||
wave image
|
||||
|
||||
string param = ""
|
||||
prompt_int_quadbg_reduction(param)
|
||||
|
||||
string data1_name = "test_data1"
|
||||
string data2_name = "test_data2"
|
||||
duplicate /o image, $data1_name, $data2_name
|
||||
wave w_data1 = $data1_name
|
||||
wave w_data2 = $data2_name
|
||||
|
||||
int_quadbg_reduction(image, w_data1, w_data2, param)
|
||||
end
|
||||
|
||||
threadsafe function int_quadbg_reduction(source, dest1, dest2, param)
|
||||
// data reduction function for adh5_load_reduced_detector
|
||||
// integrates peak area minus a quadratic backgrouind
|
||||
wave source // source wave
|
||||
// Scienta detector image, energy axis along X, angle axis along Y
|
||||
wave dest1, dest2 // destination waves
|
||||
string ¶m // parameters in a key1=value1;key2=value2;... list
|
||||
// all region parameters are relative to the image size (0...1)
|
||||
// Lcrop = size of the lower cropping region
|
||||
// Hcrop = size of the upper cropping region
|
||||
// Lsize = size of the lower background integration region
|
||||
// Hsize = size of the upper background integration region
|
||||
// Cpos = center position of the of the peak integration region
|
||||
// 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
|
||||
|
||||
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, "=", ";")
|
||||
|
||||
// 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)
|
||||
return 1 // Cpos parameter missing
|
||||
endif
|
||||
if (numtype(Csize) != 0)
|
||||
return 2 // 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)
|
||||
|
||||
adh5_setup_profile(source, dest1, 0)
|
||||
adh5_setup_profile(source, dest2, 0)
|
||||
|
||||
// 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
|
||||
|
||||
return 0 // return zero if successful, non-zero if an error occurs
|
||||
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
|
||||
Reference in New Issue
Block a user