bartels-rausch
2377e3b1f2
This commit replaces the numerically unstable Shirley background calculation in the peak fitting routine with a continuous, analytical Pseudo-Voigt Cumulative Distribution Function (CDF).
Previously, the solver struggled to converge and pushed peak parameters (area, width) into their absolute bounds due to two critical mathematical bugs in the objective function:
1. The 'flip' and Energy Axis Mutation Bug:
To handle kinetic energy scales (where the background steps up as energy decreases), the code mutated the energy axis (E = max(E)+100-E) inside the peak loop and attempted to use array flipping. This approach fundamentally decoupled the background step from the physical peak position. Array flipping mirrors data based on the length of the plot window, physically moving the Shirley step away from the true peak center. Furthermore, mutating E inside the loop corrupted the X-axis for subsequent peaks in multi-peak fits.
2. The Integration Cutoff ('idxregion') Bug:
The prior implementation used a hard Boolean cutoff (idxregion = (E >= pos)) to restrict the background calculation to only one side of the peak. Because Pseudo-Voigt peaks have long tails, this created a severe mathematical discontinuity (a sharp cliff) exactly at the peak center. Non-linear optimization solvers require smooth, continuous gradients to step correctly. The sharp cutoff destroyed the gradient, causing the solver to panic and lock parameters against their constraints.
The Fix:
- Removed all 'flip()' calls and energy array mutations.
- Replaced the hard cutoff with a continuous analytical CDF evaluated over the entire energy range.
- Handled binding vs. kinetic energy directionality cleanly by toggling between 'integral_fraction' and '(1 - integral_fraction)'.
Result:
The non-linear solver now receives perfectly smooth, mathematically continuous gradients, leading to rapid convergence. Peak parameters no longer crash into their bounds, and the plotted red fit line now perfectly aligns with the internally calculated solver residuals.
308 lines
10 KiB
Matlab
308 lines
10 KiB
Matlab
function nappAddFitToPlot(DataTable, varargin)
|
|
% nappAddFitToPlot - Plot fitted peaks with proper decomposition
|
|
% FIXED: Correct pseudo-Voigt normalization and FWHM extraction
|
|
|
|
% Parse input arguments
|
|
p = inputParser;
|
|
addOptional(p, 'tableSelection', []);
|
|
addOptional(p, 'plotHandle', gca);
|
|
addOptional(p, 'peakWidthMode', 'FWHM');
|
|
addOptional(p, 'energyScale', 'kinetic');
|
|
parse(p, varargin{:});
|
|
|
|
% Get data from table
|
|
if isempty(p.Results.tableSelection)
|
|
fitData = DataTable(1, :);
|
|
else
|
|
fitData = DataTable(p.Results.tableSelection, :);
|
|
end
|
|
|
|
% Get energy range from current figure xlim
|
|
ax = p.Results.plotHandle;
|
|
axes(ax);
|
|
xlims = xlim(ax);
|
|
energyMin = xlims(1);
|
|
energyMax = xlims(2);
|
|
|
|
% Create energy range with 0.2 eV step
|
|
energyStep = 0.2;
|
|
energyRange = (energyMin:energyStep:energyMax)';
|
|
|
|
% =====================================================================
|
|
% COUNT ONLY NON-NaN PEAKS
|
|
% =====================================================================
|
|
peakArea_vals = fitData.peakArea; % Extract from cell
|
|
validPeaks = ~isnan(peakArea_vals); % Find non-NaN values
|
|
nPeaks = sum(validPeaks); % Count valid peaks
|
|
|
|
fprintf('Found %d valid peaks (non-NaN)\n', nPeaks);
|
|
|
|
if isempty(energyRange) || nPeaks == 0
|
|
warning('No valid fit data to plot');
|
|
return;
|
|
end
|
|
|
|
% Initialize output
|
|
peakCurves = zeros(length(energyRange), nPeaks);
|
|
totalFit = zeros(size(energyRange));
|
|
peakTypes = cell(nPeaks, 1);
|
|
|
|
% Extract peak parameters (cell arrays)
|
|
peakPos_vals = fitData.peakPosition;
|
|
peakFWHM_vals = fitData.peakFWHM;
|
|
|
|
% Check for peak eta parameter (pseudo-Voigt indicator)
|
|
hasPeakEta = ismember('peakEta', fitData.Properties.VariableNames) && ...
|
|
~all(isnan(fitData.peakEta));
|
|
|
|
if hasPeakEta
|
|
peakEta_vals = fitData.peakEta;
|
|
|
|
end
|
|
|
|
% =====================================================================
|
|
% MAIN LOOP: Process each VALID peak
|
|
% =====================================================================
|
|
validIdx = find(validPeaks); % Get indices of valid peaks
|
|
for i = 1:nPeaks
|
|
iPeak = validIdx(i); % Actual peak index in array
|
|
|
|
% Extract peak parameters
|
|
pos = peakPos_vals(iPeak);
|
|
A_total = peakArea_vals(iPeak);
|
|
FWHM = peakFWHM_vals(iPeak);
|
|
|
|
% Skip if FWHM is NaN or invalid
|
|
if isnan(FWHM) || FWHM <= 0
|
|
fprintf(' WARNING: Peak %d has invalid FWHM = %g, skipping\n', iPeak, FWHM);
|
|
continue;
|
|
end
|
|
|
|
% Check if this is a pseudo-Voigt peak
|
|
if hasPeakEta && ~isnan(peakEta_vals(iPeak))
|
|
% ============================================================
|
|
% PSEUDO-VOIGT: Gaussian + Lorentzian blend
|
|
% eta = Gaussian fraction (0 = pure Lorentzian, 1 = pure Gaussian)
|
|
% ============================================================
|
|
eta = peakEta_vals(iPeak);
|
|
|
|
% === GAUSSIAN COMPONENT (FWHM-based) ===
|
|
% G(E) = A_gauss * sqrt(4*ln(2)/pi) / FWHM * exp(-4*ln(2)*(E-E0)^2/FWHM^2)
|
|
% This is normalized so integral = A_gauss
|
|
A_gauss = eta * A_total; % Gaussian gets eta fraction of total area
|
|
gauss_norm = sqrt(4*log(2)/pi) / FWHM;
|
|
peakG = A_gauss .* gauss_norm .* ...
|
|
exp(-4*log(2) * ((energyRange - pos)./FWHM).^2);
|
|
|
|
% === LORENTZIAN COMPONENT (FWHM-based) ===
|
|
% L(E) = A_lorentz * (FWHM/2) / pi / ((E-E0)^2 + (FWHM/2)^2)
|
|
% This is normalized so integral = A_lorentz
|
|
A_lorentz = (1 - eta) * A_total; % Lorentzian gets (1-eta) fraction
|
|
gamma = FWHM / 2; % Half-width at half-maximum
|
|
peakL = A_lorentz .* (gamma / pi) ./ ...
|
|
((energyRange - pos).^2 + gamma.^2);
|
|
|
|
% === COMBINED PSEUDO-VOIGT ===
|
|
peak = peakG + peakL;
|
|
|
|
peakTypes{i} = sprintf('Pseudo-Voigt (η=%.3f)', eta);
|
|
|
|
else
|
|
% ============================================================
|
|
% PURE GAUSSIAN using FWHM
|
|
% G(E) = A * sqrt(4*ln(2)/pi) / FWHM * exp(-4*ln(2)*(E-E0)^2/FWHM^2)
|
|
% ============================================================
|
|
gauss_norm = sqrt(4*log(2)/pi) / FWHM;
|
|
peak = A_total .* gauss_norm .* ...
|
|
exp(-4*log(2) * ((energyRange - pos)./FWHM).^2);
|
|
peakTypes{i} = 'Gaussian';
|
|
end
|
|
|
|
peakCurves(:, i) = peak;
|
|
totalFit = totalFit + peak;
|
|
end
|
|
|
|
% =====================================================================
|
|
% ADD BACKGROUND
|
|
% =====================================================================
|
|
% Check for Shirley background parameters
|
|
hasShirleyK = ismember('peakShirleyK', fitData.Properties.VariableNames) && ...
|
|
~all(isnan(fitData.peakShirleyK));
|
|
|
|
|
|
if hasShirleyK
|
|
% SHIRLEY BACKGROUND RECONSTRUCTION
|
|
peakShirleyK_vals = fitData.peakShirleyK;
|
|
bkgrd = zeros(size(energyRange));
|
|
|
|
for i = 1:nPeaks
|
|
iPeak = validIdx(i);
|
|
|
|
% Extract parameters
|
|
area = peakArea_vals(iPeak);
|
|
FWHM = peakFWHM_vals(iPeak);
|
|
pos = peakPos_vals(iPeak);
|
|
if ismember('peakEta', fitData.Properties.VariableNames) && ~all(isnan(fitData.peakEta))
|
|
eta = peakEta_vals(iPeak);
|
|
else
|
|
eta = 1;
|
|
end
|
|
K = peakShirleyK_vals(iPeak);
|
|
|
|
% Calculate width parameters
|
|
w = FWHM / sqrt(log(4));
|
|
sigma = w / 2;
|
|
gamma = FWHM / 2;
|
|
|
|
% Cumulative Distribution Functions (CDFs) going from 0 to 1
|
|
% This represents the integral of the peak from -infinity to E
|
|
G_cdf = 0.5 * (1 + erf((energyRange - pos) / (sqrt(2) * sigma)));
|
|
L_cdf = 0.5 + (1/pi) * atan((energyRange - pos) / gamma);
|
|
|
|
% Combined Pseudo-Voigt CDF
|
|
integral_fraction = eta * G_cdf + (1 - eta) * L_cdf;
|
|
|
|
% Shirley Contribution
|
|
if matches(p.Results.energyScale, 'kinetic')
|
|
% For Kinetic Energy: background is higher at LOWER energies.
|
|
% Step goes from 1 (at E << pos) to 0 (at E >> pos)
|
|
shirleycontrib = K * area * (1 - integral_fraction);
|
|
else
|
|
% For Binding Energy: background is higher at HIGHER energies.
|
|
% Step goes from 0 (at E << pos) to 1 (at E >> pos)
|
|
shirleycontrib = K * area * integral_fraction;
|
|
end
|
|
|
|
bkgrd = bkgrd + shirleycontrib;
|
|
end
|
|
|
|
% Add constant offset C0
|
|
if ismember('peakBackground', fitData.Properties.VariableNames)
|
|
C0 = fitData.peakBackground(2); % Second element is C0
|
|
bkgrd = bkgrd + C0;
|
|
end
|
|
|
|
backgroundType = 'Shirley';
|
|
% if hasShirleyK
|
|
% % SHIRLEY BACKGROUND RECONSTRUCTION
|
|
% peakShirleyK_vals = fitData.peakShirleyK;
|
|
%
|
|
% % Reconstruct Shirley background using the same algorithm as in fitting
|
|
% bkgrd = zeros(size(energyRange));
|
|
%
|
|
% for i = 1:nPeaks
|
|
% iPeak = validIdx(i);
|
|
%
|
|
% % Extract parameters
|
|
% area = peakArea_vals(iPeak);
|
|
% FWHM = peakFWHM_vals(iPeak);
|
|
% pos = peakPos_vals(iPeak);
|
|
% if ismember('peakEta', fitData.Properties.VariableNames) && ~all(isnan(fitData.peakEta))
|
|
% eta = peakEta_vals(iPeak);
|
|
% else
|
|
% eta = 1;
|
|
% end
|
|
% K = peakShirleyK_vals(iPeak);
|
|
%
|
|
% % Find energy points after the peak
|
|
% idxafter = energyRange >= pos;
|
|
%
|
|
% if any(idxafter)
|
|
% % Calculate width parameters
|
|
% w = FWHM / sqrt(log(4)); % Convert FWHM to Gaussian w
|
|
% sigma = w / 2;
|
|
% gamma = FWHM / 2;
|
|
%
|
|
% % Gaussian tail fraction
|
|
% Gtail = 0.5 * (1 - erf((energyRange(idxafter) - pos) / (sqrt(2) * sigma)));
|
|
%
|
|
% % Lorentzian tail fraction
|
|
% Ltail = 0.5 - (1/pi) * atan((energyRange(idxafter) - pos) / gamma);
|
|
%
|
|
% % Pseudo-Voigt tail fraction
|
|
% tailfrac = eta * Gtail + (1 - eta) * Ltail;
|
|
%
|
|
% % Shirley contribution
|
|
% shirleycontrib = K * area * (0.5 - tailfrac);
|
|
% bkgrd(idxafter) = bkgrd(idxafter) + shirleycontrib;
|
|
%
|
|
% end
|
|
%
|
|
% end
|
|
%
|
|
% % Add constant offset C0
|
|
% if ismember('peakBackground', fitData.Properties.VariableNames)
|
|
% C0 = fitData.peakBackground(2); % Second element is C0
|
|
% bkgrd = bkgrd + C0;
|
|
% end
|
|
%
|
|
% if matches(p.Results.energyScale,'kinetic')
|
|
% % This is a little work around - needs confirmation if it does what
|
|
% % it should - be careful.
|
|
% bkgrd = flip(bkgrd);
|
|
% end
|
|
%
|
|
% backgroundType = 'Shirley';
|
|
else
|
|
% LINEAR BACKGROUND (existing code)
|
|
if ismember('peakBackground', fitData.Properties.VariableNames)
|
|
bkgrd_vals = fitData.peakBackground;
|
|
if length(bkgrd_vals) >= 2
|
|
slope = bkgrd_vals(1);
|
|
intercept = bkgrd_vals(2);
|
|
else
|
|
slope = 0;
|
|
intercept = bkgrd_vals(1);
|
|
end
|
|
else
|
|
slope = 0;
|
|
intercept = 0;
|
|
end
|
|
|
|
bkgrd = slope .* energyRange + intercept;
|
|
backgroundType = sprintf('Linear (slope=%.4f, intercept=%.4f)', slope, intercept);
|
|
end
|
|
|
|
totalFit = totalFit + bkgrd;
|
|
|
|
|
|
% =====================================================================
|
|
% PLOT RESULTS
|
|
% =====================================================================
|
|
axes(ax);
|
|
hold on;
|
|
|
|
% Plot total fit (thick red line)
|
|
plot(energyRange, totalFit, '-', 'LineWidth', 2.5, 'Color', [0.8 0.2 0.2], 'DisplayName', 'Total Fit');
|
|
|
|
% Plot background (dotted gray line)
|
|
plot(energyRange, bkgrd, ':', 'LineWidth', 1.5, 'Color', [0.5 0.5 0.5], 'DisplayName', 'Background');
|
|
|
|
% Plot individual peaks (dashed colored lines)
|
|
colors = lines(nPeaks);
|
|
for i = 1:nPeaks
|
|
if any(peakCurves(:, i)) % Only plot if non-zero
|
|
plot(energyRange, peakCurves(:, i)+bkgrd, '--', 'LineWidth', 1.2, ...
|
|
'Color', colors(i,:), 'DisplayName', sprintf('Peak %d', i));
|
|
end
|
|
end
|
|
|
|
legend('show', 'Location', 'best');
|
|
hold off;
|
|
|
|
% % =====================================================================
|
|
% % DISPLAY INFORMATION
|
|
% % =====================================================================
|
|
% fprintf('\n=== FIT PLOTTED SUCCESSFULLY ===\n');
|
|
% fprintf('Valid peaks: %d\n', nPeaks);
|
|
% fprintf('Peak type(s): ');
|
|
% for i = 1:nPeaks
|
|
% fprintf('%s', peakTypes{i});
|
|
% if i < nPeaks, fprintf('; '); end
|
|
% end
|
|
% fprintf('\n');
|
|
% fprintf('Energy range: %.1f - %.1f eV\n', min(energyRange), max(energyRange));
|
|
% fprintf('=====================================\n\n');
|
|
% fprintf('Background: %s\n', backgroundType);
|
|
end
|