function plotSecondaryGaussianRange(fitResults, scanValues, varargin)
%% plotSecondaryGaussianRange
% Plots secondary Gaussian ranges (mu2 ± sigma2) in a heatmap style
% exactly like plotFitParameterPDF, with dashed lines for distribution
% edges and solid line for mean mu2.

% --- Parse optional inputs ---
p = inputParser;
addParameter(p, 'Title', '', @(x) ischar(x) || isstring(x));
addParameter(p, 'XLabel', '', @(x) ischar(x) || isstring(x));
addParameter(p, 'YLabel', '', @(x) ischar(x) || isstring(x));
addParameter(p, 'FigNum', 1, @(x) isscalar(x));
addParameter(p, 'FontName', 'Arial', @ischar);
addParameter(p, 'FontSize', 14, @isnumeric);
addParameter(p, 'SkipSaveFigures', false, @islogical);
addParameter(p, 'SaveFileName', 'SecondaryGaussianRange.fig', @ischar);
addParameter(p, 'SaveDirectory', pwd, @ischar);
addParameter(p, 'NumPoints', 200, @(x) isscalar(x));
addParameter(p, 'DataRange', [], @(x) isempty(x) || numel(x)==2);
addParameter(p, 'XLim', [], @(x) isempty(x) || numel(x)==2);
addParameter(p, 'YLim', [], @(x) isempty(x) || numel(x)==2);
addParameter(p, 'Colormap', @jet);
addParameter(p, 'PlotType', 'histogram', @(x) any(validatestring(x,{'kde','histogram'})));
addParameter(p, 'NumberOfBins', 50, @isscalar);
addParameter(p, 'NormalizeHistogram', true, @islogical);
addParameter(p, 'OverlayMeanSEM', true, @islogical);
parse(p, varargin{:});
opts = p.Results;

% --- Determine repetitions and scan parameters ---
N_params = numel(scanValues);
N_total  = numel(fitResults);
N_reps   = N_total / N_params;

% --- Extract mu2 and sigma2 ---
mu2_values    = nan(N_reps, N_params);
sigma2_values = nan(N_reps, N_params);
for k = 1:N_total
    paramIdxScan = mod(k-1, N_params) + 1;
    repIdx       = floor((k-1)/N_params) + 1;
    mu2_values(repIdx, paramIdxScan)    = fitResults(k).pFit(5);
    sigma2_values(repIdx, paramIdxScan) = fitResults(k).pFit(6);
end

% --- Mask invalid fits ---
for i = 1:N_total
    paramIdxScan = mod(i-1, N_params) + 1;
    repIdx       = floor((i-1)/N_params) + 1;
    if ~fitResults(i).isValid && all(fitResults(i).pFit == 0)
        mu2_values(repIdx, paramIdxScan)    = NaN;
        sigma2_values(repIdx, paramIdxScan) = NaN;
    end
end

% --- Compute lower and upper bounds ---
lowerBound = mu2_values - sigma2_values;
upperBound = mu2_values + sigma2_values;

% Define how much to trim (e.g., 2% on each side)
trimFraction = 0.01;  % can tune this (0.01–0.05 typical)

% Compute robust global limits
globalLower = quantile(lowerBound(:), trimFraction);
globalUpper = quantile(upperBound(:), 1 - trimFraction);

fprintf('Global lower bound (%.0f%%): %.4f\n', trimFraction*100, globalLower);
fprintf('Global upper bound (%.0f%%): %.4f\n', (1-trimFraction)*100, globalUpper);

% --- Prepare data per scan parameter ---
dataCell = cell(N_params,1);
for i = 1:N_params
    combinedData = [];
    for j = 1:N_reps
        if ~isnan(mu2_values(j,i)) && ~isnan(sigma2_values(j,i))
            combinedData = [combinedData; linspace(lowerBound(j,i), upperBound(j,i), 3)']; %#ok<AGROW>
        end
    end
    dataCell{i} = combinedData;
end

% --- Determine y-range ---
if isempty(opts.DataRange)
    allData = cell2mat(dataCell(:));
    y_min   = min(allData);
    y_max   = max(allData);
else
    y_min   = opts.DataRange(1);
    y_max   = opts.DataRange(2);
end

% --- Prepare PDF matrix ---
if strcmpi(opts.PlotType,'kde')
    y_grid = linspace(y_min, y_max, opts.NumPoints);
    pdf_matrix = zeros(numel(y_grid), N_params);
    for i = 1:N_params
        data = dataCell{i};
        if isempty(data), continue; end
        pdf_matrix(:,i) = ksdensity(data, y_grid);
    end
else
    edges = linspace(y_min, y_max, opts.NumberOfBins+1);
    binCenters = (edges(1:end-1)+edges(2:end))/2;
    pdf_matrix = zeros(numel(binCenters), N_params);
    for i = 1:N_params
        data = dataCell{i};
        if isempty(data), continue; end
        counts = histcounts(data, edges);
        if opts.NormalizeHistogram
            counts = counts / sum(counts);
        end
        pdf_matrix(:,i) = counts(:);
    end
end

% --- Mask out non-data regions ---
maskPerScan = cellfun(@(c) ~isempty(c), dataCell);
pdf_matrix(:, ~maskPerScan) = NaN;

% --- Plot heatmap ---
fig = figure(opts.FigNum); clf(fig);
set(fig, 'Color', 'w', 'Position', [100 100 950 750]);

if strcmpi(opts.PlotType,'kde')
    h = imagesc(scanValues, y_grid, pdf_matrix);
    set(h, 'AlphaData', ~isnan(pdf_matrix));
else
    h = imagesc(scanValues, binCenters, pdf_matrix);
    set(h, 'AlphaData', ~isnan(pdf_matrix));
end

colormap([1 1 1; feval(opts.Colormap)]); % white for NaNs
c = colorbar;
ylabel(c, 'Probability Density', 'Rotation', -90, 'FontName', opts.FontName, 'FontSize', opts.FontSize);

set(gca, 'YDir','normal', 'FontName', opts.FontName, 'FontSize', opts.FontSize);
xlabel(opts.XLabel, 'FontName', opts.FontName, 'FontSize', opts.FontSize);
ylabel(opts.YLabel, 'FontName', opts.FontName, 'FontSize', opts.FontSize);
title(opts.Title, 'FontName', opts.FontName, 'FontSize', opts.FontSize+2, 'FontWeight','bold');
grid on; hold on;

% --- Overlay mu2 and distribution edges ---
meanMu2  = nanmean(mu2_values,1);
minLower = min(lowerBound,[],1);
maxUpper = max(upperBound,[],1);

plot(scanValues, meanMu2, 'k-', 'LineWidth', 2);      % solid mu2
plot(scanValues, minLower, 'k--', 'LineWidth', 1.5);  % dashed lower edge
plot(scanValues, maxUpper, 'k--', 'LineWidth', 1.5);  % dashed upper edge

% --- Overlay mean ± SEM if requested ---
if opts.OverlayMeanSEM
    semMu2 = nanstd(mu2_values,0,1) ./ sqrt(sum(~isnan(mu2_values),1));
    xVec = reshape(scanValues,1,[]);
    fill([xVec, fliplr(xVec)], [meanMu2 - semMu2, fliplr(meanMu2 + semMu2)], ...
         [0.2 0.4 0.8], 'FaceAlpha',0.2, 'EdgeColor','none');
end

% --- Apply axis limits if specified ---
if ~isempty(opts.XLim)
    xlim(opts.XLim);
end
if ~isempty(opts.YLim)
    ylim(opts.YLim);
end

end