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