From 1aa130982fe094aeb8fa6b1ea7e918009f905d8d Mon Sep 17 00:00:00 2001 From: ThorBarRa Date: Thu, 9 Nov 2023 07:37:10 +0100 Subject: [PATCH] Added Gauss Peak Fitting to XPS shift determination. First implementation for positionWidth contrained fit. Not yet documented. Not yet fully tested. WIP --- beamtime/nappDetermineShiftXps.m | 149 +++++++++++++++++++------------ beamtime/nappFitGaussXps.m | 43 ++++++--- 2 files changed, 123 insertions(+), 69 deletions(-) diff --git a/beamtime/nappDetermineShiftXps.m b/beamtime/nappDetermineShiftXps.m index 1d3d1c5..b5415c0 100644 --- a/beamtime/nappDetermineShiftXps.m +++ b/beamtime/nappDetermineShiftXps.m @@ -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]; diff --git a/beamtime/nappFitGaussXps.m b/beamtime/nappFitGaussXps.m index 635598c..1f54122 100644 --- a/beamtime/nappFitGaussXps.m +++ b/beamtime/nappFitGaussXps.m @@ -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