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