New routines to compute PDF, cumulants of order parameter

Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/analyzeFolder.m

@@ -92,8 +92,16 @@ function results = analyzeFolder(options)
         bkg_img  = double(imrotate(h5read(fullFileName, append(groupList(cam), "/background")), angle));
         dark_img = double(imrotate(h5read(fullFileName, append(groupList(cam), "/dark")), angle));
         
-        refimages(:,:,k)  = subtractBackgroundOffset(cropODImage(bkg_img, center, span), fraction)';
+        if (isempty(atm_img)  && isa(atm_img, 'double')) || ...
-        absimages(:,:,k)  = subtractBackgroundOffset(cropODImage(calculateODImage(atm_img, bkg_img, dark_img, ImagingMode, PulseDuration), center, span), fraction)';
+           (isempty(bkg_img)  && isa(bkg_img, 'double')) || ...
+           (isempty(dark_img) && isa(dark_img, 'double'))
+            
+            refimages = nan(size(refimages));  % fill with NaNs
+            absimages = nan(size(absimages));
+        else
+            refimages(:,:,k)  = subtractBackgroundOffset(cropODImage(bkg_img, center, span), fraction)';
+            absimages(:,:,k)  = subtractBackgroundOffset(cropODImage(calculateODImage(atm_img, bkg_img, dark_img, ImagingMode, PulseDuration), center, span), fraction)';
+        end
     end
     
     % ===== Fringe removal =====

Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/bootstrapCumulants.m

@@ -0,0 +1,41 @@
+function [cumulants_mean, cumulants_ci, bootstrap_samples] = bootstrapCumulants(x, maxOrder, nBoot)
+% bootstrapCumulants - compute bootstrap estimates of cumulants and confidence intervals
+%
+% Syntax:
+%   [meanC, ciC, allC] = bootstrapCumulants(x, maxOrder, nBoot)
+%
+% Inputs:
+%   x        - Data vector (may contain NaNs)
+%   maxOrder - Max cumulant order (default: 6)
+%   nBoot    - Number of bootstrap samples (default: 1000)
+%
+% Outputs:
+%   cumulants_mean    - Mean of bootstrap cumulants
+%   cumulants_ci      - 95% confidence intervals [2.5th; 97.5th] percentile
+%   bootstrap_samples - All bootstrap cumulants (nBoot x maxOrder)
+
+    if nargin < 2, maxOrder = 6; end
+    if nargin < 3, nBoot = 1000; end
+    
+    x = x(:);
+    x = x(~isnan(x)); % Remove NaNs
+    
+    if isempty(x)
+        cumulants_mean = NaN(1, maxOrder);
+        cumulants_ci = NaN(2, maxOrder);
+        bootstrap_samples = NaN(nBoot, maxOrder);
+        return;
+    end
+    
+    N = numel(x);
+    bootstrap_samples = zeros(nBoot, maxOrder);
+    
+    for b = 1:nBoot
+        xb = x(randi(N, [N, 1])); % Resample with replacement
+        bootstrap_samples(b, :) = computeCumulants(xb, maxOrder);
+    end
+    
+    cumulants_mean = mean(bootstrap_samples, 1);
+    cumulants_ci = prctile(bootstrap_samples, [2.5, 97.5]);
+
+end

Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/compareAngularCorrelation.m

@@ -8,13 +8,13 @@ format long
 font          = 'Bahnschrift';
 
 % Load data
-Data                       = load('C:/Users/Karthik/Documents/GitRepositories/Calculations/Data-Analyzer/Comparison/Max_g2_DropletsToStripes.mat', 'unique_scan_parameter_values', 'mean_max_g2_values', 'std_error_g2_values');
+Data                       = load('C:/Users/Karthik/Documents/GitRepositories/Calculations/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/Max_g2_DropletsToStripes.mat', 'unique_scan_parameter_values', 'mean_max_g2_values', 'std_error_g2_values');
 
 dts_scan_parameter_values  = Data.unique_scan_parameter_values;
 dts_mean_mg2               = Data.mean_max_g2_values;
 dts_stderr_mg2             = Data.std_error_g2_values;
 
-Data                       = load('C:/Users/Karthik/Documents/GitRepositories/Calculations/Data-Analyzer/Comparison/Max_g2_StripesToDroplets.mat', 'unique_scan_parameter_values', 'mean_max_g2_values', 'std_error_g2_values');
+Data                       = load('C:/Users/Karthik/Documents/GitRepositories/Calculations/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/Max_g2_StripesToDroplets.mat', 'unique_scan_parameter_values', 'mean_max_g2_values', 'std_error_g2_values');
 
 std_scan_parameter_values  = Data.unique_scan_parameter_values;
 std_mean_mg2               = Data.mean_max_g2_values;
@@ -30,10 +30,10 @@ errorbar(std_scan_parameter_values, std_mean_mg2, std_stderr_mg2, 'o--', ...
 set(gca, 'FontSize', 14, 'YLim', [0, 1]);
 hXLabel = xlabel('$\alpha$ (degrees)', 'Interpreter', 'latex');
 hYLabel = ylabel('$\mathrm{max}[g^{(2)}_{[50,70]}(\delta\theta)]$', 'Interpreter', 'latex');
-hTitle = title('B = 2.42 G', 'Interpreter', 'tex');
+% hTitle = title('B = 2.42 G', 'Interpreter', 'tex');
 legend
 set([hXLabel, hYLabel], 'FontName', font)
 set([hXLabel, hYLabel], 'FontSize', 14)
-set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold');  % Set font and size for title
+% set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold');  % Set font and size for title
 grid on
 %%

Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/computeCumulants.m

@@ -0,0 +1,52 @@
+function cumulants = computeCumulants(x, maxOrder)
+% computeCumulants - compute cumulants up to specified order from data vector x
+%
+% Syntax: cumulants = computeCumulants(x, maxOrder)
+%
+% Inputs:
+%   x        - 1D numeric vector (may contain NaNs)
+%   maxOrder - maximum order of cumulants to compute (default: 6)
+%
+% Output:
+%   cumulants - vector [kappa_1, ..., kappa_maxOrder]
+
+    if nargin < 2
+        maxOrder = 6;
+    end
+    
+    x = x(:);
+    x = x(~isnan(x)); % Remove NaNs
+    
+    if isempty(x)
+        cumulants = NaN(1, maxOrder);
+        return;
+    end
+    
+    mu1 = mean(x, 'omitnan');
+    x_centered = x - mu1;
+    
+    cumulants = zeros(1, maxOrder);
+    cumulants(1) = mu1;
+    
+    mu = zeros(1, maxOrder);
+    for k = 2:maxOrder
+        mu(k) = mean(x_centered.^k, 'omitnan');
+    end
+    
+    if maxOrder >= 2
+        cumulants(2) = mu(2);
+    end
+    if maxOrder >= 3
+        cumulants(3) = mu(3);
+    end
+    if maxOrder >= 4
+        cumulants(4) = mu(4) - 3 * mu(2)^2;
+    end
+    if maxOrder >= 5
+        cumulants(5) = mu(5) - 10 * mu(3) * mu(2);
+    end
+    if maxOrder >= 6
+        cumulants(6) = mu(6) - 15 * mu(4) * mu(2) - 10 * mu(3)^2 + 30 * mu(2)^3;
+    end
+
+end

Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/conductSpectralAnalysis.m

@@ -3,9 +3,9 @@ groupList          = ["/images/MOT_3D_Camera/in_situ_absorption", "/images/ODT_1
                       "/images/ODT_2_Axis_Camera/in_situ_absorption", "/images/Horizontal_Axis_Camera/in_situ_absorption", ...
                       "/images/Vertical_Axis_Camera/in_situ_absorption"];
 
-folderPath         = "//DyLabNAS/Data/TwoDGas/2025/06/23/";
+folderPath         = "//DyLabNAS/Data/TwoDGas/2025/07/22/";
 
-run                = '0300';
+run                = '0021';
 
 folderPath         = strcat(folderPath, run);
 
@@ -50,32 +50,32 @@ savefileName                = 'DropletsToStripes';
 font                        = 'Bahnschrift';
 
 if strcmp(savefileName, 'DropletsToStripes')
-    scan_groups           = 0:5:45;
+    scan_groups           = 0:1:40;
     titleString           = 'Droplets to Stripes';
 elseif strcmp(savefileName, 'StripesToDroplets')
-    scan_groups           = 45:-5:0;
+    scan_groups           = 40:-1:0;
     titleString           = 'Stripes to Droplets';
 end
 
 % Flags
-skipNormalization           = false;
+skipNormalization           = true;
-skipUnshuffling             = true;
+skipUnshuffling             = false;
 skipPreprocessing           = true;
 skipMasking                 = true;
 skipIntensityThresholding   = true;
 skipBinarization            = true;
 skipMovieRender             = true;
-skipSaveFigures             = false;
+skipSaveFigures             = true;
-skipSaveOD                  = false;
+skipSaveOD                  = true;
 
 %% ===== S-D Settings =====
 groupList          = ["/images/MOT_3D_Camera/in_situ_absorption", "/images/ODT_1_Axis_Camera/in_situ_absorption", ...
                       "/images/ODT_2_Axis_Camera/in_situ_absorption", "/images/Horizontal_Axis_Camera/in_situ_absorption", ...
                       "/images/Vertical_Axis_Camera/in_situ_absorption"];
 
-folderPath         = "//DyLabNAS/Data/TwoDGas/2025/06/24/";
+folderPath         = "//DyLabNAS/Data/TwoDGas/2025/07/23/";
 
-run                = '0001';
+run                = '0055';
 
 folderPath         = strcat(folderPath, run);
 
@@ -120,10 +120,10 @@ savefileName                = 'StripesToDroplets';
 font                        = 'Bahnschrift';
 
 if strcmp(savefileName, 'DropletsToStripes')
-    scan_groups           = 0:5:45;
+    scan_groups           = 0:1:40
     titleString           = 'Droplets to Stripes';
 elseif strcmp(savefileName, 'StripesToDroplets')
-    scan_groups           = 45:-5:0;
+    scan_groups           = 40:-1:0;
     titleString           = 'Stripes to Droplets';
 end
 
@@ -135,8 +135,8 @@ skipMasking                 = true;
 skipIntensityThresholding   = true;
 skipBinarization            = true;
 skipMovieRender             = true;
-skipSaveFigures             = false;
+skipSaveFigures             = true;
-skipSaveOD                  = false;
+skipSaveOD                  = true;
 
 %% ===== Load and compute OD image, rotate and extract ROI for analysis =====
 % Get a list of all files in the folder with the desired file name pattern.
@@ -156,8 +156,16 @@ for k = 1 : length(files)
     bkg_img  = double(imrotate(h5read(fullFileName, append(groupList(cam), "/background")), angle));
     dark_img = double(imrotate(h5read(fullFileName, append(groupList(cam), "/dark")), angle));
     
-    refimages(:,:,k)  = subtractBackgroundOffset(cropODImage(bkg_img, center, span), fraction)';
+    if (isempty(atm_img)  && isa(atm_img, 'double')) || ...
-    absimages(:,:,k)  = subtractBackgroundOffset(cropODImage(calculateODImage(atm_img, bkg_img, dark_img, ImagingMode, PulseDuration), center, span), fraction)';
+       (isempty(bkg_img)  && isa(bkg_img, 'double')) || ...
+       (isempty(dark_img) && isa(dark_img, 'double'))
+        
+        refimages = nan(size(refimages));  % fill with NaNs
+        absimages = nan(size(absimages));
+    else
+        refimages(:,:,k)  = subtractBackgroundOffset(cropODImage(bkg_img, center, span), fraction)';
+        absimages(:,:,k)  = subtractBackgroundOffset(cropODImage(calculateODImage(atm_img, bkg_img, dark_img, ImagingMode, PulseDuration), center, span), fraction)';
+    end
 end
 
 % ===== Fringe removal =====

Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/extractAutocorrelation.m

@@ -86,8 +86,16 @@ for k = 1 : length(files)
     bkg_img  = double(imrotate(h5read(fullFileName, append(groupList(cam), "/background")), angle));
     dark_img = double(imrotate(h5read(fullFileName, append(groupList(cam), "/dark")), angle));
     
-    refimages(:,:,k)  = subtractBackgroundOffset(cropODImage(bkg_img, center, span), fraction)';
+    if (isempty(atm_img)  && isa(atm_img, 'double')) || ...
-    absimages(:,:,k)  = subtractBackgroundOffset(cropODImage(calculateODImage(atm_img, bkg_img, dark_img, ImagingMode, PulseDuration), center, span), fraction)';
+       (isempty(bkg_img)  && isa(bkg_img, 'double')) || ...
+       (isempty(dark_img) && isa(dark_img, 'double'))
+        
+        refimages = nan(size(refimages));  % fill with NaNs
+        absimages = nan(size(absimages));
+    else
+        refimages(:,:,k)  = subtractBackgroundOffset(cropODImage(bkg_img, center, span), fraction)';
+        absimages(:,:,k)  = subtractBackgroundOffset(cropODImage(calculateODImage(atm_img, bkg_img, dark_img, ImagingMode, PulseDuration), center, span), fraction)';
+    end
 end
 
 %% ===== Fringe removal =====
@@ -209,8 +217,8 @@ for i = 1:N_params
 
             temp(j) = num / denom;
         end
-        g2_all(i, dtheta+1)       = mean(temp);
+        g2_all(i, dtheta+1)       = mean(temp, 'omitnan');
-        g2_error_all(i, dtheta+1) = std(temp) / sqrt(length(group_idx));  % Standard error
+        g2_error_all(i, dtheta+1) = std(temp, 'omitnan') / sqrt(length(group_idx));  % Standard error
     end
 end
 
@@ -222,7 +230,7 @@ cmap = sky(nParams);  % You can also try 'jet', 'turbo', 'hot', etc.
 
 figure(1); 
 clf;
-set(gcf,'Position',[100 100 950 750])
+set(gcf, 'Color', 'w', 'Position',[100 100 950 750])
 hold on;
 legend_entries = cell(nParams, 1);
 

Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/extractCustomCorrelation.m

@@ -3,9 +3,9 @@ groupList          = ["/images/MOT_3D_Camera/in_situ_absorption", "/images/ODT_1
                       "/images/ODT_2_Axis_Camera/in_situ_absorption", "/images/Horizontal_Axis_Camera/in_situ_absorption", ...
                       "/images/Vertical_Axis_Camera/in_situ_absorption"];
 
-folderPath         = "//DyLabNAS/Data/TwoDGas/2025/06/23/";
+folderPath         = "//DyLabNAS/Data/TwoDGas/2025/07/22/";
 
-run                = '0300';
+run                = '0021';
 
 folderPath         = strcat(folderPath, run);
 
@@ -86,8 +86,16 @@ for k = 1 : length(files)
     bkg_img  = double(imrotate(h5read(fullFileName, append(groupList(cam), "/background")), angle));
     dark_img = double(imrotate(h5read(fullFileName, append(groupList(cam), "/dark")), angle));
     
-    refimages(:,:,k)  = subtractBackgroundOffset(cropODImage(bkg_img, center, span), fraction)';
+    if (isempty(atm_img)  && isa(atm_img, 'double')) || ...
-    absimages(:,:,k)  = subtractBackgroundOffset(cropODImage(calculateODImage(atm_img, bkg_img, dark_img, ImagingMode, PulseDuration), center, span), fraction)';
+       (isempty(bkg_img)  && isa(bkg_img, 'double')) || ...
+       (isempty(dark_img) && isa(dark_img, 'double'))
+        
+        refimages = nan(size(refimages));  % fill with NaNs
+        absimages = nan(size(absimages));
+    else
+        refimages(:,:,k)  = subtractBackgroundOffset(cropODImage(bkg_img, center, span), fraction)';
+        absimages(:,:,k)  = subtractBackgroundOffset(cropODImage(calculateODImage(atm_img, bkg_img, dark_img, ImagingMode, PulseDuration), center, span), fraction)';
+    end
 end
 
 %% ===== Fringe removal =====
@@ -176,7 +184,7 @@ end
 
 % Group by scan parameter values (e.g., alpha, angle, etc.)
 [unique_scan_parameter_values, ~, idx] = unique(scan_parameter_values);
-N_params = length(unique_scan_parameter_values);
+N_params                               = length(unique_scan_parameter_values);
 
 % Define angular range and conversion
 angle_range     = 180;
@@ -189,23 +197,24 @@ window_size     = 10;
 angle_threshold = 100;
 
 % Initialize containers for final results
-mean_max_g2_values       = zeros(1, N_params);
+mean_max_g2_values                       = zeros(1, N_params);
-mean_max_g2_angle_values = zeros(1, N_params);
+skew_max_g2_angle_values                 = zeros(1, N_params);
-var_max_g2_values        = zeros(1, N_params);
+var_max_g2_values                        = zeros(1, N_params);
-var_max_g2_angle_values  = zeros(1, N_params);
+kurt_max_g2_angle_values                 = zeros(1, N_params);
+fifth_order_cumulant_max_g2_angle_values = zeros(1, N_params);
+sixth_order_cumulant_max_g2_angle_values = zeros(1, N_params);
 
 % Also store raw data per group
-g2_all_per_group         = cell(1, N_params);
+max_g2_all_per_group                     = cell(1, N_params);
-angle_all_per_group      = cell(1, N_params);
+std_error_g2_values                      = zeros(1, N_params);
 
 for i = 1:N_params
-    group_idx      = find(idx == i);
+    group_idx                  = find(idx == i);
-    group_data     = delta_nkr_all(group_idx, :);
+    group_data                 = delta_nkr_all(group_idx, :);
-    N_reps         = size(group_data, 1);
+    N_reps                     = size(group_data, 1);
-
-    g2_values         = zeros(1, N_reps);
-    angle_at_max_g2   = zeros(1, N_reps);
 
+    g2_values                  = zeros(1, N_reps);
+    
     for j = 1:N_reps
         profile                = group_data(j, :);
 
@@ -225,8 +234,7 @@ for i = 1:N_params
         ref                    = profile(ref_window);
 
         correlations           = zeros(size(offsets));
-        angles                 = zeros(size(offsets));
+        
-
         for k = 1:length(offsets)
             shifted_idx        = mod(peak_idx + offsets(k) - 1, N_angular_bins) + 1;
             sec_window         = mod((shifted_idx - window_size):(shifted_idx + window_size) - 1, N_angular_bins) + 1;
@@ -237,36 +245,191 @@ for i = 1:N_params
             g2                 = num / denom;
 
             correlations(k)    = g2;
-            angles(k)          = mod((peak_idx - 1 + offsets(k)) * angle_per_bin, angle_range);
         end
 
         [max_corr, max_idx]    = max(correlations);
         g2_values(j)           = max_corr;
-        angle_at_max_g2(j)     = angles(max_idx);
     end
 
     % Store raw values
-    g2_all_per_group{i}      = g2_values;
+    max_g2_all_per_group{i}        = g2_values;
-    angle_all_per_group{i}   = angle_at_max_g2;
+    
+    % Compute cumulants
+    kappa                      = computeCumulants(g2_values(:), 6);
 
     % Final stats
-    mean_max_g2_values(i)      = mean(g2_values, 'omitnan');
+    mean_max_g2_values(i)                       = kappa(1);
-    var_max_g2_values(i)       = var(g2_values, 0, 'omitnan');
+    var_max_g2_values(i)                        = kappa(2);
-    mean_max_g2_angle_values(i)= mean(angle_at_max_g2, 'omitnan');
+    
-    var_max_g2_angle_values(i) = var(angle_at_max_g2, 0, 'omitnan');
+    N_eff                                       = sum(~isnan(g2_values));
+    std_error_g2_values(i)                      = sqrt(kappa(2)) / sqrt(N_eff);
+
+    skew_max_g2_angle_values(i)                 = kappa(3);
+    kurt_max_g2_angle_values(i)                 = kappa(4);
+    fifth_order_cumulant_max_g2_angle_values(i) = kappa(5);
+    sixth_order_cumulant_max_g2_angle_values(i) = kappa(6);
+
 end
 
-%% ── Mean ± Std vs. scan parameter ──────────────────────────────────────
+%% Plot PDF of order parameter
-% Compute standard error instead of standard deviation
+
-std_error_g2_values = zeros(1, N_params);
+if ~skipSaveFigures
-for i = 1:N_params
+    % Define folder for saving images
-    n_i = numel(g2_all_per_group{i});  % Number of repetitions for this param
+    saveFolder = [savefileName '_SavedFigures'];
-    std_error_g2_values(i) = sqrt(var_max_g2_values(i) / n_i);
+    if ~exist(saveFolder, 'dir')
+        mkdir(saveFolder);
+    end
 end
 
+figure(2);  % one persistent figure
+set(gcf, 'Color', 'w', 'Position', [100 100 950 750])
+
+for val = scan_groups
+    % Find the index i that matches this scan parameter value
+    i = find(unique_scan_parameter_values == val, 1);
+
+    % Skip if not found (sanity check)
+    if isempty(i)
+        continue;
+    end
+
+    g2_vals = g2_all_per_group{i};
+    g2_vals = g2_vals(~isnan(g2_vals));
+
+    if isempty(g2_vals)
+        continue;
+    end
+
+    % KDE estimation
+    [f, xi] = ksdensity(g2_vals, 'NumPoints', 200);
+
+    clf;
+    histogram(g2_vals, 'Normalization', 'pdf', ...
+        'NumBins', 10, ...
+        'FaceAlpha', 0.3, ...
+        'EdgeColor', 'none', ...
+        'FaceColor', [0.3 0.5 0.8]);
+
+    hold on;
+    plot(xi, f, 'LineWidth', 2, 'Color', [0 0.2 0.6]);
+
+    set(gca, 'FontSize', 16);
+    title(sprintf('%s: \\boldmath$\\alpha = %.1f^{\\circ}$', titleString, val), ...
+        'FontSize', 16, 'Interpreter', 'latex');
+
+    xlabel('$\mathrm{max}[g^{(2)}_{[50,70]}(\delta\theta)]$', 'Interpreter', 'latex', 'FontSize', 14);
+    ylabel('PDF', 'FontSize', 14);
+    xlim([0.0, 1.5]);
+    grid on;
+
+    drawnow;
+
+    % ==== Save Figure ====
+    if ~skipSaveFigures
+        % Create a filename for each averaged plot
+        fileNamePNG = fullfile(saveFolder, sprintf('max_g2_analysis_param_%03d.png', val));
+        
+        % Save current figure as PNG with high resolution
+        print(gcf, fileNamePNG, '-dpng', '-r300');  % 300 dpi for high quality
+    else
+        pause(0.5)
+    end  
+end
+
+
+%% Plot all cumulants
+figure(3)
+set(gcf, 'Color', 'w', 'Position', [100 100 950 750])
+
+scan_vals = unique_scan_parameter_values;
+
+% Define font style for consistency
+axis_fontsize  = 14;
+label_fontsize = 16;
+title_fontsize = 16;
+
+% 1. Mean with error bars
+subplot(3,2,1);
+errorbar(scan_vals, mean_max_g2_values, std_error_g2_values, 'o-', ...
+    'LineWidth', 1.5, 'MarkerSize', 6);
+title('Mean of $\mathrm{max}[g^{(2)}_{[50,70]}(\delta\theta)]$', ...
+    'Interpreter', 'latex', 'FontSize', title_fontsize);
+xlabel('$\alpha$ (degrees)', 'Interpreter', 'latex', ...
+    'FontSize', label_fontsize);
+ylabel('$\kappa_1$', 'Interpreter', 'latex', ...
+    'FontSize', label_fontsize);
+set(gca, 'FontSize', axis_fontsize);
+grid on;
+
+% 2. Variance
+subplot(3,2,2);
+plot(scan_vals, var_max_g2_values, 's-', 'LineWidth', 1.5, 'MarkerSize', 6);
+title('Variance of $\mathrm{max}[g^{(2)}_{[50,70]}(\delta\theta)]$', ...
+    'Interpreter', 'latex', 'FontSize', title_fontsize);
+xlabel('$\alpha$ (degrees)', 'Interpreter', 'latex', ...
+    'FontSize', label_fontsize);
+ylabel('$\kappa_2$', 'Interpreter', 'latex', ...
+    'FontSize', label_fontsize);
+set(gca, 'FontSize', axis_fontsize);
+grid on;
+
+% 3. Skewness
+subplot(3,2,3);
+plot(scan_vals, skew_max_g2_angle_values, 'd-', 'LineWidth', 1.5, 'MarkerSize', 6);
+title('Skewness of $\mathrm{max}[g^{(2)}_{[50,70]}(\delta\theta)]$', ...
+    'Interpreter', 'latex', 'FontSize', title_fontsize);
+xlabel('$\alpha$ (degrees)', 'Interpreter', 'latex', ...
+    'FontSize', label_fontsize);
+ylabel('$\kappa_3$', 'Interpreter', 'latex', ...
+    'FontSize', label_fontsize);
+set(gca, 'FontSize', axis_fontsize);
+grid on;
+
+% 4. Kurtosis
+subplot(3,2,4);
+plot(scan_vals, kurt_max_g2_angle_values, '^-', 'LineWidth', 1.5, 'MarkerSize', 6);
+title('Kurtosis of $\mathrm{max}[g^{(2)}_{[50,70]}(\delta\theta)]$', ...
+    'Interpreter', 'latex', 'FontSize', title_fontsize);
+xlabel('$\alpha$ (degrees)', 'Interpreter', 'latex', ...
+    'FontSize', label_fontsize);
+ylabel('$\kappa_4$', 'Interpreter', 'latex', ...
+    'FontSize', label_fontsize);
+set(gca, 'FontSize', axis_fontsize);
+grid on;
+
+% 5. 5th-order cumulant
+subplot(3,2,5);
+plot(scan_vals, fifth_order_cumulant_max_g2_angle_values, 'v-', 'LineWidth', 1.5, 'MarkerSize', 6);
+title('Fifth-order cumulant of $\mathrm{max}[g^{(2)}_{[50,70]}(\delta\theta)]$', ...
+    'Interpreter', 'latex', 'FontSize', title_fontsize);
+xlabel('$\alpha$ (degrees)', 'Interpreter', 'latex', ...
+    'FontSize', label_fontsize);
+ylabel('$\kappa_5$', 'Interpreter', 'latex', ...
+    'FontSize', label_fontsize);
+set(gca, 'FontSize', axis_fontsize);
+grid on;
+
+% 6. 6th-order cumulant
+subplot(3,2,6);
+plot(scan_vals, sixth_order_cumulant_max_g2_angle_values, '>-', 'LineWidth', 1.5, 'MarkerSize', 6);
+title('Sixth-order cumulant of $\mathrm{max}[g^{(2)}_{[50,70]}(\delta\theta)]$', ...
+    'Interpreter', 'latex', 'FontSize', title_fontsize);
+xlabel('$\alpha$ (degrees)', 'Interpreter', 'latex', ...
+    'FontSize', label_fontsize);
+ylabel('$\kappa_6$', 'Interpreter', 'latex', ...
+    'FontSize', label_fontsize);
+set(gca, 'FontSize', axis_fontsize);
+grid on;
+
+% Super title
+sgtitle(sprintf('Cumulants of Peak Offset Angular Correlation - %s', titleString), ...
+    'FontWeight', 'bold', 'FontSize', 16, 'Interpreter', 'latex');
+
+%% ── Mean ± Std vs. scan parameter ──────────────────────────────────────
+
 % Plot mean ± SEM
 figure(1);
-set(gcf,'Position',[100 100 950 750])
+set(gcf, 'Color', 'w', 'Position',[100 100 950 750])
 set(gca, 'FontSize', 14);   % For tick labels only
 errorbar(unique_scan_parameter_values, ...       % x-axis
          mean_max_g2_values, ...                 % y-axis (mean)
@@ -289,75 +452,6 @@ if ~exist(saveFolder, 'dir')
 end
 save([saveFolder savefileName '.mat'], 'unique_scan_parameter_values', 'mean_max_g2_values', 'std_error_g2_values');
 
-%{
-%% Plot histograms within 0-180 degrees only
-figure(1);
-hold on;
-
-bin_edges = 0:10:180;
-
-h1 = histogram(ref_peak_angles, 'BinEdges', bin_edges, ...
-    'FaceColor', [0.3 0.7 0.9], 'EdgeColor', 'none', 'FaceAlpha', 0.6);
-
-h2 = histogram(angle_at_max_g2, 'BinEdges', bin_edges, ...
-    'FaceColor', [0.9 0.4 0.4], 'EdgeColor', 'none', 'FaceAlpha', 0.6);
-
-h1.Normalization = 'probability';
-h2.Normalization = 'probability';
-
-xlabel('Angle (degrees)', 'FontSize', 12);
-ylabel('Probability', 'FontSize', 12);
-legend({'Reference Peak Angle', 'Angle at Max g₂'}, 'FontSize', 12);
-title('Comparison of Reference Peak and Max g₂ Angles', 'FontSize', 14);
-grid on;
-xlim([0 180]);
-hold off;
-
-% Assume ref_peak_angles and angle_at_max_g2 are row or column vectors of angles in [0,180]
-
-% Define fine angle grid for KDE evaluation
-angle_grid = linspace(0, 180, 1000);
-
-% KDE for reference peak angles
-[f_ref, xi_ref] = ksdensity(ref_peak_angles, angle_grid, 'Bandwidth', 5);
-
-% KDE for max g2 angles
-[f_g2, xi_g2] = ksdensity(angle_at_max_g2, angle_grid, 'Bandwidth', 5);
-
-% Plot KDEs
-figure(2);
-plot(xi_ref, f_ref, 'LineWidth', 2, 'DisplayName', 'Reference Peak Angles');
-hold on;
-plot(xi_g2, f_g2, 'LineWidth', 2, 'DisplayName', 'Max g_2 Angles');
-xlabel('Angle (degrees)');
-ylabel('Probability Density');
-title('KDE of Angle Distributions');
-legend;
-grid on;
-
-% Find modes (angle at max KDE value)
-[~, mode_idx_ref] = max(f_ref);
-mode_ref = xi_ref(mode_idx_ref);
-
-[~, mode_idx_g2] = max(f_g2);
-mode_g2 = xi_g2(mode_idx_g2);
-
-% Calculate difference in mode
-mode_diff = abs(mode_ref - mode_g2);
-fprintf('Mode difference between distributions: %.2f degrees\n', mode_diff);
-
-% Add vertical dashed lines at mode positions
-yl = ylim; % get y-axis limits for text positioning
-
-% Reference peak mode line and label
-xline(mode_ref, 'k--', 'LineWidth', 1.5, 'DisplayName', sprintf('Ref Mode: %.1f°', mode_ref));
-text(mode_ref, yl(2)*0.9, sprintf('%.1f°', mode_ref), 'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom', 'FontSize', 12, 'Color', 'k');
-
-% Max g2 mode line and label
-xline(mode_g2, 'r--', 'LineWidth', 1.5, 'DisplayName', sprintf('g_2 Mode: %.1f°', mode_g2));
-text(mode_g2, yl(2)*0.75, sprintf('%.1f°', mode_g2), 'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom', 'FontSize', 12, 'Color', 'r');
-%}
-
 %% Helper Functions
 function [IMGFFT, IMGPR] = computeFourierTransform(I, skipPreprocessing, skipMasking, skipIntensityThresholding, skipBinarization)
     % computeFourierSpectrum - Computes the 2D Fourier power spectrum