861 lines
51 KiB
Matlab
861 lines
51 KiB
Matlab
function DataTable = nappFitGaussXpsDevelopment(DataTable,varargin)
|
|
% nappFitGaussXpsDevelopment Playground to add Shirley background :-).
|
|
% It allows different sets of constraints for the fitting process.
|
|
%
|
|
% INPUT:
|
|
% DataTable: Table containing the spectral data to be fitted.
|
|
% varargin: Additional optional parameters for customizing the fit.
|
|
%
|
|
% OUTPUT:
|
|
% DataTable: The input table, now including the results of the peak fitting.
|
|
%
|
|
% FileList =
|
|
% nappFitGaussXps(FileList,tmpFitParameters,'tableSelection',tmpDataSelection);
|
|
%%
|
|
%--------------------------------------------------------------------------
|
|
% Variable Definitions and Purposes
|
|
%--------------------------------------------------------------------------
|
|
%
|
|
% Input Condition and Default Variables
|
|
% ------------------------------------
|
|
% defaultTableVariable - Default spectral data column name.
|
|
% expectedTableVariable - Allowed spectral data column names.
|
|
% defaultMode - Default fitting constraint mode.
|
|
% expectedMode - Allowed fitting modes.
|
|
% defaultTableSelection - Default logical array (selects all rows).
|
|
% defaultShowFigures - Default to show fit figures ('on').
|
|
% expectedShowFigures - Allowed options for figure display.
|
|
% defaultEnergyScale - Default energy axis type ('kinetic').
|
|
% expectedEnergyScale - Allowed energy axis types.
|
|
% defaultBackgroundMode - Default background subtraction mode.
|
|
% expectedBackgroundMode - Allowed background subtraction modes.
|
|
%
|
|
%
|
|
% Input Parsing and Validation
|
|
% ---------------------------
|
|
% p - inputParser object for argument handling.
|
|
% FitParameter - Structure with fit parameters and constraints.
|
|
% tableSelection - Logical array selecting table rows for fitting.
|
|
% tableVariable - Name of the column containing the spectrum.
|
|
% mode - Fitting constraint mode.
|
|
% showFigures - Whether to display fit figures ('on'/'off').
|
|
% energyScale - Energy axis type ('kinetic'/'binding').
|
|
% backgroundMode - Background subtraction type.
|
|
%
|
|
%
|
|
% Data Table and Validation
|
|
% ------------------------
|
|
% DataTable - Main input table with all spectral data.
|
|
% p.Results.tableSelection - Logical array of selected rows.
|
|
% p.Results.tableVariable - Name of spectral data column to fit.
|
|
%
|
|
%
|
|
% Color Palette Variables
|
|
% ----------------------
|
|
% PsiFarbe - Structure holding RGB color values for plotting.
|
|
% PsiFarbe.hellgrau, etc. - Individual color entries for figure styling.
|
|
%
|
|
%
|
|
% Fitting Model and Optimization Variables
|
|
% ---------------------------------------
|
|
% modelDefinitionFunction - Function defining the sum-of-Gaussians model.
|
|
% initialParameterVector - Initial guess for fit parameters.
|
|
% lowerBoundsParameterVector - Lower bounds for fit parameters.
|
|
% upperBoundsParameterVector - Upper bounds for fit parameters.
|
|
% currentParameterVector - Current/optimized parameter vector.
|
|
% modelEvaluationFunction - Function evaluating the model during fitting.
|
|
% objectiveFunction - Function to minimize (sum of squared residuals).
|
|
% optimizationProblem - Structure describing the optimization problem.
|
|
% globalSearchOptimizer - GlobalSearch optimizer object.
|
|
% standardizedFitResults - Final results vector, standardized for all fit types.
|
|
%
|
|
%
|
|
% Spectral Data Variables
|
|
% ----------------------
|
|
% fitEnergyValues - Cropped energy values used for fitting.
|
|
% fitIntensityValues - Cropped intensity values used for fitting.
|
|
% tmp1, tmp2 - Indices for cropping fit region.
|
|
%
|
|
%
|
|
% Loop and Indexing Variables
|
|
% --------------------------
|
|
% ii - Loop index over selected spectra (rows in DataTable).
|
|
%% Define Input
|
|
% % „position“ mode
|
|
% tmpFitParameters.peakArea = [abs#1 abs#2 abs#3 abs#4];
|
|
% tmpFitParameters.peakWidth = [abs#1 abs#2 abs#3 abs#4];
|
|
% tmpFitParameters.peakPos = [abs#1 rel#2 rel#3 rel#4];
|
|
%
|
|
% tmpFitParameters.peakAreaUpperBound = [#1 #2 #3 #4];
|
|
% tmpFitParameters.peakWidthUpperBound = [#1 #2 #3 #4];
|
|
% tmpFitParameters.peakPosUpperBound = #1;
|
|
% tmpFitParameters.peakAreaLowerBound = [#1 #2 #3 #4];
|
|
% tmpFitParameters.peakWidthLowerBound = [#1 #2 #3 #4];
|
|
% tmpFitParameters.peakPosLowerBound = #1;
|
|
% spectrum range
|
|
% tmpFitParameters.backgroundRange = [5 5];
|
|
%
|
|
% % „unconstrained“ mode
|
|
% tmpFitParameters.peakArea = [abs#1 abs#2 abs#3 abs#4];
|
|
% tmpFitParameters.peakWidth = [abs#1 abs#2 abs#3 abs#4];
|
|
% tmpFitParameters.peakPos = [abs#1 abs#2 abs#3 abs#4];
|
|
%
|
|
% tmpFitParameters.peakAreaUpperBound = [#1 #2 #3 #4];
|
|
% tmpFitParameters.peakWidthUpperBound = [#1 #2 #3 #4];
|
|
% tmpFitParameters.peakPosUpperBound = [#1 #2 #3 #4];
|
|
% tmpFitParameters.peakAreaLowerBound = [#1 #2 #3 #4];
|
|
% tmpFitParameters.peakWidthLowerBound = [#1 #2 #3 #4];
|
|
% tmpFitParameters.peakPosLowerBound = [#1 #2 #3 #4];
|
|
% spectrum range
|
|
% tmpFitParameters.backgroundRange = [5 5];
|
|
%
|
|
% % „widthPosition“ mode
|
|
% 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.peakAreaUpperBound = [#1 #2 #3 #4 #5];
|
|
% tmpFitParameters.peakWidthUpperBound = abs#1-#5;
|
|
% tmpFitParameters.peakPosUpperBound = #1;
|
|
% tmpFitParameters.peakAreaLowerBound = [#1 #2 #3 #4 #5];
|
|
% tmpFitParameters.peakWidthLowerBound = abs#1-#5;
|
|
% tmpFitParameters.peakPosLowerBound = #1;
|
|
% % spectrum range
|
|
% tmpFitParameters.backgroundRange = [5 5];
|
|
%
|
|
% % „ width“ mode
|
|
% tmpFitParameters.peakArea = [abs#1 abs#2 abs#3 abs#4];
|
|
% tmpFitParameters.peakWidth = [abs#1-#4];
|
|
% tmpFitParameters.peakPos = [abs#1 abs#2 abs#3 abs#4];
|
|
%
|
|
% tmpFitParameters.peakAreaUpperBound = [#1 #2 #3 #4];
|
|
% tmpFitParameters.peakWidthUpperBound = #1-#4;
|
|
% tmpFitParameters.peakPosUpperBound = [#1 #2 #3 #4];
|
|
% tmpFitParameters.peakAreaLowerBound = [#1 #2 #3 #4];
|
|
% tmpFitParameters.peakWidthLowerBound = #1-#4;
|
|
% tmpFitParameters.peakPosLowerBound = [#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];
|
|
% tmpFitParameters.peakWidth = [abs#1,2,4,5 abs#3];
|
|
% tmpFitParameters.peakPos = [rel#1 abs#2 abs#3 abs#4 abs#5];
|
|
%
|
|
% tmpFitParameters.peakAreaUpperBound = [#1 #2 #3 #4 #5];
|
|
% tmpFitParameters.peakWidthUpperBound = [#1,2,4,5 #3];
|
|
% tmpFitParameters.peakPosUpperBound = [#1 #2 #3 #4 #5];
|
|
% tmpFitParameters.peakAreaLowerBound = [#1 #2 #3 #4];
|
|
% tmpFitParameters.peakWidthLowerBound = [#1,2,4,5 #3];
|
|
% tmpFitParameters.peakPosLowerBound = [#1 #2 #3 #4];
|
|
% spectrum range
|
|
% tmpFitParameters.backgroundRange = [5 5];
|
|
%
|
|
% % „gasAndTwoAcidBasePairs" mode (novel: rel distance fit parameter)
|
|
% tmpFitParameters.peakArea = [abs#1 abs#2 abs#3 abs#4 abs#5];
|
|
% tmpFitParameters.peakWidth = [abs#1,2,4,5 abs#3];
|
|
% tmpFitParameters.peakPos = [rel#1 abs#2 abs#3 rel#4 abs#5]; % rel#4
|
|
% is realative to abs#5, rel#1 is relative to abs#2
|
|
%
|
|
% tmpFitParameters.peakAreaUpperBound = [#1 #2 #3 #4 #5];
|
|
% tmpFitParameters.peakWidthUpperBound = [#1,2,4,5 #3];
|
|
% tmpFitParameters.peakPosUpperBound = [#1 #2 #3 #4 #5];
|
|
% tmpFitParameters.peakAreaLowerBound = [#1 #2 #3 #4 #5];
|
|
% tmpFitParameters.peakWidthLowerBound = [#1,2,4,5 #3];
|
|
% tmpFitParameters.peakPosLowerBound = [#1 #2 #3 #4 #5];
|
|
% spectrum range
|
|
% tmpFitParameters.backgroundRange = [5 5];
|
|
%
|
|
%
|
|
% tableSelection (optional): logical array to filter data table
|
|
% tableVariable (optional): select source: image(default),'imageRoh'
|
|
% mode:
|
|
% 'unconstrained': Each gauss peak is fitted independent.
|
|
% 'widthPosition' (Default): With of all peaks is identical. Position of peaks
|
|
% 2-3 is set relative to peak 1.
|
|
% 'position': Position of peaks 2-3 is set relative to peak 1.
|
|
% 'width': all peaks have same width
|
|
% 'unconstrainedAndWidthPosition': One unconstrained peak #3 (for gas phase
|
|
% feature) and 4 peaks constrained in width and/or position
|
|
% energyScale: kinetic (default) or binding
|
|
% showFigures: on(default) or off: display Gauss fits in figure.
|
|
%
|
|
|
|
%% Off we go - code starts here
|
|
|
|
% Define conditions and defaults of input
|
|
defaultTableVariable = 'spectrum_countsPerSecond';
|
|
expectedTableVariable = {'spectrum_counts','spectrum_countsPerSecond','spectrum_countsNorm'};
|
|
defaultMode = 'widthPosition';
|
|
expectedMode = {'unconstrained','widthPosition','position','width','unconstrainedAndWidthPosition','gasAndTwoAcidBasePairs','widthPositionAreaRatio','widthPositionAmineAreaRatio'};
|
|
defaultTableSelection = true(size(DataTable,1),1);
|
|
defaultShowFigures = 'on';
|
|
expectedShowFigures = {'on','off'};
|
|
defaultEnergyScale = 'kinetic';
|
|
expectedEnergyScale = {'kinetic','binding'};
|
|
defaultBackgroundMode = 'linear';
|
|
expectedBackgroundMode = {'linear','shirley'};
|
|
|
|
|
|
|
|
% Parse input
|
|
p = inputParser;
|
|
addRequired(p,'FitParameter',@isstruct);
|
|
addOptional(p,'tableSelection',defaultTableSelection,@islogical);
|
|
addParameter(p,'tableVariable',defaultTableVariable,...
|
|
@(energyValues) any(validatestring(energyValues,expectedTableVariable)));
|
|
addParameter(p,'mode',defaultMode,...
|
|
@(energyValues) any(validatestring(energyValues,expectedMode)));
|
|
addParameter(p,'showFigures',defaultShowFigures,...
|
|
@(energyValues) any(validatestring(energyValues,expectedShowFigures)));
|
|
addParameter(p,'energyScale',defaultEnergyScale,...
|
|
@(energyValues) any(validatestring(energyValues,expectedEnergyScale)));
|
|
addParameter(p,'backgroundMode',defaultBackgroundMode,...
|
|
@(energyValues) any(validatestring(energyValues,expectedBackgroundMode)));
|
|
parse(p,varargin{:});
|
|
|
|
if ~all(cell2mat(cellfun(@size,DataTable(p.Results.tableSelection,:).(p.Results.tableVariable),'UniformOutput',false)),'all')
|
|
disp('The table variable appears to be empty. I quit now')
|
|
error('The table variable appears to be empty.')
|
|
end
|
|
|
|
% Typo in FitParameters: Correct any 'Bond' -> 'Bound' typos
|
|
FitParameter = fixBondToBound(p.Results.FitParameter);
|
|
|
|
|
|
|
|
|
|
%--------------------------------------------------------------------------
|
|
% Define equations,
|
|
% mostly gauss and linear to model the spectra
|
|
%--------------------------------------------------------------------------
|
|
|
|
% Define equation to fit data
|
|
% I use the area Gauss equations here, which I took from the OriginPro
|
|
% documentation. https://www.originlab.com/doc/Origin-Help/Gauss-FitFunc
|
|
%
|
|
% Gauss = area/(width*sqrt(pi/2)) * exp(-2*(energyValues-position)^2/width^2)
|
|
%
|
|
% where FWHM = width * sqrt(2*ln(2))
|
|
|
|
if matches(p.Results.mode,'position')
|
|
modelDefinitionFunction = @(currentParameterVector,energyValues)...
|
|
currentParameterVector(1)/(currentParameterVector(2)*sqrt(pi/2))*exp(-2.*((energyValues(:,1) - currentParameterVector(3)).^2/currentParameterVector(2).^2)) ... % Gauss peak 1
|
|
+ currentParameterVector(4)/(currentParameterVector(5)*sqrt(pi/2))*exp(-2.*((energyValues(:,1)- (currentParameterVector(3)+FitParameter.peakPos(1,2))).^2/currentParameterVector(5).^2))... % Gauss peak 2
|
|
+ currentParameterVector(6)/(currentParameterVector(7)*sqrt(pi/2))*exp(-2.*((energyValues(:,1)- (currentParameterVector(3)+FitParameter.peakPos(1,3))).^2/currentParameterVector(7).^2))... % Gauss peak 3
|
|
+ currentParameterVector(8)/(currentParameterVector(9)*sqrt(pi/2))*exp(-2.*((energyValues(:,1)- (currentParameterVector(3)+FitParameter.peakPos(1,4))).^2/currentParameterVector(9).^2))... % Gauss peak 4
|
|
+currentParameterVector(10).*energyValues(:,1)+currentParameterVector(11); % linear background
|
|
|
|
% initial value
|
|
initialParameterVector= [FitParameter.peakArea(1,1) FitParameter.peakWidth(1,1) FitParameter.peakPos(1,1)... % Gauss peak 1
|
|
FitParameter.peakArea(1,2) FitParameter.peakWidth(1,2)... % Gauss peak 2
|
|
FitParameter.peakArea(1,3) FitParameter.peakWidth(1,3)... % Gauss peak 3
|
|
FitParameter.peakArea(1,4) FitParameter.peakWidth(1,4)... % Gauss peak 4
|
|
1 0]; % linear background
|
|
% lower bounds
|
|
lowerBoundsParameterVector = [FitParameter.peakAreaLowerBound(1,1) FitParameter.peakWidthLowerBound(1,1) FitParameter.peakPosLowerBound(1,1) ... % Gauss peak 1
|
|
FitParameter.peakAreaLowerBound(1,2) FitParameter.peakWidthLowerBound(1,2) ... % Gauss peak 2
|
|
FitParameter.peakAreaLowerBound(1,3) FitParameter.peakWidthLowerBound(1,3) ... % Gauss peak 3
|
|
FitParameter.peakAreaLowerBound(1,4) FitParameter.peakWidthLowerBound(1,4) ... % Gauss peak 4
|
|
-Inf -Inf]; % linear background
|
|
% upper bound
|
|
upperBoundsParameterVector = [FitParameter.peakAreaUpperBound(1,1) FitParameter.peakWidthUpperBound(1,1) FitParameter.peakPosUpperBound(1,1) ... % Gauss peak 1
|
|
FitParameter.peakAreaUpperBound(1,2) FitParameter.peakWidthUpperBound(1,2) ... % Gauss peak 2
|
|
FitParameter.peakAreaUpperBound(1,3) FitParameter.peakWidthUpperBound(1,3) ... % Gauss peak 3
|
|
FitParameter.peakAreaUpperBound(1,4) FitParameter.peakWidthUpperBound(1,4) ... % Gauss peak 4
|
|
Inf Inf]; % linear background
|
|
elseif matches(p.Results.mode,'unconstrained')
|
|
% --- New July 2025: Use functions to compute equations used for fitting with variable number of peaks
|
|
% Works only with 4 peaks at the moment.
|
|
if matches(p.Results.backgroundMode,'linear') && ~exist("FitParameter.peakEta","var") % to be compatible with older function calls, we test if the new parameter peakEta is defined. If not it uses gauss only.
|
|
modelDefinitionFunction = @(currentParameterVector,energyValues) gaussianModel(currentParameterVector,energyValues, size(FitParameter.peakArea,2));
|
|
[initialParameterVector,lowerBoundsParameterVector,upperBoundsParameterVector] = gaussianModelBounds(FitParameter);
|
|
elseif matches(p.Results.backgroundMode,'linear') && exist("FitParameter.peakEta","var")
|
|
modelDefinitionFunction = @(currentParameterVector,energyValues) pseudoVoigtModel(currentParameterVector,energyValues, size(FitParameter.peakArea,2));
|
|
[initialParameterVector,lowerBoundsParameterVector,upperBoundsParameterVector] = pseudoVoigtModelBounds(FitParameter);
|
|
elseif matches(p.Results.backgroundMode,'shirley')
|
|
modelDefinitionFunction = @(currentParameterVector,energyValues) shirleyPseudoVoigtModel(currentParameterVector,energyValues, size(FitParameter.peakArea,2));
|
|
[initialParameterVector,lowerBoundsParameterVector,upperBoundsParameterVector] = shirleyPseudoVoigtModelBounds(FitParameter);
|
|
end
|
|
|
|
elseif matches(p.Results.mode,'widthPosition')
|
|
modelDefinitionFunction = @(currentParameterVector,energyValues)...
|
|
currentParameterVector(1)/(currentParameterVector(2)*sqrt(pi/2))*exp(-2.*((energyValues(:,1) - currentParameterVector(3)).^2/currentParameterVector(2).^2)) ... % Gauss peak 1
|
|
+ currentParameterVector(4)/(currentParameterVector(2)*sqrt(pi/2))*exp(-2.*((energyValues(:,1)- (currentParameterVector(3)+FitParameter.peakPos(1,2))).^2/currentParameterVector(2).^2))... % Gauss peak 2
|
|
+ currentParameterVector(5)/(currentParameterVector(2)*sqrt(pi/2))*exp(-2.*((energyValues(:,1)- (currentParameterVector(3)+FitParameter.peakPos(1,3))).^2/currentParameterVector(2).^2))... % Gauss peak 3
|
|
+ currentParameterVector(6)/(currentParameterVector(2)*sqrt(pi/2))*exp(-2.*((energyValues(:,1)- (currentParameterVector(3)+FitParameter.peakPos(1,4))).^2/currentParameterVector(2).^2))... % Gauss peak 4
|
|
+ currentParameterVector(7)/(currentParameterVector(2)*sqrt(pi/2))*exp(-2.*((energyValues(:,1)- (currentParameterVector(3)+FitParameter.peakPos(1,5))).^2/currentParameterVector(2).^2))... % Gauss peak 4
|
|
+currentParameterVector(8).*energyValues(:,1)+currentParameterVector(9); % linear background
|
|
|
|
% initial value
|
|
initialParameterVector= [FitParameter.peakArea(1,1) FitParameter.peakWidth(1,1) FitParameter.peakPos(1,1)... % Gauss peak 1
|
|
FitParameter.peakArea(1,2) ... % Gauss peak 2
|
|
FitParameter.peakArea(1,3) ... % Gauss peak 3
|
|
FitParameter.peakArea(1,4)...
|
|
FitParameter.peakArea(1,5)...% Gauss peak 5
|
|
1 0]; % linear background
|
|
% lower bounds
|
|
lowerBoundsParameterVector = [FitParameter.peakAreaLowerBound(1,1) FitParameter.peakWidthLowerBound(1,1) FitParameter.peakPosLowerBound(1,1) ... % Gauss peak 1
|
|
FitParameter.peakAreaLowerBound(1,2) ... % Gauss peak 2
|
|
FitParameter.peakAreaLowerBound(1,3) ... % Gauss peak 3
|
|
FitParameter.peakAreaLowerBound(1,4) ...
|
|
FitParameter.peakAreaLowerBound(1,5) ... % Gauss peak 5
|
|
-Inf -Inf]; % linear background
|
|
% upper bound
|
|
upperBoundsParameterVector = [FitParameter.peakAreaUpperBound(1,1) FitParameter.peakWidthUpperBound(1,1) FitParameter.peakPosUpperBound(1,1) ... % Gauss peak 1
|
|
FitParameter.peakAreaUpperBound(1,2) ... % Gauss peak 2
|
|
FitParameter.peakAreaUpperBound(1,3) ... % Gauss peak 3
|
|
FitParameter.peakAreaUpperBound(1,4) ...
|
|
FitParameter.peakAreaUpperBound(1,5) ...% Gauss peak 5
|
|
Inf Inf]; % linear background
|
|
elseif matches(p.Results.mode,'width')
|
|
modelDefinitionFunction = @(currentParameterVector,energyValues)...
|
|
currentParameterVector(1)/(currentParameterVector(2)*sqrt(pi/2))*exp(-2.*((energyValues(:,1) - currentParameterVector(3)).^2/currentParameterVector(2).^2)) ... % Gauss peak 1
|
|
+ currentParameterVector(4)/(currentParameterVector(2)*sqrt(pi/2))*exp(-2.*((energyValues(:,1)- currentParameterVector(5)).^2/currentParameterVector(2).^2))... % Gauss peak 2
|
|
+ currentParameterVector(6)/(currentParameterVector(2)*sqrt(pi/2))*exp(-2.*((energyValues(:,1)- currentParameterVector(7)).^2/currentParameterVector(2).^2))... % Gauss peak 3
|
|
+ currentParameterVector(8)/(currentParameterVector(2)*sqrt(pi/2))*exp(-2.*((energyValues(:,1)- currentParameterVector(9)).^2/currentParameterVector(2).^2))... % Gauss peak 4
|
|
+currentParameterVector(10).*energyValues(:,1)+currentParameterVector(11); % linear background
|
|
|
|
% initial value
|
|
initialParameterVector= [FitParameter.peakArea(1,1) FitParameter.peakWidth(1,1) FitParameter.peakPos(1,1)... % Gauss peak 1
|
|
FitParameter.peakArea(1,2) FitParameter.peakPos(1,2)... % Gauss peak 2
|
|
FitParameter.peakArea(1,3) FitParameter.peakPos(1,3)... % Gauss peak 3
|
|
FitParameter.peakArea(1,4) FitParameter.peakPos(1,4)... % Gauss peak 4
|
|
1 0]; % linear background
|
|
% lower bounds
|
|
lowerBoundsParameterVector = [FitParameter.peakAreaLowerBound(1,1) FitParameter.peakWidthLowerBound(1,1) FitParameter.peakPosLowerBound(1,1) ... % Gauss peak 1
|
|
FitParameter.peakAreaLowerBound(1,2) FitParameter.peakPos(1,2) ... % Gauss peak 2
|
|
FitParameter.peakAreaLowerBound(1,3) FitParameter.peakPos(1,3) ... % Gauss peak 3
|
|
FitParameter.peakAreaLowerBound(1,4) FitParameter.peakPos(1,4) ... % Gauss peak 4
|
|
-Inf -Inf]; % linear background
|
|
% upper bound
|
|
upperBoundsParameterVector = [FitParameter.peakAreaUpperBound(1,1) FitParameter.peakWidthUpperBound(1,1) FitParameter.peakPosUpperBound(1,1) ... % Gauss peak 1
|
|
FitParameter.peakAreaUpperBound(1,2) FitParameter.peakPos(1,2) ... % Gauss peak 2
|
|
FitParameter.peakAreaUpperBound(1,3) FitParameter.peakPos(1,3) ... % Gauss peak 3
|
|
FitParameter.peakAreaUpperBound(1,4) FitParameter.peakPos(1,4) ... % Gauss peak 4
|
|
Inf Inf]; % linear background
|
|
|
|
elseif matches(p.Results.mode,'unconstrainedAndWidthPosition') % Novum, also the distance between peaks is a fit parameter. [Gauss = area/(width*sqrt(pi/2)) * exp(-2*(energyValues-position)^2/width^2)]
|
|
modelDefinitionFunction = @(currentParameterVector,energyValues)...
|
|
currentParameterVector(1)/(currentParameterVector(2)*sqrt(pi/2))*exp(-2.*((energyValues(:,1) - (currentParameterVector(3)+currentParameterVector(14))).^2/currentParameterVector(2).^2)) ... % Gauss peak 1: width and position
|
|
+ currentParameterVector(4)/(currentParameterVector(2)*sqrt(pi/2))*exp(-2.*((energyValues(:,1) - currentParameterVector(3)).^2/currentParameterVector(2).^2))... % Gauss peak 2: width and position
|
|
+ currentParameterVector(5)/(currentParameterVector(6)*sqrt(pi/2))*exp(-2.*((energyValues(:,1) - currentParameterVector(7)).^2/currentParameterVector(6).^2))... % Gauss peak 3: unconstrained
|
|
+ currentParameterVector(8)/(currentParameterVector(2)*sqrt(pi/2))*exp(-2.*((energyValues(:,1) - currentParameterVector(9)).^2/currentParameterVector(2).^2))... % Gauss peak 4: width
|
|
+ currentParameterVector(10)/(currentParameterVector(2)*sqrt(pi/2))*exp(-2.*((energyValues(:,1)- currentParameterVector(11)).^2/currentParameterVector(2).^2))... % Gauss peak 5: width
|
|
+ currentParameterVector(12).*energyValues(:,1)+currentParameterVector(13); % linear background
|
|
|
|
% initial value
|
|
initialParameterVector= [FitParameter.peakArea(1,1) FitParameter.peakWidth(1,1) FitParameter.peakPos(1,2)... % Gauss peak 1
|
|
FitParameter.peakArea(1,2) ... % Gauss peak 2
|
|
FitParameter.peakArea(1,3) FitParameter.peakWidth(1,2) FitParameter.peakPos(1,3)... % Gauss peak 3
|
|
FitParameter.peakArea(1,4) FitParameter.peakPos(1,4)... % Gauss peak 4
|
|
FitParameter.peakArea(1,5) FitParameter.peakPos(1,5)... % Gauss peak 5
|
|
1 0 ... % linear background
|
|
FitParameter.peakPos(1,1)];
|
|
% lower bounds
|
|
lowerBoundsParameterVector = [FitParameter.peakAreaLowerBound(1,1) FitParameter.peakWidthLowerBound(1,1) FitParameter.peakPosLowerBound(1,2) ... % Gauss peak 1
|
|
FitParameter.peakAreaLowerBound(1,2) ... % Gauss peak 2
|
|
FitParameter.peakAreaLowerBound(1,3) FitParameter.peakWidthLowerBound(1,2) FitParameter.peakPosLowerBound(1,3) ... % Gauss peak 3
|
|
FitParameter.peakAreaLowerBound(1,4) FitParameter.peakPosLowerBound(1,4) ...
|
|
FitParameter.peakAreaLowerBound(1,5) FitParameter.peakPosLowerBound(1,5) ...% Gauss peak 4
|
|
-Inf -Inf ... % linear background
|
|
FitParameter.peakPosLowerBound(1,1)];
|
|
|
|
% upper bound
|
|
upperBoundsParameterVector = [FitParameter.peakAreaUpperBound(1,1) FitParameter.peakWidthUpperBound(1,1) FitParameter.peakPosUpperBound(1,2) ... % Gauss peak 1
|
|
FitParameter.peakAreaUpperBound(1,2) ... % Gauss peak 2
|
|
FitParameter.peakAreaUpperBound(1,3) FitParameter.peakWidthUpperBound(1,2) FitParameter.peakPosUpperBound(1,3) ... % Gauss peak 3
|
|
FitParameter.peakAreaUpperBound(1,4) FitParameter.peakPosUpperBound(1,4) ...
|
|
FitParameter.peakAreaUpperBound(1,5) FitParameter.peakPosUpperBound(1,5) ... % Gauss peak 4
|
|
Inf Inf ... % linear background
|
|
FitParameter.peakPosUpperBound(1,1)]; % shift upper bound
|
|
|
|
elseif matches(p.Results.mode,'gasAndTwoAcidBasePairs') % Novum, also the distance between peaks is a fit parameter.
|
|
% [Gauss = area/(width*sqrt(pi/2)) * exp(-2*(energyValues-position)^2/width^2)]
|
|
modelDefinitionFunction = @(currentParameterVector,energyValues)...
|
|
currentParameterVector(1)/(currentParameterVector(2)*sqrt(pi/2))*exp(-2.*((energyValues(:,1) - (currentParameterVector(3)+currentParameterVector(14))).^2/currentParameterVector(2).^2)) ... % Gauss peak 1: width and position
|
|
+ currentParameterVector(4)/(currentParameterVector(2)*sqrt(pi/2))*exp(-2.*((energyValues(:,1) - currentParameterVector(3)).^2/currentParameterVector(2).^2))... % Gauss peak 2: width and position
|
|
+ currentParameterVector(5)/(currentParameterVector(6)*sqrt(pi/2))*exp(-2.*((energyValues(:,1) - currentParameterVector(7)).^2/currentParameterVector(6).^2))... % Gauss peak 3: unconstrained
|
|
+ currentParameterVector(8)/(currentParameterVector(2)*sqrt(pi/2))*exp(-2.*((energyValues(:,1) - (currentParameterVector(11)+currentParameterVector(9))).^2/currentParameterVector(2).^2))... % Gauss peak 4: width and position
|
|
+ currentParameterVector(10)/(currentParameterVector(2)*sqrt(pi/2))*exp(-2.*((energyValues(:,1)- currentParameterVector(11)).^2/currentParameterVector(2).^2))... % Gauss peak 5: width and position
|
|
+ currentParameterVector(12).*energyValues(:,1)+currentParameterVector(13); % linear background
|
|
|
|
% initial value
|
|
initialParameterVector= [FitParameter.peakArea(1,1) FitParameter.peakWidth(1,1) FitParameter.peakPos(1,2)... % Gauss peak 1
|
|
FitParameter.peakArea(1,2) ... % Gauss peak 2
|
|
FitParameter.peakArea(1,3) FitParameter.peakWidth(1,2) FitParameter.peakPos(1,3)... % Gauss peak 3
|
|
FitParameter.peakArea(1,4) FitParameter.peakPos(1,4)... % Gauss peak 4
|
|
FitParameter.peakArea(1,5) FitParameter.peakPos(1,5)... % Gauss peak 5
|
|
1 0 ... % linear background
|
|
FitParameter.peakPos(1,1)];
|
|
% lower bounds
|
|
lowerBoundsParameterVector = [FitParameter.peakAreaLowerBound(1,1) FitParameter.peakWidthLowerBound(1,1) FitParameter.peakPosLowerBound(1,2) ... % Gauss peak 1
|
|
FitParameter.peakAreaLowerBound(1,2) ... % Gauss peak 2
|
|
FitParameter.peakAreaLowerBound(1,3) FitParameter.peakWidthLowerBound(1,2) FitParameter.peakPosLowerBound(1,3) ... % Gauss peak 3
|
|
FitParameter.peakAreaLowerBound(1,4) FitParameter.peakPosLowerBound(1,4) ...
|
|
FitParameter.peakAreaLowerBound(1,5) FitParameter.peakPosLowerBound(1,5) ...% Gauss peak 4
|
|
-Inf -Inf ... % linear background
|
|
FitParameter.peakPosLowerBound(1,1)];
|
|
|
|
% upper bound
|
|
upperBoundsParameterVector = [FitParameter.peakAreaUpperBound(1,1) FitParameter.peakWidthUpperBound(1,1) FitParameter.peakPosUpperBound(1,2) ... % Gauss peak 1
|
|
FitParameter.peakAreaUpperBound(1,2) ... % Gauss peak 2
|
|
FitParameter.peakAreaUpperBound(1,3) FitParameter.peakWidthUpperBound(1,2) FitParameter.peakPosUpperBound(1,3) ... % Gauss peak 3
|
|
FitParameter.peakAreaUpperBound(1,4) FitParameter.peakPosUpperBound(1,4) ...
|
|
FitParameter.peakAreaUpperBound(1,5) FitParameter.peakPosUpperBound(1,5) ... % Gauss peak 4
|
|
Inf Inf ... % linear background
|
|
FitParameter.peakPosUpperBound(1,1)]; % shift upper bound
|
|
|
|
elseif matches(p.Results.mode,'widthPositionAreaRatio')
|
|
modelDefinitionFunction = @(currentParameterVector,energyValues)...
|
|
currentParameterVector(1)/(currentParameterVector(2)*sqrt(pi/2))*exp(-2.*((energyValues(:,1) - currentParameterVector(3)).^2/currentParameterVector(2).^2)) ... % Gauss peak 1
|
|
+ currentParameterVector(1)*FitParameter.peakArea(1,2)/(currentParameterVector(2)*sqrt(pi/2))*exp(-2.*((energyValues(:,1)- (currentParameterVector(3)+FitParameter.peakPos(1,2))).^2/currentParameterVector(2).^2))... % Gauss peak 2
|
|
+ currentParameterVector(1)*FitParameter.peakArea(1,3)/(currentParameterVector(2)*sqrt(pi/2))*exp(-2.*((energyValues(:,1)- (currentParameterVector(3)+FitParameter.peakPos(1,3))).^2/currentParameterVector(2).^2))... % Gauss peak 3
|
|
+ currentParameterVector(1)*FitParameter.peakArea(1,4)/(currentParameterVector(2)*sqrt(pi/2))*exp(-2.*((energyValues(:,1)- (currentParameterVector(3)+FitParameter.peakPos(1,4))).^2/currentParameterVector(2).^2))... % Gauss peak 4
|
|
+ currentParameterVector(1)*FitParameter.peakArea(1,5)/(currentParameterVector(2)*sqrt(pi/2))*exp(-2.*((energyValues(:,1)- (currentParameterVector(3)+FitParameter.peakPos(1,5))).^2/currentParameterVector(2).^2))... % Gauss peak 4
|
|
+currentParameterVector(4).*energyValues(:,1)+currentParameterVector(5); % linear background
|
|
|
|
% initial value
|
|
initialParameterVector= [FitParameter.peakArea(1,1) FitParameter.peakWidth(1,1) FitParameter.peakPos(1,1)... % Gauss peak 1
|
|
1 0]; % linear background
|
|
% lower bounds
|
|
lowerBoundsParameterVector = [FitParameter.peakAreaLowerBound(1,1) FitParameter.peakWidthLowerBound(1,1) FitParameter.peakPosLowerBound(1,1) ... % Gauss peak 1
|
|
-Inf -Inf]; % linear background
|
|
% upper bound
|
|
upperBoundsParameterVector = [FitParameter.peakAreaUpperBound(1,1) FitParameter.peakWidthUpperBound(1,1) FitParameter.peakPosUpperBound(1,1) ... % Gauss peak 1
|
|
Inf Inf]; % linear background
|
|
|
|
elseif matches(p.Results.mode,'widthPositionAmineAreaRatio')
|
|
modelDefinitionFunction = @(currentParameterVector,energyValues)...
|
|
currentParameterVector(1)/(currentParameterVector(2)*sqrt(pi/2))*exp(-2.*((energyValues(:,1) - currentParameterVector(3)).^2/currentParameterVector(2).^2)) ... % Gauss peak 1
|
|
+ currentParameterVector(4)/(currentParameterVector(2)*sqrt(pi/2))*exp(-2.*((energyValues(:,1)- (currentParameterVector(3)+FitParameter.peakPos(1,2))).^2/currentParameterVector(2).^2))... % Gauss peak 2
|
|
+ currentParameterVector(4)*FitParameter.peakArea(1,3)/(currentParameterVector(2)*sqrt(pi/2))*exp(-2.*((energyValues(:,1)- (currentParameterVector(3)+FitParameter.peakPos(1,3))).^2/currentParameterVector(2).^2))... % Gauss peak 3
|
|
+ currentParameterVector(5)/(currentParameterVector(2)*sqrt(pi/2))*exp(-2.*((energyValues(:,1)- (currentParameterVector(3)+FitParameter.peakPos(1,4))).^2/currentParameterVector(2).^2))... % Gauss peak 4
|
|
+ currentParameterVector(6)/(currentParameterVector(2)*sqrt(pi/2))*exp(-2.*((energyValues(:,1)- (currentParameterVector(3)+FitParameter.peakPos(1,5))).^2/currentParameterVector(2).^2))... % Gauss peak 4
|
|
+currentParameterVector(7).*energyValues(:,1)+currentParameterVector(8); % linear background
|
|
|
|
% initial value
|
|
initialParameterVector= [FitParameter.peakArea(1,1) FitParameter.peakWidth(1,1) FitParameter.peakPos(1,1)... % Gauss peak 1
|
|
FitParameter.peakArea(1,2) ... % Gauss peak 3
|
|
FitParameter.peakArea(1,4)...
|
|
FitParameter.peakArea(1,5)...% Gauss peak 5
|
|
1 0]; % linear background
|
|
% lower bounds
|
|
lowerBoundsParameterVector = [FitParameter.peakAreaLowerBound(1,1) FitParameter.peakWidthLowerBound(1,1) FitParameter.peakPosLowerBound(1,1) ... % Gauss peak 1
|
|
FitParameter.peakAreaLowerBound(1,2) ... % Gauss peak 3
|
|
FitParameter.peakAreaLowerBound(1,4) ...
|
|
FitParameter.peakAreaLowerBound(1,5) ... % Gauss peak 5
|
|
-Inf -Inf]; % linear background
|
|
% upper bound
|
|
upperBoundsParameterVector = [FitParameter.peakAreaUpperBound(1,1) FitParameter.peakWidthUpperBound(1,1) FitParameter.peakPosUpperBound(1,1) ... % Gauss peak 1
|
|
FitParameter.peakAreaUpperBound(1,2) ... % Gauss peak 3
|
|
FitParameter.peakAreaUpperBound(1,4) ...
|
|
FitParameter.peakAreaUpperBound(1,5) ...% Gauss peak 5
|
|
Inf Inf]; % linear background
|
|
end
|
|
|
|
|
|
|
|
for ii = 1:size(DataTable(p.Results.tableSelection,:),1)
|
|
|
|
if contains(p.Results.energyScale,'kinetic')
|
|
%--------------------------------------------------------------------------
|
|
% Prepare data (kinetic energy scale)
|
|
% optionally cut to a section of relevance
|
|
%--------------------------------------------------------------------------
|
|
|
|
[fitEnergyValues, fitIntensityValues] = prepareCurveData(DataTable(p.Results.tableSelection,:).scientaEkin_eV{ii},...
|
|
DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii});
|
|
% make data slimmer so that background fit improved
|
|
if FitParameter.peakPos(1,1) >= 10
|
|
tmp1 = find(fitEnergyValues>=FitParameter.peakPos(1,1)-FitParameter.backgroundRange(1,1),1);
|
|
tmp2 = find(fitEnergyValues>=FitParameter.peakPos(1,1)+FitParameter.backgroundRange(1,2),1);
|
|
else
|
|
tmp1 = find(fitEnergyValues>=FitParameter.peakPos(1,2)-FitParameter.backgroundRange(1,1),1);
|
|
tmp2 = find(fitEnergyValues>=FitParameter.peakPos(1,2)+FitParameter.backgroundRange(1,2),1);
|
|
end
|
|
if isempty(tmp2)
|
|
fitIntensityValues = fitIntensityValues(tmp1:end);
|
|
fitEnergyValues = fitEnergyValues(tmp1:end);
|
|
else
|
|
fitIntensityValues = fitIntensityValues(tmp1:tmp2);
|
|
fitEnergyValues = fitEnergyValues(tmp1:tmp2);
|
|
end
|
|
|
|
%--------------------------------------------------------------------------
|
|
% Core part,
|
|
% the fitting is done here (kinetic energy scale)
|
|
%--------------------------------------------------------------------------
|
|
|
|
% Fit equations: Gauss eq. 1 with data
|
|
modelEvaluationFunction = @(currentParameterVector,energyValues) modelDefinitionFunction(currentParameterVector,energyValues);
|
|
% We find an optimal fit by minimising the residuals:
|
|
objectiveFunction = @(currentParameterVector) sum((fitIntensityValues - modelEvaluationFunction(currentParameterVector,fitEnergyValues)).^2);
|
|
optimizationProblem = createOptimProblem('fmincon', 'x0',initialParameterVector, 'lb',lowerBoundsParameterVector,'ub',upperBoundsParameterVector, 'objective',objectiveFunction);
|
|
globalSearchOptimizer = GlobalSearch('PlotFcns',@gsplotbestf);
|
|
currentParameterVector = run(globalSearchOptimizer,optimizationProblem);
|
|
elseif contains(p.Results.energyScale,'binding')
|
|
%--------------------------------------------------------------------------
|
|
% Prepare data (binding energy scale)
|
|
% optionally cut to a section of relevance
|
|
%--------------------------------------------------------------------------
|
|
|
|
[fitEnergyValues, fitIntensityValues] = prepareCurveData(DataTable(p.Results.tableSelection,:).bindingEnergy_eV{ii},...
|
|
DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii});
|
|
% make data slimmer so that background fit improves
|
|
if FitParameter.peakPos(1,1) >= 10
|
|
tmp2 = find(fitEnergyValues<=FitParameter.peakPos(1,1)-FitParameter.backgroundRange(1,1),1);
|
|
tmp1 = find(fitEnergyValues<=FitParameter.peakPos(1,1)+FitParameter.backgroundRange(1,2),1);
|
|
else
|
|
tmp2 = find(fitEnergyValues<=FitParameter.peakPos(1,2)-FitParameter.backgroundRange(1,1),1);
|
|
tmp1 = find(fitEnergyValues<=FitParameter.peakPos(1,2)+FitParameter.backgroundRange(1,2),1);
|
|
end
|
|
if isempty(tmp2)
|
|
fitIntensityValues = fitIntensityValues(tmp1:end);
|
|
fitEnergyValues = fitEnergyValues(tmp1:end);
|
|
else
|
|
fitIntensityValues = fitIntensityValues(tmp1:tmp2);
|
|
fitEnergyValues = fitEnergyValues(tmp1:tmp2);
|
|
end
|
|
%--------------------------------------------------------------------------
|
|
% Core part,
|
|
% the fitting is done here (binding energy scale)
|
|
%--------------------------------------------------------------------------
|
|
|
|
% Fit equations: Gauss eq. 1 with data
|
|
modelEvaluationFunction = @(currentParameterVector,energyValues) modelDefinitionFunction(currentParameterVector,energyValues);
|
|
% We find an optimal fit by minimising the residuals:
|
|
objectiveFunction = @(currentParameterVector) sum((fitIntensityValues - modelEvaluationFunction(currentParameterVector,fitEnergyValues)).^2);
|
|
optimizationProblem = createOptimProblem('fmincon', 'x0',initialParameterVector, 'lowerBoundsParameterVector',lowerBoundsParameterVector,'upperBoundsParameterVector',upperBoundsParameterVector, 'objective',objectiveFunction);
|
|
globalSearchOptimizer = GlobalSearch('PlotFcns',@gsplotbestf);
|
|
currentParameterVector = run(globalSearchOptimizer,optimizationProblem);
|
|
end
|
|
|
|
|
|
|
|
%--------------------------------------------------------------------------
|
|
% Export results
|
|
% prepare for data table
|
|
%--------------------------------------------------------------------------
|
|
tmpPeakPos = [];
|
|
tmpPeakFWHW = [];
|
|
tmpPeakArea = [];
|
|
tmpBkgrd = [];
|
|
tmpPeakEta = [];
|
|
tmpPeakShirleyK = [];
|
|
|
|
if matches(p.Results.backgroundMode,'linear') && ~exist("FitParameter.peakEta","var")
|
|
|
|
nPeaks = size(FitParameter.peakArea,2);
|
|
% separate pos, width, and area into sep. variables for any number of peaks.
|
|
for ij = 1:nPeaks
|
|
idx = (ij - 1)*3;
|
|
tmpPeakArea = [tmpPeakArea, currentParameterVector(idx+1)];
|
|
tmpPeakFWHW = [tmpPeakFWHW, currentParameterVector(idx+2).*sqrt( 2 * log(2))];
|
|
tmpPeakPos = [tmpPeakPos, currentParameterVector(idx+3)] ;
|
|
end
|
|
% add a line for each spectrum.
|
|
resultsPeakPos(ii,:) = tmpPeakPos;
|
|
resultsPeakFWHW(ii,:)= tmpPeakFWHW;
|
|
resultsPeakArea(ii,:)= tmpPeakArea;
|
|
resultsBkgrd = [currentParameterVector(end-1), currentParameterVector(end)] ;
|
|
|
|
elseif matches(p.Results.backgroundMode,'linear') && exist("FitParameter.peakEta","var")
|
|
|
|
elseif matches(p.Results.backgroundMode,'shirley')
|
|
|
|
else
|
|
resultsPeakPos(ii,:) = [standardizedFitResults(3), standardizedFitResults(6), standardizedFitResults(9), standardizedFitResults(12), standardizedFitResults(16)] ; %#ok<*AGROW>
|
|
resultsPeakFWHW(ii,:) = [standardizedFitResults(2).*sqrt( 2 * log(2)), standardizedFitResults(5).*sqrt( 2 * log(2)), standardizedFitResults(8).*sqrt( 2 * log(2)), standardizedFitResults(11).*sqrt( 2 * log(2)), standardizedFitResults(17).*sqrt( 2 * log(2))] ;
|
|
resultsPeakArea(ii,:) = [standardizedFitResults(1), standardizedFitResults(4), standardizedFitResults(7), standardizedFitResults(10) standardizedFitResults(15)] ;
|
|
resultsBkgrd(ii,:) = [standardizedFitResults(13), standardizedFitResults(14)] ;
|
|
end
|
|
% close the fiting progess window
|
|
fh = findobj( 'Type', 'Figure', 'Name', 'GlobalSearch');
|
|
close(fh)
|
|
end
|
|
|
|
|
|
% write results to data table
|
|
try
|
|
existColumn = strcmp('peakPosition',DataTable.Properties.VariableNames);
|
|
if ~any(existColumn)
|
|
DataTable.peakPosition = NaN(size(DataTable,1),size(resultsPeakPos,2));
|
|
end
|
|
existColumn = strcmp('peakArea',DataTable.Properties.VariableNames);
|
|
if ~any(existColumn)
|
|
DataTable.peakArea = NaN(size(DataTable,1),size(resultsPeakArea,2));
|
|
end
|
|
existColumn = strcmp('peakFWHM',DataTable.Properties.VariableNames);
|
|
if ~any(existColumn)
|
|
DataTable.peakFWHM = NaN(size(DataTable,1),size(resultsPeakFWHW,2));
|
|
end
|
|
existColumn = strcmp('peakBackground',DataTable.Properties.VariableNames);
|
|
if ~any(existColumn)
|
|
DataTable.peakBackground = NaN(size(DataTable,1),size(resultsBkgrd,2));
|
|
end
|
|
|
|
n=max(numel(resultsPeakPos),numel(DataTable(p.Results.tableSelection,:).peakPosition));
|
|
|
|
resultsPeakPos(end+1:n)=nan;
|
|
resultsPeakFWHW(end+1:n)=nan;
|
|
resultsPeakArea(end+1:n)=nan;
|
|
|
|
%if contains(p.Results.energyScale,'binding')
|
|
DataTable(p.Results.tableSelection,:).peakPosition = resultsPeakPos(:,:);
|
|
%end
|
|
DataTable(p.Results.tableSelection,:).peakFWHM = resultsPeakFWHW(:,:);
|
|
DataTable(p.Results.tableSelection,:).peakArea = resultsPeakArea(:,:);
|
|
DataTable(p.Results.tableSelection,:).peakBackground = resultsBkgrd(:,:);
|
|
catch
|
|
|
|
disp('Could not write results to table')
|
|
|
|
end
|
|
end
|
|
|
|
%% ------------------ Functions ---------------
|
|
|
|
function correctedStruct = fixBondToBound(inputStruct)
|
|
% Recursively corrects all field names in a struct: 'Bond' -> 'Bound'
|
|
|
|
correctedStruct = inputStruct;
|
|
fieldNames = fieldnames(inputStruct);
|
|
|
|
for i = 1:numel(fieldNames)
|
|
thisField = fieldNames{i};
|
|
value = inputStruct.(thisField);
|
|
|
|
% Recursively fix nested structs
|
|
if isstruct(value)
|
|
value = fixBondToBound(value);
|
|
end
|
|
|
|
% If 'Bond' appears in the field name, rename it
|
|
if contains(thisField, 'Bond')
|
|
newField = strrep(thisField, 'Bond', 'Bound');
|
|
correctedStruct.(newField) = value;
|
|
correctedStruct = rmfield(correctedStruct, thisField);
|
|
else
|
|
correctedStruct.(thisField) = value;
|
|
end
|
|
end
|
|
end
|
|
|
|
function y = gaussianModel(params, x, nPeaks)
|
|
%
|
|
% construct equation with variable number of Gaussian peaks and linear
|
|
% background
|
|
y = zeros(size(x(:,1)));
|
|
|
|
|
|
% ------------------- model function ---------
|
|
% Add all Gaussian peaks
|
|
for i = 1:nPeaks
|
|
idx = (i-1)*3 + 1;
|
|
area = params(idx);
|
|
width = params(idx+1);
|
|
pos = params(idx+2);
|
|
|
|
y = y + area/(width*sqrt(pi/2)) * exp(-2*((x(:,1) - pos).^2/width^2));
|
|
end
|
|
% Add linear background
|
|
y = y + params(end-1)*x(:,1) + params(end);
|
|
|
|
end
|
|
|
|
function [initialParameterVector,lowerBoundsParameterVector,upperBoundsParameterVector] = gaussianModelBounds(FitParameter)
|
|
%
|
|
% Construct initial, upper bounds, and lower bounds arrays for gaussian
|
|
% model with variable peak number
|
|
|
|
% ------------------- initial parameters ---------
|
|
initialParameterVector = [];
|
|
for i = 1:size(FitParameter.peakArea,2)
|
|
initialParameterVector = [initialParameterVector, FitParameter.peakArea(1,i)];
|
|
initialParameterVector = [initialParameterVector, FitParameter.peakWidth(1,i)];
|
|
initialParameterVector = [initialParameterVector, FitParameter.peakPos(1,i)];
|
|
end
|
|
initialParameterVector = [initialParameterVector, 1, 0]; % linear background
|
|
|
|
% Build lower bounds vector with proper interleaving
|
|
lowerBoundsParameterVector = [];
|
|
for i = 1:size(FitParameter.peakArea,2)
|
|
lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakAreaLowerBound(1,i)];
|
|
lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakWidthLowerBound(1,i)];
|
|
lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakPosLowerBound(1,i)];
|
|
end
|
|
lowerBoundsParameterVector = [lowerBoundsParameterVector, -Inf, -Inf]; % linear background
|
|
|
|
% Build upper bounds vector with proper interleaving
|
|
upperBoundsParameterVector = [];
|
|
for i = 1:size(FitParameter.peakArea,2)
|
|
upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakAreaUpperBound(1,i)];
|
|
upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakWidthUpperBound(1,i)];
|
|
upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakPosUpperBound(1,i)];
|
|
end
|
|
upperBoundsParameterVector = [upperBoundsParameterVector, Inf, Inf]; % linear background
|
|
end
|
|
|
|
function y = pseudoVoigtModel(params, x, nPeaks)
|
|
y = zeros(size(x(:,1)));
|
|
|
|
% Add all pseudo-Voigt peaks
|
|
for i = 1:nPeaks
|
|
idx = (i-1)*4 + 1;
|
|
area = params(idx);
|
|
width = params(idx+1); % Gaussian width (sigma)
|
|
pos = params(idx+2);
|
|
eta = params(idx+3); % Mixing parameter
|
|
|
|
% Gaussian component
|
|
sigma = width;
|
|
gaussianComponent = (1./(sigma.*sqrt(2.*pi))) .* exp(-((x(:,1) - pos).^2)./(2.*sigma.^2));
|
|
|
|
% Lorentzian component
|
|
FWHM = sigma * 2 * sqrt(2*log(2));
|
|
gamma = FWHM/2;
|
|
lorentzianComponent = (1/pi) * (gamma ./ ((x(:,1) - pos).^2 + gamma.^2));
|
|
|
|
% Pseudo-Voigt peak
|
|
pseudoVoigtPeak = area .* (eta .* gaussianComponent + (1 - eta) .* lorentzianComponent);
|
|
y = y + pseudoVoigtPeak;
|
|
end
|
|
|
|
% Add linear background
|
|
y = y + params(end-1)*x(:,1) + params(end);
|
|
end
|
|
|
|
function [initialParameterVector,lowerBoundsParameterVector,upperBoundsParameterVector] = pseudoVoigtModelBounds(FitParameter)
|
|
%
|
|
% Construct initial, upper bounds, and lower bounds arrays for gaussian
|
|
% model with variable peak number
|
|
|
|
% ------------------- initial parameters ---------
|
|
initialParameterVector = [];
|
|
for i = 1:size(FitParameter.peakArea,2)
|
|
initialParameterVector = [initialParameterVector, FitParameter.peakArea(1,i)];
|
|
initialParameterVector = [initialParameterVector, FitParameter.peakWidth(1,i)];
|
|
initialParameterVector = [initialParameterVector, FitParameter.peakPos(1,i)];
|
|
initialParameterVector = [initialParameterVector, FitParameter.peakEta(1,i)];
|
|
end
|
|
initialParameterVector = [initialParameterVector, 1, 0]; % linear background
|
|
|
|
% Build lower bounds vector with proper interleaving
|
|
|
|
for i = 1:size(FitParameter.peakArea,2)
|
|
lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakAreaLowerBound(1,i)];
|
|
lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakWidthLowerBound(1,i)];
|
|
lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakPosLowerBound(1,i)];
|
|
lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakEtaLowerBound(1,i)];
|
|
end
|
|
lowerBoundsParameterVector = [lowerBoundsParameterVector, -Inf, -Inf]; % linear background
|
|
|
|
% Build upper bounds vector with proper interleaving
|
|
upperBoundsParameterVector = [];
|
|
for i = 1:size(FitParameter.peakArea,2)
|
|
upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakAreaUpperBound(1,i)];
|
|
upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakWidthUpperBound(1,i)];
|
|
upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakPosUpperBound(1,i)];
|
|
lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakEtaLowerBound(1,i)];
|
|
end
|
|
upperBoundsParameterVector = [upperBoundsParameterVector, Inf, Inf]; % linear background
|
|
end
|
|
|
|
|
|
|
|
function y = shirleyPseudoVoigtModel(params, x, nPeaks)
|
|
% Initialize model output
|
|
y = zeros(size(x(:,1)));
|
|
|
|
% Constant background level
|
|
bgConstant = params(end);
|
|
shirleyBackground = bgConstant * ones(size(x(:,1)));
|
|
|
|
% Process each peak
|
|
for i = 1:nPeaks
|
|
idx = (i-1)*5 + 1;
|
|
area = params(idx);
|
|
width = params(idx+1); % Gaussian width (sigma)
|
|
pos = params(idx+2);
|
|
eta = params(idx+3); % Mixing parameter (0=Lorentzian, 1=Gaussian)
|
|
k = params(idx+4); % Shirley scattering factor
|
|
|
|
% Calculate Gaussian component (using sigma directly)
|
|
sigma = width;
|
|
gaussianComponent = (1/(sigma*sqrt(2*pi))) * exp(-((x(:,1) - pos).^2)/(2*sigma^2));
|
|
|
|
% Calculate Lorentzian component
|
|
FWHM = sigma * 2 * sqrt(2*log(2)); % Convert sigma to FWHM
|
|
gamma = FWHM/2; % Lorentzian HWHM
|
|
lorentzianComponent = (1/pi) * (gamma ./ ((x(:,1) - pos).^2 + gamma^2));
|
|
|
|
% Pseudo-Voigt peak (weighted sum)
|
|
pseudoVoigtPeak = area * (eta * gaussianComponent + (1 - eta) * lorentzianComponent);
|
|
y = y + pseudoVoigtPeak;
|
|
|
|
% Add Shirley background contribution (if k > 0)
|
|
if k > 0
|
|
% Shirley background rises after peak (higher binding energy)
|
|
idxAfter = (x(:,1) >= pos);
|
|
|
|
if any(idxAfter)
|
|
% Calculate cumulative tail fractions
|
|
gaussianTail = 0.5 * (1 - erf((x(idxAfter,1) - pos)/(sqrt(2)*sigma)));
|
|
lorentzianTail = 0.5 - (1/pi)*atan((x(idxAfter,1) - pos)/gamma);
|
|
|
|
% Combined tail fraction for pseudo-Voigt
|
|
combinedTail = eta * gaussianTail + (1 - eta) * lorentzianTail;
|
|
|
|
% Shirley contribution for this peak
|
|
shirleyContrib = k * area * (0.5 - combinedTail);
|
|
|
|
% Add to background
|
|
shirleyTemp = zeros(size(x(:,1)));
|
|
shirleyTemp(idxAfter) = shirleyContrib;
|
|
shirleyBackground = shirleyBackground + shirleyTemp;
|
|
end
|
|
end
|
|
end
|
|
|
|
% Final model = peaks + Shirley background
|
|
y = y + shirleyBackground;
|
|
end
|
|
function [initialParameterVector,lowerBoundsParameterVector,upperBoundsParameterVector] = shirleyPseudoVoigtModelBounds(FitParameter)
|
|
%
|
|
% Construct initial, upper bounds, and lower bounds arrays for gaussian
|
|
% model with variable peak number
|
|
|
|
% ------------------- initial parameters ---------
|
|
initialParameterVector = [];
|
|
for i = 1:size(FitParameter.peakArea,2)
|
|
initialParameterVector = [initialParameterVector, FitParameter.peakArea(1,i)];
|
|
initialParameterVector = [initialParameterVector, FitParameter.peakWidth(1,i)];
|
|
initialParameterVector = [initialParameterVector, FitParameter.peakPos(1,i)];
|
|
initialParameterVector = [initialParameterVector, FitParameter.peakEta(1,i)];
|
|
initialParameterVector = [initialParameterVector, FitParameter.peakShirleyK(1,i)];
|
|
end
|
|
initialParameterVector = [initialParameterVector, 1, 0]; % linear background
|
|
|
|
% Build lower bounds vector with proper interleaving
|
|
lowerBoundsParameterVector = [];
|
|
for i = 1:size(FitParameter.peakArea,2)
|
|
lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakAreaLowerBound(1,i)];
|
|
lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakWidthLowerBound(1,i)];
|
|
lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakPosLowerBound(1,i)];
|
|
lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakEtaLowerBound(1,i)];
|
|
lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakShirleyKLowerBound(1,i)];
|
|
end
|
|
lowerBoundsParameterVector = [lowerBoundsParameterVector, -Inf, -Inf]; % linear background
|
|
|
|
% Build upper bounds vector with proper interleaving
|
|
upperBoundsParameterVector = [];
|
|
for i = 1:size(FitParameter.peakArea,2)
|
|
upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakAreaUpperBound(1,i)];
|
|
upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakWidthUpperBound(1,i)];
|
|
upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakPosUpperBound(1,i)];
|
|
upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakEtaUpperBound(1,i)];
|
|
upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakShirleyKUpperBound(1,i)];
|
|
end
|
|
|
|
end |