From 0e07d2972b049978f4a007fb8a025fbae04ecf26 Mon Sep 17 00:00:00 2001 From: ThorBarRa Date: Thu, 5 Mar 2026 11:11:37 +0100 Subject: [PATCH] Finally, several spectra can be fitted with same constraints again...function now loops through any number of spectra --- beamtime/nappFitGaussXps.m | 1688 ++++++++++++++++++------------------ 1 file changed, 843 insertions(+), 845 deletions(-) diff --git a/beamtime/nappFitGaussXps.m b/beamtime/nappFitGaussXps.m index e89fec0..b547cce 100644 --- a/beamtime/nappFitGaussXps.m +++ b/beamtime/nappFitGaussXps.m @@ -6,7 +6,6 @@ function DataTable = nappFitGaussXps(DataTable,varargin) % % 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. % % DataTable = NAPPFITGAUSSXPSDEVELOPMENT(..., 'Name', Value) specifies % additional options using name-value pair arguments: @@ -621,13 +620,6 @@ end % 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. -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 @@ -1653,16 +1645,16 @@ for iuz = 1 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. -% ========================================================================= + + %% 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. + % ========================================================================= @@ -1804,7 +1796,7 @@ for iuz = 1 initialParameterVector = [initialParameterVector, FitParameter.peakArea(1, 3)]; lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakAreaLowerBound(1, 3)]; upperBoundsParameterVector = [upperBoundsParameterVector, FitParameter.peakAreaUpperBound(1, 3)]; - + % #5: Width 3+4 initialParameterVector = [initialParameterVector, FitParameter.peakWidth(1, 2)]; lowerBoundsParameterVector = [lowerBoundsParameterVector, FitParameter.peakWidthLowerBound(1, 2)]; @@ -1968,6 +1960,8 @@ end %% Core part: Run the fit + +% loop through data for ii = 1:size(DataTable(p.Results.tableSelection,:),1) if contains(p.Results.energyScale,'kinetic') @@ -1998,7 +1992,7 @@ for ii = 1:size(DataTable(p.Results.tableSelection,:),1) %% Core fitting part, % the fitting is done here (kinetic energy scale) %-------------------------------------------------------------------------- - + clearvars currentParameterVector % Fit equations: Gauss eq. 1 with data modelEvaluationFunction = @(currentParameterVector,energyValues) modelDefinitionFunction(currentParameterVector,energyValues); % We find an optimal fit by minimising the residuals: @@ -2033,7 +2027,7 @@ for ii = 1:size(DataTable(p.Results.tableSelection,:),1) %% Core fitting part % the fitting is done here (binding energy scale) %-------------------------------------------------------------------------- - + clearvars currentParameterVector % Fit equations: Gauss eq. 1 with data modelEvaluationFunction = @(currentParameterVector,energyValues) modelDefinitionFunction(currentParameterVector,energyValues); % We find an optimal fit by minimising the residuals: @@ -2045,889 +2039,893 @@ for ii = 1:size(DataTable(p.Results.tableSelection,:),1) % 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. -resultsPeakPos = []; -resultsPeakArea = []; -resultsPeakFWHM = []; -resultsBkgrd = []; -resultsPeakEta = []; -resultsPeakShirleyK = []; -nPeaks = size(FitParameter.peakArea, 2); + %% Results Extraction + % + % Extract fitting results from currentParameterVector for EACH MODEL TYPE + % and populate result arrays: resultsPeakPos, resultsPeakArea, resultsPeakFWHM, resultsBkgrd, etc. -for iii= 1 + resultsPeakPos = []; + resultsPeakArea = []; + resultsPeakFWHM = []; + resultsBkgrd = []; + resultsPeakEta = []; + resultsPeakShirleyK = []; + nPeaks = size(FitParameter.peakArea, 2); - %% Generic Models + for iii= 1 - if matches(p.Results.mode, 'unconstrained') + %% Generic Models - %% MODE: UNCONSTRAINED + if matches(p.Results.mode, 'unconstrained') - if matches(p.Results.backgroundMode, 'linear') && ~isfield(FitParameter,'peakEta') - % ===== GAUSSIAN + LINEAR (UNCONSTRAINED) ===== + %% MODE: UNCONSTRAINED - for iPeak = 1:nPeaks - areaIdx = 3*iPeak - 2; - widthIdx = 3*iPeak - 1; - posIdx = 3*iPeak; + if matches(p.Results.backgroundMode, 'linear') && ~isfield(FitParameter,'peakEta') + % ===== GAUSSIAN + LINEAR (UNCONSTRAINED) ===== - resultsPeakArea = [resultsPeakArea, currentParameterVector(areaIdx)]; - resultsPeakFWHM = [resultsPeakFWHM, currentParameterVector(widthIdx) ]; - resultsPeakPos = [resultsPeakPos, currentParameterVector(posIdx)]; - end - resultsBkgrd = [currentParameterVector(end-1), currentParameterVector(end)]; + for iPeak = 1:nPeaks + areaIdx = 3*iPeak - 2; + widthIdx = 3*iPeak - 1; + posIdx = 3*iPeak; - elseif matches(p.Results.backgroundMode, 'linear') && isfield(FitParameter,'peakEta') - % ===== PSEUDO-VOIGT + LINEAR (UNCONSTRAINED) ===== + resultsPeakArea = [resultsPeakArea, currentParameterVector(areaIdx)]; + resultsPeakFWHM = [resultsPeakFWHM, currentParameterVector(widthIdx) ]; + resultsPeakPos = [resultsPeakPos, currentParameterVector(posIdx)]; + end + resultsBkgrd = [currentParameterVector(end-1), currentParameterVector(end)]; - for iPeak = 1:nPeaks - areaIdx = 4*iPeak - 3; - widthIdx = 4*iPeak - 2; - posIdx = 4*iPeak - 1; - etaIdx = 4*iPeak; + elseif matches(p.Results.backgroundMode, 'linear') && isfield(FitParameter,'peakEta') + % ===== PSEUDO-VOIGT + LINEAR (UNCONSTRAINED) ===== - 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)]; + for iPeak = 1:nPeaks + areaIdx = 4*iPeak - 3; + widthIdx = 4*iPeak - 2; + posIdx = 4*iPeak - 1; + etaIdx = 4*iPeak; - 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)]; %#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)]; - 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.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); - 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; + 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)]; - resultsPeakArea = [resultsPeakArea, currentParameterVector(areaIdx)]; - resultsPeakFWHM = [resultsPeakFWHM, currentParameterVector(fwhmIdx)]; - resultsPeakPos = [resultsPeakPos, currentParameterVector(posIdx)]; + 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 - for iPeak = 1:nPeaks - resultsPeakShirleyK = [resultsPeakShirleyK, currentParameterVector(3*nPeaks+iPeak)]; - end - % background (Shirley intercept only) - resultsBkgrd = [0,currentParameterVector(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 + %% Custom Models - 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)]; + if matches(p.Results.mode, 'unconstrainedAndWidthPosition') && matches(p.Results.backgroundMode, 'linear') + %% unconstrainedAndWidthPosition + % ===== GAUSSIAN + LINEAR ===== - 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; + % 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); - resultsPeakArea = [resultsPeakArea, currentParameterVector(areaIdx)]; - resultsPeakFWHM = [resultsPeakFWHM, currentParameterVector(widthIdx)]; - resultsPeakEta = [resultsPeakEta, currentParameterVector(etaIdx)]; - end - pos1_fit = currentParameterVector(3*nPeaks + 1); + % Calculate absolute positions + pos1 = pos2_rel + pos1_shift; + pos2 = pos2_rel; + pos4 = pos4_offset; + pos5 = pos5_offset; - 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)]; + % Results: areas, positions, FWHM + resultsPeakArea = [area1, area2, area3, area4, area5]; + resultsPeakPos = [pos1, pos2, pos3, pos4, pos5]; + resultsPeakFWHM = [width1, width1, width3, width1, width1]; - 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; + % Background + resultsBkgrd = [slope, intercept]; - resultsPeakArea = [resultsPeakArea, currentParameterVector(areaIdx)]; - resultsPeakFWHM = [resultsPeakFWHM, currentParameterVector(widthIdx) ]; - resultsPeakEta = [resultsPeakEta, currentParameterVector(etaIdx)]; - resultsPeakShirleyK = [resultsPeakShirleyK, currentParameterVector(kIdx)]; - end - pos1_fit = currentParameterVector(3*nPeaks + 1); + elseif matches(p.Results.mode, 'TwoSpinOrbitDoublet') + %% TwoSpinOrbitDoublet + % ===== GAUSSIAN + LINEAR ===== - 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; + % 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); - resultsPeakArea = [resultsPeakArea, currentParameterVector(areaIdx)]; - resultsPeakFWHM = [resultsPeakFWHM, currentParameterVector(widthIdx) ]; - resultsPeakShirleyK = [resultsPeakShirleyK, currentParameterVector(kIdx)]; - end - pos1_fit = currentParameterVector(2*nPeaks+1); + % Results: areas, positions, FWHM + resultsPeakArea = [area1, area2, area3, area4]; + resultsPeakPos = [pos1, pos2, pos3, pos4]; + resultsPeakFWHM = [width1, width2, width3, width4]; - 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 - %% 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]; + % Background + resultsBkgrd = [slope, intercept]; elseif matches(p.Results.mode, 'GasAndSurfaceSpinOrbitDoublet') - %% GasAndSurfaceSpinOrbitDoublet - % ===== GAUSSIAN + LINEAR ===== + %% 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); + % 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]; + % Results: areas, positions, FWHM + resultsPeakArea = [area1, area2, area3, area4]; + resultsPeakPos = [pos1, pos2, pos3, pos4]; + resultsPeakFWHM = [width1, width2, width3, width4]; - % Background - resultsBkgrd = [slope, intercept]; + % Background + resultsBkgrd = [slope, intercept]; - elseif matches(p.Results.mode, 'TwoSpinOrbitDoubletPseudoVoigtShirley') - %% TwoSpinOrbitDoubletPseudoVoigtShirley - % ===== PSEUDO-VOIGT + SHIRLEY ===== + 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); + % 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]; + % 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)]; + % 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 ===== + elseif matches(p.Results.mode,'gasAndTwoAcidBasePairs') + %% gasAndTwoAcidBasePairs + % ===== GAUSSIAN + LINEAR ===== - % Extract parameters - area1 = currentParameterVector(1); - width1 = currentParameterVector(2); - pos1 = currentParameterVector(3)+currentParameterVector(14); + % 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); + width2 = currentParameterVector(2); + pos2 = currentParameterVector(3); - resultsPeakArea = [resultsPeakArea, area2, area3, area4, area5]; + area3 = currentParameterVector(5); + width3 = currentParameterVector(6); + pos3 = currentParameterVector(7); - for iPeak = 2:5 - resultsPeakPos = [resultsPeakPos, pos1 + FitParameter.peakPos(1, iPeak)]; + 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 - 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); + end + end - resultsPeakArea = [resultsPeakArea, area2, area3, area4, area5]; + %% Write Results to table - for iPeak = 2:5 - resultsPeakPos = [resultsPeakPos, pos1 + FitParameter.peakPos(1, iPeak)]; + % 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 + + %% Print Results + if matches(p.Results.showFigures,'on') + disp(DataTable(p.Results.tableSelection,["name","peakPosition"])) + disp(DataTable(p.Results.tableSelection,["name","peakArea"])) + disp(DataTable(p.Results.tableSelection,["name","peakFWHM"])) + disp(DataTable(p.Results.tableSelection,["name","peakBackground"])) + try + if ~isnan(DataTable(p.Results.tableSelection,:).peakEta(1)) + disp(DataTable(p.Results.tableSelection,["name","peakEta",])) end - - resultsPeakFWHM = ones(1, 5) * width_shared ; - resultsBkgrd = [currentParameterVector(7), currentParameterVector(8)]; + catch end + try + if ~isnan(DataTable(p.Results.tableSelection,:).peakShirleyK(1)) + disp(DataTable(p.Results.tableSelection,["name","peakShirleyK"])) + end + catch - 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.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"])) - disp(DataTable(p.Results.tableSelection,["name","peakArea"])) - disp(DataTable(p.Results.tableSelection,["name","peakFWHM"])) - disp(DataTable(p.Results.tableSelection,["name","peakBackground"])) - try - if ~isnan(DataTable(p.Results.tableSelection,:).peakEta(1)) - disp(DataTable(p.Results.tableSelection,["name","peakEta",])) end - catch end - try - if ~isnan(DataTable(p.Results.tableSelection,:).peakShirleyK(1)) - disp(DataTable(p.Results.tableSelection,["name","peakShirleyK"])) - end - catch + %% 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 + + 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