diff --git a/beamtime/commit.md b/beamtime/commit.md new file mode 100644 index 0000000..7af1f81 --- /dev/null +++ b/beamtime/commit.md @@ -0,0 +1,42 @@ +feat: Major refactor - generic model architecture with full Shirley background support + +BREAKING CHANGES: +- Single-spectrum fitting only (enforces exactly one row selected) +- Gauss fitting now uses FWHM parameter instead of Gaussian width sigma +- Removed multi-row batch processing capability + +NEW FEATURES: +- Generic model architecture: arbitrary number of peaks (determined by length of FitParameter.peakArea) +- Full Shirley background support for both Gaussian and Pseudo-Voigt peak shapes +- Comprehensive parameter validation system with detailed error messages +- Support for shared vs per-peak Shirley K parameters +- New fitting modes: 'positionShirley', 'widthShirley', 'widthPositionShirley' +- New custom mode: 'TwoSpinOrbitDoubletPseudoVoigtShirley' +- Automatic field name typo correction ('Bond' → 'Bound') +- Enhanced results export: now includes peakEta and peakShirleyK columns + +REFACTORING: +- Centralized peak and background equations in helper functions + - gaussPeak(), lorentzPeak(), pseudoVoigtPeak() + - linearBackground(), activeShirleyBackgroundFwhm() + - fwhm2w() conversion utility +- Extracted model generation to genericModelGenerator() +- Modular validation helpers: checkFieldDimension(), checkBoundsPair() +- Standardized parameter vector layouts documented for all mode/background/peak-type combinations +- Improved column resizing logic with resizeColumnHelper() +- Enhanced plotting integration with nappPlot.m and nappAddFitToPlot.m + +ARCHITECTURE: +- Unified parameter extraction for all fitting modes in centralized switch-case structure +- Consistent FWHM-based width parameterization across all peak types +- Mode-specific Shirley K sizing (scalar for shared modes, array for per-peak modes) +- Comprehensive inline documentation of parameter vector layouts (50+ combinations) + +IMPROVEMENTS: +- Backward compatibility for 'gauss' width mode (auto-converts to FWHM) +- Robust bounds checking with detailed diagnostic messages +- Enhanced figure generation: residuals plot + fit visualization in tiledlayout +- Better error handling for empty/malformed data structures +- NaN-padded result arrays for variable peak counts + +FILE SIZE: ~73KB → ~152KB (107% increase due to comprehensive documentation and generic architecture) diff --git a/beamtime/documentedExampleScripts/FileList.mat b/beamtime/documentedExampleScripts/FileList.mat index 226c049..e354ff1 100644 Binary files a/beamtime/documentedExampleScripts/FileList.mat and b/beamtime/documentedExampleScripts/FileList.mat differ diff --git a/beamtime/documentedExampleScripts/generateSpectra.m b/beamtime/documentedExampleScripts/generateSpectra.m index 2667838..d578d73 100644 --- a/beamtime/documentedExampleScripts/generateSpectra.m +++ b/beamtime/documentedExampleScripts/generateSpectra.m @@ -302,4 +302,11 @@ tmpDataSelection = month(FileList.lastModifiedDatenum) == 11 ... FileList = nappGenerateSpectra(FileList,tmpDataSelection); nappPlot(FileList,tmpDataSelection,'tableVariable','spectrum_countsPerSecond','mode','AllInOne'); disp(FileList(tmpDataSelection,{'name','lastModifiedDatestr','scientaRegionIterations','scientaSequenceIterations','scientaDwellTime_ms','scientaPassEnergy_eV'})) +end + +%% helper functions +% Check if the environment is MATLAB or Octave +if exist('OCTAVE_VERSION', 'builtin') +% Anonymous function mimicking contains +contains = @(str, pattern) ~cellfun('isempty', strfind(str, pattern)); end \ No newline at end of file diff --git a/beamtime/documentedExampleScripts/plotExportedData.m b/beamtime/documentedExampleScripts/plotExportedData.m new file mode 100644 index 0000000..7826b3d --- /dev/null +++ b/beamtime/documentedExampleScripts/plotExportedData.m @@ -0,0 +1,151 @@ +function plotExportedData(dataDirectory) +%plotExportedData Plot all exported ASCII data files in one figure +% +% This function reads all .txt files from the exported data (created by +% exportFigureData) and plots them in a single figure for visual comparison. +% File names are parsed to extract axis labels and legend entries. +% +% SYNTAX: +% plotExportedData() % Uses current working directory +% plotExportedData(dataDirectory) % Uses specified directory path +% +% INPUTS: +% dataDirectory : char or string (optional) +% Path to directory containing exported .txt files. +% If omitted, uses current working directory (pwd). +% +% OUTPUT: +% A new figure with all data plotted as lines. Axis labels and legend +% entries are reconstructed from the file names. +% +% NOTES: +% - Reads all .txt files in the specified directory +% - File name format expected: _[_]_vs__.txt +% - Plots are created with underscores converted back to spaces for display +% - X and Y axis labels are set based on parsed filename +% - Legend shows line labels from filenames +% - All data is plotted in a single axes for direct visual comparison +% +% EXAMPLE: +% % After running exportFigureData(): +% plotExportedData() +% % or +% plotExportedData('../exported_data/') + +% Handle input: use current directory if none provided +if nargin == 0 || isempty(dataDirectory) + dataDirectory = pwd; +end + +% Get all .txt files in the directory +fileStruct = dir(fullfile(dataDirectory, '*.txt')); +if isempty(fileStruct) + warning('plotExportedData:NoFiles', 'No .txt files found in %s', dataDirectory); + return; +end + +% Create a new figure +fig = figure('Name', 'Exported Data Comparison', 'NumberTitle', 'off'); +ax = axes(fig); +hold(ax, 'on'); +grid(ax, 'on'); + +% Storage for axis labels (to ensure consistency) +xLabelStr = ''; +yLabelStr = ''; +legendEntries = {}; + +% Loop over all .txt files +for ii = 1:length(fileStruct) + fileName = fileStruct(ii).name; + filePath = fullfile(dataDirectory, fileName); + + % Read the data + try + data = readmatrix(filePath); + catch ME + warning('plotExportedData:ReadFailed', 'Could not read file: %s\n%s', fileName, ME.message); + continue; + end + + % Check if data has at least 2 columns + if size(data, 2) < 2 || size(data, 1) < 1 + warning('plotExportedData:InvalidData', 'File %s does not have valid X,Y columns', fileName); + continue; + end + + xData = data(:, 1); + yData = data(:, 2); + + % Parse filename to extract labels + % Expected format: _[_]_vs__.txt + + % Remove .txt extension + nameNoExt = fileName(1:end-4); + + % Split by '_vs_' to separate y-axis label from the rest + parts = strsplit(nameNoExt, '_vs_'); + + if length(parts) >= 2 + % Extract line label (after _vs__) + rightPart = parts{2}; % _ + + % Find the last underscore to separate XLabel from LineLabel + underscores = find(rightPart == '_'); + if ~isempty(underscores) + lastUnderscore = underscores(end); + xLabelExtracted = rightPart(1:lastUnderscore-1); + lineLabel = rightPart(lastUnderscore+1:end); + else + xLabelExtracted = rightPart; + lineLabel = 'data'; + end + + % Extract y-axis label (before _vs_) + yLabelExtracted = parts{1}; + + % Remove figure title (first part before first underscore) + % This is tricky, so we'll keep what we have and note it's approximate + + % Convert underscores back to spaces for display + xLabelDisplay = strrep(xLabelExtracted, '_', ' '); + yLabelDisplay = strrep(yLabelExtracted, '_', ' '); + lineLabel = strrep(lineLabel, '_', ' '); + + % Store axis labels (use the first valid ones found) + if isempty(xLabelStr) + xLabelStr = xLabelDisplay; + end + if isempty(yLabelStr) + yLabelStr = yLabelDisplay; + end + + else + % Fallback if filename parsing fails + lineLabel = nameNoExt; + end + + % Plot the data + plot(ax, xData, yData, 'DisplayName', lineLabel, 'LineWidth', 1.5, 'Marker', 'o', 'MarkerSize', 4); + + % Store legend entry + legendEntries{end+1} = lineLabel; +end + +% Set axis labels +xlabel(ax, xLabelStr, 'Interpreter', 'none'); +ylabel(ax, yLabelStr, 'Interpreter', 'none'); +title(ax, 'Exported Data Comparison', 'Interpreter', 'none'); + +% Add legend +if ~isempty(legendEntries) + legend(ax, legendEntries, 'Interpreter', 'none', 'Location', 'best'); +end + +% Optimize figure appearance +grid(ax, 'on'); +hold(ax, 'off'); + +fprintf('Plotted %d data files from: %s\n', length(fileStruct), dataDirectory); + +end diff --git a/beamtime/nappAddFitToPlot.m b/beamtime/nappAddFitToPlot.m new file mode 100644 index 0000000..c2fc45c --- /dev/null +++ b/beamtime/nappAddFitToPlot.m @@ -0,0 +1,253 @@ +function nappAddFitToPlot(DataTable, varargin) +% nappAddFitToPlot - Plot fitted peaks with proper decomposition +% FIXED: Correct pseudo-Voigt normalization and FWHM extraction + +% Parse input arguments +p = inputParser; +addOptional(p, 'tableSelection', []); +addOptional(p, 'plotHandle', gca); +addOptional(p, 'peakWidthMode', 'FWHM'); +addOptional(p, 'energyScale', 'kinetic'); +parse(p, varargin{:}); + +% Get data from table +if isempty(p.Results.tableSelection) + fitData = DataTable(1, :); +else + fitData = DataTable(p.Results.tableSelection, :); +end + +% Get energy range from current figure xlim +ax = p.Results.plotHandle; +axes(ax); +xlims = xlim(ax); +energyMin = xlims(1); +energyMax = xlims(2); + +% Create energy range with 0.2 eV step +energyStep = 0.2; +energyRange = (energyMin:energyStep:energyMax)'; + +% ===================================================================== +% COUNT ONLY NON-NaN PEAKS +% ===================================================================== +peakArea_vals = fitData.peakArea; % Extract from cell +validPeaks = ~isnan(peakArea_vals); % Find non-NaN values +nPeaks = sum(validPeaks); % Count valid peaks + +fprintf('Found %d valid peaks (non-NaN)\n', nPeaks); + +if isempty(energyRange) || nPeaks == 0 + warning('No valid fit data to plot'); + return; +end + +% Initialize output +peakCurves = zeros(length(energyRange), nPeaks); +totalFit = zeros(size(energyRange)); +peakTypes = cell(nPeaks, 1); + +% Extract peak parameters (cell arrays) +peakPos_vals = fitData.peakPosition; +peakFWHM_vals = fitData.peakFWHM; + +% Check for peak eta parameter (pseudo-Voigt indicator) +hasPeakEta = ismember('peakEta', fitData.Properties.VariableNames) && ... + ~all(isnan(fitData.peakEta)); + +if hasPeakEta + peakEta_vals = fitData.peakEta; + +end + +% ===================================================================== +% MAIN LOOP: Process each VALID peak +% ===================================================================== +validIdx = find(validPeaks); % Get indices of valid peaks +for i = 1:nPeaks + iPeak = validIdx(i); % Actual peak index in array + + % Extract peak parameters + pos = peakPos_vals(iPeak); + A_total = peakArea_vals(iPeak); + FWHM = peakFWHM_vals(iPeak); + + % Skip if FWHM is NaN or invalid + if isnan(FWHM) || FWHM <= 0 + fprintf(' WARNING: Peak %d has invalid FWHM = %g, skipping\n', iPeak, FWHM); + continue; + end + + % Check if this is a pseudo-Voigt peak + if hasPeakEta && ~isnan(peakEta_vals(iPeak)) + % ============================================================ + % PSEUDO-VOIGT: Gaussian + Lorentzian blend + % eta = Gaussian fraction (0 = pure Lorentzian, 1 = pure Gaussian) + % ============================================================ + eta = peakEta_vals(iPeak); + + % === GAUSSIAN COMPONENT (FWHM-based) === + % G(E) = A_gauss * sqrt(4*ln(2)/pi) / FWHM * exp(-4*ln(2)*(E-E0)^2/FWHM^2) + % This is normalized so integral = A_gauss + A_gauss = eta * A_total; % Gaussian gets eta fraction of total area + gauss_norm = sqrt(4*log(2)/pi) / FWHM; + peakG = A_gauss .* gauss_norm .* ... + exp(-4*log(2) * ((energyRange - pos)./FWHM).^2); + + % === LORENTZIAN COMPONENT (FWHM-based) === + % L(E) = A_lorentz * (FWHM/2) / pi / ((E-E0)^2 + (FWHM/2)^2) + % This is normalized so integral = A_lorentz + A_lorentz = (1 - eta) * A_total; % Lorentzian gets (1-eta) fraction + gamma = FWHM / 2; % Half-width at half-maximum + peakL = A_lorentz .* (gamma / pi) ./ ... + ((energyRange - pos).^2 + gamma.^2); + + % === COMBINED PSEUDO-VOIGT === + peak = peakG + peakL; + + peakTypes{i} = sprintf('Pseudo-Voigt (η=%.3f)', eta); + + else + % ============================================================ + % PURE GAUSSIAN using FWHM + % G(E) = A * sqrt(4*ln(2)/pi) / FWHM * exp(-4*ln(2)*(E-E0)^2/FWHM^2) + % ============================================================ + gauss_norm = sqrt(4*log(2)/pi) / FWHM; + peak = A_total .* gauss_norm .* ... + exp(-4*log(2) * ((energyRange - pos)./FWHM).^2); + peakTypes{i} = 'Gaussian'; + end + + peakCurves(:, i) = peak; + totalFit = totalFit + peak; +end + +% ===================================================================== +% ADD BACKGROUND +% ===================================================================== +% Check for Shirley background parameters +hasShirleyK = ismember('peakShirleyK', fitData.Properties.VariableNames) && ... + ~all(isnan(fitData.peakShirleyK)); + +if hasShirleyK + % SHIRLEY BACKGROUND RECONSTRUCTION + peakShirleyK_vals = fitData.peakShirleyK; + + % Reconstruct Shirley background using the same algorithm as in fitting + bkgrd = zeros(size(energyRange)); + + for i = 1:nPeaks + iPeak = validIdx(i); + + % Extract parameters + area = peakArea_vals(iPeak); + FWHM = peakFWHM_vals(iPeak); + pos = peakPos_vals(iPeak); + if ismember('peakEta', fitData.Properties.VariableNames) && ~all(isnan(fitData.peakEta)) + eta = peakEta_vals(iPeak); + else + eta = 1; + end + K = peakShirleyK_vals(iPeak); + + % Find energy points after the peak + idxafter = energyRange >= pos; + + if any(idxafter) + % Calculate width parameters + w = FWHM / sqrt(log(4)); % Convert FWHM to Gaussian w + sigma = w / 2; + gamma = FWHM / 2; + + % Gaussian tail fraction + Gtail = 0.5 * (1 - erf((energyRange(idxafter) - pos) / (sqrt(2) * sigma))); + + % Lorentzian tail fraction + Ltail = 0.5 - (1/pi) * atan((energyRange(idxafter) - pos) / gamma); + + % Pseudo-Voigt tail fraction + tailfrac = eta * Gtail + (1 - eta) * Ltail; + + % Shirley contribution + shirleycontrib = K * area * (0.5 - tailfrac); + bkgrd(idxafter) = bkgrd(idxafter) + shirleycontrib; + + end + + end + + % Add constant offset C0 + if ismember('peakBackground', fitData.Properties.VariableNames) + C0 = fitData.peakBackground(2); % Second element is C0 + bkgrd = bkgrd + C0; + end + + if matches(p.Results.energyScale,'kinetic') + % This is a little work around - needs confirmation if it does what + % it should - be careful. + bkgrd = flip(bkgrd); + end + + backgroundType = 'Shirley'; +else + % LINEAR BACKGROUND (existing code) + if ismember('peakBackground', fitData.Properties.VariableNames) + bkgrd_vals = fitData.peakBackground; + if length(bkgrd_vals) >= 2 + slope = bkgrd_vals(1); + intercept = bkgrd_vals(2); + else + slope = 0; + intercept = bkgrd_vals(1); + end + else + slope = 0; + intercept = 0; + end + + bkgrd = slope .* energyRange + intercept; + backgroundType = sprintf('Linear (slope=%.4f, intercept=%.4f)', slope, intercept); +end + +totalFit = totalFit + bkgrd; + + +% ===================================================================== +% PLOT RESULTS +% ===================================================================== +axes(ax); +hold on; + +% Plot total fit (thick red line) +plot(energyRange, totalFit, '-', 'LineWidth', 2.5, 'Color', [0.8 0.2 0.2], 'DisplayName', 'Total Fit'); + +% Plot background (dotted gray line) +plot(energyRange, bkgrd, ':', 'LineWidth', 1.5, 'Color', [0.5 0.5 0.5], 'DisplayName', 'Background'); + +% Plot individual peaks (dashed colored lines) +colors = lines(nPeaks); +for i = 1:nPeaks + if any(peakCurves(:, i)) % Only plot if non-zero + plot(energyRange, peakCurves(:, i)+bkgrd, '--', 'LineWidth', 1.2, ... + 'Color', colors(i,:), 'DisplayName', sprintf('Peak %d', i)); + end +end + +legend('show', 'Location', 'best'); +hold off; + +% % ===================================================================== +% % DISPLAY INFORMATION +% % ===================================================================== +% fprintf('\n=== FIT PLOTTED SUCCESSFULLY ===\n'); +% fprintf('Valid peaks: %d\n', nPeaks); +% fprintf('Peak type(s): '); +% for i = 1:nPeaks +% fprintf('%s', peakTypes{i}); +% if i < nPeaks, fprintf('; '); end +% end +% fprintf('\n'); +% fprintf('Energy range: %.1f - %.1f eV\n', min(energyRange), max(energyRange)); +% fprintf('=====================================\n\n'); +% fprintf('Background: %s\n', backgroundType); +end diff --git a/beamtime/nappCombineSpectra.m b/beamtime/nappCombineSpectra.m index 9e93d59..7c117f7 100644 --- a/beamtime/nappCombineSpectra.m +++ b/beamtime/nappCombineSpectra.m @@ -201,16 +201,16 @@ else end % check if x has different scaling or size -for ii = 1:size(p.Results.DataTable(p.Results.tableSelection,:),1) - testX(:,ii) = tmpX{ii}; -end +% for ii = 1:size(p.Results.DataTable(p.Results.tableSelection,:),1) +% testX(:,ii) = tmpX{ii}; +% end -if all(all(testX == testX(1,:))) - xq = tmpX{1}; - for ii = 1:size(p.Results.DataTable(p.Results.tableSelection,:),1) - vq(:,ii) = p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii}; - end -else +%if all(all(testX == testX(1,:))) +% xq = tmpX{1}; +% for ii = 1:size(p.Results.DataTable(p.Results.tableSelection,:),1) +% vq(:,ii) = p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii}; +% end +%else xq = tmpX{p.Results.tableRow}(tmpStartIndex:tmpEndIndex); % repeat this for all columns in input. for ii = 1:size(p.Results.DataTable(p.Results.tableSelection,:),1) @@ -219,7 +219,7 @@ else vq(:,ii) = interp1(x,v,xq); end -end +% end if contains(p.Results.output,'simple') combinedSpectrum(:,1) = xq; if contains(p.Results.mode,'sum') diff --git a/beamtime/nappFitGaussXps.m b/beamtime/nappFitGaussXps.m index b28757a..ec505a7 100644 --- a/beamtime/nappFitGaussXps.m +++ b/beamtime/nappFitGaussXps.m @@ -1,176 +1,583 @@ function DataTable = nappFitGaussXps(DataTable,varargin) -% nappFitGaussXps fits photoelectron (PE) spectra with up to 4 or 5 Gaussian peaks. -% It allows different sets of constraints for the fitting process. +%NAPPFITGAUSSXPSDEVELOPMENT Fit XPS spectra with Gaussian/Pseudo-Voigt peaks. +% NAPPFITGAUSSXPS fits XPS data with Gauss or pseudo-voigt peaks and +% linear or Shirley background using initial guesses and bounds defined +% in FitParameter. % -% INPUT: -% DataTable: Table containing the spectral data to be fitted. -% varargin: Additional optional parameters for customizing the fit. +% DataTable = NAPPFITGAUSSXPSDEVELOPMENT(DataTable, FitParameter, tableSelection) +% fits the single row specified by the logical vector tableSelection. +% The function enforces that exactly one row is selected at a time. % -% OUTPUT: -% DataTable: The input table, now including the results of the peak fitting. +% DataTable = NAPPFITGAUSSXPSDEVELOPMENT(..., 'Name', Value) specifies +% additional options using name-value pair arguments: % -% FileList = -% nappFitGaussXps(FileList,tmpFitParameters,'tableSelection',tmpDataSelection); -%% -%-------------------------------------------------------------------------- -% Variable Definitions and Purposes -%-------------------------------------------------------------------------- +% Name-Value Pair Arguments: +% -------------------------- +% 'tableVariable' - Column name containing intensity data. +% Options: 'spectrum_countsPerSecond' (default), +% 'spectrum_counts', 'spectrum_countsNorm'. % -% 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. +% 'mode' - Fitting constraint mode controlling parameter sharing. +% Options: 'widthPosition' (default), 'unconstrained', +% 'position', 'width', 'widthPositionShirley', +% 'positionShirley','widthShirley', +% 'unconstrainedAndWidthPosition', +% 'gasAndTwoAcidBasePairs', +% 'widthPositionAreaRatio', +% 'widthPositionAmineAreaRatio'. +% +% 'energyScale' - Energy axis type for X-data. +% Options: 'kinetic' (default), 'binding'. +% +% 'backgroundMode' - Background subtraction model. +% Options: 'linear' (default), 'shirley'. +% +% 'showFigures' - Display fit figures and results. +% Options: 'on' (default), 'off'. +% +% 'peakWidthMode' - Input format for peak width parameter. +% Options: 'gauss' (default), 'FWHM'. +% +% NEW in 2026: +% ------------ +% Major upgrade to the code +% - generic models run with any number of peaks, length of +% FitParams.peakArea defines their number +% - Gauss and PV both work with Shirley and linear background +% - Gauss fitting now based on FWHM, no longer Gauss width +% - Plotting uses nappPlot.m and addFitToPlot.m function, no longer +% hard-coded +% - Functions and background equations are coded in functions at the end of the +% document. One place to modify or validate them. +% +% INPUT DATA REQUIREMENTS: +% ------------------------ +% DataTable must be a MATLAB table containing: +% - Energy column: 'scientaEkin_eV' or 'bindingEnergy_eV' +% - Intensity column: as specified by 'tableVariable' +% +% FitParameter must be a struct with fields (1xnPeaks row vectors): +% - peakArea, peakPos, peakWidth : Initial guesses +% - peakAreaLowerBound, peakAreaUpperBound : Area constraints +% - peakPosLowerBound, peakPosUpperBound : Position constraints +% - peakWidthLowerBound, peakWidthUpperBound : Width constraints +% - backgroundRange : [low high] energy range for background +% +% Optional FitParameter fields for Pseudo-Voigt: +% - peakEta, peakEtaLowerBound, peakEtaUpperBound +% +% Optional FitParameter fields for Shirley background: +% - peakShirleyK, peakShirleyKLowerBound, peakShirleyKUpperBound +% - peakShirleyOffset +% +% OUTPUT: +% ------- +% DataTable - Input table with appended fit result columns. +% +% FITTING MODES: +% -------------- +% 'unconstrained' - All peaks independent +% 'position' - Independent widths; linked positions +% 'widthPosition' - Shared width; linked positions +% 'width' - Shared width; independent positions +% +% EXAMPLE: +% -------- +% selection = false(height(MyData),1); selection(5) = true; +% FitParams.peakArea = [1000 500]; +% FitParams.peakPos = [710 725]; +% FitParams.peakWidth = 2.0; +% FitParams.peakAreaLowerBound = [0 0]; +% FitParams.peakAreaUpperBound = [Inf Inf]; +% FitParams.peakWidthLowerBound = 0.5; +% FitParams.peakWidthUpperBound = 5.0; +% FitParams.peakPosLowerBound = 700; +% FitParams.peakPosUpperBound = 720; +% FitParams.backgroundRange = [10 10]; +% Result = nappFitGaussXpsDevelopment(MyData, FitParams, selection, ... +% 'mode', 'widthPosition', 'energyScale', 'binding'); +% +% See also INPUTPARSER, GLOBALSEARCH, FMINCON. +% ========================================================================= +% DESIGN PHILOSOPHY: Physics-Driven Model Customization +% ========================================================================= +% +% The core design principle of nappFitGaussXps is to enable researchers to embed +% physical knowledge directly into the fitting model. Instead of forcing all analyses +% into predefined constraint patterns, users can construct custom fitting modes that +% reflect their specific physical interpretation—such as linking peak positions via +% fixed energy separations (acid-base pair splitting), constraining area ratios +% (spin-orbit doublet intensities), or imposing width sharing based on mechanistic +% understanding. This flexibility is achieved through straightforward parameter +% vector definitions and simple switch-case model definitions at the end of the +% function (see CUSTOM MODES section below). For users who prefer standard +% approaches without customization, the function also provides several well-tested +% generic models (unconstrained, position, width, widthPosition) that support any +% number of peaks with commonly-used constraint strategies. Both approaches—physics- +% based customization and generic flexibility—can leverage the same underlying +% infrastructure (Gaussian, Pseudo-Voigt, linear, and Shirley background models). +% +% =========================================================================% ========================================================================= +% PARAMETER VECTOR LAYOUTS FOR ALL MODE / BACKGROUND / PEAK-TYPE COMBINATIONS +% ========================================================================= +% +% NOTE: +% - All layouts below use FWHM as the width parameter (not Gaussian sigma). +% - "pos1_abs" means the absolute position of peak 1 (fit parameter). +% - "pos2_rel", "pos4_rel", etc. are offsets relative to a reference peak. +% - "width_shared" / "FWHM_shared" means a single FWHM shared across peaks. +% - "K_i" are Shirley K parameters per peak, "C0" is the Shirley intercept. +% - "slope" and "intercept" are linear background parameters: BG(E) = slope*E + intercept. +% - nPeaks is determined from numel(FitParameter.peakArea). +% - All generic modes (unconstrained, position, width, widthPosition) support any nPeaks >= 1. +% - Custom modes (unconstrainedAndWidthPosition, gasAndTwoAcidBasePairs, etc.) are hard-coded for specific peak counts. +% +% ========================================================================= +% GENERIC MODES (support arbitrary number of peaks) +% ========================================================================= +% +% ------------------------------------------------------------------------- +% UNCONSTRAINED MODE +% - All peaks: independent FWHM, independent position. +% - All FWHM values are free parameters. +% - All peak positions are free parameters. +% ------------------------------------------------------------------------- +% +% ===== GAUSSIAN + LINEAR (UNCONSTRAINED, FWHM as parameter) ===== +% PARAMETER VECTOR: +% [area1, FWHM1, pos1, area2, FWHM2, pos2, ..., areaN, FWHMN, posN, ... +% slope, intercept] +% +% ===== PSEUDO-VOIGT + LINEAR (UNCONSTRAINED, FWHM as parameter) ===== +% PARAMETER VECTOR: +% [area1, FWHM1, pos1, eta1, ... +% area2, FWHM2, pos2, eta2, ... +% ... +% areaN, FWHMN, posN, etaN, ... +% slope, intercept] +% +% ===== GAUSSIAN + SHIRLEY (UNCONSTRAINED, FWHM as parameter) ===== +% PARAMETER VECTOR: +% [area1, FWHM1, pos1, ... +% area2, FWHM2, pos2, ... +% ... +% areaN, FWHMN, posN, ... +% K1, K2, ..., KN, ... +% C0] +% +% ===== PSEUDO-VOIGT + SHIRLEY (UNCONSTRAINED, FWHM as parameter) ===== +% PARAMETER VECTOR: +% [area1, FWHM1, pos1, eta1, ... +% area2, FWHM2, pos2, eta2, ... +% ... +% areaN, FWHMN, posN, etaN, ... +% K1, K2, ..., KN, ... +% C0] % % -% 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. +% ------------------------------------------------------------------------- +% POSITION MODE +% - Peak 1: absolute position (fit parameter), independent FWHM. +% - Peaks 2..N: positions relative to peak 1 via FitParameter.peakPos(1,i), +% independent FWHM. +% - All FWHM values are independent (no sharing). +% - Shirley mode: K1, K2, ..., KN per-peak (NOT shared). +% ------------------------------------------------------------------------- +% +% ===== GAUSSIAN + LINEAR (POSITION MODE, FWHM as parameter) ===== +% PARAMETER VECTOR: +% [area1, FWHM1, area2, FWHM2, ..., areaN, FWHMN, ... +% pos1_abs, ... +% slope, intercept] +% +% ===== PSEUDO-VOIGT + LINEAR (POSITION MODE, FWHM as parameter) ===== +% PARAMETER VECTOR: +% [area1, FWHM1, eta1, ... +% area2, FWHM2, eta2, ... +% ... +% areaN, FWHMN, etaN, ... +% pos1_abs, ... +% slope, intercept] +% +% ===== GAUSSIAN + SHIRLEY (POSITION MODE, FWHM as parameter) ===== +% PARAMETER VECTOR: +% [area1, FWHM1, area2, FWHM2, ..., areaN, FWHMN, ... +% pos1_abs, ... +% K1, K2, ..., KN, ... +% C0] +% +% ===== PSEUDO-VOIGT + SHIRLEY (POSITION MODE, FWHM as parameter) ===== +% PARAMETER VECTOR: +% [area1, FWHM1, eta1, ... +% area2, FWHM2, eta2, ... +% ... +% areaN, FWHMN, etaN, ... +% pos1_abs, ... +% K1, K2, ..., KN, ... +% C0] % % -% 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. +% ------------------------------------------------------------------------- +% POSITION-SHIRLEY MODE +% - Peak 1: absolute position (fit parameter). +% - Peaks 2..N: positions relative to peak 1. +% - All peaks: independent FWHM. +% - Shirley background ONLY (no linear option). +% - K value is SHARED (single K for all peaks, NOT per-peak). +% ------------------------------------------------------------------------- +% +% ===== GAUSSIAN + SHIRLEY (POSITION-SHIRLEY MODE, Shirley shared) ===== +% PARAMETER VECTOR: +% [area1, FWHM1, area2, FWHM2, ..., areaN, FWHMN, ... +% pos1_abs, ... +% K, ... +% C0] +% +% ===== PSEUDO-VOIGT + SHIRLEY (POSITION-SHIRLEY MODE, Shirley shared) ===== +% PARAMETER VECTOR: +% [area1, FWHM1, eta1, ... +% area2, FWHM2, eta2, ... +% ... +% areaN, FWHMN, etaN, ... +% pos1_abs, ... +% K, ... +% C0] % % -% Color Palette Variables -% ---------------------- -% PsiFarbe - Structure holding RGB color values for plotting. -% PsiFarbe.hellgrau, etc. - Individual color entries for figure styling. +% ------------------------------------------------------------------------- +% WIDTH-POSITION MODE +% - All peaks: shared FWHM (single FWHM_shared for all). +% - Peak 1: absolute position (fit parameter). +% - Peaks 2..N: positions relative to peak 1 via FitParameter.peakPos(1,i). +% - Shirley mode: K1, K2, ..., KN per-peak (NOT shared). +% ------------------------------------------------------------------------- +% +% ===== GAUSSIAN + LINEAR (WIDTH-POSITION MODE, FWHM as parameter) ===== +% PARAMETER VECTOR: +% [area1, area2, ..., areaN, ... +% FWHM_shared, ... +% pos1_abs, ... +% slope, intercept] +% +% ===== PSEUDO-VOIGT + LINEAR (WIDTH-POSITION MODE, FWHM as parameter) ===== +% PARAMETER VECTOR: +% [area1, eta1, ... +% area2, eta2, ... +% ... +% areaN, etaN, ... +% FWHM_shared, ... +% pos1_abs, ... +% slope, intercept] +% +% ===== GAUSSIAN + SHIRLEY (WIDTH-POSITION MODE, FWHM as parameter) ===== +% PARAMETER VECTOR: +% [area1, area2, ..., areaN, ... +% FWHM_shared, ... +% pos1_abs, ... +% K1, K2, ..., KN, ... +% C0] +% +% ===== PSEUDO-VOIGT + SHIRLEY (WIDTH-POSITION MODE, FWHM as parameter) ===== +% PARAMETER VECTOR: +% [area1, eta1, ... +% area2, eta2, ... +% ... +% areaN, etaN, ... +% FWHM_shared, ... +% pos1_abs, ... +% K1, K2, ..., KN, ... +% C0] % % -% 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. +% ------------------------------------------------------------------------- +% WIDTH-POSITION-SHIRLEY MODE +% - All peaks: shared FWHM (single FWHM_shared for all). +% - Peak 1: absolute position (fit parameter). +% - Peaks 2..N: positions relative to peak 1. +% - Shirley background ONLY (no linear option). +% - K value is SHARED (single K for all peaks, NOT per-peak). +% ------------------------------------------------------------------------- +% +% ===== GAUSSIAN + SHIRLEY (WIDTH-POSITION-SHIRLEY MODE, Shirley shared) ===== +% PARAMETER VECTOR: +% [area1, area2, ..., areaN, ... +% FWHM_shared, ... +% pos1_abs, ... +% K, ... +% C0] +% +% ===== PSEUDO-VOIGT + SHIRLEY (WIDTH-POSITION-SHIRLEY MODE, Shirley shared) ===== +% PARAMETER VECTOR: +% [area1, eta1, ... +% area2, eta2, ... +% ... +% areaN, etaN, ... +% FWHM_shared, ... +% pos1_abs, ... +% K, ... +% C0] % % -% Spectral Data Variables -% ---------------------- -% fitEnergyValues - Cropped energy values used for fitting. -% fitIntensityValues - Cropped intensity values used for fitting. -% tmp1, tmp2 - Indices for cropping fit region. +% ------------------------------------------------------------------------- +% WIDTH MODE +% - All peaks: shared FWHM (single FWHM_shared for all). +% - All peaks: independent absolute positions. +% - Shirley mode: K1, K2, ..., KN per-peak (NOT shared). +% ------------------------------------------------------------------------- +% +% ===== GAUSSIAN + LINEAR (WIDTH MODE, FWHM as parameter) ===== +% PARAMETER VECTOR: +% [area1, pos1_abs, ... +% area2, pos2_abs, ... +% ... +% areaN, posN_abs, ... +% FWHM_shared, ... +% slope, intercept] +% +% ===== PSEUDO-VOIGT + LINEAR (WIDTH MODE, FWHM as parameter) ===== +% PARAMETER VECTOR: +% [area1, pos1_abs, eta1, ... +% area2, pos2_abs, eta2, ... +% ... +% areaN, posN_abs, etaN, ... +% FWHM_shared, ... +% slope, intercept] +% +% ===== GAUSSIAN + SHIRLEY (WIDTH MODE, FWHM as parameter) ===== +% PARAMETER VECTOR: +% [area1, pos1_abs, ... +% area2, pos2_abs, ... +% ... +% areaN, posN_abs, ... +% FWHM_shared, ... +% K1, K2, ..., KN, ... +% C0] +% +% ===== PSEUDO-VOIGT + SHIRLEY (WIDTH MODE, FWHM as parameter) ===== +% PARAMETER VECTOR: +% [area1, pos1_abs, eta1, ... +% area2, pos2_abs, eta2, ... +% ... +% areaN, posN_abs, etaN, ... +% FWHM_shared, ... +% K1, K2, ..., KN, ... +% C0] % % -% Loop and Indexing Variables -% -------------------------- -% ii - Loop index over selected spectra (rows in DataTable). +% ------------------------------------------------------------------------- +% WIDTH-SHIRLEY MODE +% - All peaks: shared FWHM (single FWHM_shared for all). +% - All peaks: independent absolute positions. +% - Shirley background ONLY (no linear option). +% - K value is SHARED (single K for all peaks, NOT per-peak). +% ------------------------------------------------------------------------- % -% Fit Parameters -% ------------- -% tmpFitParameters.backgroundRange = [5 5]; cuts the spectra to -5 +5 from -% the first peak position for better background fitting. -% tmpFitParameters.weight = [647.5 648 0.1]; allows to define a -% region in energy between which a lower weight is used for fitting. Third -% parameter defines the weight (1 is maximum and default). +% ===== GAUSSIAN + SHIRLEY (WIDTH-SHIRLEY MODE, Shirley shared) ===== +% PARAMETER VECTOR: +% [area1, pos1_abs, ... +% area2, pos2_abs, ... +% ... +% areaN, posN_abs, ... +% FWHM_shared, ... +% K, ... +% C0] % +% ===== PSEUDO-VOIGT + SHIRLEY (WIDTH-SHIRLEY MODE, Shirley shared) ===== +% PARAMETER VECTOR: +% [area1, pos1_abs, eta1, ... +% area2, pos2_abs, eta2, ... +% ... +% areaN, posN_abs, etaN, ... +% FWHM_shared, ... +% K, ... +% C0] +% +% +% ========================================================================= +% CUSTOM MODES (hard-coded constraint patterns) +% ========================================================================= +% +% ------------------------------------------------------------------------- +% CUSTOM MODE: unconstrainedAndWidthPosition +% - 5-peak hard-coded constraint pattern. +% - Peak 1 and 2: share width (FWHM1), positions coupled by shift parameter. +% - Peak 3: independent width (FWHM3), independent position. +% - Peaks 4 and 5: share width with peak 1, independent relative positions. +% ------------------------------------------------------------------------- +% +% ===== GAUSSIAN + LINEAR (UNCONSTRAINED+WIDTH-POSITION HYBRID, 5 PEAKS) ===== +% PARAMETER VECTOR: +% [area1, FWHM1, pos2_rel, ... +% area2, ... +% area3, FWHM3, pos3_abs, ... +% area4, pos4_rel, ... +% area5, pos5_rel, ... +% slope, intercept, ... +% pos1_abs_shift] +% +% ===== PSEUDO-VOIGT + LINEAR (UNCONSTRAINED+WIDTH-POSITION HYBRID, 5 PEAKS) ===== +% PARAMETER VECTOR: +% [area1, FWHM1, pos2_rel, eta1, ... +% area2, eta2, ... +% area3, FWHM3, pos3_abs, eta3, ... +% area4, pos4_rel, eta4, ... +% area5, pos5_rel, eta5, ... +% slope, intercept, ... +% pos1_abs_shift] +% +% ===== PSEUDO-VOIGT + SHIRLEY (UNCONSTRAINED+WIDTH-POSITION HYBRID, 5 PEAKS) ===== +% PARAMETER VECTOR: +% [area1, FWHM1, pos2_rel, eta1, K1, ... +% area2, eta2, K2, ... +% area3, FWHM3, pos3_abs, eta3, K3, ... +% area4, pos4_rel, eta4, K4, ... +% area5, pos5_rel, eta5, K5, ... +% C0, ... +% pos1_abs_shift] +% +% +% ------------------------------------------------------------------------- +% CUSTOM MODE: gasAndTwoAcidBasePairs +% - 5-peak hard-coded pattern: one gas peak + two acid-base pairs. +% - Peak 1: gas peak, independent width, positions of peaks 2/4 coupled to it via shift. +% - Peak 2: first acid (bound to peak 1). +% - Peak 3: base pair, independent width and position. +% - Peaks 4 and 5: second acid-base pair (relative positions). +% ------------------------------------------------------------------------- +% +% ===== GAUSSIAN + LINEAR (GAS + TWO ACID-BASE PAIRS, 5 PEAKS) ===== +% PARAMETER VECTOR: +% [area1, FWHM1, pos2_rel, ... +% area2, ... +% area3, FWHM3, pos3_abs, ... +% area4, pos4_rel, ... +% area5, pos5_rel, ... +% slope, intercept, ... +% pos1_abs_shift] +% +% ===== PSEUDO-VOIGT + LINEAR (GAS + TWO ACID-BASE PAIRS, 5 PEAKS) ===== +% PARAMETER VECTOR: +% [area1, FWHM1, pos2_rel, eta1, ... +% area2, eta2, ... +% area3, FWHM3, pos3_abs, eta3, ... +% area4, pos4_rel, eta4, ... +% area5, pos5_rel, eta5, ... +% slope, intercept, ... +% pos1_abs_shift] +% +% ===== PSEUDO-VOIGT + SHIRLEY (GAS + TWO ACID-BASE PAIRS, 5 PEAKS) ===== +% PARAMETER VECTOR: +% [area1, FWHM1, pos2_rel, eta1, K1, ... +% area2, eta2, K2, ... +% area3, FWHM3, pos3_abs, eta3, K3, ... +% area4, pos4_rel, eta4, K4, ... +% area5, pos5_rel, eta5, K5, ... +% C0, ... +% pos1_abs_shift] +% +% +% ------------------------------------------------------------------------- +% CUSTOM MODE: TwoSpinOrbitDoublet +% - 4-peak model: two spin-orbit doublets (each doublet has a fixed area ratio +% and energy separation, representing J=1/2 and J=3/2 spin-orbit splitting). +% - Peak 1: first doublet, intensity I (area1). +% - Peak 2: first doublet satellite, intensity I*R (fixed ratio), energy shift fixed. +% - Peak 3: second doublet, intensity I3 (area3). +% - Peak 4: second doublet satellite, intensity I3*R (fixed ratio), energy shift fixed. +% - All peaks share the same FWHM. +% ------------------------------------------------------------------------- +% +% ===== GAUSSIAN + LINEAR (TWO SPIN-ORBIT DOUBLET, 4 PEAKS) ===== +% PARAMETER VECTOR: +% [area1, FWHM_shared, pos1_abs, ... +% area3, pos3_abs, ... +% slope, intercept] +% +% ===== PSEUDO-VOIGT + SHIRLEY (TWO SPIN-ORBIT DOUBLET, 4 PEAKS) ===== +% PARAMETER VECTOR: +% [area1, FWHM_shared, pos1_abs, eta_shared, ... +% area3, pos3_abs, ... +% K_shared, ... +% C0] +% +% +% ------------------------------------------------------------------------- +% CUSTOM MODE: widthPositionAreaRatio +% - 2-peak model: second peak area and position are constrained by fixed ratios. +% - Peak 1: free area, free position, shared width. +% - Peak 2: area = area1 * FitParameter.peakArea(1,2) (fixed ratio), +% position = pos1 + FitParameter.peakPos(1,2) (fixed offset). +% ------------------------------------------------------------------------- +% +% ===== GAUSSIAN + LINEAR (2 PEAKS, AREA RATIO, SHARED WIDTH) ===== +% PARAMETER VECTOR: +% [area1, FWHM_shared, pos1_abs, ... +% slope, intercept] +% +% ===== PSEUDO-VOIGT + LINEAR (2 PEAKS, AREA RATIO, SHARED WIDTH) ===== +% PARAMETER VECTOR: +% [area1, FWHM_shared, eta_shared, ... +% pos1_abs, ... +% slope, intercept] +% +% ===== PSEUDO-VOIGT + SHIRLEY (2 PEAKS, AREA RATIO, SHARED WIDTH) ===== +% PARAMETER VECTOR: +% [area1, FWHM_shared, eta_shared, K_shared, ... +% pos1_abs, ... +% C0] +% +% +% ------------------------------------------------------------------------- +% CUSTOM MODE: widthPositionAmineAreaRatio +% - 5-peak model: specialized for amine deconvolution with constrained areas. +% - Peak 1: primary amine (free area, free position), shared width. +% - Peak 2: secondary amine (area constrained), relative position. +% - Peak 3: secondary amine satellite (area constrained), relative position. +% - Peaks 4 and 5: tertiary amines, shared width. +% ------------------------------------------------------------------------- +% +% ===== GAUSSIAN + LINEAR (AMINE AREA RATIO, 5 PEAKS) ===== +% PARAMETER VECTOR: +% [area1, FWHM_shared, pos1_abs, ... +% area2, ... +% area2*ratio, ... +% area4, ... +% slope, intercept] +% +% +% ========================================================================= +% SYMBOLS AND SHORT-HANDS +% ========================================================================= +% +% area_i : area (intensity integral) of peak i +% FWHM_i : full width at half maximum of peak i +% FWHM_shared : single shared FWHM used by all (or multiple) peaks +% pos_i : center energy of peak i +% pos1_abs : absolute fitted position of peak 1 (fit parameter) +% pos*_rel : relative energy offset from a reference peak +% (calculated from FitParameter.peakPos, NOT fit parameter) +% pos1_abs_shift : additional shift applied to peak 1 position in custom modes +% eta_i : Pseudo-Voigt mixing parameter of peak i (0 = pure Gaussian, +% 1 = pure Lorentzian) +% eta_shared : single shared eta for all peaks +% K_i : Shirley background K parameter for peak i +% K : shared Shirley K parameter (used in shared modes) +% C0 : Shirley background intercept constant +% slope : linear background slope (coefficient of energy) +% intercept : linear background y-intercept +% +% ========================================================================= -%% 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; -% -% % „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]; -% -% % „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; -% -% % „ 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]; -% -% % „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]; -% -% % „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]; -%% Off we go - code starts here +%% Define function input -% Define conditions and defaults of input +% Define input pairs and defaults defaultTableVariable = 'spectrum_countsPerSecond'; expectedTableVariable = {'spectrum_counts','spectrum_countsPerSecond','spectrum_countsNorm'}; defaultMode = 'widthPosition'; -expectedMode = {'unconstrained','widthPosition','position','width','unconstrainedAndWidthPosition','gasAndTwoAcidBasePairs','widthPositionAreaRatio','widthPositionAmineAreaRatio','TwoSpinOrbitDoublet'}; +expectedMode = {'unconstrained','widthPosition','position','width','widthPositionShirley', ... + 'positionShirley','widthShirley','unconstrainedAndWidthPosition','gasAndTwoAcidBasePairs',... + 'widthPositionAreaRatio','widthPositionAmineAreaRatio','TwoSpinOrbitDoublet','TwoSpinOrbitDoubletPseudoVoigtShirley'}; defaultTableSelection = true(size(DataTable,1),1); defaultShowFigures = 'on'; expectedShowFigures = {'on','off'}; @@ -178,7 +585,8 @@ defaultEnergyScale = 'kinetic'; expectedEnergyScale = {'kinetic','binding'}; defaultBackgroundMode = 'linear'; expectedBackgroundMode = {'linear','shirley'}; - +defaultPeakWidthMode = 'gauss'; +expectedPeakWidthMode = {'gauss','FWHM'}; % Parse input @@ -195,281 +603,1328 @@ addParameter(p,'energyScale',defaultEnergyScale,... @(energyValues) any(validatestring(energyValues,expectedEnergyScale))); addParameter(p,'backgroundMode',defaultBackgroundMode,... @(energyValues) any(validatestring(energyValues,expectedBackgroundMode))); +addParameter(p,'peakWidthMode',defaultPeakWidthMode,... + @(energyValues) any(validatestring(energyValues,expectedPeakWidthMode))); 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); +% backward comparability +if matches(p.Results.peakWidthMode,'gauss') + % input is in gauss-width, this function requires FWHM: + FitParameter.peakWidth = FitParameter.peakWidth .* sqrt(log(4)); +end +% ========================================================================= +%% VALIDATION +% ========================================================================= +% This code section has been written by Perplexity. Use with care. Its +% purpose is to give detailed error messages if the input should be false. -% load color scheme -PsiFarbe = getPsiColors(); +selectedRowCount = sum(p.Results.tableSelection); +if selectedRowCount ~= 1 + %% single row + error('ERROR: tableSelection must contain EXACTLY ONE row. Found %d rows selected. This function fits only single spectra at a time.', selectedRowCount); +end +selectedRowIdx = find(p.Results.tableSelection, 1); +fprintf('Fitting row %d of DataTable...\n', selectedRowIdx); + +if matches(p.Results.backgroundMode, 'shirley') + %% background shirley + if ~isfield(FitParameter, 'peakShirleyK') + error('ERROR: backgroundMode is set to ''shirley'' but FitParameter does not have peakShirleyK field. '); + end +end + +if matches(p.Results.backgroundMode, 'linear') + %% background linear + if isfield(FitParameter, 'peakShirleyK') + error('ERROR: backgroundMode is set to ''linear'' but FitParameter does specify peakShirleyK field. '); + end +end +if isfield(FitParameter, 'peakShirleyK') && ~matches(p.Results.backgroundMode, 'shirley') + warning('WARNING: FitParameter has peakShirleyK but backgroundMode=''%s''. peakShirleyK will be ignored. Use backgroundMode=''shirley'' to enable Shirley background.', p.Results.backgroundMode); +end + +for iaz = 1 + %% FitParameter Structure + fprintf('Validating FitParameter structure...\n'); + % Determine number of peaks + if ~isfield(FitParameter, 'peakArea') + error('FitParameter must contain ''peakArea'' field.'); + end + peakArea = FitParameter.peakArea; + if ~isvector(peakArea) || size(peakArea, 1) ~= 1 + error('peakArea must be a 1×nPeaks row vector, got size [%d, %d].', ... + size(peakArea, 1), size(peakArea, 2)); + end + nPeaks = size(peakArea, 2); + if nPeaks < 1 + error('peakArea must contain at least 1 peak.'); + end + fprintf(' Number of peaks detected: %d\n', nPeaks); -%-------------------------------------------------------------------------- -% 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)) + % Validate core parameters required for ALL modes + fprintf(' Checking core parameters (required for all modes)...\n'); -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 + checkFieldDimension(FitParameter, 'peakAreaLowerBound', nPeaks, p.Results.mode, 'lower area bounds'); + checkFieldDimension(FitParameter, 'peakAreaUpperBound', nPeaks, p.Results.mode, 'upper area bounds'); + checkBoundsPair(FitParameter, 'peakAreaLowerBound', 'peakAreaUpperBound', 'peakArea'); - % 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); + + % Mode-specific validation + fprintf(' Checking mode-specific parameters for mode ''%s''...\n', p.Results.mode); + + switch p.Results.mode + case 'unconstrained' + % All fields are arrays with size 1×nPeaks + checkFieldDimension(FitParameter, 'peakWidth', nPeaks, p.Results.mode, 'peak widths (FWHM)'); + checkFieldDimension(FitParameter, 'peakPos', nPeaks, p.Results.mode, 'peak positions'); + checkFieldDimension(FitParameter, 'peakWidthLowerBound', nPeaks, p.Results.mode, 'lower width bounds'); + checkFieldDimension(FitParameter, 'peakWidthUpperBound', nPeaks, p.Results.mode, 'upper width bounds'); + checkFieldDimension(FitParameter, 'peakPosLowerBound', nPeaks, p.Results.mode, 'lower position bounds'); + checkFieldDimension(FitParameter, 'peakPosUpperBound', nPeaks, p.Results.mode, 'upper position bounds'); + checkBoundsPair(FitParameter, 'peakWidthLowerBound', 'peakWidthUpperBound', 'peakWidth'); + checkBoundsPair(FitParameter, 'peakPosLowerBound', 'peakPosUpperBound', 'peakPos'); + case 'position' + % peakWidth and peakPos are arrays; position bounds are SCALAR + checkFieldDimension(FitParameter, 'peakWidth', nPeaks, p.Results.mode, 'peak widths (independent)'); + checkFieldDimension(FitParameter, 'peakPos', nPeaks, p.Results.mode, 'peak positions (pos[1] abs, others rel)'); + checkFieldDimension(FitParameter, 'peakWidthLowerBound', nPeaks, p.Results.mode, 'lower width bounds'); + checkFieldDimension(FitParameter, 'peakWidthUpperBound', nPeaks, p.Results.mode, 'upper width bounds'); + % Position bounds MUST be SCALAR + checkFieldDimension(FitParameter, 'peakPosLowerBound', 0, p.Results.mode, 'lower position bound (SCALAR for peak 1)'); + checkFieldDimension(FitParameter, 'peakPosUpperBound', 0, p.Results.mode, 'upper position bound (SCALAR for peak 1)'); + checkBoundsPair(FitParameter, 'peakWidthLowerBound', 'peakWidthUpperBound', 'peakWidth'); + checkBoundsPair(FitParameter, 'peakPosLowerBound', 'peakPosUpperBound', 'peakPos'); + case 'widthPosition' + % peakWidth is SCALAR; peakPos is array; all bounds are SCALAR + checkFieldDimension(FitParameter, 'peakWidth', 0, p.Results.mode, 'peak width (SHARED, scalar)'); + checkFieldDimension(FitParameter, 'peakPos', nPeaks, p.Results.mode, 'peak positions (pos[1] abs, others rel)'); + checkFieldDimension(FitParameter, 'peakWidthLowerBound', 0, p.Results.mode, 'lower width bound (SCALAR, shared)'); + checkFieldDimension(FitParameter, 'peakWidthUpperBound', 0, p.Results.mode, 'upper width bound (SCALAR, shared)'); + checkFieldDimension(FitParameter, 'peakPosLowerBound', 0, p.Results.mode, 'lower position bound (SCALAR for peak 1)'); + checkFieldDimension(FitParameter, 'peakPosUpperBound', 0, p.Results.mode, 'upper position bound (SCALAR for peak 1)'); + checkBoundsPair(FitParameter, 'peakWidthLowerBound', 'peakWidthUpperBound', 'peakWidth'); + checkBoundsPair(FitParameter, 'peakPosLowerBound', 'peakPosUpperBound', 'peakPos'); + case 'width' + % peakWidth is SCALAR; peakPos and pos bounds are arrays + checkFieldDimension(FitParameter, 'peakWidth', 0, p.Results.mode, 'peak width (SHARED, scalar)'); + checkFieldDimension(FitParameter, 'peakPos', nPeaks, p.Results.mode, 'peak positions (all absolute, independent)'); + checkFieldDimension(FitParameter, 'peakWidthLowerBound', 0, p.Results.mode, 'lower width bound (SCALAR, shared)'); + checkFieldDimension(FitParameter, 'peakWidthUpperBound', 0, p.Results.mode, 'upper width bound (SCALAR, shared)'); + checkFieldDimension(FitParameter, 'peakPosLowerBound', nPeaks, p.Results.mode, 'lower position bounds (array, per-peak)'); + checkFieldDimension(FitParameter, 'peakPosUpperBound', nPeaks, p.Results.mode, 'upper position bounds (array, per-peak)'); + checkBoundsPair(FitParameter, 'peakWidthLowerBound', 'peakWidthUpperBound', 'peakWidth'); + checkBoundsPair(FitParameter, 'peakPosLowerBound', 'peakPosUpperBound', 'peakPos'); + otherwise + % For custom modes, do basic checks + warning('Custom mode ''%s'' detected. Performing basic checks only.', p.Results.mode); 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 + % Background model validation + fprintf(' Checking background model parameters for ''%s''...\n', p.Results.backgroundMode); - % 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 + if contains(p.Results.mode, 'Shirley') % Conditional K validation: + shirleyKExpectedSize = 0; % scalar for shared modes + else + shirleyKExpectedSize = nPeaks; % array for per-peak modes + end - % 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 + switch p.Results.backgroundMode + case 'shirley' + % Shirley background requires 4 fields + checkFieldDimension(FitParameter, 'peakShirleyK', shirleyKExpectedSize, p.Results.backgroundMode, 'Shirley K values'); + checkFieldDimension(FitParameter, 'peakShirleyKLowerBound', shirleyKExpectedSize, p.Results.backgroundMode, 'lower Shirley K bounds'); + checkFieldDimension(FitParameter, 'peakShirleyKUpperBound', shirleyKExpectedSize, p.Results.backgroundMode, 'upper Shirley K bounds'); + checkFieldDimension(FitParameter, 'peakShirleyOffset', 0, p.Results.backgroundMode, 'Shirley offset (C0, scalar)'); + checkBoundsPair(FitParameter, 'peakShirleyKLowerBound', 'peakShirleyKUpperBound', 'peakShirleyK'); + case 'linear' + % Linear background has no additional parameters + % But warn if Shirley parameters are present (will be ignored) + if isfield(FitParameter, 'peakShirleyK') + warning('backgroundMode is ''linear'' but peakShirleyK is present. It will be ignored.'); + end + if isfield(FitParameter, 'peakShirleyOffset') + warning('backgroundMode is ''linear'' but peakShirleyOffset is present. It will be ignored.'); + end + end -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 + % Pseudo-Voigt extension check + if isfield(FitParameter, 'peakEta') + fprintf(' Checking Pseudo-Voigt extension parameters...\n'); - % 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)]; + % eta can have fewer or equal entries than nPeaks + nEta = numel(FitParameter.peakEta); + if nEta > nPeaks + error('peakEta has more entries (%d) than number of peaks (%d).', nEta, nPeaks); + end - % 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 + checkFieldDimension(FitParameter, 'peakEta', nEta, 'pseudo-Voigt', 'eta parameters'); + checkFieldDimension(FitParameter, 'peakEtaLowerBound', nEta, 'pseudo-Voigt', 'lower eta bounds'); + checkFieldDimension(FitParameter, 'peakEtaUpperBound', nEta, 'pseudo-Voigt', 'upper eta bounds'); + checkBoundsPair(FitParameter, 'peakEtaLowerBound', 'peakEtaUpperBound', 'peakEta'); -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 + % Check eta range + if any(FitParameter.peakEta < 0 | FitParameter.peakEta > 1) + warning('peakEta values should be in range [0, 1] (0=Gaussian, 1=Lorentzian). Found values outside this range.'); + end + end + fprintf('✓ FitParameter validation successful for %d peak(s) in ''%s'' mode with ''%s'' background.\n\n', ... + nPeaks, p.Results.mode, p.Results.backgroundMode); - % 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)]; +end +% ========================================================================= +%% Define Models and Generate Parameter Structures +% ========================================================================= +% +% This section builds the parameter vectors for the optimization: +% +% 1. Determine number of peaks from input FitParameter.peakArea array size +% 2. Check if Pseudo-Voigt eta parameter exists (optional shape control) +% 3. Select appropriate model based on: +% - Background mode: 'linear' or 'shirley' +% - Peak shape: Gaussian (no eta) or Pseudo-Voigt (with eta) +% +% For each mode, three vectors are constructed: +% +% initialParameterVector : Starting guess values for optimization +% lowerBoundsParameterVector : Minimum allowed values (constraints) +% upperBoundsParameterVector : Maximum allowed values (constraints) +% +% CRITICAL: Each parameter vector must: +% • Contain ONLY valid numbers (no NaN or empty fields) +% • Have EXACTLY as many elements as there are fitting parameters +% • Match the specific model's parameter count +% +% Since different models have different fitting parameters, the vectors +% MUST be defined separately for EACH model combination: +% +% Example for UNCONSTRAINED + LINEAR + Gaussian: +% [area₁, width₁, pos₁, area₂, width₂, pos₂, ..., slope, intercept] +% +% Example for UNCONSTRAINED + LINEAR + Pseudo-Voigt: +% [area₁, width₁, pos₁, eta₁, area₂, width₂, pos₂, eta₂, ..., slope, intercept] +% +% Example for UNCONSTRAINED + SHIRLEY + Gaussian: +% [area₁, width₁, pos₁, ..., K₁, K₂, ..., C₀] +% +% The vectors are flattened into 1D arrays and passed to the optimizer +% along with modelDefinitionFunction (how to calculate spectrum from params) +% +% ========================================================================= - % 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 +for iuz = 1 + % initialize Parameter Vectors + initialParameterVector = []; + lowerBoundsParameterVector = []; + upperBoundsParameterVector = []; + %% Generic Models + % + % The modelDefinitionFunction is a function handle that computes the + % theoretical spectrum from a parameter vector and energy values. Rather + % than hard coding specific equations for each model combination, a generic + % model generator is used: + % + % @(currentParameterVector, energyValues) ... + % genericModelGenerator(currentParameterVector, energyValues, FitParameter, p) + % + % This function dynamically constructs the fit based on: + % • Current parameter values (area, width, position, eta, background, etc.) + % • Energy scale points + % • User-selected mode (unconstrained, position, width, etc.) + % • Background type (linear, shirley) + % • Peak shape (gaussian, pseudo-voigt) + % + % The actual mathematical equations (Gaussian, Pseudo-Voigt, Shirley background) + % are generated inside genericModelGenerator based on the input configuration. + % This avoids code duplication and makes adding new models more straightforward. + % + % ========================================================================= -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 + if matches(p.Results.mode, 'unconstrained') + %% unconstrained + nPeaks = size(FitParameter.peakArea, 2); + peakEta_exists = isfield(FitParameter,'peakEta'); + + if matches(p.Results.backgroundMode, 'linear') && ~peakEta_exists + % ===== GAUSSIAN + LINEAR (UNCONSTRAINED) ===== + modelDefinitionFunction = @(currentParameterVector,energyValues) ... + genericModelGenerator(currentParameterVector,energyValues,FitParameter,p); + + for iPeak = 1:nPeaks + initialParameterVector = [initialParameterVector, ... + FitParameter.peakArea(1,iPeak), ... + FitParameter.peakWidth(1,iPeak), ... + FitParameter.peakPos(1,iPeak)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, ... + FitParameter.peakAreaLowerBound(1,iPeak), ... + FitParameter.peakWidthLowerBound(1,iPeak), ... + FitParameter.peakPosLowerBound(1,iPeak)]; + upperBoundsParameterVector = [upperBoundsParameterVector, ... + FitParameter.peakAreaUpperBound(1,iPeak), ... + FitParameter.peakWidthUpperBound(1,iPeak), ... + FitParameter.peakPosUpperBound(1,iPeak)]; + end + % Append linear background (slope, intercept) + initialParameterVector = [initialParameterVector, 0, 0]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, -Inf, -Inf]; + upperBoundsParameterVector = [upperBoundsParameterVector, Inf, Inf]; + + elseif matches(p.Results.backgroundMode, 'linear') && peakEta_exists + % ===== PSEUDO-VOIGT + LINEAR (UNCONSTRAINED) ===== + modelDefinitionFunction = @(currentParameterVector,energyValues) ... + genericModelGenerator(currentParameterVector,energyValues,FitParameter,p); + + initialParameterVector = []; + lowerBoundsParameterVector = []; + upperBoundsParameterVector = []; + for iPeak = 1:nPeaks + initialParameterVector = [initialParameterVector, ... + FitParameter.peakArea(1,iPeak), ... + FitParameter.peakWidth(1,iPeak), ... + FitParameter.peakPos(1,iPeak), ... + FitParameter.peakEta(1,iPeak)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, ... + FitParameter.peakAreaLowerBound(1,iPeak), ... + FitParameter.peakWidthLowerBound(1,iPeak), ... + FitParameter.peakPosLowerBound(1,iPeak), ... + FitParameter.peakEtaLowerBound(1,iPeak)]; + upperBoundsParameterVector = [upperBoundsParameterVector, ... + FitParameter.peakAreaUpperBound(1,iPeak), ... + FitParameter.peakWidthUpperBound(1,iPeak), ... + FitParameter.peakPosUpperBound(1,iPeak), ... + FitParameter.peakEtaUpperBound(1,iPeak)]; + end + % Append linear background (slope, intercept) + initialParameterVector = [initialParameterVector, 0, 0]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, -Inf, -Inf]; + upperBoundsParameterVector = [upperBoundsParameterVector, Inf, Inf]; + + elseif contains(p.Results.backgroundMode, 'shirley') && ~peakEta_exists + % ===== GAUSSIAN + SHIRLEY (UNCONSTRAINED) ===== + modelDefinitionFunction = @(currentParameterVector,energyValues) ... + genericModelGenerator(currentParameterVector,energyValues,FitParameter,p); + + initialParameterVector = []; + lowerBoundsParameterVector = []; + upperBoundsParameterVector = []; + + for iPeak = 1:nPeaks + initialParameterVector = [initialParameterVector, ... + FitParameter.peakArea(1,iPeak), ... + FitParameter.peakWidth(1,iPeak), ... + FitParameter.peakPos(1,iPeak)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, ... + FitParameter.peakAreaLowerBound(1,iPeak), ... + FitParameter.peakWidthLowerBound(1,iPeak), ... + FitParameter.peakPosLowerBound(1,iPeak)]; + upperBoundsParameterVector = [upperBoundsParameterVector, ... + FitParameter.peakAreaUpperBound(1,iPeak), ... + FitParameter.peakWidthUpperBound(1,iPeak), ... + FitParameter.peakPosUpperBound(1,iPeak)]; + end + + % Append kappa + for iPeak = 1:nPeaks + initialParameterVector = [initialParameterVector, FitParameter.peakShirleyK(1,iPeak)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakShirleyKLowerBound(1,iPeak)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakShirleyKUpperBound(1,iPeak)]; + end + + % Append C0 + initialParameterVector = [initialParameterVector, FitParameter.peakShirleyOffset]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, 0]; + upperBoundsParameterVector = [upperBoundsParameterVector, Inf]; + + elseif matches(p.Results.backgroundMode, 'shirley') && peakEta_exists + % ===== PSEUDO-VOIGT + SHIRLEY (UNCONSTRAINED) ===== + modelDefinitionFunction = @(currentParameterVector,energyValues) ... + genericModelGenerator(currentParameterVector,energyValues,FitParameter,p); + + for iPeak = 1:nPeaks + initialParameterVector = [initialParameterVector, ... + FitParameter.peakArea(1,iPeak), ... + FitParameter.peakWidth(1,iPeak), ... + FitParameter.peakPos(1,iPeak), ... + FitParameter.peakEta(1,iPeak)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, ... + FitParameter.peakAreaLowerBound(1,iPeak), ... + FitParameter.peakWidthLowerBound(1,iPeak), ... + FitParameter.peakPosLowerBound(1,iPeak), ... + FitParameter.peakEtaLowerBound(1,iPeak) ]; + upperBoundsParameterVector = [upperBoundsParameterVector, ... + FitParameter.peakAreaUpperBound(1,iPeak), ... + FitParameter.peakWidthUpperBound(1,iPeak), ... + FitParameter.peakPosUpperBound(1,iPeak), ... + FitParameter.peakEtaUpperBound(1,iPeak) ]; + end + + % Append kappa + for iPeak = 1:nPeaks + initialParameterVector = [initialParameterVector, FitParameter.peakShirleyK(1,iPeak)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakShirleyKLowerBound(1,iPeak)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakShirleyKUpperBound(1,iPeak)]; + end + + % Append C0 + initialParameterVector = [initialParameterVector, FitParameter.peakShirleyOffset]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, 0]; + upperBoundsParameterVector = [upperBoundsParameterVector, Inf]; + + else + error('Unconstrained mode: invalid combination of backgroundMode and peakEta.'); + end + + elseif matches(p.Results.mode, 'position') + %% position constrained + + nPeaks = size(FitParameter.peakArea, 2); + peakEtaexists = isfield(FitParameter, 'peakEta'); + + if matches(p.Results.backgroundMode, 'linear') + if ~peakEtaexists % GAUSSIAN LINEAR + modelDefinitionFunction = @(currentParameterVector,energyValues) ... + genericModelGenerator(currentParameterVector,energyValues,FitParameter,p); + + for iPeak = 1:nPeaks + initialParameterVector = [initialParameterVector, ... + FitParameter.peakArea(1, iPeak), ... + FitParameter.peakWidth(1, iPeak)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, ... + FitParameter.peakAreaLowerBound(1, iPeak), ... + FitParameter.peakWidthLowerBound(1, iPeak)]; + upperBoundsParameterVector = [upperBoundsParameterVector, ... + FitParameter.peakAreaUpperBound(1, iPeak), ... + FitParameter.peakWidthUpperBound(1, iPeak)]; + end + + % Peak 1 absolute position + initialParameterVector = [initialParameterVector, FitParameter.peakPos(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakPosLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakPosUpperBound(1, 1)]; + + % Linear background: slope, intercept + initialParameterVector = [initialParameterVector, 0, 0]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, -Inf, -Inf]; + upperBoundsParameterVector = [upperBoundsParameterVector, Inf, Inf]; + + elseif peakEtaexists + modelDefinitionFunction = @(currentParameterVector,energyValues) ... + genericModelGenerator(currentParameterVector,energyValues,FitParameter,p); + + for iPeak = 1:nPeaks + initialParameterVector = [initialParameterVector, ... + FitParameter.peakArea(1, iPeak), ... + FitParameter.peakWidth(1, iPeak), ... + FitParameter.peakEta(1, iPeak)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, ... + FitParameter.peakAreaLowerBound(1, iPeak), ... + FitParameter.peakWidthLowerBound(1, iPeak), ... + FitParameter.peakEtaLowerBound(1, iPeak)]; + upperBoundsParameterVector = [upperBoundsParameterVector, ... + FitParameter.peakAreaUpperBound(1, iPeak), ... + FitParameter.peakWidthUpperBound(1, iPeak), ... + FitParameter.peakEtaUpperBound(1, iPeak)]; + end + + % Peak 1 absolute position + initialParameterVector = [initialParameterVector, FitParameter.peakPos(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakPosLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakPosUpperBound(1, 1)]; + + % Linear background + initialParameterVector = [initialParameterVector, 0, 0]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, -Inf, -Inf]; + upperBoundsParameterVector = [upperBoundsParameterVector, Inf, Inf]; + end + + elseif matches(p.Results.backgroundMode, 'shirley') + if ~peakEtaexists + modelDefinitionFunction = @(currentParameterVector,energyValues) ... + genericModelGenerator(currentParameterVector,energyValues,FitParameter,p); + + for iPeak = 1:nPeaks + initialParameterVector = [initialParameterVector, ... + FitParameter.peakArea(1, iPeak), ... + FitParameter.peakWidth(1, iPeak)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, ... + FitParameter.peakAreaLowerBound(1, iPeak), ... + FitParameter.peakWidthLowerBound(1, iPeak)]; + upperBoundsParameterVector = [upperBoundsParameterVector, ... + FitParameter.peakAreaUpperBound(1, iPeak), ... + FitParameter.peakWidthUpperBound(1, iPeak)]; + end + + % Peak 1 absolute position + initialParameterVector = [initialParameterVector, FitParameter.peakPos(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakPosLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakPosUpperBound(1, 1)]; + + % Append kappa + for iPeak = 1:nPeaks + initialParameterVector = [initialParameterVector, FitParameter.peakShirleyK(1,iPeak)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakShirleyKLowerBound(1,iPeak)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakShirleyKUpperBound(1,iPeak)]; + end + + % Append C0 + initialParameterVector = [initialParameterVector, FitParameter.peakShirleyOffset]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, 0]; + upperBoundsParameterVector = [upperBoundsParameterVector, Inf]; + + elseif peakEtaexists + modelDefinitionFunction = @(currentParameterVector,energyValues) ... + genericModelGenerator(currentParameterVector,energyValues,FitParameter,p); + + for iPeak = 1:nPeaks + initialParameterVector = [initialParameterVector, ... + FitParameter.peakArea(1, iPeak), ... + FitParameter.peakWidth(1, iPeak), ... + FitParameter.peakEta(1, iPeak)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, ... + FitParameter.peakAreaLowerBound(1, iPeak), ... + FitParameter.peakWidthLowerBound(1, iPeak), ... + FitParameter.peakEtaLowerBound(1, iPeak)]; + upperBoundsParameterVector = [upperBoundsParameterVector, ... + FitParameter.peakAreaUpperBound(1, iPeak), ... + FitParameter.peakWidthUpperBound(1, iPeak), ... + FitParameter.peakEtaUpperBound(1, iPeak)]; + end + + initialParameterVector = [initialParameterVector, FitParameter.peakPos(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakPosLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakPosUpperBound(1, 1)]; + + % Append kappa + for iPeak = 1:nPeaks + initialParameterVector = [initialParameterVector, FitParameter.peakShirleyK(1,iPeak)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakShirleyKLowerBound(1,iPeak)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakShirleyKUpperBound(1,iPeak)]; + end + + % Append C0 + initialParameterVector = [initialParameterVector, FitParameter.peakShirleyOffset]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, 0]; + upperBoundsParameterVector = [upperBoundsParameterVector, Inf]; + + else + error('Position mode: invalid combination of backgroundMode and peakEta.'); + end + else + error('Position mode: unsupported backgroundMode ''%s''.', p.Results.backgroundMode); + end + + elseif matches(p.Results.mode, 'positionShirley') + %% position & Shirley constrained + + nPeaks = size(FitParameter.peakArea, 2); + peakEtaexists = isfield(FitParameter, 'peakEta'); + + if matches(p.Results.backgroundMode, 'shirley') + if ~peakEtaexists + modelDefinitionFunction = @(currentParameterVector,energyValues) ... + genericModelGenerator(currentParameterVector,energyValues,FitParameter,p); + + for iPeak = 1:nPeaks + initialParameterVector = [initialParameterVector, ... + FitParameter.peakArea(1, iPeak), ... + FitParameter.peakWidth(1, iPeak)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, ... + FitParameter.peakAreaLowerBound(1, iPeak), ... + FitParameter.peakWidthLowerBound(1, iPeak)]; + upperBoundsParameterVector = [upperBoundsParameterVector, ... + FitParameter.peakAreaUpperBound(1, iPeak), ... + FitParameter.peakWidthUpperBound(1, iPeak)]; + end + + % Peak 1 absolute position + initialParameterVector = [initialParameterVector, FitParameter.peakPos(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakPosLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakPosUpperBound(1, 1)]; + + % Append kappa + initialParameterVector = [initialParameterVector, FitParameter.peakShirleyK(1,1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakShirleyKLowerBound(1,1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakShirleyKUpperBound(1,1)]; + + % Append C0 + initialParameterVector = [initialParameterVector, FitParameter.peakShirleyOffset]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, 0]; + upperBoundsParameterVector = [upperBoundsParameterVector, Inf]; + + elseif peakEtaexists + modelDefinitionFunction = @(currentParameterVector,energyValues) ... + genericModelGenerator(currentParameterVector,energyValues,FitParameter,p); + + for iPeak = 1:nPeaks + initialParameterVector = [initialParameterVector, ... + FitParameter.peakArea(1, iPeak), ... + FitParameter.peakWidth(1, iPeak), ... + FitParameter.peakEta(1, iPeak)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, ... + FitParameter.peakAreaLowerBound(1, iPeak), ... + FitParameter.peakWidthLowerBound(1, iPeak), ... + FitParameter.peakEtaLowerBound(1, iPeak)]; + upperBoundsParameterVector = [upperBoundsParameterVector, ... + FitParameter.peakAreaUpperBound(1, iPeak), ... + FitParameter.peakWidthUpperBound(1, iPeak), ... + FitParameter.peakEtaUpperBound(1, iPeak)]; + end + + % Peak 1 absolute position + initialParameterVector = [initialParameterVector, FitParameter.peakPos(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakPosLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakPosUpperBound(1, 1)]; + + % Append kappa + initialParameterVector = [initialParameterVector, FitParameter.peakShirleyK(1,1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakShirleyKLowerBound(1,1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakShirleyKUpperBound(1,1)]; + + % Append C0 + initialParameterVector = [initialParameterVector, FitParameter.peakShirleyOffset]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, 0]; + upperBoundsParameterVector = [upperBoundsParameterVector, Inf]; + end + end + + elseif matches(p.Results.mode,'widthPosition') + %% position & width constrained + + nPeaks = size(FitParameter.peakArea, 2); + peakEtaexists = isfield(FitParameter, 'peakEta'); + + if matches(p.Results.backgroundMode, 'linear') + + if ~peakEtaexists + modelDefinitionFunction = @(currentParameterVector,energyValues) ... + genericModelGenerator(currentParameterVector,energyValues,FitParameter,p); + + for iPeak = 1:nPeaks + initialParameterVector = [initialParameterVector, ... + FitParameter.peakArea(1, iPeak)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, ... + FitParameter.peakAreaLowerBound(1, iPeak)]; + upperBoundsParameterVector = [upperBoundsParameterVector, ... + FitParameter.peakAreaUpperBound(1, iPeak)]; + end + + % Shared width + initialParameterVector = [initialParameterVector, FitParameter.peakWidth(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakWidthLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakWidthUpperBound(1, 1)]; + + % Peak 1 absolute position + initialParameterVector = [initialParameterVector, FitParameter.peakPos(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakPosLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakPosUpperBound(1, 1)]; + + % Linear background: slope, intercept + initialParameterVector = [initialParameterVector, 0, 0]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, -Inf, -Inf]; + upperBoundsParameterVector = [upperBoundsParameterVector, Inf, Inf]; + + elseif peakEtaexists + modelDefinitionFunction = @(currentParameterVector,energyValues) ... + genericModelGenerator(currentParameterVector,energyValues,FitParameter,p); + + for iPeak = 1:nPeaks + initialParameterVector = [initialParameterVector, ... + FitParameter.peakArea(1, iPeak), ... + FitParameter.peakEta(1, iPeak)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, ... + FitParameter.peakAreaLowerBound(1, iPeak), ... + FitParameter.peakEtaLowerBound(1, iPeak)]; + upperBoundsParameterVector = [upperBoundsParameterVector, ... + FitParameter.peakAreaUpperBound(1, iPeak), ... + FitParameter.peakEtaUpperBound(1, iPeak)]; + end + + % Shared width + initialParameterVector = [initialParameterVector, FitParameter.peakWidth(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakWidthLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakWidthUpperBound(1, 1)]; + + % Peak 1 absolute position + initialParameterVector = [initialParameterVector, FitParameter.peakPos(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakPosLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakPosUpperBound(1, 1)]; + + % Linear background: slope, intercept + initialParameterVector = [initialParameterVector, 0, 0]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, -Inf, -Inf]; + upperBoundsParameterVector = [upperBoundsParameterVector, Inf, Inf]; + end + + elseif matches(p.Results.backgroundMode, 'shirley') + + if ~peakEtaexists + modelDefinitionFunction = @(currentParameterVector,energyValues) ... + genericModelGenerator(currentParameterVector,energyValues,FitParameter,p); + + for iPeak = 1:nPeaks + initialParameterVector = [initialParameterVector, ... + FitParameter.peakArea(1, iPeak)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, ... + FitParameter.peakAreaLowerBound(1, iPeak)]; + upperBoundsParameterVector = [upperBoundsParameterVector, ... + FitParameter.peakAreaUpperBound(1, iPeak)]; + end + + % Shared width + initialParameterVector = [initialParameterVector, FitParameter.peakWidth(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakWidthLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakWidthUpperBound(1, 1)]; + + % Peak 1 absolute position (pos1) + initialParameterVector = [initialParameterVector, FitParameter.peakPos(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakPosLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakPosUpperBound(1, 1)]; + + % Grouped kappas: K1, K2, ..., KN + for iPeak = 1:nPeaks + initialParameterVector = [initialParameterVector, FitParameter.peakShirleyK(1, iPeak)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakShirleyKLowerBound(1, iPeak)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakShirleyKUpperBound(1, iPeak)]; + end + + % C0 (Shirley offset) + initialParameterVector = [initialParameterVector, FitParameter.peakShirleyOffset]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, 0]; + upperBoundsParameterVector = [upperBoundsParameterVector, Inf]; + + elseif peakEtaexists + + modelDefinitionFunction = @(currentParameterVector,energyValues) ... + genericModelGenerator(currentParameterVector,energyValues,FitParameter,p); + + for iPeak = 1:nPeaks + initialParameterVector = [initialParameterVector, ... + FitParameter.peakArea(1, iPeak), ... + FitParameter.peakEta(1, iPeak)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, ... + FitParameter.peakAreaLowerBound(1, iPeak), ... + FitParameter.peakEtaLowerBound(1, iPeak)]; + upperBoundsParameterVector = [upperBoundsParameterVector, ... + FitParameter.peakAreaUpperBound(1, iPeak), ... + FitParameter.peakEtaUpperBound(1, iPeak)]; + end + + % Shared width + initialParameterVector = [initialParameterVector, FitParameter.peakWidth(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakWidthLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakWidthUpperBound(1, 1)]; + + % Peak 1 absolute position (pos1) + initialParameterVector = [initialParameterVector, FitParameter.peakPos(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakPosLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakPosUpperBound(1, 1)]; + + % Grouped kappas: K1, K2, ..., KN + for iPeak = 1:nPeaks + initialParameterVector = [initialParameterVector, FitParameter.peakShirleyK(1, iPeak)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakShirleyKLowerBound(1, iPeak)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakShirleyKUpperBound(1, iPeak)]; + end + + % C0 (Shirley offset) + initialParameterVector = [initialParameterVector, FitParameter.peakShirleyOffset]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, 0]; + upperBoundsParameterVector = [upperBoundsParameterVector, Inf]; + + else + error('widthPosition mode: invalid combination of backgroundMode and peakEta.'); + end + + else + error('widthPosition mode: unsupported backgroundMode ''%s''.', p.Results.backgroundMode); + end + + elseif matches(p.Results.mode,'widthPositionShirley') + %% position & width & Shirley constrained + + nPeaks = size(FitParameter.peakArea, 2); + peakEtaexists = isfield(FitParameter, 'peakEta'); + + if matches(p.Results.backgroundMode, 'shirley') + + if ~peakEtaexists + modelDefinitionFunction = @(currentParameterVector,energyValues) ... + genericModelGenerator(currentParameterVector,energyValues,FitParameter,p); + + for iPeak = 1:nPeaks + initialParameterVector = [initialParameterVector, ... + FitParameter.peakArea(1, iPeak)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, ... + FitParameter.peakAreaLowerBound(1, iPeak)]; + upperBoundsParameterVector = [upperBoundsParameterVector, ... + FitParameter.peakAreaUpperBound(1, iPeak)]; + end + + % Shared width + initialParameterVector = [initialParameterVector, FitParameter.peakWidth(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakWidthLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakWidthUpperBound(1, 1)]; + + % Peak 1 absolute position (pos1) + initialParameterVector = [initialParameterVector, FitParameter.peakPos(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakPosLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakPosUpperBound(1, 1)]; + + % Kappa + initialParameterVector = [initialParameterVector, FitParameter.peakShirleyK(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakShirleyKLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakShirleyKUpperBound(1, 1)]; + + % C0 (Shirley offset) + initialParameterVector = [initialParameterVector, FitParameter.peakShirleyOffset]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, 0]; + upperBoundsParameterVector = [upperBoundsParameterVector, Inf]; + + elseif peakEtaexists + + modelDefinitionFunction = @(currentParameterVector,energyValues) ... + genericModelGenerator(currentParameterVector,energyValues,FitParameter,p); + + for iPeak = 1:nPeaks + initialParameterVector = [initialParameterVector, ... + FitParameter.peakArea(1, iPeak), ... + FitParameter.peakEta(1, iPeak)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, ... + FitParameter.peakAreaLowerBound(1, iPeak), ... + FitParameter.peakEtaLowerBound(1, iPeak)]; + upperBoundsParameterVector = [upperBoundsParameterVector, ... + FitParameter.peakAreaUpperBound(1, iPeak), ... + FitParameter.peakEtaUpperBound(1, iPeak)]; + end + + % Shared width + initialParameterVector = [initialParameterVector, FitParameter.peakWidth(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakWidthLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakWidthUpperBound(1, 1)]; + + % Peak 1 absolute position (pos1) + initialParameterVector = [initialParameterVector, FitParameter.peakPos(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakPosLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakPosUpperBound(1, 1)]; + + % Kappa + initialParameterVector = [initialParameterVector, FitParameter.peakShirleyK(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakShirleyKLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakShirleyKUpperBound(1, 1)]; + + % C0 (Shirley offset) + initialParameterVector = [initialParameterVector, FitParameter.peakShirleyOffset]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, 0]; + upperBoundsParameterVector = [upperBoundsParameterVector, Inf]; + + end + end + + elseif matches(p.Results.mode,'width') + %% width constrained + + nPeaks = size(FitParameter.peakArea, 2); + peakEtaexists = isfield(FitParameter, 'peakEta'); + + if matches(p.Results.backgroundMode, 'linear') + + if ~peakEtaexists + modelDefinitionFunction = @(currentParameterVector,energyValues) ... + genericModelGenerator(currentParameterVector,energyValues,FitParameter,p); + + for iPeak = 1:nPeaks + initialParameterVector = [initialParameterVector, ... + FitParameter.peakArea(1, iPeak), ... + FitParameter.peakPos(1, iPeak)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, ... + FitParameter.peakAreaLowerBound(1, iPeak), ... + FitParameter.peakPosLowerBound(1, iPeak)]; + upperBoundsParameterVector = [upperBoundsParameterVector, ... + FitParameter.peakAreaUpperBound(1, iPeak), ... + FitParameter.peakPosUpperBound(1, iPeak)]; + end + + % Shared FWHM + initialParameterVector = [initialParameterVector, FitParameter.peakWidth(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakWidthLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakWidthUpperBound(1, 1)]; + + % Linear background: slope, intercept + initialParameterVector = [initialParameterVector, 0, 0]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, -Inf, -Inf]; + upperBoundsParameterVector = [upperBoundsParameterVector, Inf, Inf]; + + elseif peakEtaexists + modelDefinitionFunction = @(currentParameterVector,energyValues) ... + genericModelGenerator(currentParameterVector,energyValues,FitParameter,p); + + for iPeak = 1:nPeaks + initialParameterVector = [initialParameterVector, ... + FitParameter.peakArea(1, iPeak), ... + FitParameter.peakPos(1, iPeak), ... + FitParameter.peakEta(1, iPeak)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, ... + FitParameter.peakAreaLowerBound(1, iPeak), ... + FitParameter.peakPosLowerBound(1, iPeak), ... + FitParameter.peakEtaLowerBound(1, iPeak)]; + upperBoundsParameterVector = [upperBoundsParameterVector, ... + FitParameter.peakAreaUpperBound(1, iPeak), ... + FitParameter.peakPosUpperBound(1, iPeak), ... + FitParameter.peakEtaUpperBound(1, iPeak)]; + end + + % Shared FWHM + initialParameterVector = [initialParameterVector, FitParameter.peakWidth(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakWidthLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakWidthUpperBound(1, 1)]; + + % Linear background: slope, intercept + initialParameterVector = [initialParameterVector, 0, 0]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, -Inf, -Inf]; + upperBoundsParameterVector = [upperBoundsParameterVector, Inf, Inf]; + end + + elseif matches(p.Results.backgroundMode, 'shirley') + + if ~peakEtaexists + + modelDefinitionFunction = @(currentParameterVector,energyValues) ... + genericModelGenerator(currentParameterVector,energyValues,FitParameter,p); + + for iPeak = 1:nPeaks + initialParameterVector = [initialParameterVector, ... + FitParameter.peakArea(1, iPeak), ... + FitParameter.peakPos(1, iPeak)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, ... + FitParameter.peakAreaLowerBound(1, iPeak), ... + FitParameter.peakPosLowerBound(1, iPeak)]; + upperBoundsParameterVector = [upperBoundsParameterVector, ... + FitParameter.peakAreaUpperBound(1, iPeak), ... + FitParameter.peakPosUpperBound(1, iPeak)]; + end + + % Shared FWHM + initialParameterVector = [initialParameterVector, FitParameter.peakWidth(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakWidthLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakWidthUpperBound(1, 1)]; + + % Grouped kappas: K1..KN + for iPeak = 1:nPeaks + initialParameterVector = [initialParameterVector, FitParameter.peakShirleyK(1, iPeak)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakShirleyKLowerBound(1, iPeak)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakShirleyKUpperBound(1, iPeak)]; + end + + % C0 (Shirley offset) + initialParameterVector = [initialParameterVector, FitParameter.peakShirleyOffset]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, 0]; + upperBoundsParameterVector = [upperBoundsParameterVector, Inf]; + + elseif peakEtaexists + + modelDefinitionFunction = @(currentParameterVector,energyValues) ... + genericModelGenerator(currentParameterVector,energyValues,FitParameter,p); + + for iPeak = 1:nPeaks + initialParameterVector = [initialParameterVector, ... + FitParameter.peakArea(1, iPeak), ... + FitParameter.peakPos(1, iPeak), ... + FitParameter.peakEta(1, iPeak)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, ... + FitParameter.peakAreaLowerBound(1, iPeak), ... + FitParameter.peakPosLowerBound(1, iPeak), ... + FitParameter.peakEtaLowerBound(1, iPeak)]; + upperBoundsParameterVector = [upperBoundsParameterVector, ... + FitParameter.peakAreaUpperBound(1, iPeak), ... + FitParameter.peakPosUpperBound(1, iPeak), ... + FitParameter.peakEtaUpperBound(1, iPeak)]; + end + + % Shared FWHM + initialParameterVector = [initialParameterVector, FitParameter.peakWidth(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakWidthLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakWidthUpperBound(1, 1)]; + + % Grouped kappas: K1..KN + for iPeak = 1:nPeaks + initialParameterVector = [initialParameterVector, FitParameter.peakShirleyK(1, iPeak)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakShirleyKLowerBound(1, iPeak)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakShirleyKUpperBound(1, iPeak)]; + end + + % C0 (Shirley offset) + initialParameterVector = [initialParameterVector, FitParameter.peakShirleyOffset]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, 0]; + upperBoundsParameterVector = [upperBoundsParameterVector, Inf]; + + else + error('width mode: invalid combination of backgroundMode and peakEta.'); + end + + else + error('width mode: unsupported backgroundMode ''%s''.', p.Results.backgroundMode); + end + + elseif matches(p.Results.mode,'widthShirley') + %% width & Shirley constrained + + nPeaks = size(FitParameter.peakArea, 2); + peakEtaexists = isfield(FitParameter, 'peakEta'); + if matches(p.Results.backgroundMode, 'shirley') + + if ~peakEtaexists + + modelDefinitionFunction = @(currentParameterVector,energyValues) ... + genericModelGenerator(currentParameterVector,energyValues,FitParameter,p); + + for iPeak = 1:nPeaks + initialParameterVector = [initialParameterVector, ... + FitParameter.peakArea(1, iPeak), ... + FitParameter.peakPos(1, iPeak)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, ... + FitParameter.peakAreaLowerBound(1, iPeak), ... + FitParameter.peakPosLowerBound(1, iPeak)]; + upperBoundsParameterVector = [upperBoundsParameterVector, ... + FitParameter.peakAreaUpperBound(1, iPeak), ... + FitParameter.peakPosUpperBound(1, iPeak)]; + end + + % Shared FWHM + initialParameterVector = [initialParameterVector, FitParameter.peakWidth(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakWidthLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakWidthUpperBound(1, 1)]; + + % Kappa + initialParameterVector = [initialParameterVector, FitParameter.peakShirleyK(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakShirleyKLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakShirleyKUpperBound(1, 1)]; + + % C0 (Shirley offset) + initialParameterVector = [initialParameterVector, FitParameter.peakShirleyOffset]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, 0]; + upperBoundsParameterVector = [upperBoundsParameterVector, Inf]; + + elseif peakEtaexists + + modelDefinitionFunction = @(currentParameterVector,energyValues) ... + genericModelGenerator(currentParameterVector,energyValues,FitParameter,p); + + for iPeak = 1:nPeaks + initialParameterVector = [initialParameterVector, ... + FitParameter.peakArea(1, iPeak), ... + FitParameter.peakPos(1, iPeak), ... + FitParameter.peakEta(1, iPeak)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, ... + FitParameter.peakAreaLowerBound(1, iPeak), ... + FitParameter.peakPosLowerBound(1, iPeak), ... + FitParameter.peakEtaLowerBound(1, iPeak)]; + upperBoundsParameterVector = [upperBoundsParameterVector, ... + FitParameter.peakAreaUpperBound(1, iPeak), ... + FitParameter.peakPosUpperBound(1, iPeak), ... + FitParameter.peakEtaUpperBound(1, iPeak)]; + end + + % Shared FWHM + initialParameterVector = [initialParameterVector, FitParameter.peakWidth(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakWidthLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakWidthUpperBound(1, 1)]; + + % Kappa + initialParameterVector = [initialParameterVector, FitParameter.peakShirleyK(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakShirleyKLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakShirleyKUpperBound(1, 1)]; + + % C0 (Shirley offset) + initialParameterVector = [initialParameterVector, FitParameter.peakShirleyOffset]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, 0]; + upperBoundsParameterVector = [upperBoundsParameterVector, Inf]; + + else + error('width mode: invalid combination of backgroundMode and peakEta.'); + end + + else + error('width mode: unsupported backgroundMode ''%s''.', p.Results.backgroundMode); + end + end + + +%% Custom Models +% +% Beyond generic peak fitting, custom models are available for specific +% chemical systems where peaks follow known relationships or constraints. +% Custom models are defined and hard-coded in the mode selection and +% generate their own parameter vectors. Simply add one if needed with a +% specific name. Add the name to the mode options. Also add a section to +% the RESULTS EXTRACTION part and the shirley model function. +% ========================================================================= + + + + if matches(p.Results.mode,'unconstrainedAndWidthPosition') && ... + matches(p.Results.backgroundMode,'linear') + %% unconstrainedAndWidthPosition + modelDefinitionFunction = @(currentParameterVector, energyValues) ... + gaussPeak(energyValues,currentParameterVector(1), currentParameterVector(2), currentParameterVector(3)+currentParameterVector(14)) + ... % Peak 1: a1, w1, (pos2+shift) + gaussPeak(energyValues,currentParameterVector(4), currentParameterVector(2), currentParameterVector(3)) + ... % Peak 2: a2, w1, pos2 + gaussPeak(energyValues,currentParameterVector(5), currentParameterVector(6), currentParameterVector(7)) + ... % Peak 3: a3, w3, pos3 + gaussPeak(energyValues,currentParameterVector(8), currentParameterVector(2), currentParameterVector(9)) + ... % Peak 4: a4, w1, pos4 + gaussPeak(energyValues,currentParameterVector(10), currentParameterVector(2), currentParameterVector(11)) + ... % Peak 5: a5, w1, pos5 + linearBackground(energyValues,currentParameterVector(12), currentParameterVector(13)); % Background: slope, intercept + + initialParameterVector = []; + lowerBoundsParameterVector = []; + upperBoundsParameterVector = []; + + % #1: Area 1 + initialParameterVector = [initialParameterVector, FitParameter.peakArea(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakAreaLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakAreaUpperBound(1, 1)]; + + % #2: Width 1 (shared) + initialParameterVector = [initialParameterVector, FitParameter.peakWidth(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakWidthLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakWidthUpperBound(1, 1)]; + + % #3: Pos 2 (reference) + initialParameterVector = [initialParameterVector, FitParameter.peakPos(1, 2)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakPosLowerBound(1, 2)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakPosUpperBound(1, 2)]; + + % #4: Area 2 + initialParameterVector = [initialParameterVector, FitParameter.peakArea(1, 2)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakAreaLowerBound(1, 2)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakAreaUpperBound(1, 2)]; + + % #5,6,7: Area 3, Width 3, Pos 3 + initialParameterVector = [initialParameterVector, FitParameter.peakArea(1, 3), FitParameter.peakWidth(1, 2), FitParameter.peakPos(1, 3)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakAreaLowerBound(1, 3), FitParameter.peakWidthLowerBound(1, 2), FitParameter.peakPosLowerBound(1, 3)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakAreaUpperBound(1, 3), FitParameter.peakWidthUpperBound(1, 2), FitParameter.peakPosUpperBound(1, 3)]; + + % #8,9 Area 4, Pos 4 + initialParameterVector = [initialParameterVector, FitParameter.peakArea(1, 4), FitParameter.peakPos(1, 4)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakAreaLowerBound(1, 4), FitParameter.peakPosLowerBound(1, 4)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakAreaUpperBound(1, 4) FitParameter.peakPosUpperBound(1, 4)]; + + % #11,12 Area 5, Pos 5 + initialParameterVector = [initialParameterVector, FitParameter.peakArea(1, 5), FitParameter.peakPos(1, 5)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakAreaLowerBound(1, 5), FitParameter.peakPosLowerBound(1, 5)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakAreaUpperBound(1, 5), FitParameter.peakPosUpperBound(1, 5)]; + + % #13,14 Background: slope, intercept + initialParameterVector = [initialParameterVector, 0, 0]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, -Inf, -Inf]; + upperBoundsParameterVector = [upperBoundsParameterVector, Inf, Inf]; + + + % #14: Pos 1 (shift) + initialParameterVector = [initialParameterVector, FitParameter.peakPos(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakPosLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakPosUpperBound(1, 1)]; - % 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,'TwoSpinOrbitDoublet') - % 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)).^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(4)/(currentParameterVector(2)*sqrt(pi/2))*exp(-2.*((energyValues(:,1)- (currentParameterVector(3)+currentParameterVector(5))).^2/currentParameterVector(2).^2))... % Gauss peak 3 - + currentParameterVector(4)*FitParameter.peakArea(1,2)/(currentParameterVector(2)*sqrt(pi/2))*exp(-2.*((energyValues(:,1)- (currentParameterVector(3)+currentParameterVector(5)+FitParameter.peakPos(1,2))).^2/currentParameterVector(2).^2))... % Gauss peak 4 - +currentParameterVector(6).*energyValues(:,1)+currentParameterVector(7); % linear background + %% TwoSpinOrbitDoublet - % initial value - initialParameterVector= [FitParameter.peakArea(1,1) FitParameter.peakWidth(1,1) FitParameter.peakPos(1,1)... % Gauss peak 1 - FitParameter.peakArea(1,3) FitParameter.peakPos(1,3)... - 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,3) FitParameter.peakPosLowerBound(1,3)... - -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,3) FitParameter.peakPosUpperBound(1,3)... - Inf Inf]; % linear background + modelDefinitionFunction = @(currentParameterVector,energyValues)... + gaussPeak(energyValues, currentParameterVector(1),currentParameterVector(2), currentParameterVector(3)) ... % Gauss peak 1 + + gaussPeak(energyValues, currentParameterVector(1)*FitParameter.peakArea(1,2), currentParameterVector(2), currentParameterVector(3)+FitParameter.peakPos(1,2))... % Gauss peak 2 + + gaussPeak(energyValues, currentParameterVector(4),currentParameterVector(2), currentParameterVector(3)+currentParameterVector(5))... % Gauss peak 3 + + gaussPeak(energyValues, currentParameterVector(4)*FitParameter.peakArea(1,2),currentParameterVector(2),currentParameterVector(3)+currentParameterVector(5)+FitParameter.peakPos(1,2)) ... % Gauss peak 4 + + linearBackground(energyValues,currentParameterVector(6), currentParameterVector(7)); + + initialParameterVector = []; + lowerBoundsParameterVector = []; + upperBoundsParameterVector = []; + + % #1: Area 1 + initialParameterVector = [initialParameterVector, FitParameter.peakArea(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakAreaLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakAreaUpperBound(1, 1)]; + + % #2: Width 1 (shared) + initialParameterVector = [initialParameterVector, FitParameter.peakWidth(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakWidthLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakWidthUpperBound(1, 1)]; + + % #3: Pos 1 + initialParameterVector = [initialParameterVector, FitParameter.peakPos(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakPosLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakPosUpperBound(1, 1)]; + + % #4: Area 3 + initialParameterVector = [initialParameterVector, FitParameter.peakArea(1, 3)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakAreaLowerBound(1, 3)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakAreaUpperBound(1, 3)]; + + % #5 Pos 3 + initialParameterVector = [initialParameterVector, FitParameter.peakPos(1, 3)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakPosLowerBound(1, 3)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakPosUpperBound(1, 3)]; + + % #6 & #7 Background: slope, intercept + initialParameterVector = [initialParameterVector, 0, 0]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, -Inf, -Inf]; + upperBoundsParameterVector = [upperBoundsParameterVector, Inf, Inf]; + + elseif matches(p.Results.mode,'TwoSpinOrbitDoubletPseudoVoigtShirley') + %% TwoSpinOrbitDoubletPseudoVoigtShirley + + modelDefinitionFunction = @(currentParameterVector,energyValues)... + pseudoVoigtPeak(energyValues, currentParameterVector(1),currentParameterVector(2), currentParameterVector(3),currentParameterVector(6)) ... % Gauss peak 1 + + pseudoVoigtPeak(energyValues, currentParameterVector(1)*FitParameter.peakArea(1,2), currentParameterVector(2), currentParameterVector(3)+FitParameter.peakPos(1,2),currentParameterVector(6))... % Gauss peak 2 + + pseudoVoigtPeak(energyValues, currentParameterVector(4),currentParameterVector(2), currentParameterVector(3)+currentParameterVector(5),currentParameterVector(6))... % Gauss peak 3 + + pseudoVoigtPeak(energyValues, currentParameterVector(4)*FitParameter.peakArea(1,2),currentParameterVector(2),currentParameterVector(3)+currentParameterVector(5)+FitParameter.peakPos(1,2),currentParameterVector(6)) ... % Gauss peak 4 + + activeShirleyBackgroundFwhm(energyValues,currentParameterVector,FitParameter,p,nPeaks)+currentParameterVector(end); + + initialParameterVector = []; + lowerBoundsParameterVector = []; + upperBoundsParameterVector = []; + + % #1: Area 1 + initialParameterVector = [initialParameterVector, FitParameter.peakArea(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakAreaLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakAreaUpperBound(1, 1)]; + + % #2: Width 1 (shared) + initialParameterVector = [initialParameterVector, FitParameter.peakWidth(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakWidthLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakWidthUpperBound(1, 1)]; + + % #3: Pos 1 + initialParameterVector = [initialParameterVector, FitParameter.peakPos(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakPosLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakPosUpperBound(1, 1)]; + + % #4: Area 3 + initialParameterVector = [initialParameterVector, FitParameter.peakArea(1, 3)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakAreaLowerBound(1, 3)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakAreaUpperBound(1, 3)]; + + % #5 Pos 3, rel to Pos 1 + initialParameterVector = [initialParameterVector, FitParameter.peakPos(1, 3)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakPosLowerBound(1, 3)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakPosUpperBound(1, 3)]; + + % #6 eta + initialParameterVector = [initialParameterVector, FitParameter.peakEta(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakEtaLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector,FitParameter.peakEtaLowerBound(1, 1)]; + + % #7 Kappa + initialParameterVector = [initialParameterVector, FitParameter.peakShirleyK(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakShirleyKLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector,FitParameter.peakShirleyKUpperBound(1, 1)]; + + % #7 C0 + initialParameterVector = [initialParameterVector, FitParameter.peakShirleyOffset(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, -Inf]; + upperBoundsParameterVector = [upperBoundsParameterVector,Inf]; + + elseif matches(p.Results.mode,'gasAndTwoAcidBasePairs') + %% gasAndTwoAcidBasePairs + + modelDefinitionFunction = @(currentParameterVector,energyValues)... + gaussPeak(energyValues, currentParameterVector(1), currentParameterVector(2),currentParameterVector(3)+currentParameterVector(14)) ... % Gauss peak 1: width and position + + gaussPeak(energyValues, currentParameterVector(4),currentParameterVector(2),currentParameterVector(3))... % Gauss peak 2: width and position + + gaussPeak(energyValues, currentParameterVector(5),currentParameterVector(6),currentParameterVector(7))... % Gauss peak 3: unconstrained + + gaussPeak(energyValues, currentParameterVector(8),currentParameterVector(2),currentParameterVector(11)+currentParameterVector(9))... % Gauss peak 4: width and position + + gaussPeak(energyValues, currentParameterVector(10),currentParameterVector(2),currentParameterVector(11))... % Gauss peak 5: width and position + + linearBackground(energyValues,currentParameterVector(12),currentParameterVector(13)); % 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,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)]; - % 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 + % 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') + %% widthPositionAreaRatio + modelDefinitionFunction = @(currentParameterVector,energyValues)... + gaussPeak(energyValues,currentParameterVector(1),currentParameterVector(2),currentParameterVector(3)) ... % Gauss peak 1 + + gaussPeak(energyValues,currentParameterVector(1)*FitParameter.peakArea(1,2),currentParameterVector(2),currentParameterVector(3)+FitParameter.peakPos(1,2))... % Gauss peak 2 + + gaussPeak(energyValues,currentParameterVector(1)*FitParameter.peakArea(1,3),currentParameterVector(2),currentParameterVector(3)+FitParameter.peakPos(1,3))... % Gauss peak 3 + + gaussPeak(energyValues,currentParameterVector(1)*FitParameter.peakArea(1,4),currentParameterVector(2),currentParameterVector(3)+FitParameter.peakPos(1,4))... % Gauss peak 4 + + gaussPeak(energyValues,currentParameterVector(1)*FitParameter.peakArea(1,5),currentParameterVector(2),currentParameterVector(3)+FitParameter.peakPos(1,5))... % Gauss peak 5 + +linearBackground(energyValues,currentParameterVector(4),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') + %% widthPositionAmineAreaRatio + + modelDefinitionFunction = @(currentParameterVector,energyValues)... + gaussPeak(energyValues, currentParameterVector(1),currentParameterVector(2), currentParameterVector(3)) ... % Gauss peak 1 + + gaussPeak(energyValues, currentParameterVector(4),currentParameterVector(2), currentParameterVector(3)+FitParameter.peakPos(1,2))... % Gauss peak 2 + + gaussPeak(energyValues, currentParameterVector(4)*FitParameter.peakArea(1,3),currentParameterVector(2), currentParameterVector(3)+FitParameter.peakPos(1,3))... % Gauss peak 3 + + gaussPeak(energyValues, currentParameterVector(5),currentParameterVector(2),currentParameterVector(3)+FitParameter.peakPos(1,4)) ... % Gauss peak 4 + + gaussPeak(energyValues, currentParameterVector(6),currentParameterVector(2),currentParameterVector(3)+FitParameter.peakPos(1,5)) ... % Gauss peak 5 + + linearBackground(energyValues,currentParameterVector(7), currentParameterVector(8)); + + % 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 end +%% Core part: Run the fit for ii = 1:size(DataTable(p.Results.tableSelection,:),1) if contains(p.Results.energyScale,'kinetic') %-------------------------------------------------------------------------- - % Prepare data (kinetic energy scale) + %% Prepare data (kinetic energy scale) % optionally cut to a section of relevance %-------------------------------------------------------------------------- @@ -492,39 +1947,20 @@ for ii = 1:size(DataTable(p.Results.tableSelection,:),1) end %-------------------------------------------------------------------------- - % Test to implement regions with less weight - % - %-------------------------------------------------------------------------- - if isfield(FitParameter, 'weight') - % Define the energy range you want to weight less - - lowWeightRangeMin = FitParameter.weight(1); % Adjust to your range - lowWeightRangeMax = FitParameter.weight(2); % Adjust to your range - lowWeight = FitParameter.weight(3); % Weight factor (0.1 = 10% of normal weight) - - % Create weight vector - weights = ones(size(fitIntensityValues)); - lowWeightIndices = (fitEnergyValues >= lowWeightRangeMin) & (fitEnergyValues <= lowWeightRangeMax); - weights(lowWeightIndices) = lowWeight; - else - weights = ones(size(fitIntensityValues)); - end - - %-------------------------------------------------------------------------- - % Core part, + %% Core fitting 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(weights .*(fitIntensityValues - modelEvaluationFunction(currentParameterVector,fitEnergyValues)).^2); + 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) + %% Prepare data (binding energy scale) % optionally cut to a section of relevance %-------------------------------------------------------------------------- @@ -545,378 +1981,1280 @@ for ii = 1:size(DataTable(p.Results.tableSelection,:),1) fitIntensityValues = fitIntensityValues(tmp1:tmp2); fitEnergyValues = fitEnergyValues(tmp1:tmp2); end - - if isfield(FitParameter, 'weight') - % Define the energy range you want to weight less - - lowWeightRangeMin = FitParameter.weight(1); % Adjust to your range - lowWeightRangeMax = FitParameter.weight(2); % Adjust to your range - lowWeight = FitParameter.weight(3); % Weight factor (0.1 = 10% of normal weight) - - % Create weight vector - weights = ones(size(fitIntensityValues)); - lowWeightIndices = (fitEnergyValues <= lowWeightRangeMin) & (fitEnergyValues >= lowWeightRangeMax); - weights(lowWeightIndices) = lowWeight; - else - weights = ones(size(fitIntensityValues)); - end %-------------------------------------------------------------------------- - % Core part, + %% Core fitting 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(weights .*(fitIntensityValues - modelEvaluationFunction(currentParameterVector,fitEnergyValues)).^2); + 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); end - %-------------------------------------------------------------------------- - % Compile results, into variable standardizedFitResults with same dimension - % for all fit-types (Area Width Position) - %-------------------------------------------------------------------------- - try - - if matches(p.Results.mode,'position') - standardizedFitResults = [currentParameterVector(1) currentParameterVector(2) currentParameterVector(3)... - currentParameterVector(4) currentParameterVector(5) currentParameterVector(3)+FitParameter.peakPos(1,2) ... - currentParameterVector(6) currentParameterVector(7) currentParameterVector(3)+FitParameter.peakPos(1,3)... - currentParameterVector(8) currentParameterVector(9) currentParameterVector(3)+FitParameter.peakPos(1,4)... - currentParameterVector(10) currentParameterVector(11)... - NaN NaN NaN]; - elseif matches(p.Results.mode,'unconstrained') - standardizedFitResults = [currentParameterVector(1) currentParameterVector(2) currentParameterVector(3)... - currentParameterVector(4) currentParameterVector(5) currentParameterVector(6) ... - currentParameterVector(7) currentParameterVector(8) currentParameterVector(9) ... - currentParameterVector(10) currentParameterVector(11) currentParameterVector(12) ... - currentParameterVector(13) currentParameterVector(14)... - NaN NaN NaN]; - - - elseif matches(p.Results.mode,'widthPosition') - standardizedFitResults = [currentParameterVector(1) currentParameterVector(2) currentParameterVector(3)... - currentParameterVector(4) currentParameterVector(2) currentParameterVector(3)+FitParameter.peakPos(1,2)... - currentParameterVector(5) currentParameterVector(2) currentParameterVector(3)+FitParameter.peakPos(1,3)... - currentParameterVector(6) currentParameterVector(2) currentParameterVector(3)+FitParameter.peakPos(1,4)... - currentParameterVector(8) currentParameterVector(9)... - currentParameterVector(7) currentParameterVector(2) currentParameterVector(3)+FitParameter.peakPos(1,5)]; - elseif matches(p.Results.mode,'width') - standardizedFitResults = [currentParameterVector(1) currentParameterVector(2) currentParameterVector(3)... - currentParameterVector(4) currentParameterVector(2) currentParameterVector(5)... - currentParameterVector(6) currentParameterVector(2) currentParameterVector(7)... - currentParameterVector(8) currentParameterVector(2) currentParameterVector(9)... - currentParameterVector(10) currentParameterVector(11)... - NaN NaN NaN]; - elseif matches(p.Results.mode,'unconstrainedAndWidthPosition') - standardizedFitResults = [currentParameterVector(1) currentParameterVector(2) currentParameterVector(3)+currentParameterVector(14)... - currentParameterVector(4) currentParameterVector(2) currentParameterVector(3)... - currentParameterVector(5) currentParameterVector(6) currentParameterVector(7)... - currentParameterVector(8) currentParameterVector(2) currentParameterVector(9)... - currentParameterVector(12) currentParameterVector(13) ... - currentParameterVector(10) currentParameterVector(2) currentParameterVector(11)]; - elseif matches(p.Results.mode,'gasAndTwoAcidBasePairs') - standardizedFitResults = [currentParameterVector(1) currentParameterVector(2) currentParameterVector(3)+currentParameterVector(14)... - currentParameterVector(4) currentParameterVector(2) currentParameterVector(3)... - currentParameterVector(5) currentParameterVector(6) currentParameterVector(7)... - currentParameterVector(8) currentParameterVector(2) currentParameterVector(11)+currentParameterVector(9)... - currentParameterVector(12) currentParameterVector(13) ... - currentParameterVector(10) currentParameterVector(2) currentParameterVector(11)]; - elseif matches(p.Results.mode,'widthPositionAreaRatio') - standardizedFitResults = [currentParameterVector(1) currentParameterVector(2) currentParameterVector(3)... - currentParameterVector(1)*FitParameter.peakArea(1,2) currentParameterVector(2) currentParameterVector(3)+FitParameter.peakPos(1,2)... - currentParameterVector(1)*FitParameter.peakArea(1,3) currentParameterVector(2) currentParameterVector(3)+FitParameter.peakPos(1,3)... - currentParameterVector(1)*FitParameter.peakArea(1,4) currentParameterVector(2) currentParameterVector(3)+FitParameter.peakPos(1,4)... - currentParameterVector(4) currentParameterVector(5)... - currentParameterVector(1)*FitParameter.peakArea(1,5) currentParameterVector(2) currentParameterVector(3)+FitParameter.peakPos(1,5)]; - elseif matches(p.Results.mode,'TwoSpinOrbitDoublet') - standardizedFitResults = [currentParameterVector(1) currentParameterVector(2) currentParameterVector(3)... - currentParameterVector(1)*FitParameter.peakArea(1,2) currentParameterVector(2) currentParameterVector(3)+FitParameter.peakPos(1,2)... - currentParameterVector(4) currentParameterVector(2) currentParameterVector(3)+currentParameterVector(5)... - currentParameterVector(4)*FitParameter.peakArea(1,2) currentParameterVector(2) currentParameterVector(3)+currentParameterVector(5)+FitParameter.peakPos(1,2)... - currentParameterVector(6) currentParameterVector(7)... - NaN NaN NaN]; - - elseif matches(p.Results.mode,'widthPositionAmineAreaRatio') - standardizedFitResults = [currentParameterVector(1) currentParameterVector(2) currentParameterVector(3)... - currentParameterVector(4) currentParameterVector(2) currentParameterVector(3)+FitParameter.peakPos(1,2)... - currentParameterVector(4)*FitParameter.peakArea(1,3) currentParameterVector(2) currentParameterVector(3)+FitParameter.peakPos(1,3)... - currentParameterVector(5) currentParameterVector(2) currentParameterVector(3)+FitParameter.peakPos(1,4)... - currentParameterVector(7) currentParameterVector(8)... - currentParameterVector(6) currentParameterVector(2) currentParameterVector(3)+FitParameter.peakPos(1,5)]; - end - - %-------------------------------------------------------------------------- - % Print results - %-------------------------------------------------------------------------- - - if matches(p.Results.showFigures,'on') - % Display results in table - disp(' ' ); - disp(' ' ); - disp('Fit results') - disp(DataTable(p.Results.tableSelection,:).name{ii}); - disp(' ' ); - if (standardizedFitResults(1) <= 1 || standardizedFitResults(4) <= 1) - disp(cell2table({sprintf('%.1f',standardizedFitResults(3)),... - sprintf('%.1f',standardizedFitResults(2)*sqrt( 2 * log(2))),... - sprintf('%.1f',standardizedFitResults(2)),... - sprintf('%.3f',standardizedFitResults(1));... - sprintf('%.1f',standardizedFitResults(6)),... - sprintf('%.1f',standardizedFitResults(5)*sqrt( 2 * log(2))),... - sprintf('%.1f',standardizedFitResults(5)),... - sprintf('%.3f',standardizedFitResults(4));... - sprintf('%.1f',standardizedFitResults(9)),... - sprintf('%.1f',standardizedFitResults(8)*sqrt( 2 * log(2))),... - sprintf('%.1f',standardizedFitResults(8)),... - sprintf('%.3f',standardizedFitResults(7));... - sprintf('%.1f',standardizedFitResults(12)),... - sprintf('%.1f',standardizedFitResults(11)*sqrt( 2 * log(2))),... - sprintf('%.1f',standardizedFitResults(11)),... - sprintf('%.3f',standardizedFitResults(10)); ... - sprintf('%.1f',standardizedFitResults(17)),... - sprintf('%.1f',standardizedFitResults(16)*sqrt( 2 * log(2))),... - sprintf('%.1f',standardizedFitResults(16)),... - sprintf('%.3f',standardizedFitResults(15))},... - 'VariableNames',{'position';'FWHW';'width';'area'},... - 'RowNames',{'peak1','peak2','peak3','peak4','peak5'})); % 'RowNames',{'peak1','peak2','peak3','peak4','peak5'})); - else - disp(cell2table({sprintf('%.1f',standardizedFitResults(3)),... - sprintf('%.1f',standardizedFitResults(2)*sqrt( 2 * log(2))),... - sprintf('%.1f',standardizedFitResults(2)),... - sprintf('%7.0f',standardizedFitResults(1));... - sprintf('%.1f',standardizedFitResults(6)),... - sprintf('%.1f',standardizedFitResults(5)*sqrt( 2 * log(2))),... - sprintf('%.1f',standardizedFitResults(5)),... - sprintf('%7.0f',standardizedFitResults(4));... - sprintf('%.1f',standardizedFitResults(9)),... - sprintf('%.1f',standardizedFitResults(8)*sqrt( 2 * log(2))),... - sprintf('%.1f',standardizedFitResults(8)),... - sprintf('%7.0f',standardizedFitResults(7));... - sprintf('%.1f',standardizedFitResults(12)),... - sprintf('%.1f',standardizedFitResults(11)*sqrt( 2 * log(2))),... - sprintf('%.1f',standardizedFitResults(11)),... - sprintf('%7.0f',standardizedFitResults(10)); ... - sprintf('%.1f',standardizedFitResults(17)),... - sprintf('%.1f',standardizedFitResults(16)*sqrt( 2 * log(2))),... - sprintf('%.1f',standardizedFitResults(16)),... - sprintf('%7.0f',standardizedFitResults(15))},... - 'VariableNames',{'position';'FWHW';'width';'area'},... - 'RowNames',{'peak1','peak2','peak3','peak4','peak5'})); % 'RowNames',{'peak1','peak2','peak3','peak4','peak5'})); - end - %-------------------------------------------------------------------------- - % Plot results - %-------------------------------------------------------------------------- - - % eliminate NaN in standardizedFitResults - standardizedFitResults(isnan(standardizedFitResults))=0.01; - % define equations to be plotted - Gauss = @(energyValues, A, Pos, Wdt) A./(Wdt.*sqrt(pi/2))*exp(-2.*((energyValues - Pos).^2/Wdt.^2)); - linear = @(energyValues,slp,intrcpt) slp*energyValues+intrcpt; - figure( 'Name', DataTable(p.Results.tableSelection,:).name{ii}) - t = tiledlayout(3,1); - - ax(1) = nexttile(2,[2,1]); - % define high-resolution energyValues for plotting - if fitEnergyValues(1)<= fitEnergyValues(end) - plotEnergyValues = fitEnergyValues(1):0.01:fitEnergyValues(end); - else - plotEnergyValues = fitEnergyValues(1):-0.01:fitEnergyValues(end); - end - - %Define gauss peaks - yData_pk1= Gauss(plotEnergyValues ,standardizedFitResults(1), standardizedFitResults(3), standardizedFitResults(2)); - yData_pk2= Gauss(plotEnergyValues ,standardizedFitResults(4), standardizedFitResults(6), standardizedFitResults(5)); - yData_pk3= Gauss(plotEnergyValues ,standardizedFitResults(7), standardizedFitResults(9), standardizedFitResults(8)); - yData_pk4= Gauss(plotEnergyValues ,standardizedFitResults(10), standardizedFitResults(12), standardizedFitResults(11)); - yData_pk5= Gauss(plotEnergyValues ,standardizedFitResults(15), standardizedFitResults(17), standardizedFitResults(16)); - bkgrd = linear(plotEnergyValues,standardizedFitResults(13),standardizedFitResults(14)); - - - % plot fit - h = plot(plotEnergyValues,yData_pk1+yData_pk2+yData_pk3+yData_pk4+yData_pk5+bkgrd); - hold on - h.LineWidth=2; - h.Color = PsiFarbe.dunkelblau; - % plot each gauss peak - h = plot(plotEnergyValues,yData_pk1+bkgrd); - h.LineWidth=1.5; - h.Color = PsiFarbe.gelb; - h = plot(plotEnergyValues,yData_pk2+bkgrd); - h.LineWidth=1.5; - h.Color = PsiFarbe.gelb; - h = plot(plotEnergyValues,yData_pk3+bkgrd); - h.LineWidth=1.5; - h.Color = PsiFarbe.gelb; - h = plot(plotEnergyValues,yData_pk4+bkgrd); - h.LineWidth=1.5; - h.Color = PsiFarbe.gelb; - h = plot(plotEnergyValues,yData_pk5+bkgrd); - h.LineWidth=1.5; - h.Color = PsiFarbe.gelb; - plot(plotEnergyValues, bkgrd); - if fitEnergyValues(1)<= fitEnergyValues(end) - xlabel('Kinetic Energy [eV]'); - else - set(ax(1),'Xdir','reverse'); - xlabel('Binding Energy [eV]'); - end - % plot data - plot(fitEnergyValues,fitIntensityValues,'o',"MarkerFaceColor","white") - legend OFF - if matches(p.Results.tableVariable,'spectrum_countsPerSecond') - ylabel('PE Intensity [counts per second]'); - else - ylabel('PE intensity [-]'); - end - grid on - - -% plot residuals -ax(2) = nexttile(1,[1,1]); -%Define gauss peaks with data scaling - yData_pk1= Gauss(fitEnergyValues ,standardizedFitResults(1), standardizedFitResults(3), standardizedFitResults(2)); - yData_pk2= Gauss(fitEnergyValues ,standardizedFitResults(4), standardizedFitResults(6), standardizedFitResults(5)); - yData_pk3= Gauss(fitEnergyValues ,standardizedFitResults(7), standardizedFitResults(9), standardizedFitResults(8)); - yData_pk4= Gauss(fitEnergyValues ,standardizedFitResults(10), standardizedFitResults(12), standardizedFitResults(11)); - yData_pk5= Gauss(fitEnergyValues ,standardizedFitResults(15), standardizedFitResults(17), standardizedFitResults(16)); - bkgrd = linear(fitEnergyValues,standardizedFitResults(13),standardizedFitResults(14)); - -residuals = fitIntensityValues - (yData_pk1+yData_pk2+yData_pk3+yData_pk4+yData_pk5+bkgrd); % Calculate residuals -plot(fitEnergyValues, residuals, 'b.'); -ylabel('Residuals (Data - Fit)'); -% title('Residuals Plot'); -grid on; - -if fitEnergyValues(1)<= fitEnergyValues(end) - % xlabel('Kinetic Energy [eV]'); -else - set(ax(2),'Xdir','reverse'); - % xlabel('Binding Energy [eV]'); -end - -title(t,DataTable(p.Results.tableSelection,:).name{ii}, 'Interpreter', 'none'); - end - - catch - disp("Plotting and listing results did not work") - end - %-------------------------------------------------------------------------- - % Export results - % prepare for data table - %-------------------------------------------------------------------------- - resultsPeakPos = []; -resultsPeakFWHW = []; -resultsPeakArea = []; -resultsBkgrd = []; -resultsPeakEta = []; -resultsPeakShirleyK = []; - - if matches(p.Results.backgroundMode,'linear') && ~exist("FitParameter.peakEta","var") - - for i = 1:size(FitParameter.peakArea,2) - 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)] ; - end - resultsBkgrd(ii,:) = [standardizedFitResults(13), standardizedFitResults(14)] ; - - 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 +%% Results Extraction +% +% Extract fitting results from currentParameterVector for EACH MODEL TYPE +% and populate result arrays: resultsPeakPos, resultsPeakArea, resultsPeakFWHM, resultsBkgrd, etc. -% write results to data table -try - existColumn = strcmp('peakPosition',DataTable.Properties.VariableNames); - if ~any(existColumn) - DataTable.peakPosition = NaN(size(DataTable,1),size(resultsPeakPos,2)); +resultsPeakPos = []; +resultsPeakArea = []; +resultsPeakFWHM = []; +resultsBkgrd = []; +resultsPeakEta = []; +resultsPeakShirleyK = []; +nPeaks = size(FitParameter.peakArea, 2); + +for iii= 1 + + %% Generic Models + + if matches(p.Results.mode, 'unconstrained') + + %% MODE: UNCONSTRAINED + + if matches(p.Results.backgroundMode, 'linear') && ~isfield(FitParameter,'peakEta') + % ===== GAUSSIAN + LINEAR (UNCONSTRAINED) ===== + + for iPeak = 1:nPeaks + areaIdx = 3*iPeak - 2; + widthIdx = 3*iPeak - 1; + posIdx = 3*iPeak; + + resultsPeakArea = [resultsPeakArea, currentParameterVector(areaIdx)]; + resultsPeakFWHM = [resultsPeakFWHM, currentParameterVector(widthIdx) ]; + resultsPeakPos = [resultsPeakPos, currentParameterVector(posIdx)]; + end + resultsBkgrd = [currentParameterVector(end-1), currentParameterVector(end)]; + + elseif matches(p.Results.backgroundMode, 'linear') && isfield(FitParameter,'peakEta') + % ===== PSEUDO-VOIGT + LINEAR (UNCONSTRAINED) ===== + + for iPeak = 1:nPeaks + areaIdx = 4*iPeak - 3; + widthIdx = 4*iPeak - 2; + posIdx = 4*iPeak - 1; + etaIdx = 4*iPeak; + + resultsPeakArea = [resultsPeakArea, currentParameterVector(areaIdx)]; %#ok<*AGROW> + resultsPeakFWHM = [resultsPeakFWHM, currentParameterVector(widthIdx) ]; + resultsPeakPos = [resultsPeakPos, currentParameterVector(posIdx)]; + resultsPeakEta = [resultsPeakEta, currentParameterVector(etaIdx)]; + end + resultsBkgrd = [currentParameterVector(end-1), currentParameterVector(end)]; + + elseif matches(p.Results.backgroundMode, 'shirley') && isfield(FitParameter,'peakEta') + % ===== PSEUDO-VOIGT + SHIRLEY (UNCONSTRAINED) ===== + for iPeak = 1:nPeaks + areaIdx = 4*iPeak - 3; + fwhmIdx = 4*iPeak - 2; + posIdx = 4*iPeak - 1; + etaIdx = 4*iPeak; + + resultsPeakArea = [resultsPeakArea, currentParameterVector(areaIdx)]; + resultsPeakFWHM = [resultsPeakFWHM, currentParameterVector(fwhmIdx)]; + resultsPeakPos = [resultsPeakPos, currentParameterVector(posIdx)]; + resultsPeakEta = [resultsPeakEta, currentParameterVector(etaIdx)]; + end + + for iPeak = 1:nPeaks + resultsPeakShirleyK = [resultsPeakShirleyK, currentParameterVector(4*nPeaks+iPeak)]; + end + % background (Shirley intercept only) + resultsBkgrd = [0,currentParameterVector(end)]; + + elseif matches(p.Results.backgroundMode, 'shirley') && ~isfield(FitParameter,'peakEta') + % ===== GAUSSIAN + SHIRLEY (UNCONSTRAINED) ===== + for iPeak = 1:nPeaks + areaIdx = 3*iPeak - 2; + fwhmIdx = 3*iPeak - 1; + posIdx = 3*iPeak; + + resultsPeakArea = [resultsPeakArea, currentParameterVector(areaIdx)]; + resultsPeakFWHM = [resultsPeakFWHM, currentParameterVector(fwhmIdx)]; + resultsPeakPos = [resultsPeakPos, currentParameterVector(posIdx)]; + end + + for iPeak = 1:nPeaks + resultsPeakShirleyK = [resultsPeakShirleyK, currentParameterVector(3*nPeaks+iPeak)]; + end + % background (Shirley intercept only) + resultsBkgrd = [0,currentParameterVector(end)]; + + end + + elseif matches(p.Results.mode, 'position') + %% MODE: Position + if matches(p.Results.backgroundMode, 'linear') && ~isfield(FitParameter,'peakEta') + % ===== GAUSSIAN + LINEAR (position) ===== + for iPeak = 1:nPeaks + areaIdx = 2*iPeak - 1; + widthIdx = 2*iPeak; + resultsPeakArea = [resultsPeakArea, currentParameterVector(areaIdx)]; + resultsPeakFWHM = [resultsPeakFWHM, currentParameterVector(widthIdx)]; + end + pos1_fit = currentParameterVector(2*nPeaks + 1); + + for iPeak = 1:nPeaks + if iPeak == 1 + peakPos = pos1_fit; + else + peakPos = pos1_fit + FitParameter.peakPos(1, iPeak); + end + resultsPeakPos = [resultsPeakPos, peakPos]; + end + resultsBkgrd = [currentParameterVector(end-1), currentParameterVector(end)]; + + elseif matches(p.Results.backgroundMode, 'linear') && isfield(FitParameter,'peakEta') + % ===== PSEUDO-VOIGT + LINEAR (position) ===== + for iPeak = 1:nPeaks + areaIdx = 3*iPeak - 2; + widthIdx = 3*iPeak - 1; + etaIdx = 3*iPeak; + + resultsPeakArea = [resultsPeakArea, currentParameterVector(areaIdx)]; + resultsPeakFWHM = [resultsPeakFWHM, currentParameterVector(widthIdx)]; + resultsPeakEta = [resultsPeakEta, currentParameterVector(etaIdx)]; + end + pos1_fit = currentParameterVector(3*nPeaks + 1); + + for iPeak = 1:nPeaks + if iPeak == 1 + peakPos = pos1_fit; + else + peakPos = pos1_fit + FitParameter.peakPos(1, iPeak); + end + resultsPeakPos = [resultsPeakPos, peakPos]; + end + resultsBkgrd = [currentParameterVector(end-1), currentParameterVector(end)]; + + elseif matches(p.Results.backgroundMode, 'shirley') && isfield(FitParameter,'peakEta') + % ===== PSEUDO-VOIGT + SHIRLEY (position) ===== + for iPeak = 1:nPeaks + areaIdx = 3*iPeak - 2; + widthIdx = 3*iPeak - 1; + etaIdx = 3*iPeak; + kIdx = 3*nPeaks +1 +iPeak; + + resultsPeakArea = [resultsPeakArea, currentParameterVector(areaIdx)]; + resultsPeakFWHM = [resultsPeakFWHM, currentParameterVector(widthIdx) ]; + resultsPeakEta = [resultsPeakEta, currentParameterVector(etaIdx)]; + resultsPeakShirleyK = [resultsPeakShirleyK, currentParameterVector(kIdx)]; + end + pos1_fit = currentParameterVector(3*nPeaks + 1); + + for iPeak = 1:nPeaks + if iPeak == 1 + peakPos = pos1_fit; + else + peakPos = pos1_fit + FitParameter.peakPos(1, iPeak); + end + resultsPeakPos = [resultsPeakPos, peakPos]; + end + % background (Shirley intercept only) + resultsBkgrd = [0,currentParameterVector(end)]; + elseif matches(p.Results.backgroundMode, 'shirley') && ~isfield(FitParameter,'peakEta') + % ===== GAUSSIAN + SHIRLEY (position) ===== + for iPeak = 1:nPeaks + areaIdx = 2*iPeak - 1; + widthIdx = 2*iPeak; + kIdx = 2*nPeaks+1+iPeak; + + resultsPeakArea = [resultsPeakArea, currentParameterVector(areaIdx)]; + resultsPeakFWHM = [resultsPeakFWHM, currentParameterVector(widthIdx) ]; + resultsPeakShirleyK = [resultsPeakShirleyK, currentParameterVector(kIdx)]; + end + pos1_fit = currentParameterVector(2*nPeaks+1); + + for iPeak = 1:nPeaks + if iPeak == 1 + peakPos = pos1_fit; + else + peakPos = pos1_fit + FitParameter.peakPos(1, iPeak); + end + resultsPeakPos = [resultsPeakPos, peakPos]; + end + % background (Shirley intercept only) + resultsBkgrd = [0,currentParameterVector(end)]; + end + + elseif matches(p.Results.mode, 'positionShirley') + %% MODE: Position-Shirley + + if matches(p.Results.backgroundMode, 'shirley') && isfield(FitParameter,'peakEta') + % ===== PSEUDO-VOIGT + SHIRLEY (positionShirley) ===== + for iPeak = 1:nPeaks + areaIdx = 3*iPeak - 2; + widthIdx = 3*iPeak - 1; + etaIdx = 3*iPeak; + + resultsPeakArea = [resultsPeakArea, currentParameterVector(areaIdx)]; + resultsPeakFWHM = [resultsPeakFWHM, currentParameterVector(widthIdx) ]; + resultsPeakEta = [resultsPeakEta, currentParameterVector(etaIdx)]; + resultsPeakShirleyK = [resultsPeakShirleyK, currentParameterVector(3*nPeaks +2)]; + end + pos1_fit = currentParameterVector(3*nPeaks + 1); + + for iPeak = 1:nPeaks + if iPeak == 1 + peakPos = pos1_fit; + else + peakPos = pos1_fit + FitParameter.peakPos(1, iPeak); + end + resultsPeakPos = [resultsPeakPos, peakPos]; + end + % background (Shirley intercept only) + resultsBkgrd = [0,currentParameterVector(end)]; + elseif matches(p.Results.backgroundMode, 'shirley') && ~isfield(FitParameter,'peakEta') + + for iPeak = 1:nPeaks + areaIdx = 2*iPeak - 1; + widthIdx = 2*iPeak; + kIdx = 2*nPeaks+2; + + resultsPeakArea = [resultsPeakArea, currentParameterVector(areaIdx)]; + resultsPeakFWHM = [resultsPeakFWHM, currentParameterVector(widthIdx) ]; + resultsPeakShirleyK = [resultsPeakShirleyK, currentParameterVector(kIdx)]; + end + pos1_fit = currentParameterVector(2*nPeaks+1); + + for iPeak = 1:nPeaks + if iPeak == 1 + peakPos = pos1_fit; + else + peakPos = pos1_fit + FitParameter.peakPos(1, iPeak); + end + resultsPeakPos = [resultsPeakPos, peakPos]; + end + % background (Shirley intercept only) + resultsBkgrd = [0,currentParameterVector(end)]; + end + %% MODE: WIDTH + + elseif matches(p.Results.mode, 'width') + + if matches(p.Results.backgroundMode, 'linear') && ~isfield(FitParameter,'peakEta') + % ===== GAUSSIAN + LINEAR (width) ===== + for iPeak = 1:nPeaks + areaIdx = 2*iPeak - 1; + posIdx = 2*iPeak; + + resultsPeakArea = [resultsPeakArea, currentParameterVector(areaIdx)]; + resultsPeakPos = [resultsPeakPos, currentParameterVector(posIdx)]; + end + width_shared = currentParameterVector(2*nPeaks + 1); + resultsPeakFWHM = ones(1, nPeaks) * width_shared ; + resultsBkgrd = [currentParameterVector(end-1), currentParameterVector(end)]; + + elseif matches(p.Results.backgroundMode, 'linear') && isfield(FitParameter,'peakEta') + % ===== PSEUDO-VOIGT + LINEAR (width) ===== + for iPeak = 1:nPeaks + areaIdx = 3*iPeak - 2; + posIdx = 3*iPeak - 1; + etaIdx = 3*iPeak; + + resultsPeakArea = [resultsPeakArea, currentParameterVector(areaIdx)]; + resultsPeakPos = [resultsPeakPos, currentParameterVector(posIdx)]; + resultsPeakEta = [resultsPeakEta, currentParameterVector(etaIdx)]; + end + width_shared = currentParameterVector(3*nPeaks + 1); + resultsPeakFWHM = ones(1, nPeaks) * width_shared ; + resultsBkgrd = [currentParameterVector(end-1), currentParameterVector(end)]; + + elseif matches(p.Results.backgroundMode, 'shirley') && ~isfield(FitParameter,'peakEta') + % ===== GAUSSIAN + SHIRLEY (width) ===== + for iPeak = 1:nPeaks + areaIdx = 2*iPeak - 1; + posIdx = 2*iPeak; + kIdx = 2*nPeaks +1 +iPeak; + + resultsPeakArea = [resultsPeakArea, currentParameterVector(areaIdx)]; + resultsPeakPos = [resultsPeakPos, currentParameterVector(posIdx)]; + resultsPeakShirleyK = [resultsPeakShirleyK, currentParameterVector(kIdx)]; + end + width_shared = currentParameterVector(2*nPeaks +1); + resultsPeakFWHM = ones(1, nPeaks) * width_shared ; + % background (Shirley intercept only) + resultsBkgrd = [0,currentParameterVector(end)]; + + elseif matches(p.Results.backgroundMode, 'shirley') && isfield(FitParameter,'peakEta') + % ===== PSEUDO-VOIGT + LINEAR (width) ===== + for iPeak = 1:nPeaks + areaIdx = 3*iPeak - 2; + posIdx = 3*iPeak - 1; + etaIdx = 3*iPeak; + kIdx = 3*nPeaks +1 +iPeak; + + resultsPeakArea = [resultsPeakArea, currentParameterVector(areaIdx)]; + resultsPeakPos = [resultsPeakPos, currentParameterVector(posIdx)]; + resultsPeakEta = [resultsPeakEta, currentParameterVector(etaIdx)]; + resultsPeakShirleyK = [resultsPeakShirleyK, currentParameterVector(kIdx)]; + end + width_shared = currentParameterVector(3*nPeaks + 1); + resultsPeakFWHM = ones(1, nPeaks) * width_shared ; + % background (Shirley intercept only) + resultsBkgrd = [0,currentParameterVector(end)]; + end + + elseif matches(p.Results.mode, 'widthShirley') + %% MODE: widthShirley + if matches(p.Results.backgroundMode, 'shirley') && isfield(FitParameter,'peakEta') + % ===== PSEUDO-VOIGT + SHIRLEY (widthShirley) ===== + for iPeak = 1:nPeaks + areaIdx = 3*iPeak - 2; + posIdx = 3*iPeak - 1; + etaIdx = 3*iPeak; + kIdx = 3*nPeaks +2; + + resultsPeakArea = [resultsPeakArea, currentParameterVector(areaIdx)]; + resultsPeakPos = [resultsPeakPos, currentParameterVector(posIdx)]; + resultsPeakEta = [resultsPeakEta, currentParameterVector(etaIdx)]; + resultsPeakShirleyK = [resultsPeakShirleyK, currentParameterVector(kIdx)]; + end + width_shared = currentParameterVector(3*nPeaks + 1); + resultsPeakFWHM = ones(1, nPeaks) * width_shared ; + % background (Shirley intercept only) + resultsBkgrd = [0,currentParameterVector(end)]; + + elseif matches(p.Results.backgroundMode, 'shirley') && ~isfield(FitParameter,'peakEta') + % ===== GAUSSIAN + SHIRLEY (widthShirley) ===== + for iPeak = 1:nPeaks + areaIdx = 2*iPeak - 1; + posIdx = 2*iPeak; + kIdx = 2*nPeaks +2; + + resultsPeakArea = [resultsPeakArea, currentParameterVector(areaIdx)]; + resultsPeakPos = [resultsPeakPos, currentParameterVector(posIdx)]; + resultsPeakShirleyK = [resultsPeakShirleyK, currentParameterVector(kIdx)]; + end + width_shared = currentParameterVector(2*nPeaks +1); + resultsPeakFWHM = ones(1, nPeaks) * width_shared ; + % background (Shirley intercept only) + resultsBkgrd = [0,currentParameterVector(end)]; + end + %% MODE: WIDTH-POSITION + + elseif matches(p.Results.mode, 'widthPosition') + % ===== GAUSSIAN + LINEAR (widthPosition) ===== + if matches(p.Results.backgroundMode, 'linear') && ~isfield(FitParameter,'peakEta') + for iPeak = 1:nPeaks + resultsPeakArea = [resultsPeakArea, currentParameterVector(iPeak)]; + end + pos1_fit = currentParameterVector(nPeaks + 2); + width_shared = currentParameterVector(nPeaks + 1); + + for iPeak = 1:nPeaks + if iPeak == 1 + peakPos = pos1_fit; + else + peakPos = pos1_fit + FitParameter.peakPos(1, iPeak); + end + resultsPeakPos = [resultsPeakPos, peakPos]; + end + resultsPeakFWHM = ones(1, nPeaks) * width_shared ; + resultsBkgrd = [currentParameterVector(end-1), currentParameterVector(end)]; + + elseif matches(p.Results.backgroundMode, 'linear') && isfield(FitParameter,'peakEta') + % ===== PSEUDO-VOIGT + LINEAR (widthPosition) ===== + for iPeak = 1:nPeaks + areaIdx = 2*iPeak - 1; + etaIdx = 2*iPeak; + + resultsPeakArea = [resultsPeakArea, currentParameterVector(areaIdx)]; + resultsPeakEta = [resultsPeakEta, currentParameterVector(etaIdx)]; + end + pos1_fit = currentParameterVector(2*nPeaks + 2); + width_shared = currentParameterVector(2*nPeaks + 1); + + for iPeak = 1:nPeaks + if iPeak == 1 + peakPos = pos1_fit; + else + peakPos = pos1_fit + FitParameter.peakPos(1, iPeak); + end + resultsPeakPos = [resultsPeakPos, peakPos]; + end + resultsPeakFWHM = ones(1, nPeaks) * width_shared ; + resultsBkgrd = [currentParameterVector(end-1), currentParameterVector(end)]; + + elseif matches(p.Results.backgroundMode, 'shirley') && ~isfield(FitParameter,'peakEta') + % ===== GAUSSIAN + SHIRLEY (widthPosition) ===== + for iPeak = 1:nPeaks + areaIdx = iPeak; + kIdx = nPeaks + 2 +iPeak; + + resultsPeakArea = [resultsPeakArea, currentParameterVector(areaIdx)]; + resultsPeakShirleyK = [resultsPeakShirleyK, currentParameterVector(kIdx)]; + end + pos1_fit = currentParameterVector(nPeaks + 2); + width_shared = currentParameterVector(nPeaks + 1); + + for iPeak = 1:nPeaks + if iPeak == 1 + peakPos = pos1_fit; + else + peakPos = pos1_fit + FitParameter.peakPos(1, iPeak); + end + resultsPeakPos = [resultsPeakPos, peakPos]; + end + resultsPeakFWHM = ones(1, nPeaks) * width_shared ; + % background (Shirley intercept only) + resultsBkgrd = [0,currentParameterVector(end)]; + + elseif matches(p.Results.backgroundMode, 'shirley') && isfield(FitParameter,'peakEta') + % ===== PSEUDO-VOIGT + LINEAR (widthPosition) ===== + for iPeak = 1:nPeaks + areaIdx = 2*iPeak - 1; + etaIdx = 2*iPeak; + kIdx = 2*nPeaks + 2 +iPeak; + + resultsPeakArea = [resultsPeakArea, currentParameterVector(areaIdx)]; + resultsPeakEta = [resultsPeakEta, currentParameterVector(etaIdx)]; + resultsPeakShirleyK = [resultsPeakShirleyK, currentParameterVector(kIdx)]; + end + pos1_fit = currentParameterVector(2*nPeaks + 2); + width_shared = currentParameterVector(2*nPeaks + 1); + + for iPeak = 1:nPeaks + if iPeak == 1 + peakPos = pos1_fit; + else + peakPos = pos1_fit + FitParameter.peakPos(1, iPeak); + end + resultsPeakPos = [resultsPeakPos, peakPos]; + end + resultsPeakFWHM = ones(1, nPeaks) * width_shared ; + % background (Shirley intercept only) + resultsBkgrd = [0,currentParameterVector(end)]; + end + + elseif matches(p.Results.mode, 'widthPositionShirley') + %% MODE: widthPositionShirley + + if matches(p.Results.backgroundMode, 'shirley') && ~isfield(FitParameter,'peakEta') + % ===== GAUSSIAN + SHIRLEY (widthPositionShirley) ===== + for iPeak = 1:nPeaks + areaIdx = iPeak; + kIdx = nPeaks + 3; + + resultsPeakArea = [resultsPeakArea, currentParameterVector(areaIdx)]; + resultsPeakShirleyK = [resultsPeakShirleyK, currentParameterVector(kIdx)]; + end + pos1_fit = currentParameterVector(nPeaks + 2); + width_shared = currentParameterVector(nPeaks + 1); + + for iPeak = 1:nPeaks + if iPeak == 1 + peakPos = pos1_fit; + else + peakPos = pos1_fit + FitParameter.peakPos(1, iPeak); + end + resultsPeakPos = [resultsPeakPos, peakPos]; + end + resultsPeakFWHM = ones(1, nPeaks) * width_shared ; + % background (Shirley intercept only) + resultsBkgrd = [0,currentParameterVector(end)]; + + elseif matches(p.Results.backgroundMode, 'shirley') && isfield(FitParameter,'peakEta') + % ===== PSEUDO-VOIGT + SHIRLEY (widthPositionShirley) ===== + for iPeak = 1:nPeaks + areaIdx = 2*iPeak - 1; + etaIdx = 2*iPeak; + kIdx = 2*nPeaks + 3; + + resultsPeakArea = [resultsPeakArea, currentParameterVector(areaIdx)]; + resultsPeakEta = [resultsPeakEta, currentParameterVector(etaIdx)]; + resultsPeakShirleyK = [resultsPeakShirleyK, currentParameterVector(kIdx)]; + end + pos1_fit = currentParameterVector(2*nPeaks + 2); + width_shared = currentParameterVector(2*nPeaks + 1); + + for iPeak = 1:nPeaks + if iPeak == 1 + peakPos = pos1_fit; + else + peakPos = pos1_fit + FitParameter.peakPos(1, iPeak); + end + resultsPeakPos = [resultsPeakPos, peakPos]; + end + resultsPeakFWHM = ones(1, nPeaks) * width_shared ; + % background (Shirley intercept only) + resultsBkgrd = [0,currentParameterVector(end)]; + + end end - existColumn = strcmp('peakArea',DataTable.Properties.VariableNames); - if ~any(existColumn) - DataTable.peakArea = NaN(size(DataTable,1),size(resultsPeakArea,2)); + %% Custom Models + + + if matches(p.Results.mode, 'unconstrainedAndWidthPosition') && matches(p.Results.backgroundMode, 'linear') + %% unconstrainedAndWidthPosition + % ===== GAUSSIAN + LINEAR ===== + + % Extract parameters + area1 = currentParameterVector(1); + width1 = currentParameterVector(2); + pos2_rel = currentParameterVector(3); + area2 = currentParameterVector(4); + area3 = currentParameterVector(5); + width3 = currentParameterVector(6); + pos3 = currentParameterVector(7); + area4 = currentParameterVector(8); + pos4_offset = currentParameterVector(9); + area5 = currentParameterVector(10); + pos5_offset = currentParameterVector(11); + slope = currentParameterVector(12); + intercept = currentParameterVector(13); + pos1_shift = currentParameterVector(14); + + % Calculate absolute positions + pos1 = pos2_rel + pos1_shift; + pos2 = pos2_rel; + pos4 = pos4_offset; + pos5 = pos5_offset; + + % Results: areas, positions, FWHM + resultsPeakArea = [area1, area2, area3, area4, area5]; + resultsPeakPos = [pos1, pos2, pos3, pos4, pos5]; + resultsPeakFWHM = [width1, width1, width3, width1, width1]; + + % Background + resultsBkgrd = [slope, intercept]; + + elseif matches(p.Results.mode, 'TwoSpinOrbitDoublet') + %% TwoSpinOrbitDoublet + % ===== GAUSSIAN + LINEAR ===== + + % Extract parameters + area1 = currentParameterVector(1); + width1 = currentParameterVector(2); + pos1 = currentParameterVector(3); + area2 = currentParameterVector(1)*FitParameter.peakArea(1,2); + width2 = currentParameterVector(2); + pos2 = currentParameterVector(3)+FitParameter.peakPos(1,2); + area3 = currentParameterVector(4); + width3 = currentParameterVector(2); + pos3 = currentParameterVector(3)+currentParameterVector(5); + area4 = currentParameterVector(4)*FitParameter.peakArea(1,2); + width4 = currentParameterVector(2); + pos4 = currentParameterVector(3)+currentParameterVector(5)++FitParameter.peakPos(1,2); + slope = currentParameterVector(6); + intercept = currentParameterVector(7); + + % Results: areas, positions, FWHM + resultsPeakArea = [area1, area2, area3, area4]; + resultsPeakPos = [pos1, pos2, pos3, pos4]; + resultsPeakFWHM = [width1, width2, width3, width4]; + + % Background + resultsBkgrd = [slope, intercept]; + + elseif matches(p.Results.mode, 'TwoSpinOrbitDoubletPseudoVoigtShirley') + %% TwoSpinOrbitDoubletPseudoVoigtShirley + % ===== PSEUDO-VOIGT + SHIRLEY ===== + + % Extract parameters + area1 = currentParameterVector(1); + width1 = currentParameterVector(2); + pos1 = currentParameterVector(3); + area2 = currentParameterVector(1)*FitParameter.peakArea(1,2); + width2 = currentParameterVector(2); + pos2 = currentParameterVector(3)+FitParameter.peakPos(1,2); + area3 = currentParameterVector(4); + width3 = currentParameterVector(2); + pos3 = currentParameterVector(3)+currentParameterVector(5); + area4 = currentParameterVector(4)*FitParameter.peakArea(1,2); + width4 = currentParameterVector(2); + pos4 = currentParameterVector(3)+currentParameterVector(5)++FitParameter.peakPos(1,2); + + % Results: areas, positions, FWHM + resultsPeakArea = [area1, area2, area3, area4]; + resultsPeakPos = [pos1, pos2, pos3, pos4]; + resultsPeakFWHM = [width1, width2, width3, width4]; + + % Eta + resultsPeakEta = [currentParameterVector(6), currentParameterVector(6), currentParameterVector(6), currentParameterVector(6)]; + % Kappa + resultsPeakShirleyK = [currentParameterVector(7), currentParameterVector(7), currentParameterVector(7), currentParameterVector(7)]; + % background (Shirley intercept only) + resultsBkgrd = [0,currentParameterVector(end)]; + + elseif matches(p.Results.mode,'gasAndTwoAcidBasePairs') + %% gasAndTwoAcidBasePairs + % ===== GAUSSIAN + LINEAR ===== + + % Extract parameters + area1 = currentParameterVector(1); + width1 = currentParameterVector(2); + pos1 = currentParameterVector(3)+currentParameterVector(14); + + area2 = currentParameterVector(4); + width2 = currentParameterVector(2); + pos2 = currentParameterVector(3); + + area3 = currentParameterVector(5); + width3 = currentParameterVector(6); + pos3 = currentParameterVector(7); + + area4 = currentParameterVector(8); + width4 = currentParameterVector(2); + pos4 = currentParameterVector(9)+currentParameterVector(11); + + area5 = currentParameterVector(10); + width5 = currentParameterVector(2); + pos5 = currentParameterVector(11); + + slope = currentParameterVector(12); + intercept = currentParameterVector(13); + + % Results: areas, positions, FWHM + resultsPeakArea = [area1, area2, area3, area4, area5]; + resultsPeakPos = [pos1, pos2, pos3, pos4, pos5]; + resultsPeakFWHM = [width1, width2, width3, width4, width5]; + + % Background + resultsBkgrd = [slope, intercept]; + + elseif matches(p.Results.mode,'widthPositionAreaRatio') + %% widthPositionAreaRatio + % ===== GAUSSIAN + LINEAR ===== + + % Extract parameters + area1 = currentParameterVector(1); + width1 = currentParameterVector(2); + pos1 = currentParameterVector(3); + + area2 = currentParameterVector(1)*FitParameter.peakArea(1,2); + width2 = currentParameterVector(2); + pos2 = currentParameterVector(3)+FitParameter.peakPos(1,2); + + area3 = currentParameterVector(1)*FitParameter.peakArea(1,3); + width3 = currentParameterVector(2); + pos3 = currentParameterVector(3)+FitParameter.peakPos(1,3); + + area4 = currentParameterVector(1)*FitParameter.peakArea(1,4); + width4 = currentParameterVector(2); + pos4 = currentParameterVector(3)+FitParameter.peakPos(1,4); + + area5 = currentParameterVector(1)*FitParameter.peakArea(1,5); + width5 = currentParameterVector(2); + pos5 = currentParameterVector(3)+FitParameter.peakPos(1,5); + + slope = currentParameterVector(4); + intercept = currentParameterVector(5); + + % Results: areas, positions, FWHM + resultsPeakArea = [area1, area2, area3, area4, area5]; + resultsPeakPos = [pos1, pos2, pos3, pos4, pos5]; + resultsPeakFWHM = [width1, width2, width3, width4, width5]; + + % Background + resultsBkgrd = [slope, intercept]; + + elseif matches(p.Results.mode, 'widthPositionAmineAreaRatio') + %% MODE: WIDTH-POSITION-AMINE-AREA-RATIO + % ===== GAUSSIAN + LINEAR ===== + + area1 = currentParameterVector(1); + width_shared = currentParameterVector(2); + pos1 = currentParameterVector(3); + + resultsPeakPos = pos1; + resultsPeakArea = area1; + + % Remaining peaks: area2, area4, area5 from parameters + % Peak 3 = area2 * ratio + if isfield(FitParameter, 'peakEta') + % With eta: area1, width, pos1, area2, eta values for each peak, area4, area5, slope, intercept + area2 = currentParameterVector(4); + area4 = currentParameterVector(5); + area5 = currentParameterVector(6); + area3 = area2 * FitParameter.peakArea(1, 3); + + resultsPeakArea = [resultsPeakArea, area2, area3, area4, area5]; + + for iPeak = 2:5 + resultsPeakPos = [resultsPeakPos, pos1 + FitParameter.peakPos(1, iPeak)]; + end + + resultsPeakFWHM = ones(1, 5) * width_shared ; + resultsBkgrd = [currentParameterVector(7), currentParameterVector(8)]; + else + area2 = currentParameterVector(4); + area4 = currentParameterVector(5); + area5 = currentParameterVector(6); + area3 = area2 * FitParameter.peakArea(1, 3); + + resultsPeakArea = [resultsPeakArea, area2, area3, area4, area5]; + + for iPeak = 2:5 + resultsPeakPos = [resultsPeakPos, pos1 + FitParameter.peakPos(1, iPeak)]; + end + + resultsPeakFWHM = ones(1, 5) * width_shared ; + resultsBkgrd = [currentParameterVector(7), currentParameterVector(8)]; + end + end - existColumn = strcmp('peakFWHM',DataTable.Properties.VariableNames); - if ~any(existColumn) - DataTable.peakFWHM = NaN(size(DataTable,1),size(resultsPeakFWHW,2)); +end + +%% Write Results to table + +% Ensure all columns exist with appropriate sizing +if ~ismember('peakArea', DataTable.Properties.VariableNames) + DataTable.peakArea = NaN(height(DataTable), length(resultsPeakArea)); +else + oldWidth = size(DataTable.peakArea, 2); + newWidth = length(resultsPeakArea); + if newWidth > oldWidth + DataTable = resizeColumnHelper(DataTable, 'peakArea', newWidth); end - existColumn = strcmp('peakBackground',DataTable.Properties.VariableNames); - if ~any(existColumn) - DataTable.peakBackground = NaN(size(DataTable,1),size(resultsBkgrd,2)); +end +if ~ismember('peakPosition', DataTable.Properties.VariableNames) + DataTable.peakPosition = NaN(height(DataTable), length(resultsPeakPos)); +else + oldWidth = size(DataTable.peakPosition, 2); + newWidth = length(resultsPeakPos); + if newWidth > oldWidth + DataTable = resizeColumnHelper(DataTable, 'peakPosition', newWidth); end - %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 +end +if ~ismember('peakFWHM', DataTable.Properties.VariableNames) + DataTable.peakFWHM = NaN(height(DataTable), length(resultsPeakFWHM)); +else + oldWidth = size(DataTable.peakFWHM, 2); + newWidth = length(resultsPeakFWHM); + if newWidth > oldWidth + DataTable = resizeColumnHelper(DataTable, 'peakFWHM', newWidth); + end +end +if ~ismember('peakBackground', DataTable.Properties.VariableNames) + DataTable.peakBackground = NaN(height(DataTable), length(resultsBkgrd)); +else + oldWidth = size(DataTable.peakBackground, 2); + newWidth = length(resultsBkgrd); + if newWidth > oldWidth + DataTable = resizeColumnHelper(DataTable, 'peakBackground', newWidth); + end +end +% Handle optional peakEta (write if populated, clear to NaN if not) +if ~ismember('peakEta', DataTable.Properties.VariableNames) + DataTable.peakEta = NaN(height(DataTable), nPeaks); +else + oldWidth = size(DataTable.peakEta, 2); + newWidth = max(nPeaks, length(resultsPeakEta)); + if newWidth > oldWidth + DataTable = resizeColumnHelper(DataTable, 'peakEta', newWidth); + end +end +% Handle optional peakShirleyK (write if populated, clear to NaN if not) +if ~ismember('peakShirleyK', DataTable.Properties.VariableNames) + DataTable.peakShirleyK = NaN(height(DataTable), nPeaks); +else + oldWidth = size(DataTable.peakShirleyK, 2); + newWidth = max(nPeaks, length(resultsPeakShirleyK)); + if newWidth > oldWidth + DataTable = resizeColumnHelper(DataTable, 'peakShirleyK', newWidth); + end +end +% Pad results to match column width and write results +if ~isempty(resultsPeakArea) + resultsPeakArea = [resultsPeakArea, NaN(1, max(0, size(DataTable.peakArea, 2) - length(resultsPeakArea)))]; + DataTable.peakArea(selectedRowIdx, :) = resultsPeakArea; +end +if ~isempty(resultsPeakPos) + resultsPeakPos = [resultsPeakPos, NaN(1, max(0, size(DataTable.peakPosition, 2) - length(resultsPeakPos)))]; + DataTable.peakPosition(selectedRowIdx, :) = resultsPeakPos; +end +if ~isempty(resultsPeakFWHM) + resultsPeakFWHM = [resultsPeakFWHM, NaN(1, max(0, size(DataTable.peakFWHM, 2) - length(resultsPeakFWHM)))]; + DataTable.peakFWHM(selectedRowIdx, :) = resultsPeakFWHM; +end +if ~isempty(resultsBkgrd) + resultsBkgrd = [resultsBkgrd, NaN(1, max(0, size(DataTable.peakBackground, 2) - length(resultsBkgrd)))]; + DataTable.peakBackground(selectedRowIdx, :) = resultsBkgrd; +end +if ~isempty(resultsPeakEta) + % Pad eta to match column width + resultsPeakEta = [resultsPeakEta, NaN(1, max(0, size(DataTable.peakEta, 2) - length(resultsPeakEta)))]; + DataTable.peakEta(selectedRowIdx, :) = resultsPeakEta; +else + DataTable.peakEta(selectedRowIdx, :) = NaN; +end +if ~isempty(resultsPeakShirleyK) + % Pad K to match column width + resultsPeakShirleyK = [resultsPeakShirleyK, NaN(1, max(0, size(DataTable.peakShirleyK, 2) - length(resultsPeakShirleyK)))]; + DataTable.peakShirleyK(selectedRowIdx, :) = resultsPeakShirleyK; +else + DataTable.peakShirleyK(selectedRowIdx, :) = NaN; +end + +%% Print Results +if matches(p.Results.showFigures,'on') + disp(DataTable(p.Results.tableSelection,["name","peakPosition","peakArea","peakFWHM"])) + disp(DataTable(p.Results.tableSelection,["name","peakBackground"])) try - disp('Could not write results to table. Writing only 4 entries') - 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(resultsPeakFWHW,2)); - end - %if contains(p.Results.energyScale,'binding') - DataTable(p.Results.tableSelection,:).peakPosition = resultsPeakPos(:,1:4); - %end - DataTable(p.Results.tableSelection,:).peakFWHM = resultsPeakFWHW(:,1:4); - DataTable(p.Results.tableSelection,:).peakArea = resultsPeakArea(:,1:4); - DataTable(p.Results.tableSelection,:).peakBackground = resultsBkgrd(:,:); + disp(DataTable(p.Results.tableSelection,["name","peakEta",])) + catch + end + try + disp(DataTable(p.Results.tableSelection,["name","peakShirleyK"])) catch - disp('Could not write results to table. Writing only 2 entries') - 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(resultsPeakFWHW,2)); - end - %if contains(p.Results.energyScale,'binding') - DataTable(p.Results.tableSelection,:).peakPosition = resultsPeakPos(:,1:2); - %end - DataTable(p.Results.tableSelection,:).peakFWHM = resultsPeakFWHW(:,1:2); - DataTable(p.Results.tableSelection,:).peakArea = resultsPeakArea(:,1:2); - DataTable(p.Results.tableSelection,:).peakBackground = resultsBkgrd(:,1:2); end end +%% Plot Results +% +% This uses nappPlot and then moves the figure in a tiledlayout to also +% plot the residuals in a second tile. + +if matches(p.Results.showFigures,'on') + close all + + % Let nappPlot + nappAddFitToPlot do their thing + nappPlot(DataTable,p.Results.tableSelection, ... + 'tableVariable',p.Results.tableVariable, ... + 'energyScale',p.Results.energyScale); + axOld = gca; % axes created by nappPlot + lines = findobj(axOld,'-regexp','Type','line'); + set(lines,'LineStyle','none','Marker','o','MarkerSize',5); + nappAddFitToPlot(DataTable,p.Results.tableSelection, ... + 'energyScale',p.Results.energyScale); + % Create a NEW figure with tiledlayout(3,1) + fNew = figure; + t = tiledlayout(fNew,3,1); + % Tile 1: residuals (top tile) + ax1 = nexttile; % tile 1 + optimalModel = modelEvaluationFunction(currentParameterVector, fitEnergyValues); + residuals = fitIntensityValues - optimalModel; + plot(fitEnergyValues, residuals, 'o-'); + ylabel('residuals') + % Copy original axes and make it span the lower 2 tiles + axCopy = copyobj(axOld,t); + axCopy.Layout.Tile = 2; + axCopy.Layout.TileSpan = [2 1]; + % link axes + linkaxes([ax1 axCopy],'x'); + ax1.XDir = axCopy.XDir; + % Close the original nappPlot figure + close(ancestor(axOld,'figure')); +end + +end + + +%% MATHEMATICAL DEFINITIONS: Gaussian and Lorentzian Width Parameters +% +% GAUSSIAN PEAK (https://www.originlab.com/doc/Origin-Help/Gauss-FitFunc): +% G = A / (w * sqrt(pi/2)) * exp(-2 * ((E - E0) / w)^2) +% +% The width parameter 'w' in the Gaussian definition above is related to +% sigma (σ), the standard deviation of the gaussian, by: +% w = 2 * σ +% +% The Full Width at Half Maximum (FWHM) of the Gaussian, which is +% frequently reported in XPS, is: +% FWHM_Gauss = w * sqrt(log(4)) ≈ 1.17741 * w ≈ 2.35482 * σ +% +% +% LORENTZIAN PEAK: +% L = (A/pi) * (0.5*w) / ((E - E0)^2 + (0.5*w)^2) +% +% For the Lorentzian, the width parameter 'w' is directly related to the FWHM: +% FWHM_Lorentz = w (exactly!) +% The gamma parameter (half-width at half-maximum) is: gamma = 0.5*w +% +% +% RELATIONSHIP BETWEEN GAUSSIAN AND LORENTZIAN WIDTH: +% +% When using the SAME width parameter 'w' for both Gaussian and Lorentzian components: +% FWHM_Lorentz / FWHM_Gauss = w / (w * sqrt(log(4))) = 1 / sqrt(log(4)) ≈ 0.8493 +% +% FITTING APPROACH: Use a common FWHM parameter for both components. +% During fitting, FWHM is the fit parameter. Internally converted to: +% w_Gauss = FWHM / sqrt(log(4)) ≈ FWHM * 0.8493 +% w_Lorentz = FWHM (no conversion) +% +% RESULT: Both components have the SAME FWHM. FWHM_peak = FWHM (constant, independent of eta) + + +function w = fwhm2w(FWHM) +% Convert Gaussian FWHM to width parameter w +w = FWHM ./ sqrt(log(4)); +end + +function G = gaussPeak(x, area, FWHM, pos) +% Gaussian using width w (not sigma) +G = area ./ (fwhm2w(FWHM) .* sqrt(pi/2)) .* exp(-2 * ((x - pos)./fwhm2w(FWHM)).^2); +end + +function L = lorentzPeak(x, area, FWHM, pos) +L = area/pi .* (0.5*FWHM) ./ ((x - pos).^2 + (0.5*FWHM).^2); +end + +function peak = pseudoVoigtPeak(x, A, FWHM, pos, eta) +G = gaussPeak(x, A, FWHM, pos); +L = lorentzPeak(x, A, FWHM, pos); +peak = eta .* G + (1 - eta) .* L; +end + +function bkg = linearBackground(x, slope, intercept) +% Simple linear background +bkg = slope .* x + intercept; +end + +function B = activeShirleyBackgroundFwhm(E, peak_params,FitParameter,p,nPeaks) +% activeShirleyBackgroundFwhm - Compute Shirley background +% +% INPUT: +% E - Energy values (column vector) +% peak_params - Matrix with rows: [area, FWHM, pos, eta, K] for each peak +% +% OUTPUT: +% B - Shirley background curve (offset-free, added to model) + + +peakEta_exists = isfield(FitParameter,'peakEta'); +% totalParams = length(peak_params); +% nPeaks = (totalParams - 1) / 5; +B = zeros(size(E)); + +for iPeak = 1:nPeaks + if matches(p.Results.mode,"unconstrained") && ~peakEta_exists + area = peak_params(3*iPeak - 2); + FWHM = peak_params(3*iPeak - 1); + pos = peak_params(3*iPeak); + eta = 1; % Gauss only + K = peak_params(3*nPeaks +iPeak); + elseif matches(p.Results.mode,"unconstrained") && peakEta_exists + area = peak_params(4*iPeak - 3); + FWHM = peak_params(4*iPeak - 2); + pos = peak_params(4*iPeak - 1); + eta = peak_params(4*iPeak); + K = peak_params(4*nPeaks +iPeak); + elseif matches(p.Results.mode,"position") && ~peakEta_exists + area = peak_params(2*iPeak - 1); + FWHM = peak_params(2*iPeak); + pos1 = peak_params(2*nPeaks + 1); + if iPeak == 1 + pos = pos1; + else + relPos = FitParameter.peakPos(1, iPeak); + pos = pos1 + relPos; + end + eta = 1; % Gauss only + K = peak_params(2*nPeaks +1 +iPeak); + elseif matches(p.Results.mode,"position") && peakEta_exists + pos1 = peak_params(3*nPeaks + 1); + area = peak_params(3*iPeak - 2); + FWHM = peak_params(3*iPeak - 1); + eta = peak_params(3*iPeak); + if iPeak == 1 + pos = pos1; + else + relPos = FitParameter.peakPos(1, iPeak); + pos = pos1 + relPos; + end + K = peak_params(3*nPeaks +1 +iPeak); + + elseif matches(p.Results.mode,"widthPosition") && ~peakEta_exists + area = peak_params(iPeak); + FWHM = peak_params(nPeaks + 1); + pos1 = peak_params(nPeaks + 2); + if iPeak == 1 + pos = pos1; + else + relPos = FitParameter.peakPos(1, iPeak); + pos = pos1 + relPos; + end + eta = 1; % Gauss only + K = peak_params(nPeaks + 2 +iPeak); + + elseif matches(p.Results.mode,"widthPosition") && peakEta_exists + area = peak_params(2*iPeak - 1); + eta = peak_params(2*iPeak); + FWHM = peak_params(2*nPeaks + 1); + pos1 = peak_params(2*nPeaks + 2); + + if iPeak == 1 + pos = pos1; + else + relPos = FitParameter.peakPos(1, iPeak); + pos = pos1 + relPos; + end + K = peak_params(2*nPeaks +2 +iPeak); + + elseif matches(p.Results.mode,"width") && ~peakEta_exists + area = peak_params(2*iPeak - 1); + pos = peak_params(2*iPeak); + FWHM = peak_params(2*nPeaks + 1); + eta = 1; % Gauss only + K = peak_params(2*nPeaks +1 +iPeak); + + elseif matches(p.Results.mode,"width") && peakEta_exists + area = peak_params(3*iPeak - 2); + pos = peak_params(3*iPeak - 1); + eta = peak_params(3*iPeak); + FWHM = peak_params(3*nPeaks + 1); + K = peak_params(3*nPeaks +1 +iPeak); + + elseif matches(p.Results.mode,"positionShirley") && ~peakEta_exists + area = peak_params(2*iPeak - 1); + FWHM = peak_params(2*iPeak); + pos1 = peak_params(2*nPeaks + 1); + if iPeak == 1 + pos = pos1; + else + relPos = FitParameter.peakPos(1, iPeak); + pos = pos1 + relPos; + end + eta = 1; % Gauss only + K = peak_params(2*nPeaks +2); + + elseif matches(p.Results.mode,"positionShirley") && peakEta_exists + pos1 = peak_params(3*nPeaks + 1); + area = peak_params(3*iPeak - 2); + FWHM = peak_params(3*iPeak - 1); + eta = peak_params(3*iPeak); + if iPeak == 1 + pos = pos1; + else + relPos = FitParameter.peakPos(1, iPeak); + pos = pos1 + relPos; + end + K = peak_params(3*nPeaks +2); + + elseif matches(p.Results.mode,"widthPositionShirley") && ~peakEta_exists + area = peak_params(iPeak); + FWHM = peak_params(nPeaks + 1); + pos1 = peak_params(nPeaks + 2); + if iPeak == 1 + pos = pos1; + else + relPos = FitParameter.peakPos(1, iPeak); + pos = pos1 + relPos; + end + eta = 1; % Gauss only + K = peak_params(nPeaks + 3); + + elseif matches(p.Results.mode,"widthPositionShirley") && peakEta_exists + area = peak_params(2*iPeak - 1); + eta = peak_params(2*iPeak); + FWHM = peak_params(2*nPeaks + 1); + pos1 = peak_params(2*nPeaks + 2); + + if iPeak == 1 + pos = pos1; + else + relPos = FitParameter.peakPos(1, iPeak); + pos = pos1 + relPos; + end + K = peak_params(2*nPeaks +3); + + elseif matches(p.Results.mode,"widthShirley") && ~peakEta_exists + area = peak_params(2*iPeak - 1); + pos = peak_params(2*iPeak); + FWHM = peak_params(2*nPeaks + 1); + eta = 1; % Gauss only + K = peak_params(2*nPeaks +2); + + elseif matches(p.Results.mode,"widthShirley") && peakEta_exists + area = peak_params(3*iPeak - 2); + pos = peak_params(3*iPeak - 1); + eta = peak_params(3*iPeak); + FWHM = peak_params(3*nPeaks + 1); + K = peak_params(3*nPeaks +2); + + elseif matches(p.Results.mode,"TwoSpinOrbitDoubletPseudoVoigtShirley") + if iPeak == 1 + area = peak_params(1); + FWHM = peak_params(2); + pos = peak_params(3); + eta = peak_params(6); + K = peak_params(7); + elseif iPeak == 2 + area = peak_params(1)*FitParameter.peakArea(1,2); + FWHM = peak_params(2); + pos = peak_params(3)+FitParameter.peakPos(1,2); + eta = peak_params(6); + K = peak_params(7); + elseif iPeak == 3 + area = peak_params(4); + FWHM = peak_params(2); + pos = peak_params(3)+peak_params(5); + eta = peak_params(6); + K = peak_params(7); + elseif iPeak == 4 + area = peak_params(4)*FitParameter.peakArea(1,2); + FWHM = peak_params(2); + pos = peak_params(3)+peak_params(5)+FitParameter.peakPos(1,2); + eta = peak_params(6); + K = peak_params(7); + end + + end + + + % Calculate sigma and gamma from FWHM + % These match the definitions used in your helper functions + w = fwhm2w(FWHM); % Convert FWHM to width parameter + sigma = w / 2; % Gaussian sigma + gamma = FWHM / 2; % Lorentzian gamma + + if matches(p.Results.energyScale, 'binding') + % Background on HIGH binding energy side (E >= E0, right) + idxregion = (E >= pos); + else % kinetic + % error('ERROR: Shirley works only in binding energy mode') + E = max(E)+100-E; + pos = max(E)+100-pos; + idxregion = (E >= pos); + end + + % Use complementary CDF + Gtail = 0.5 * (1 - erf((E(idxregion) - pos) / (sqrt(2) * sigma))); + Ltail = 0.5 - (1/pi) * atan((E(idxregion) - pos) / gamma); + tailfrac = eta * Gtail + (1 - eta) * Ltail; + shirleycontrib = K * area * (0.5 - tailfrac); + Btemp = zeros(size(E)); + Btemp(idxregion) = shirleycontrib; + B = B + Btemp; +end +if matches(p.Results.energyScale, 'kinetic') + %B = flip(B); +end end -%% ------------------ Functions --------------- + +%% Create Generic Models +% +% This needs to be a function to allow variable number of peaks. + +function y = genericModelGenerator(params, energyValues, FitParameter, p ) +x = energyValues(:,1); +y = zeros(size(x)); +nPeaks = size(FitParameter.peakArea, 2); +peakEta_exists = isfield(FitParameter,'peakEta'); + +% peaks +for iPeak = 1:nPeaks + if matches(p.Results.mode,"unconstrained") && ~peakEta_exists + area = params(3*iPeak - 2); + FWHM = params(3*iPeak - 1); + pos = params(3*iPeak); + y = y + gaussPeak(x, area, FWHM, pos); + + elseif matches(p.Results.mode,"unconstrained") && peakEta_exists + area = params(4*iPeak - 3); + FWHM = params(4*iPeak - 2); + pos = params(4*iPeak - 1); + eta = params(4*iPeak); + y = y + pseudoVoigtPeak(x, area, FWHM, pos, eta); + + elseif startsWith(p.Results.mode,"position") && ~peakEta_exists + pos1 = params(2*nPeaks + 1); + area = params(2*iPeak - 1); + FWHM = params(2*iPeak); + + if iPeak == 1 + pos = pos1; + else + relPos = FitParameter.peakPos(1, iPeak); + pos = pos1 + relPos; + end + y = y + gaussPeak(x, area, FWHM, pos); + + elseif startsWith(p.Results.mode,"position") && peakEta_exists + pos1 = params(3*nPeaks + 1); + area = params(3*iPeak - 2); + FWHM = params(3*iPeak - 1); + eta = params(3*iPeak); + if iPeak == 1 + pos = pos1; + else + relPos = FitParameter.peakPos(1, iPeak); + pos = pos1 + relPos; + end + y = y + pseudoVoigtPeak(x, area, FWHM, pos, eta); + + elseif startsWith(p.Results.mode,"widthPosition") && ~peakEta_exists + pos1 = params(nPeaks + 2); + FWHM = params(nPeaks + 1); + area = params(iPeak); + + if iPeak == 1 + pos = pos1; + else + relPos = FitParameter.peakPos(1, iPeak); + pos = pos1 + relPos; + + end + y = y + gaussPeak(x, area, FWHM, pos); + + elseif startsWith(p.Results.mode,"widthPosition") && peakEta_exists + pos1 = params(2*nPeaks + 2); + FWHM = params(2*nPeaks + 1); + area = params(2*iPeak-1); + eta = params(2*iPeak); + if iPeak == 1 + pos = pos1; + else + relPos = FitParameter.peakPos(1, iPeak); + pos = pos1 + relPos; + end + y = y + pseudoVoigtPeak(x, area, FWHM, pos, eta); + + + elseif startsWith(p.Results.mode,"width") && ~peakEta_exists + FWHM = params(2*nPeaks + 1); + area = params(2*iPeak - 1); + pos = params(2*iPeak); + + y = y + gaussPeak(x, area, FWHM, pos); + + elseif startsWith(p.Results.mode,"width") && peakEta_exists + FWHM = params(3*nPeaks + 1); + area = params(3*iPeak - 2); + pos = params(3*iPeak-1); + eta = params(3*iPeak); + + y = y + pseudoVoigtPeak(x, area, FWHM, pos, eta); + end + +end + +% background +if matches(p.Results.backgroundMode, 'linear') + slope = params(end-1); + intercept = params(end); + y = y + linearBackground(x, slope, intercept); +elseif matches(p.Results.backgroundMode, 'shirley') + C0 = params(end); + y = y + activeShirleyBackgroundFwhm(x, params,FitParameter,p,nPeaks)+ C0; +end + + +end + + +%% HELPER FUNCTION +% Call this function to safely resize a column while preserving all existing data +function DataTable = resizeColumnHelper(DataTable, colName, newWidth) +oldWidth = size(DataTable.(colName), 2); +if newWidth > oldWidth + % Create new column with NaN, width = newWidth + oldData = DataTable.(colName); + DataTable.(colName) = NaN(height(DataTable), newWidth); + % Restore old data in first oldWidth columns, other rows unchanged + DataTable.(colName)(:, 1:oldWidth) = oldData; +end +end function correctedStruct = fixBondToBound(inputStruct) % Recursively corrects all field names in a struct: 'Bond' -> 'Bound' @@ -944,220 +3282,76 @@ for i = 1:numel(fieldNames) end end -function y = gaussianModel(params, x, nPeaks) -% -% construct equation with variable number of Gaussian peaks and linear -% background -y = zeros(size(x(:,1))); +% Helper function: Check bounds consistency +function checkBoundsPair(FitParam, lowerName, upperName, fieldLabel) +if isfield(FitParam, lowerName) && isfield(FitParam, upperName) + lower = FitParam.(lowerName); + upper = FitParam.(upperName); - -% ------------------- 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; + violatingIdx = find(lower > upper); + if ~isempty(violatingIdx) + if numel(violatingIdx) <= 5 + error('%s bounds inconsistency: lower > upper at index(es) %s', ... + fieldLabel, num2str(violatingIdx)); + else + error('%s bounds inconsistency: lower > upper at %d indices (first: %d)', ... + fieldLabel, length(violatingIdx), violatingIdx(1)); 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)]; +% Helper function: Check field dimensions +function checkFieldDimension(FitParam, fieldName, expectedSize, modeStr, fieldLabel) +if ~isfield(FitParam, fieldName) + error('Mode ''%s'': Missing required field ''%s'' (%s).', ... + modeStr, fieldName, fieldLabel); 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)]; +fieldValue = FitParam.(fieldName); +% Scalar case (expectedSize == 0): allow non‑scalar, but warn and ignore extras +if expectedSize == 0 + sz = size(fieldValue); + if ~isscalar(fieldValue) + warning(['Mode ''%s'': Field ''%s'' is expected to be scalar, ' ... + 'but has size [%d, %d]. Only the first element will be used; ' ... + 'extra values are ignored. (%s)'], ... + modeStr, fieldName, sz(1), sz(2), fieldLabel); + % Reduce to scalar (first element, preserving NaN if that is first) + fieldValue = fieldValue(1); + FitParam.(fieldName) = fieldValue; + end + + % Optional warning if the (used) scalar is NaN + if isnan(fieldValue) + warning(['Mode ''%s'': Field ''%s'' is NaN. This is allowed, but may ' ... + 'lead to ineffective constraints. (%s)'], ... + modeStr, fieldName, fieldLabel); + end + return; 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)]; +% Non‑scalar: allow extra NaNs beyond expectedSize +if ~isvector(fieldValue) || size(fieldValue,1) ~= 1 + error('Mode ''%s'': Field ''%s'' must be a 1×N row vector, but got size [%d, %d]. (%s)', ... + modeStr, fieldName, size(fieldValue,1), size(fieldValue,2), fieldLabel); +end + +nanMask = isnan(fieldValue); +nNonNaN = sum(~nanMask); + +if nNonNaN ~= expectedSize + error(['Mode ''%s'': Field ''%s'' must contain %d non-NaN value(s) ' ... + '(one per peak), but found %d non-NaN value(s). (%s)'], ... + modeStr, fieldName, expectedSize, nNonNaN, fieldLabel); +end + +% Optional warning if vector longer than needed (extra NaNs or values) +if numel(fieldValue) > expectedSize + warning(['Mode ''%s'': Field ''%s'' has %d element(s), but only ' ... + '%d non-NaN value(s) are used for %d peak(s). ' ... + 'Extra entries are ignored (typically NaN). (%s)'], ... + modeStr, fieldName, numel(fieldValue), nNonNaN, expectedSize, fieldLabel); end end \ No newline at end of file diff --git a/beamtime/nappFitGaussXpsDevelopment.m b/beamtime/nappFitGaussXpsDevelopment.m deleted file mode 100644 index 24e5777..0000000 --- a/beamtime/nappFitGaussXpsDevelopment.m +++ /dev/null @@ -1,861 +0,0 @@ -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 \ No newline at end of file