Added Gauss Peak Fitting to XPS shift determination. First implementation for positionWidth contrained fit. Not yet documented. Not yet fully tested. WIP

This commit is contained in:
bartels-rausch committed 2023-11-09 07:37:10 +01:00
1 parent b23894e06c
commit 1aa130982f
2 files changed
+123 -69

No files matched your search

+90 -59
View File
@@ -1,6 +1,6 @@
function DataTable = nappDetermineShiftXps(DataTable,varargin)
function DataTable = nappDetermineShiftXps(DataTable, varargin)
%nappDetermineShiftXps determines the relative shift of each sequence
%iteration in images or the absolute shift of spectra to a binding energy reference
%iteration in images or the absolute shift of spectra to a binding energy reference
%
% FileList = nappDetermineShiftXps(FileList,eKinStart,eKinEnd,'tableSelection',tmpDataSelection,...
% 'showFigures','off','tableVariable','spectrum_countsPerSecond','binEnReference',535);
@@ -9,12 +9,12 @@ function DataTable = nappDetermineShiftXps(DataTable,varargin)
% 'showFigures','on','tableVariable','image');
%
% *input*
% FileList: data table
% FileList: data table
% eKinStart (optional): smallest value of energy range to use for
% peak fitting. Peak fitting is used to get position of peak
% eKinEnd (optional): largest value of energy energy range to use for
% peak fitting. Peak fitting is used to get position of peak
% tableSelection (optional): logical array to filter data table
% tableSelection (optional): logical array to filter data table
% tableVariable: choose to process either
% 'spectrum_counts','spectrum_countsPerSecond' (default), => absolute binding energy shift will be determined or
% ,'image','imageRoh'=> relative kinetic energy shift will be determined.
@@ -25,22 +25,19 @@ function DataTable = nappDetermineShiftXps(DataTable,varargin)
% minPeakHeight: findpeak parameter, default 0.1
% movMeanWindow: findpeak parameter, default 1
% showFigures, on(default): show fit in figure, off do not show fit in figure
%
%
%
% *output*
% new table variable kineticEnergyShift for images
% new table variable kineticEnergyShift for images
% new table variable bindingEnergyShift for spectra
%
%
% *ToDo*
%
% FIXME: make use of O1s fit more general => use nappFitGaussXps.m function
% FIXME: Do not use peakfind for image fit => use nappFitGaussXps.m function
%
% November 2023: Use nappFitGaussXps.m function to find peak position
%
% See also nappnappAlignImages.m nappBindingEnergyXps nappFitGaussXps.m nappCopyShift.m
defaultKinEnReference = 1;
defaultKinEnReference = 1;
defaultTableVariable = 'spectrum_countsPerSecond';
expectedTableVariable = {'spectrum_counts','spectrum_countsPerSecond','image','imageRoh'};
defaultTableSelection = true(size(DataTable,1));
@@ -49,6 +46,10 @@ defaultMinPeakHeight = 0.1;
defaultMovMeanWindow = 1;
defaultShowFigures = 'off';
expectedShowFigures = {'on','off'};
defaultMode = 'gaussFit';
expectedMode = {'gaussFit','findPeak'};
defaultFitMode = 'widthPosition';
expectedFitMode = {'unconstrained','widthPosition','position','width','unconstrainedAndWidthPosition'};
%%
@@ -58,6 +59,7 @@ p = inputParser;
addRequired(p,'DataTable',@istable);
addOptional(p,'eKinStart',@isnumeric)
addOptional(p,'eKinEnd',@isnumeric)
addParameter(p,'FitParameter',@isstruct);
addParameter(p,'binEnReference',@isnumeric)
addParameter(p,'kinEnReference',defaultKinEnReference,@isnumeric)
addParameter(p,'tableSelection',defaultTableSelection,@islogical);
@@ -68,8 +70,11 @@ addParameter(p,'tableVariable',defaultTableVariable,...
@(x) any(validatestring(x,expectedTableVariable)));
addParameter(p,'showFigures',defaultShowFigures,...
@(x) any(validatestring(x,expectedShowFigures)));
addParameter(p,'mode',defaultMode, @(x) any(validatestring(x,expectedMode)));
addParameter(p,'fitMode',defaultFitMode, @(x) any(validatestring(x,expectedFitMode)));
parse(p,DataTable,varargin{:});
% check input (all elements of input variable need to have non-zero size)
if ~all(cell2mat(cellfun(@size,p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable),'UniformOutput',false)),'all')
disp('The table variable appears to be empty. I quit now')
@@ -131,7 +136,7 @@ if matches(p.Results.tableVariable,'image')
warning('There is something wrong with the eKinStart and/or eKinEnd parameters');
disp('There is something wrong with the eKinStart and/or eKinEnd parameters');
end
for ij = 1:size(p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii},2)
[B,x1,y1] = GaussFitOnePeak([p.Results.DataTable(p.Results.tableSelection,:).scientaEkin_eV{ii}(start:stop)',...
p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii}(start:stop,ij)]);
@@ -162,14 +167,14 @@ if matches(p.Results.tableVariable,'image')
h = plot(x1,yData_pk2+bkgrd);
h.LineWidth=1.5;
h.Color = PsiFarbe.gelb;
legend OFF
xlabel('Kinetic Energy [eV]');
ylabel('PE intensity [-]');
grid on
% List results in table
disp('Results');
disp(' ' );
disp(cell2table({sprintf('%.1f',B(3)),...
@@ -181,8 +186,8 @@ if matches(p.Results.tableVariable,'image')
end
end
%%
DataTable(p.Results.tableSelection,:).kineticEnergyShift{ii} = tmp_1-tmp_1(p.Results.kinEnReference);
clearvars tmp_1
fh = findobj( 'Type', 'Figure', 'Name', 'Gauss Multi Fit ');
@@ -193,7 +198,7 @@ if matches(p.Results.tableVariable,'image')
elseif matches(p.Results.tableVariable,'imageRoh')
% repeat this for all columns in input.
for ii = 1:size(p.Results.DataTable(p.Results.tableSelection,:),1)
%%
try
if isnumeric(p.Results.eKinStart) && isempty(DataTable(p.Results.tableSelection,:).scientaEkinRoh_eV{ii})
@@ -225,7 +230,7 @@ elseif matches(p.Results.tableVariable,'imageRoh')
[B,x1,y1] = GaussFitOnePeak([DataTable(p.Results.tableSelection,:).scientaEkinRoh_eV{ii}(start:stop)',...
p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii}(start:stop,ij)]);
end
tmp_1(ij) = p.Results.DataTable(p.Results.tableSelection,:).xRayEkin_eV{ii}-B(3);
% plot result
if matches(p.Results.showFigures,'on')
@@ -253,14 +258,14 @@ elseif matches(p.Results.tableVariable,'imageRoh')
h = plot(x1,yData_pk2+bkgrd);
h.LineWidth=1.5;
h.Color = PsiFarbe.gelb;
legend OFF
xlabel('Kinetic Energy [eV]');
ylabel('PE intensity [-]');
grid on
% List results in table
disp('Results');
disp(' ' );
disp(cell2table({sprintf('%.1f',B(3)),...
@@ -272,7 +277,7 @@ elseif matches(p.Results.tableVariable,'imageRoh')
end
end
%%
DataTable(p.Results.tableSelection,:).kineticEnergyShift{ii} = tmp_1-tmp_1(p.Results.kinEnReference);
clearvars tmp_1
fh = findobj( 'Type', 'Figure', 'Name', 'Gauss Multi Fit ');
@@ -282,43 +287,69 @@ elseif matches(p.Results.tableVariable,'imageRoh')
end
else
% repeat this for all columns in input.
for ii = 1:size(p.Results.DataTable(p.Results.tableSelection,:),1)
%%
try
if isnumeric(p.Results.eKinStart)
start = find(p.Results.DataTable(p.Results.tableSelection,:).scientaEkin_eV{ii}...
>= p.Results.eKinStart,1);
stop = find(p.Results.DataTable(p.Results.tableSelection,:).scientaEkin_eV{ii}...
>= p.Results.eKinEnd,1);
else
start = 1;
stop = length(p.Results.DataTable(p.Results.tableSelection,:).scientaEkin_eV{ii});
if matches(p.Results.mode,'findPeak')
for ii = 1:size(p.Results.DataTable(p.Results.tableSelection,:),1)
%%
try
if isnumeric(p.Results.eKinStart)
start = find(p.Results.DataTable(p.Results.tableSelection,:).scientaEkin_eV{ii}...
>= p.Results.eKinStart,1);
stop = find(p.Results.DataTable(p.Results.tableSelection,:).scientaEkin_eV{ii}...
>= p.Results.eKinEnd,1);
else
start = 1;
stop = length(p.Results.DataTable(p.Results.tableSelection,:).scientaEkin_eV{ii});
end
catch
warning('There is something wrong with the eKinStart and/or eKinEnd parameters');
disp('There is something wrong with the eKinStart and/or eKinEnd parameters');
end
catch
warning('There is something wrong with the eKinStart and/or eKinEnd parameters');
disp('There is something wrong with the eKinStart and/or eKinEnd parameters');
end
if contains(p.Results.showFigures,'on')
figure
findpeaks(movmean(p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii}(start:stop),...
if contains(p.Results.showFigures,'on')
figure
findpeaks(movmean(p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii}(start:stop),...
p.Results.movMeanWindow),...
p.Results.DataTable(p.Results.tableSelection,:).scientaEkin_eV{ii}(start:stop),...
'MinPeakWidth',p.Results.minPeakWidth,'MinPeakHeight',p.Results.minPeakHeight)
end
[pks,locs]=findpeaks(movmean(p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii}(start:stop),...
p.Results.movMeanWindow),...
p.Results.DataTable(p.Results.tableSelection,:).scientaEkin_eV{ii}(start:stop),...
'MinPeakWidth',p.Results.minPeakWidth,'MinPeakHeight',p.Results.minPeakHeight)
'MinPeakWidth',p.Results.minPeakWidth,'MinPeakHeight',p.Results.minPeakHeight);
%%
[~,I] = max(pks);
tmp_1(ii,1) = p.Results.DataTable(p.Results.tableSelection,:).xRayEkin_eV{ii}-locs(I)-p.Results.binEnReference;
clearvars tmp_2
end
[pks,locs]=findpeaks(movmean(p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii}(start:stop),...
p.Results.movMeanWindow),...
p.Results.DataTable(p.Results.tableSelection,:).scientaEkin_eV{ii}(start:stop),...
'MinPeakWidth',p.Results.minPeakWidth,'MinPeakHeight',p.Results.minPeakHeight);
DataTable.bindingEnergyShift(p.Results.tableSelection) = tmp_1;
elseif matches(p.Results.mode,'gaussFit')
tmpFitResults = nappFitGaussXps(DataTable,p.Results.tableSelection,p.Results.FitParameter,'tableSelection',p.Results.tableSelection,...
'energyScale','kinetic','mode',p.Results.fitMode,'showFigures','on');
DataTable.bindingEnergyShift(p.Results.tableSelection) = tmpFitResults(p.Results.tableSelection,:).xRayEkinRange_eV(:,1) - tmpFitResults(p.Results.tableSelection,:).peakPosition(:,1) - repmat(p.Results.binEnReference,size(DataTable.bindingEnergyShift(p.Results.tableSelection)));
%%
[~,I] = max(pks);
tmp_1(ii,1) = p.Results.DataTable(p.Results.tableSelection,:).xRayEkin_eV{ii}-locs(I)-p.Results.binEnReference;
clearvars tmp_2
end
DataTable.bindingEnergyShift(p.Results.tableSelection) = tmp_1;
end
try
fh = findobj( 'Type', 'Figure', 'Name', 'Gauss Multi Fit ');
close(fh)
catch
end
try
fh = findobj( 'Type', 'Figure', 'Name', 'GlobalSearch');
close(fh)
catch
end
end
@@ -329,7 +360,7 @@ function [B,x1,y1] = GaussFitOnePeak(data1)
% make data slimmer so that background makes sense
[~,I] = max(y1);
tmp1 = find(x1>=x1(I)-5,1);
tmp1 = find(x1>=x1(I)-8,1);
tmp2 = find(x1>=x1(I)+5,1);
if isempty(tmp2)
y1 = y1(tmp1:end);
@@ -354,20 +385,20 @@ yfcn1 = @(b,x)...
+ b(4)/(b(5)*sqrt(pi/2))*exp(-2.*((x(:,1)- b(6)).^2/b(5).^2))... % Gauss peak 2: C-OOH, Ac-
+b(7).*x(:,1)+b(8); % linear background
B0 = [80000 1.5 x1(I) ...
10 1.0 x1(I)-4 ...
0 0];
B0 = [1e5 1.5 x1(I) ...
1e5 1.0 x1(I)-4 ...
1 0];
% lower bonds
lb = [1000 0.5 x1(I)-1. ...
10 1.0 x1(I)-5 ...
lb = [100 0.5 x1(I)-1. ...
100 1.0 x1(I)-5 ...
-Inf -Inf];
% upper bond
ub = [Inf 3 x1(I)+1 ...
10 1.0 x1(I)-3 ...
Inf 5 x1(I)-1 ...
Inf Inf];
+33 -10
View File
@@ -21,7 +21,8 @@ function DataTable = nappFitGaussXps(DataTable,~,varargin)
% tmpFitParameters.peakAreaLowerBond = [#1 #2 #3 #4];
% tmpFitParameters.peakWidthLowerBond = [#1 #2 #3 #4];
% tmpFitParameters.peakPosLowerBond = #1;
%
% spectrum range
% tmpFitParameters.backgroundRange = [5 5];
%
% % „unconstrained“ mode
% tmpFitParameters.peakArea = [abs#1 abs#2 abs#3 abs#4];
@@ -34,10 +35,22 @@ function DataTable = nappFitGaussXps(DataTable,~,varargin)
% tmpFitParameters.peakAreaLowerBond = [#1 #2 #3 #4];
% tmpFitParameters.peakWidthLowerBond = [#1 #2 #3 #4];
% tmpFitParameters.peakPosLowerBond = [#1 #2 #3 #4];
%
% spectrum range
% tmpFitParameters.backgroundRange = [5 5];
%
% % „ widthPosition“ mode
% 1.5
% tmpFitParameters.peakArea = [abs#1 abs#2 abs#3 abs#4 abs#5];
% tmpFitParameters.peakWidth = [abs#1-#5];
% tmpFitParameters.peakPos = [abs#1 rel#2 rel#3 rel#4 rel#5];
%
% tmpFitParameters.peakAreaUpperBond = [#1 #2 #3 #4 #5];
% tmpFitParameters.peakWidthUpperBond = abs#1-#5;
% tmpFitParameters.peakPosUpperBond = #1;
% tmpFitParameters.peakAreaLowerBond = [#1 #2 #3 #4 #5];
% tmpFitParameters.peakWidthLowerBond = abs#1-#5;
% tmpFitParameters.peakPosLowerBond = #1;
% % spectrum range
% tmpFitParameters.backgroundRange = [5 5];
%
% % „ width“ mode
% tmpFitParameters.peakArea = [abs#1 abs#2 abs#3 abs#4];
@@ -50,7 +63,8 @@ function DataTable = nappFitGaussXps(DataTable,~,varargin)
% tmpFitParameters.peakAreaLowerBond = [#1 #2 #3 #4];
% tmpFitParameters.peakWidthLowerBond = #1-#4;
% tmpFitParameters.peakPosLowerBond = [#1 #2 #3 #4];
%
% spectrum range
% tmpFitParameters.backgroundRange = [5 5];
%
% % „ unconstrainedAndWidthPosition" mode (novel: rel distance fit parameter)
% tmpFitParameters.peakArea = [abs#1 abs#2 abs#3 abs#4 abs#5];
@@ -63,7 +77,8 @@ function DataTable = nappFitGaussXps(DataTable,~,varargin)
% tmpFitParameters.peakAreaLowerBond = [#1 #2 #3 #4];
% tmpFitParameters.peakWidthLowerBond = [#1,2,4,5 #3];
% tmpFitParameters.peakPosLowerBond = [#1 #2 #3 #4];
%
% spectrum range
% tmpFitParameters.backgroundRange = [5 5];
%
%
% tableSelection (optional): logical array to filter data table
@@ -209,25 +224,29 @@ elseif matches(p.Results.mode,'widthPosition')
+ tmpB(4)/(tmpB(2)*sqrt(pi/2))*exp(-2.*((x(:,1)- (tmpB(3)+p.Results.FitParameter.peakPos(1,2))).^2/tmpB(2).^2))... % Gauss peak 2
+ tmpB(5)/(tmpB(2)*sqrt(pi/2))*exp(-2.*((x(:,1)- (tmpB(3)+p.Results.FitParameter.peakPos(1,3))).^2/tmpB(2).^2))... % Gauss peak 3
+ tmpB(6)/(tmpB(2)*sqrt(pi/2))*exp(-2.*((x(:,1)- (tmpB(3)+p.Results.FitParameter.peakPos(1,4))).^2/tmpB(2).^2))... % Gauss peak 4
+tmpB(7).*x(:,1)+tmpB(8); % linear background
+ tmpB(7)/(tmpB(2)*sqrt(pi/2))*exp(-2.*((x(:,1)- (tmpB(3)+p.Results.FitParameter.peakPos(1,5))).^2/tmpB(2).^2))... % Gauss peak 4
+tmpB(8).*x(:,1)+tmpB(9); % linear background
% initial value
B0= [p.Results.FitParameter.peakArea(1,1) p.Results.FitParameter.peakWidth(1,1) p.Results.FitParameter.peakPos(1,1)... % Gauss peak 1
p.Results.FitParameter.peakArea(1,2) ... % Gauss peak 2
p.Results.FitParameter.peakArea(1,3) ... % Gauss peak 3
p.Results.FitParameter.peakArea(1,4)... % Gauss peak 4
p.Results.FitParameter.peakArea(1,4)...
p.Results.FitParameter.peakArea(1,5)...% Gauss peak 5
0 0]; % linear background
% lower bonds
lb = [p.Results.FitParameter.peakAreaLowerBond(1,1) p.Results.FitParameter.peakWidthLowerBond(1,1) p.Results.FitParameter.peakPosLowerBond(1,1) ... % Gauss peak 1
p.Results.FitParameter.peakAreaLowerBond(1,2) ... % Gauss peak 2
p.Results.FitParameter.peakAreaLowerBond(1,3) ... % Gauss peak 3
p.Results.FitParameter.peakAreaLowerBond(1,4) ... % Gauss peak 4
p.Results.FitParameter.peakAreaLowerBond(1,4) ...
p.Results.FitParameter.peakAreaLowerBond(1,5) ... % Gauss peak 5
0 0]; % linear background
% upper bond
ub = [p.Results.FitParameter.peakAreaUpperBond(1,1) p.Results.FitParameter.peakWidthUpperBond(1,1) p.Results.FitParameter.peakPosUpperBond(1,1) ... % Gauss peak 1
p.Results.FitParameter.peakAreaUpperBond(1,2) ... % Gauss peak 2
p.Results.FitParameter.peakAreaUpperBond(1,3) ... % Gauss peak 3
p.Results.FitParameter.peakAreaUpperBond(1,4) ... % Gauss peak 4
p.Results.FitParameter.peakAreaUpperBond(1,4) ...
p.Results.FitParameter.peakAreaUpperBond(1,5) ...% Gauss peak 5
Inf Inf]; % linear background
elseif matches(p.Results.mode,'width')
yfcn1 = @(tmpB,x)...
@@ -365,7 +384,8 @@ for ii = 1:size(p.Results.DataTable(p.Results.tableSelection,:),1)
tmpB(4) tmpB(2) tmpB(3)+p.Results.FitParameter.peakPos(1,2)...
tmpB(5) tmpB(2) tmpB(3)+p.Results.FitParameter.peakPos(1,3)...
tmpB(6) tmpB(2) tmpB(3)+p.Results.FitParameter.peakPos(1,4)...
tmpB(7) tmpB(8)...
tmpB(8) tmpB(9)...
tmpB(7) tmpB(2) tmpB(3)+p.Results.FitParameter.peakPos(1,5)...
NaN NaN NaN];
elseif matches(p.Results.mode,'width')
B = [tmpB(1) tmpB(2) tmpB(3)...
@@ -475,6 +495,9 @@ for ii = 1:size(p.Results.DataTable(p.Results.tableSelection,:),1)
resultsPeakFWHW(ii,:) = [B(2).*sqrt(log(4)), B(5).*sqrt(log(4)), B(8).*sqrt(log(4)), B(11).*sqrt(log(4)), B(17).*sqrt(log(4))] ;
resultsPeakArea(ii,:) = [B(1), B(4), B(7), B(10) B(15)] ;
% resultsBkgrd(ii,:) = [B(13), B(14)] ;
% close the fiting progess window
fh = findobj( 'Type', 'Figure', 'Name', 'GlobalSearch');
close(fh)
end