249 lines
8.4 KiB
Matlab
249 lines
8.4 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.01 eV step
|
|
energyStep = 0.01;
|
|
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';
|
|
|
|
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
|