function results = analyzeAutocorrelation(autocorrresults)
%%   Computes features from g2(theta) curves and aggregates results.
%   Assumes theta_values are given in radians.

    feature_list = table();
    per_group    = cell(numel(autocorrresults.scan_parameter_values),1);

    % --- Loop over groups ---
    for i = 1:numel(autocorrresults.scan_parameter_values)
        curves  = autocorrresults.g2_curves{i};
        thetas  = autocorrresults.theta_values;

        group_tbl = table();
        for j = 1:size(curves,1)
            g2curve      = curves(j,:);
            feats        = extractFeatures(thetas, g2curve);
            group_tbl    = [group_tbl; feats]; %#ok<AGROW>
            feature_list = [feature_list; feats]; %#ok<AGROW>
        end
        per_group{i} = group_tbl;
    end

    %% Compute group summaries (mean & SEM of features)
    summary_table = table();
    N_params = numel(autocorrresults.scan_parameter_values);

    for i = 1:N_params
        FT             = per_group{i};

        % Scalars
        vars           = {'A2','A4','A6','S2','Q4','H6','Contrast','NumberOfPeaks','MaxProm','MaxHeight'};
        
        % --- Handle spacing ---
        spacing_means  = cellfun(@mean, FT.Spacing, 'UniformOutput', false);
        spacing_stds   = cellfun(@std,  FT.Spacing, 'UniformOutput', false);
        FT.SpacingMean = cell2mat(spacing_means);
        FT.SpacingStd  = cell2mat(spacing_stds);

        % --- Handle widths ---
        width_means    = cellfun(@mean, FT.Widths, 'UniformOutput', false);
        width_stds     = cellfun(@std,  FT.Widths, 'UniformOutput', false);
        FT.WidthsMean  = cell2mat(width_means);
        FT.WidthsStd   = cell2mat(width_stds);

        vars           = [vars, {'SpacingMean','SpacingStd','WidthsMean','WidthsStd'}];

        % Group averages
        meanVals       = varfun(@mean, FT, 'InputVariables', vars);
        semVals        = varfun(@(x) std(x,0,1,'omitnan')/sqrt(size(FT,1)), ...
                          FT, 'InputVariables', vars);

        row = table(autocorrresults.scan_parameter_values(i), ...
            'VariableNames', {'scanParamVal'});
        summary_table = [summary_table; [row meanVals semVals]]; %#ok<AGROW>
    end

    %% --- Compute cumulants from full distribution of curves ---
    N_theta  = numel(autocorrresults.theta_values);
    g2_cums  = cell(N_params,1);

    for i = 1:N_params
        curves = autocorrresults.g2_curves{i}; % [N_reps × N_theta]
        cumul_mat = zeros(N_theta,4); % θ × first 4 cumulants
        for th = 1:N_theta
            g2vals = curves(:,th); % all repetitions at this theta
            cumul_mat(th,:) = Calculator.computeCumulants(g2vals,4);
        end
        g2_cums{i} = cumul_mat;
    end

    %% Package results
    results = struct();
    results.features_all   = feature_list;   % all repetitions pooled
    results.features_group = per_group;      % per scan parameter
    results.summaries      = summary_table;  % group-level stats
    results.g2_cumulants   = g2_cums;        % cumulants from full distribution
    results.theta_values   = autocorrresults.theta_values;
end

function features = extractFeatures(thetas, g2curve)
    %% --- Step 1: Restrict theta to [0, π] ---
    mask           = (thetas >= 0 & thetas <= pi);
    thetasWithinPi = thetas(mask);
    g2WithinPi     = g2curve(mask);

    %% --- Step 2: DC removal ---
    g2cWithinPi    = g2WithinPi - mean(g2WithinPi);

    %% --- Step 3: Fourier projections (n=2,4,6) ---
    nlist = [2 4 6];
    An = zeros(size(nlist));
    for k = 1:numel(nlist)
        n = nlist(k);
        An(k) = abs(mean(g2cWithinPi .* exp(-1i*n*thetasWithinPi)));
    end

    %% --- Step 3b: Symmetry fractions ---
    totalAmp = sum(An);
    if totalAmp > 0
        S2 = An(1)/totalAmp;
        Q4 = An(2)/totalAmp;
        H6 = An(3)/totalAmp;
    else
        S2 = 0; Q4 = 0; H6 = 0;
    end

    %% --- Step 4: RMS fluctuation (contrast) ---
    RMSval = sqrt(mean(g2cWithinPi.^2));

    %% --- Step 5: Peak finding with dynamic threshold ---
    alpha        = 0.25;                      
    dynamicProm  = alpha * range(g2WithinPi);     
    minDist      = round(numel(thetasWithinPi)/12); 

    [pks,locs,w,prom] = findpeaks(g2WithinPi, ...
        'MinPeakProminence', dynamicProm, ...
        'MinPeakDistance', minDist);

    peakAngles = thetasWithinPi(locs);

    %% --- Step 6: Spacing calculation (wrap within π) ---
    if numel(peakAngles) > 1
        d = diff([peakAngles, peakAngles(1)+pi]);
    else
        d = NaN;
    end

    %% --- Step 7: Safely handle empty peaks ---
    if isempty(prom), MaxProm = NaN; else, MaxProm = max(prom(:),[],'omitnan'); end
    if isempty(pks), MaxHeight = NaN; else, MaxHeight = max(pks(:),[],'omitnan'); end

    %% --- Step 8: Package as one row table ---
    features = table( ...
        An(1), An(2), An(3), RMSval, numel(pks), ...
        MaxProm, MaxHeight, ...
        {d(:)'}, {w(:)'}, ...
        S2, Q4, H6, ...
        'VariableNames', {'A2','A4','A6','Contrast','NumberOfPeaks','MaxProm','MaxHeight','Spacing','Widths','S2','Q4','H6'});
end