From d72e8b4ddc8c3540d67bd9b7e4d5e6e991c59936 Mon Sep 17 00:00:00 2001 From: ThorBarRa Date: Fri, 8 May 2026 20:43:17 +0200 Subject: [PATCH] Added a model where Shirley and PV eta are shared for peaks with shared width and relative position. This is usefull for adv. Carbon. --- .../TestAndDevelopShirley.m | 65 + beamtime/nappFitGaussXps.m | 1845 +++++++++-------- 2 files changed, 1036 insertions(+), 874 deletions(-) diff --git a/beamtime/documentedExampleScripts/TestAndDevelopShirley.m b/beamtime/documentedExampleScripts/TestAndDevelopShirley.m index 93aecd0..7507c2e 100644 --- a/beamtime/documentedExampleScripts/TestAndDevelopShirley.m +++ b/beamtime/documentedExampleScripts/TestAndDevelopShirley.m @@ -1026,6 +1026,71 @@ FitList = nappFitGaussXps(FitList,tmpFitParameters,tmpDataSelection,... quickPlotBinding(FitList,tmpDataSelection) xlim([34,39]) +%% TEST 3.13: widthPosition kinetic PV Shirley & Eta Shared +% Pseudo-Voigt + Shirley, shared width, first peak position free, other peaks offset fixed, shared Shirley K, shared PV eta +clearvars tmpFitParameters; +clc +tmpFitParameters.peakArea = [5e5 100]; +tmpFitParameters.peakWidth = 1.6; +tmpFitParameters.peakPos = [113 -2.32]; +% Bounds for areas (lower and upper) +tmpFitParameters.peakAreaLowerBound = [1e3 1e3]; +tmpFitParameters.peakAreaUpperBound = [Inf Inf]; +% Bounds for widths +tmpFitParameters.peakWidthLowerBound = 1.0; +tmpFitParameters.peakWidthUpperBound = 2.0; +% Bounds for positions +tmpFitParameters.peakPosLowerBound = 112; +tmpFitParameters.peakPosUpperBound = 114; +tmpFitParameters.backgroundRange = [12 6]; +% PSEUDO-VOIGT: + tmpFitParameters.peakEta = 1; + tmpFitParameters.peakEtaLowerBound = 0.0; + tmpFitParameters.peakEtaUpperBound = 1.0; +% SHIRLEY variant: + tmpFitParameters.peakShirleyK = 0.1; + tmpFitParameters.peakShirleyKLowerBound = 0.001; + tmpFitParameters.peakShirleyKUpperBound = 10; + tmpFitParameters.peakShirleyOffset = 1E4; + +FitList = nappFitGaussXps(FitList,tmpFitParameters,tmpDataSelection,... + 'tableVariable','spectrum_countsNorm','energyScale','kinetic',... + 'mode','widthPositionShirleyEta','backgroundMode','shirley','showFigures','on'); +%% TEST 3.14: widthPosition binding PV Shirley & Eta Shared +% Pseudo-Voigt + Shirley, shared width, first peak position free, other +% peaks offset fixed, shared Shirley K and PV eta + +clearvars tmpFitParameters; +clc +tmpFitParameters.peakArea = [1e5 100]; +tmpFitParameters.peakWidth = 1.6; +tmpFitParameters.peakPos = [38 -1.8]; +% Bounds for areas (lower and upper) +tmpFitParameters.peakAreaLowerBound = [0.01 0.01]; +tmpFitParameters.peakAreaUpperBound = [Inf Inf]; +% Bounds for widths +tmpFitParameters.peakWidthLowerBound = 0.2; +tmpFitParameters.peakWidthUpperBound = 2.0; +% Bounds for positions +tmpFitParameters.peakPosLowerBound = 37; +tmpFitParameters.peakPosUpperBound = 39; +tmpFitParameters.backgroundRange = [12 6]; +% PSEUDO-VOIGT + tmpFitParameters.peakEta = 1; + tmpFitParameters.peakEtaLowerBound = 0; + tmpFitParameters.peakEtaUpperBound = 1.0; +% SHIRLEY variant: + tmpFitParameters.peakShirleyK = 0.01; + tmpFitParameters.peakShirleyKLowerBound = 0.001; + tmpFitParameters.peakShirleyKUpperBound = 100; + tmpFitParameters.peakShirleyOffset = 80; + +FitList = nappFitGaussXps(FitList,tmpFitParameters,tmpDataSelection,... + 'tableVariable','spectrum_counts','energyScale','binding',... + 'mode','widthPositionShirleyEta','backgroundMode','shirley','showFigures','off'); +quickPlotBinding(FitList,tmpDataSelection) +xlim([34,39]) + %% TEST 4.1: Width kinetic Gauss linear % Gaussian + Linear, shared width, individual areas and positions diff --git a/beamtime/nappFitGaussXps.m b/beamtime/nappFitGaussXps.m index efe194e..390947b 100644 --- a/beamtime/nappFitGaussXps.m +++ b/beamtime/nappFitGaussXps.m @@ -575,7 +575,7 @@ function DataTable = nappFitGaussXps(DataTable,varargin) defaultTableVariable = 'spectrum_countsPerSecond'; expectedTableVariable = {'spectrum_counts','spectrum_countsPerSecond','spectrum_countsNorm'}; defaultMode = 'widthPosition'; -expectedMode = {'unconstrained','widthPosition','position','width','widthPositionShirley', ... +expectedMode = {'unconstrained','widthPosition','position','width','widthPositionShirley', 'widthPositionShirleyEta',... 'positionShirley','widthShirley','unconstrainedAndWidthPosition','gasAndTwoAcidBasePairs',... 'widthPositionAreaRatio','widthPositionAmineAreaRatio','TwoSpinOrbitDoublet','TwoSpinOrbitDoubletPseudoVoigtShirley','GasAndSurfaceSpinOrbitDoublet','TwoUnconstrainedAndWidthPositionPseudoVoigtPeak'}; defaultTableSelection = true(size(DataTable,1),1); @@ -662,15 +662,15 @@ for iaz = 1 % % Validate core parameters required for ALL modes % fprintf(' Checking core parameters (required for all modes)...\n'); - % + % % 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'); - % - % + % + % % % 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 @@ -717,16 +717,16 @@ for iaz = 1 % % For custom modes, do basic checks % warning('Custom mode ''%s'' detected. Performing basic checks only.', p.Results.mode); % end - % + % % % Background model validation % fprintf(' Checking background model parameters for ''%s''...\n', p.Results.backgroundMode); - % + % % if contains(p.Results.mode, 'Shirley') % Conditional K validation: % shirleyKExpectedSize = 0; % scalar for shared modes % else % shirleyKExpectedSize = nPeaks; % array for per-peak modes % end - % + % % switch p.Results.backgroundMode % case 'shirley' % % Shirley background requires 4 fields @@ -745,25 +745,25 @@ for iaz = 1 % warning('backgroundMode is ''linear'' but peakShirleyOffset is present. It will be ignored.'); % end % end - % + % % % Pseudo-Voigt extension check % if isfield(FitParameter, 'peakEta') % fprintf(' Checking Pseudo-Voigt extension parameters...\n'); - % + % % % 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 - % + % % 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'); - % + % % % 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.'); + % warning('peakEta values should be in range [0, 1] (0=Lorentzian, 1=Gaussian). Found values outside this range.'); % end % end % fprintf('✓ FitParameter validation successful for %d peak(s) in ''%s'' mode with ''%s'' background.\n\n', ... @@ -1418,6 +1418,50 @@ for iuz = 1 end end + elseif matches(p.Results.mode, 'widthPositionShirleyEta') + %% position & width & Shirley & eta — shared/constrained + + nPeaks = size(FitParameter.peakArea, 2); + + if matches(p.Results.backgroundMode, 'shirley') + modelDefinitionFunction = @(currentParameterVector,energyValues) ... + genericModelGenerator(currentParameterVector,energyValues,FitParameter,p); + + % Areas (one per peak) + for iPeak = 1:nPeaks + initialParameterVector = [initialParameterVector, FitParameter.peakArea(1, iPeak)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakAreaLowerBound(1, iPeak)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakAreaUpperBound(1, iPeak)]; + end + + % Shared FWHM + initialParameterVector = [initialParameterVector, FitParameter.peakWidth(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakWidthLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakWidthUpperBound(1, 1)]; + + % Shared eta + initialParameterVector = [initialParameterVector, FitParameter.peakEta(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakEtaLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakEtaUpperBound(1, 1)]; + + % Peak 1 absolute position + initialParameterVector = [initialParameterVector, FitParameter.peakPos(1, 1)]; + lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakPosLowerBound(1, 1)]; + upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakPosUpperBound(1, 1)]; + + % Shared 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('widthPositionShirleyEta mode requires backgroundMode ''shirley''.'); + end + elseif matches(p.Results.mode,'width') %% width constrained @@ -1720,7 +1764,7 @@ for iuz = 1 lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakPosLowerBound(1, 1)]; upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakPosUpperBound(1, 1)]; -%% TwoUnconstrainedAndWidthPositionPseudoVoigtPeak + %% TwoUnconstrainedAndWidthPositionPseudoVoigtPeak elseif matches(p.Results.mode,'TwoUnconstrainedAndWidthPositionPseudoVoigtPeak') && ... matches(p.Results.backgroundMode,'shirley') modelDefinitionFunction = @(currentParameterVector, energyValues) ... @@ -1785,7 +1829,7 @@ for iuz = 1 upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakShirleyOffsetUpperBound]; - + elseif matches(p.Results.mode,'TwoSpinOrbitDoublet') %% TwoSpinOrbitDoublet @@ -2056,7 +2100,7 @@ for ii = 1:size(DataTable(p.Results.tableSelection,:),1) end %-------------------------------------------------------------------------- - %% Core fitting part, + %% Core fitting part, % the fitting is done here (kinetic energy scale) %-------------------------------------------------------------------------- clearvars currentParameterVector @@ -2069,18 +2113,18 @@ for ii = 1:size(DataTable(p.Results.tableSelection,:),1) % ------------------------------------------------------------------------- % PUBLISH MODE DETECTION - % If MATLAB is publishing a document, completely silence the optimization + % If MATLAB is publishing a document, completely silence the optimization % output and live plots so they don't flood the final PDF/HTML. stackInfo = dbstack; if any(contains({stackInfo.name}, 'publish')) || any(contains({stackInfo.name}, 'evalmxdom')) globalSearchOptimizer.Display = 'off'; % Turns off the scrolling text % globalSearchOptimizer.PlotFcns = []; % (Commented out: leaves the live fitting plot active) - + % Silence the local fmincon solver by attaching options directly to the optimizationProblem - optimizationProblem.options = optimoptions('fmincon', 'Display', 'off'); + optimizationProblem.options = optimoptions('fmincon', 'Display', 'off'); end % ------------------------------------------------------------------------- - + currentParameterVector = run(globalSearchOptimizer,optimizationProblem); elseif contains(p.Results.energyScale,'binding') %-------------------------------------------------------------------------- @@ -2105,7 +2149,7 @@ for ii = 1:size(DataTable(p.Results.tableSelection,:),1) fitIntensityValues = fitIntensityValues(tmp1:tmp2); fitEnergyValues = fitEnergyValues(tmp1:tmp2); end - %-------------------------------------------------------------------------- + %-------------------------------------------------------------------------- %% Core fitting part % the fitting is done here (binding energy scale) %-------------------------------------------------------------------------- @@ -2613,445 +2657,472 @@ for ii = 1:size(DataTable(p.Results.tableSelection,:),1) resultsBkgrd = [0,currentParameterVector(end)]; end - end - %% Custom Models - if matches(p.Results.mode, 'unconstrainedAndWidthPosition') && matches(p.Results.backgroundMode, 'linear') - %% unconstrainedAndWidthPosition - % ===== GAUSSIAN + LINEAR ===== + elseif matches(p.Results.mode,"widthPositionShirleyEta") + %% MODE: widthPositionShirleyEta - % 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); + % ===== PSEUDO-VOIGT + SHIRLEY (shared FWHM, shared eta, shared K) ===== - % Calculate absolute positions - pos1 = pos2_rel + pos1_shift; - pos2 = pos2_rel; - pos4 = pos4_offset; - pos5 = pos5_offset; + fwhm_shared = currentParameterVector(nPeaks + 1); + eta_shared = currentParameterVector(nPeaks + 2); + pos1_fit = currentParameterVector(nPeaks + 3); + K_shared = currentParameterVector(nPeaks + 4); - % Results: areas, positions, FWHM - resultsPeakArea = [area1, area2, area3, area4, area5]; - resultsPeakPos = [pos1, pos2, pos3, pos4, pos5]; - resultsPeakFWHM = [width1, width1, width3, width1, width1]; + for iPeak = 1:nPeaks + resultsPeakArea = [resultsPeakArea, currentParameterVector(iPeak)]; + resultsPeakEta = [resultsPeakEta, eta_shared]; + resultsPeakShirleyK = [resultsPeakShirleyK, K_shared]; - % Background - resultsBkgrd = [slope, intercept]; -elseif matches(p.Results.mode, 'TwoUnconstrainedAndWidthPositionPseudoVoigtPeak') && matches(p.Results.backgroundMode,'shirley') - %% TwoSpinOrbitDoubletPseudoVoigtShirley - % ===== PSEUDO-VOIGT + SHIRLEY ===== - - % Extract parameters - area1 = currentParameterVector(1); - width1 = currentParameterVector(2); - pos1 = currentParameterVector(3)+currentParameterVector(11); - area2 = currentParameterVector(4); - width2 = currentParameterVector(2); - pos2 = currentParameterVector(3); - area3 = currentParameterVector(5); - width3 = currentParameterVector(6); - pos3 = currentParameterVector(7); - area4 = currentParameterVector(8); - width4 = currentParameterVector(9); - pos4 = currentParameterVector(10); - - % Results: areas, positions, FWHM - resultsPeakArea = [area1, area2, area3, area4]; - resultsPeakPos = [pos1, pos2, pos3, pos4]; - resultsPeakFWHM = [width1, width2, width3, width4]; - - % Eta - resultsPeakEta = [currentParameterVector(13), currentParameterVector(13), currentParameterVector(14), currentParameterVector(15)]; - % Kappa - resultsPeakShirleyK = [currentParameterVector(16), currentParameterVector(16), currentParameterVector(17), currentParameterVector(18)]; - % background (Shirley intercept only) - resultsBkgrd = [0,currentParameterVector(end)]; - - 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, 'GasAndSurfaceSpinOrbitDoublet') - %% GasAndSurfaceSpinOrbitDoublet - % ===== 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(5); - pos3 = currentParameterVector(6); - area4 = currentParameterVector(4)*FitParameter.peakArea(1,4); - width4 = currentParameterVector(5); - pos4 = currentParameterVector(6)+FitParameter.peakPos(1,4); - slope = currentParameterVector(7); - intercept = currentParameterVector(8); - - % 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)]; + if iPeak == 1 + peakPos = pos1_fit; + else + peakPos = pos1_fit + 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)]; + resultsPeakPos = [resultsPeakPos, peakPos]; end + resultsPeakFWHM = ones(1, nPeaks) * fwhm_shared; + % background (Shirley intercept only) + resultsBkgrd = [0, currentParameterVector(end)]; + end + %% Custom Models - end - end - %% Write Results to table + if matches(p.Results.mode, 'unconstrainedAndWidthPosition') && matches(p.Results.backgroundMode, 'linear') + %% unconstrainedAndWidthPosition + % ===== GAUSSIAN + LINEAR ===== - % 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 - 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 - 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(p.Results.tableSelection,:).peakArea(ii, :) = resultsPeakArea; - end - if ~isempty(resultsPeakPos) - resultsPeakPos = [resultsPeakPos, NaN(1, max(0, size(DataTable.peakPosition, 2) - length(resultsPeakPos)))]; - DataTable(p.Results.tableSelection,:).peakPosition(ii, :) = resultsPeakPos; - end - if ~isempty(resultsPeakFWHM) - resultsPeakFWHM = [resultsPeakFWHM, NaN(1, max(0, size(DataTable.peakFWHM, 2) - length(resultsPeakFWHM)))]; - DataTable(p.Results.tableSelection,:).peakFWHM(ii, :) = resultsPeakFWHM; - end - if ~isempty(resultsBkgrd) - resultsBkgrd = [resultsBkgrd, NaN(1, max(0, size(DataTable.peakBackground, 2) - length(resultsBkgrd)))]; - DataTable(p.Results.tableSelection,:).peakBackground(ii, :) = resultsBkgrd; - end - if ~isempty(resultsPeakEta) - % Pad eta to match column width - resultsPeakEta = [resultsPeakEta, NaN(1, max(0, size(DataTable.peakEta, 2) - length(resultsPeakEta)))]; - DataTable(p.Results.tableSelection,:).peakEta(ii, :) = resultsPeakEta; - else - DataTable(p.Results.tableSelection,:).peakEta(ii, :) = NaN; - end - if ~isempty(resultsPeakShirleyK) - % Pad K to match column width - resultsPeakShirleyK = [resultsPeakShirleyK, NaN(1, max(0, size(DataTable.peakShirleyK, 2) - length(resultsPeakShirleyK)))]; - DataTable(p.Results.tableSelection,:).peakShirleyK(ii, :) = resultsPeakShirleyK; - else - DataTable(p.Results.tableSelection,:).peakShirleyK(ii, :) = NaN; - end + % 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); - - %% Plot Results - % - % This uses nappPlot and then moves the figure in a tiledlayout to also - % plot the residuals in a second tile. + % Calculate absolute positions + pos1 = pos2_rel + pos1_shift; + pos2 = pos2_rel; + pos4 = pos4_offset; + pos5 = pos5_offset; - if matches(p.Results.showFigures,'on') - - % define data to plot - TmpData2Plot = DataTable(p.Results.tableSelection,:); - % Let nappPlot + nappAddFitToPlot do their thing - nappPlot(TmpData2Plot(ii,:), ... - '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); - title(TmpData2Plot(ii,:).name); - nappAddFitToPlot(TmpData2Plot(ii,:), ... - '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 - - % --- Define variables for the Universal Logging Block --- - numRowsTable = height(DataTable); - processIdx = find(p.Results.tableSelection); - % ------------------------------------------------------------------- + % Results: areas, positions, FWHM + resultsPeakArea = [area1, area2, area3, area4, area5]; + resultsPeakPos = [pos1, pos2, pos3, pos4, pos5]; + resultsPeakFWHM = [width1, width1, width3, width1, width1]; - - % 1. Dynamically get this function's name - thisFuncName = mfilename; - colDate = sprintf('log%sDate', thisFuncName); - colCaller = sprintf('log%sCaller', thisFuncName); - - - % 2. Probe for column existence, create empty string arrays if missing - if ~ismember(colDate, DataTable.Properties.VariableNames) - DataTable.(colDate) = repmat("", numRowsTable, 1); + % Background + resultsBkgrd = [slope, intercept]; + elseif matches(p.Results.mode, 'TwoUnconstrainedAndWidthPositionPseudoVoigtPeak') && matches(p.Results.backgroundMode,'shirley') + %% TwoSpinOrbitDoubletPseudoVoigtShirley + % ===== PSEUDO-VOIGT + SHIRLEY ===== + + % Extract parameters + area1 = currentParameterVector(1); + width1 = currentParameterVector(2); + pos1 = currentParameterVector(3)+currentParameterVector(11); + area2 = currentParameterVector(4); + width2 = currentParameterVector(2); + pos2 = currentParameterVector(3); + area3 = currentParameterVector(5); + width3 = currentParameterVector(6); + pos3 = currentParameterVector(7); + area4 = currentParameterVector(8); + width4 = currentParameterVector(9); + pos4 = currentParameterVector(10); + + % Results: areas, positions, FWHM + resultsPeakArea = [area1, area2, area3, area4]; + resultsPeakPos = [pos1, pos2, pos3, pos4]; + resultsPeakFWHM = [width1, width2, width3, width4]; + + % Eta + resultsPeakEta = [currentParameterVector(13), currentParameterVector(13), currentParameterVector(14), currentParameterVector(15)]; + % Kappa + resultsPeakShirleyK = [currentParameterVector(16), currentParameterVector(16), currentParameterVector(17), currentParameterVector(18)]; + % background (Shirley intercept only) + resultsBkgrd = [0,currentParameterVector(end)]; + + 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, 'GasAndSurfaceSpinOrbitDoublet') + %% GasAndSurfaceSpinOrbitDoublet + % ===== 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(5); + pos3 = currentParameterVector(6); + area4 = currentParameterVector(4)*FitParameter.peakArea(1,4); + width4 = currentParameterVector(5); + pos4 = currentParameterVector(6)+FitParameter.peakPos(1,4); + slope = currentParameterVector(7); + intercept = currentParameterVector(8); + + % 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 + 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 + 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 + 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(p.Results.tableSelection,:).peakArea(ii, :) = resultsPeakArea; + end + if ~isempty(resultsPeakPos) + resultsPeakPos = [resultsPeakPos, NaN(1, max(0, size(DataTable.peakPosition, 2) - length(resultsPeakPos)))]; + DataTable(p.Results.tableSelection,:).peakPosition(ii, :) = resultsPeakPos; + end + if ~isempty(resultsPeakFWHM) + resultsPeakFWHM = [resultsPeakFWHM, NaN(1, max(0, size(DataTable.peakFWHM, 2) - length(resultsPeakFWHM)))]; + DataTable(p.Results.tableSelection,:).peakFWHM(ii, :) = resultsPeakFWHM; + end + if ~isempty(resultsBkgrd) + resultsBkgrd = [resultsBkgrd, NaN(1, max(0, size(DataTable.peakBackground, 2) - length(resultsBkgrd)))]; + DataTable(p.Results.tableSelection,:).peakBackground(ii, :) = resultsBkgrd; + end + if ~isempty(resultsPeakEta) + % Pad eta to match column width + resultsPeakEta = [resultsPeakEta, NaN(1, max(0, size(DataTable.peakEta, 2) - length(resultsPeakEta)))]; + DataTable(p.Results.tableSelection,:).peakEta(ii, :) = resultsPeakEta; + else + DataTable(p.Results.tableSelection,:).peakEta(ii, :) = NaN; + end + if ~isempty(resultsPeakShirleyK) + % Pad K to match column width + resultsPeakShirleyK = [resultsPeakShirleyK, NaN(1, max(0, size(DataTable.peakShirleyK, 2) - length(resultsPeakShirleyK)))]; + DataTable(p.Results.tableSelection,:).peakShirleyK(ii, :) = resultsPeakShirleyK; + else + DataTable(p.Results.tableSelection,:).peakShirleyK(ii, :) = NaN; + 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') + + % define data to plot + TmpData2Plot = DataTable(p.Results.tableSelection,:); + % Let nappPlot + nappAddFitToPlot do their thing + nappPlot(TmpData2Plot(ii,:), ... + '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); + title(TmpData2Plot(ii,:).name); + nappAddFitToPlot(TmpData2Plot(ii,:), ... + '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 + + % --- Define variables for the Universal Logging Block --- + numRowsTable = height(DataTable); + processIdx = find(p.Results.tableSelection); + % ------------------------------------------------------------------- + + + % 1. Dynamically get this function's name + thisFuncName = mfilename; + colDate = sprintf('log%sDate', thisFuncName); + colCaller = sprintf('log%sCaller', thisFuncName); + + + % 2. Probe for column existence, create empty string arrays if missing + if ~ismember(colDate, DataTable.Properties.VariableNames) + DataTable.(colDate) = repmat("", numRowsTable, 1); + end + if ~ismember(colCaller, DataTable.Properties.VariableNames) + DataTable.(colCaller) = repmat("", numRowsTable, 1); + end + + % 3. Get current time as a string + currentTimeStr = string(datetime('now', 'Format', 'yyyy-MM-dd HH:mm:ss')); + + % 4. Determine the calling script using the call stack + stackInfo = dbstack; + if length(stackInfo) > 1 + callerScript = string(stackInfo(2).file); % The script that called this function + else + callerScript = " "; % Called directly from the command window + end + + % 5. Write the log data only to the rows being processed + DataTable.(colDate)(processIdx) = currentTimeStr; + DataTable.(colCaller)(processIdx) = callerScript; end - if ~ismember(colCaller, DataTable.Properties.VariableNames) - DataTable.(colCaller) = repmat("", numRowsTable, 1); - end - - % 3. Get current time as a string - currentTimeStr = string(datetime('now', 'Format', 'yyyy-MM-dd HH:mm:ss')); - - % 4. Determine the calling script using the call stack - stackInfo = dbstack; - if length(stackInfo) > 1 - callerScript = string(stackInfo(2).file); % The script that called this function - else - callerScript = " "; % Called directly from the command window - end - - % 5. Write the log data only to the rows being processed - DataTable.(colDate)(processIdx) = currentTimeStr; - DataTable.(colCaller)(processIdx) = callerScript; -end - %% Print Results + %% Print Results if matches(p.Results.showFigures,'on') disp(DataTable(p.Results.tableSelection,["name","peakPosition"])) disp(DataTable(p.Results.tableSelection,["name","peakArea"])) @@ -3110,501 +3181,527 @@ end % 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 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 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 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 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 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) + 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)); + 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); + 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(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); + 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); + 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(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,"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 + 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,"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(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); + 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); + 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(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,"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); + elseif matches(p.Results.mode,"widthPositionShirleyEta") + area = peak_params(iPeak); + FWHM = peak_params(nPeaks + 1); + eta = peak_params(nPeaks + 2); + pos1 = peak_params(nPeaks + 3); + if iPeak == 1 + pos = pos1; + else + relPos = FitParameter.peakPos(1, iPeak); + pos = pos1 + relPos; + end + K = peak_params(nPeaks + 4); + + 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 + + elseif matches(p.Results.mode,"TwoUnconstrainedAndWidthPositionPseudoVoigtPeak") && peakEta_exists + if iPeak == 1 + area = peak_params(1); + FWHM = peak_params(2); + pos = peak_params(3)+ peak_params(11); + eta = peak_params(13); + K = peak_params(16); + elseif iPeak == 2 + area = peak_params(4); + FWHM = peak_params(2); + pos = peak_params(3); + eta = peak_params(13); + K = peak_params(16); + elseif iPeak == 3 + area = peak_params(5); + FWHM = peak_params(6); + pos = peak_params(7); + eta = peak_params(14); + K = peak_params(17); + elseif iPeak == 4 + area = peak_params(8); + FWHM = peak_params(9); + pos = peak_params(10); + eta = peak_params(15); + K = peak_params(18); + 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 + + % Calculate the continuous CDF over the ENTIRE energy range + % This gives the solver a smooth gradient. + G_cdf = 0.5 * (1 + erf((E - pos) / (sqrt(2) * sigma))); + L_cdf = 0.5 + (1/pi) * atan((E - pos) / gamma); + + % Combined Pseudo-Voigt integral fraction (0 to 1) + integral_fraction = eta * G_cdf + (1 - eta) * L_cdf; + + % Shirley contribution direction depends on the energy scale + if matches(p.Results.energyScale, 'kinetic') + % Kinetic Energy: Background steps UP as kinetic energy goes DOWN + shirleycontrib = K * area * (1 - integral_fraction); + else + % Binding Energy: Background steps UP as binding energy goes UP + shirleycontrib = K * area * integral_fraction; + end + + % Add continuously to the total background + B = B + shirleycontrib; end - elseif matches(p.Results.mode,"TwoUnconstrainedAndWidthPositionPseudoVoigtPeak") && peakEta_exists - if iPeak == 1 - area = peak_params(1); - FWHM = peak_params(2); - pos = peak_params(3)+ peak_params(11); - eta = peak_params(13); - K = peak_params(16); - elseif iPeak == 2 - area = peak_params(4); - FWHM = peak_params(2); - pos = peak_params(3); - eta = peak_params(13); - K = peak_params(16); - elseif iPeak == 3 - area = peak_params(5); - FWHM = peak_params(6); - pos = peak_params(7); - eta = peak_params(14); - K = peak_params(17); - elseif iPeak == 4 - area = peak_params(8); - FWHM = peak_params(9); - pos = peak_params(10); - eta = peak_params(15); - K = peak_params(18); - 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 - - % Calculate the continuous CDF over the ENTIRE energy range - % This gives the solver a smooth gradient. - G_cdf = 0.5 * (1 + erf((E - pos) / (sqrt(2) * sigma))); - L_cdf = 0.5 + (1/pi) * atan((E - pos) / gamma); - - % Combined Pseudo-Voigt integral fraction (0 to 1) - integral_fraction = eta * G_cdf + (1 - eta) * L_cdf; - - % Shirley contribution direction depends on the energy scale - if matches(p.Results.energyScale, 'kinetic') - % Kinetic Energy: Background steps UP as kinetic energy goes DOWN - shirleycontrib = K * area * (1 - integral_fraction); - else - % Binding Energy: Background steps UP as binding energy goes UP - shirleycontrib = K * area * integral_fraction; - end - - % Add continuously to the total background - B = B + shirleycontrib; -end - - -end %% 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'); + 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); + % 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 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); + 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); + 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,"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); + elseif matches(p.Results.mode,"widthPositionShirleyEta") + area = params(iPeak); + FWHM = params(nPeaks + 1); + eta = params(nPeaks + 2); + pos1 = params(nPeaks + 3); + if iPeak == 1 + pos = pos1; + else + relPos = FitParameter.peakPos(1, iPeak); + pos = pos1 + relPos; + end + y = y + pseudoVoigtPeak(x, area, FWHM, pos, eta); - if iPeak == 1 - pos = pos1; - else - relPos = FitParameter.peakPos(1, iPeak); - pos = pos1 + relPos; + 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 - 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; + % 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 - 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 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' + function correctedStruct = fixBondToBound(inputStruct) + % Recursively corrects all field names in a struct: 'Bond' -> 'Bound' -correctedStruct = inputStruct; -fieldNames = fieldnames(inputStruct); + correctedStruct = inputStruct; + fieldNames = fieldnames(inputStruct); -for i = 1:numel(fieldNames) - thisField = fieldNames{i}; - value = inputStruct.(thisField); + for i = 1:numel(fieldNames) + thisField = fieldNames{i}; + value = inputStruct.(thisField); - % Recursively fix nested structs - if isstruct(value) - value = fixBondToBound(value); - end + % 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; + % 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 -end -end % Helper function: Check bounds consistency -function checkBoundsPair(FitParam, lowerName, upperName, fieldLabel) -if isfield(FitParam, lowerName) && isfield(FitParam, upperName) - lower = FitParam.(lowerName); - upper = FitParam.(upperName); + function checkBoundsPair(FitParam, lowerName, upperName, fieldLabel) + if isfield(FitParam, lowerName) && isfield(FitParam, upperName) + lower = FitParam.(lowerName); + upper = FitParam.(upperName); - 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)); + 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 end -end -end % 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 + function checkFieldDimension(FitParam, fieldName, expectedSize, modeStr, fieldLabel) + if ~isfield(FitParam, fieldName) + error('Mode ''%s'': Missing required field ''%s'' (%s).', ... + modeStr, fieldName, fieldLabel); + end -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 + 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 + % 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 -% 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 + % 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); + 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 + 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 + % 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 + end \ No newline at end of file