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
FermiFuncLinDOS2D_corr
threadsafe variable FermiFuncLinDOS2D_corr(variable w, threadsafe x, wave y)
Definition: fermi-edge-analysis.ipf:216
slit_shift
threadsafe variable slit_shift(variable ypix, variable epass)
Definition: fermi-edge-analysis.ipf:267
mcp_radius_pix
static const variable mcp_radius_pix
MCP radius seen by the camera in pixels.
Definition: fermi-edge-analysis.ipf:259
mcp_radius_mm
static const variable mcp_radius_mm
physical size (radius) of the MCP in mm
Definition: fermi-edge-analysis.ipf:261
ad_profile_y
threadsafe wave ad_profile_y(wave dataset, variable p1, variable p2, string destname, variable noavg=defaultValue)
1D cut through 2D dataset along Y dimension, new destination wave.
Definition: pearl-area-profiles.ipf:542
analyse_curved_edge
variable analyse_curved_edge(wave data)
Definition: fermi-edge-analysis.ipf:27
hemi_radius_mm
static const variable hemi_radius_mm
physical size (radius) of the hemisphere in mm
Definition: fermi-edge-analysis.ipf:263
analyser_energy_resolution
variable analyser_energy_resolution(variable epass, variable slit)
calculate the energy resolution of the analyser
Definition: fermi-edge-analysis.ipf:306
slit_correction
variable slit_correction(wave data, wave data_out, variable epass)
Definition: fermi-edge-analysis.ipf:173
mcp_radius_epass
static const variable mcp_radius_epass
energy range imaged on MCP relative to the pass energy
Definition: fermi-edge-analysis.ipf:265
integrate_curved_edge
variable integrate_curved_edge(wave data, wave data_sig)
Definition: fermi-edge-analysis.ipf:116
show_shift
variable show_shift(wave data)
Definition: fermi-edge-analysis.ipf:284
ad_profile_x
threadsafe wave ad_profile_x(wave dataset, variable q1, variable q2, string destname, variable noavg=defaultValue)
1D cut through 2D dataset along X dimension, new destination wave.
Definition: pearl-area-profiles.ipf:480
record_results
variable record_results(variable index)
Definition: fermi-edge-analysis.ipf:98