fermi-edge-analysis.ipf

Go to the documentation of this file.

#pragma rtGlobals=3     // Use modern global access method and strict wave access.
#include "pearl-area-profiles"

// $Id$



function analyse_curved_edge(data)
    wave data   // 2D counts data, x-scale = kinetic energy, y-scale = analyser angle
    variable guess_EF   // initial guess of the Fermi level
    
    dfref savedf = GetDataFolderDFR()
    dfref datadf = GetWavesDataFolderDFR(data)
    setdatafolder datadf
    
    make /n=5 /d /o w_coef_int
    make /n=7 /d /o w_coef_curved

    // 1) integrate data
    wave xint = ad_profile_x(data, -inf, inf, "", noavg=1)
    duplicate /free xint, xint_sig
    xint_sig = sqrt(xint * 4) / 4
    variable xmin = wavemin(xint)
    variable xmax = wavemax(xint)
    variable xmean = mean(xint)

    // 2) fit integrated data for reference
    w_coef_int[0] = xmin
    w_coef_int[1] = 0
    w_coef_int[2] = xmax - xmin
    w_coef_int[3] = dimoffset(xint, 0) + dimdelta(xint, 0) * dimsize(xint, 0) / 2
    w_coef_int[4] = 100
    FuncFit /NTHR=0 FermiFuncLinDOS w_coef_int xint /D /I=1 /W=xint_sig
    wave w_sigma
    duplicate /o w_sigma, w_sigma_int
    
    // 3) normalize data
    wave yavg = ad_profile_y(data, -inf, inf, "")
    variable ymean = mean(yavg)

    duplicate /o data, data_norm, data_sig
    data_sig = sqrt(data * 4) / 4
    data_norm = data_norm * ymean / yavg[q]
    data_sig = data_sig * ymean / yavg[q]
    
    // 4) fit normalized data
    wave xavg = ad_profile_x(data, -inf, inf, "")
    xmin = wavemin(xavg)
    xmax = wavemax(xavg)
    xmean = mean(xavg)
    
    w_coef_curved[0] = {0, 0, 1, 95.5, 100, 0, -0.0001}
    w_coef_curved[0] = xmin
    w_coef_curved[2] = xmax - xmin
    w_coef_curved[3] = w_coef_int[3]
    w_coef_curved[4] = w_coef_int[4]
    
    //variable box = min(11, numpnts(xavg) / 5)
    //FindLevel /B=(box) /EDGE=2 /Q xavg, (xmin + xmax) / 2
    //if (v_flag == 0)
    //  w_coef_curved[3] = v_levelx
    //else
    //  w_coef_curved[3] =  dimoffset(data, 0) + dimdelta(data, 0) * dimsize(data, 0) / 2
    //endif
    
    duplicate /o data_norm, fit_data_norm
    FuncFitMD /X=1 /NTHR=0 FermiFuncLinDOS_2Dcorr w_coef_curved  data_norm /D=fit_data_norm /I=1 /W=data_sig
    wave w_sigma
    duplicate /o w_sigma, w_sigma_curved

    display /k=1; appendimage data_norm
    ModifyImage data_norm ctab= {xmin,xmax,Grays,0}
    AppendMatrixContour fit_data_norm
    
    setdatafolder savedf
    return 0
end

function record_results(index)
    variable index
    
    dfref savedf = GetDataFolderDFR()
    wave /sdfr=root: Tint, Tint_sig, Tcurv, Tcurv_sig, EFcurv, EFcurv_sig
    wave w_coef_int, w_coef_curved
    wave w_sigma_int, w_sigma_curved
    
    Tint[index] = w_coef_int[4]
    Tint_sig[index] = w_sigma_int[4]
    Tcurv[index] = w_coef_curved[4]
    Tcurv_sig[index] = w_sigma_curved[4]
    EFcurv[index] = w_coef_curved[3]
    EFcurv_sig[index] = w_sigma_curved[3]
    
    setdatafolder savedf
end

function integrate_curved_edge(data, data_sig)
    wave data
    wave /z data_sig
    //wave coef
    
    string name = nameofwave(data)
    
    duplicate /o data, data_out
    redimension /n=(dimsize(data,0)) data_out
    data_out = 0
    
    if (waveexists(data_sig))
        duplicate /o data_sig, sig_out
        redimension /n=(dimsize(data,0)) sig_out
        sig_out = 0
    endif

    wave ywgt = ad_profile_y(data, -inf, inf, "", noavg=1)
    ywgt = 1 / ywgt / 4
    //ywgt = 1 / 4
    variable sum_ywgt = sum(ywgt)
    
    variable nx = dimsize(data, 0)
    variable ny = dimsize(data, 1)
    variable iy
    variable yy
    variable dx
    variable dy
    variable dp
    variable dp_min = 0
    
    wave PassEnergy = :attr:PassEnergy
    wave NumSlices = :attr:NumSlices
    
    for (iy = 0; iy < ny; iy += 1)
        dy = dimoffset(data, 1) + dimdelta(data, 1) * iy
        dy = dy / dimdelta(data, 1)
        dx = slit_shift(dy * 902 / NumSlices[0], PassEnergy[0])
        dp = round(dx / dimdelta(data, 0))  // <= 0
        dp_min = min(dp_min, dp)
        
        data_out[] = p+dp >= 0 ? data_out + data[p+dp][iy] * ywgt[iy] : data_out

        if (waveexists(data_sig))
            sig_out = p+dp >= 0 ? sig_out + data_sig[p+dp][iy] * ywgt[iy] : sig_out
        endif
    endfor
    
    data_out /= sum_ywgt
    data_out[0, -dp_min][] = nan

    if (waveexists(sig_out))
        sig_out /= sum_ywgt
        //sig_out[0, -dp_min] = nan
    endif
end

function slit_correction(data, data_out, epass)
    wave data   // must be image with original dimensions (no cropping!)
    wave /z data_out    // 2D or 1D wave to receive the result
        // X dimension must be identical to the one of data
        // if 2D, Y dimension must be either identical to the one of data
        // if 1D, the result will be the sum of the corrected slices
    variable epass  // pass energy
    
    if (!WaveExists(data_out))
        string name = nameofwave(data) + "_corr"
        duplicate /o data, $name
        wave data_out = $name
        //redimension /n=(dimsize(data,0)) data_out
    endif
    data_out = 0
    
    variable nx = dimsize(data, 0)
    variable ny = dimsize(data, 1)
    variable iy
    variable yy
    variable dx
    variable dy
    variable dp
    variable dp_min = 0
    
    for (iy = 0; iy < ny; iy += 1)
        dy = dimoffset(data, 1) + dimdelta(data, 1) * iy
        dy = dy / dimdelta(data, 1)
        dx = slit_shift(dy * 902 / ny, epass)
        dp = round(dx / dimdelta(data, 0))  // <= 0
        dp_min = min(dp_min, dp)
        
        if (wavedims(data_out) >= 2)
            data_out[][iy] = p+dp >= 0 ? data_out + data[p+dp][iy] : data_out
        else
            data_out = p+dp >= 0 ? data_out + data[p+dp][iy] : data_out
        endif
    endfor
    
    data_out[0, -dp_min][] = nan
end

//------------------------------------------------------------------------------
threadsafe Function FermiFuncLinDOS2D_corr(w,x,y) : FitFunc
// linear density of states below Fermi level
// 2D data with corrections:
// - straight slit (slit shift)
// - transmission function (polynomial)
//------------------------------------------------------------------------------
    Wave w; Variable x; variable y
    // w[0] = background far above the fermi level
    // w[1] = slope of the linear background
    // w[2] = amplitude
    // w[3] = fermi level in eV
    // w[4] = temperature in K
    // w[5] = transmission - linear term
    // w[6] = transmission - quadratic term
    // w[7] = transmission - cubic term
    // w[8] = pass energy (hold this value)
    // w[9] = y scale: pixels / unit of y

    variable pos = w[3] + slit_shift(y * w[9], w[8])
    variable transm = 1 + w[5] * y + w[6] * y^2 + w[7] * y^3
    variable fermi = (w[1] * min(x - pos, 0) + w[2]) / ( exp( (x - pos)  /  (kBoltzmann * w[4]) ) + 1.0  )
    return transm * (fermi + w[0])
end

//------------------------------------------------------------------------------
threadsafe Function FermiFuncLinDOS_2Dcorr_old(w,x,y) : FitFunc
// linear density of states below Fermi level
// 2D data with a polynomial shift of the Fermi level in the second dimension
//------------------------------------------------------------------------------
    Wave w; Variable x; variable y
    // w[0] = background far above the fermi level
    // w[1] = slope of the linear background
    // w[2] = amplitude
    // w[3] = fermi level in eV
    // w[4] = temperature in K
    // w[5] = shift - linear term
    // w[6] = shift - quadratic term

    variable pos = w[3] + w[5] * y + w[6] * y^2
    return w[0] + (w[1] * min(x - pos, 0) + w[2]) / ( exp( (x - pos)  /  (kBoltzmann * w[4]) ) + 1.0  )
end

static constant mcp_radius_pix = 555
static constant mcp_radius_mm = 20
static constant hemi_radius_mm = 200
static constant mcp_radius_epass = 0.04

threadsafe function slit_shift(ypix, epass)
    // calculates the energy shift due to straight slit
    // the radius of the curve is 1/2 of the hemisphere radius
    variable ypix   // vertical (angle/position) pixels distance from center
        // = slice coordinate * 902 / number of slices
    variable epass // pass energy
    
    variable rpix = mcp_radius_pix * hemi_radius_mm / mcp_radius_mm
    //variable rpix = mcp_radius_pix * hemi_radius_mm / 2 / mcp_radius_mm
    variable rene = epass * mcp_radius_epass * hemi_radius_mm / 2 / mcp_radius_mm
    
    variable isin = asin(ypix / rpix)
    variable dene = rene * (cos(isin) - 1)
    
    return dene
end

function show_shift(data)
    wave data
    variable epass
    
    variable ny = dimsize(data, 1)
    make /o /n=(ny) shift_x
    setscale /i x -ny/2, ny/2, "", shift_x
    wave PassEnergy = :attr:PassEnergy
    wave NumSlices = :attr:NumSlices
    
    shift_x = slit_shift(x * 902 / NumSlices[0], PassEnergy[0])
    
    variable ecenter = dimoffset(data, 0) + dimdelta(data, 0) * dimsize(data, 0) / 2
    shift_x += ecenter
end

function analyser_energy_resolution(epass, slit)
    variable epass
    variable slit

    variable respow
    if (epass < 4)
        respow = 1110
    elseif (epass < 8)
        respow = 1400
    else
        respow = 1750
    endif
    
    return epass * max(0.2, slit) / 0.2 / respow    
end