function []=figureS2b()
close all
figure()
hold on
%% Data parameters
conversion_factor=1/1.00663684;
short_distance1=3.5;
short_distance2=5.5;
max_fitting_distance1=15;
max_fitting_distance2=15;
Sampledividing=204;
sample_indices=[200,204,205,209,214,219,223,227];
NS=size(sample_indices,2);
xi6s=zeros(NS,1);
xi6_errors=zeros(NS,1);

%%
Ts=zeros(NS,1);
xi_guess=1;
for sample_index_dummy=1:NS
    sample_index=sample_indices(sample_index_dummy);
    % color
    colorB=(sample_index_dummy-1)/(NS-1);
    colorR=1-colorB;
    color=[colorR,0,colorB];
    %% Reading
    % gr
    data_name=sprintf('gr_%d.dat',sample_index);
    fid = fopen(data_name,'r'); numbers=fscanf(fid, '%f'); fclose('all');
    T=numbers(1);
    Nr=numbers(2);
    numbers=numbers(3:end);
    numbers=reshape(numbers,2,Nr).';
    Ts(sample_index_dummy)=T;
    r=numbers(:,1)*conversion_factor;
    gr=numbers(:,2);
    % gl(with l=angular momentum)
    data_name=sprintf('g6_%d.dat',sample_index);
    fid = fopen(data_name,'r'); numbers=fscanf(fid, '%f'); fclose('all');
    T=numbers(1);
    Nr=numbers(2);
    numbers=numbers(3:end);
    numbers=reshape(numbers,2,Nr).';
    Ts(sample_index_dummy)=T;
    gBO=numbers(:,2);
    bx=r;
    by=gBO./gr;
    plot(bx,by,'Color',color)
    
    %% Fitting
    if sample_index<=Sampledividing
        short_distance=short_distance1;
        max_fitting_distance=max_fitting_distance1;
    else
        short_distance=short_distance2;
        max_fitting_distance=max_fitting_distance2;
    end
    x=short_distance:0.01:max_fitting_distance;
    [f,xi_fit,xi_error,bxpeak,bypeak]=exponential_fit_length(bx,by,x,short_distance,max_fitting_distance,xi_guess);
    xi_guess=xi_fit;
    plot(x,f,'Color',color)
    scatter(bxpeak,bypeak,100,color,'filled')
    xi6s(sample_index_dummy)=xi_fit;
    xi6_errors(sample_index_dummy)=1.96*xi_error;

end
set(gca,'yscale','log')
xlim([0,15])
ylim([10^(-3),0.2])
xlabel('$r$','Interpreter','latex');
ylabel('$g_6(r)/g(r)$','Interpreter','latex');
print(gcf, '-dpdf', 'figureS2b.pdf');

end

function [y_fit,xi_fit,xi_error,bxpeak,bypeak]=exponential_fit_length(bx,by,x_interest,short_distance,max_fitting_distance,xi_guess)
surrounding=0.5;
peak_threshold=10^(-3);
Nbx=size(bx,1);
%% First selection of local maxima
I_localmaxima=zeros(Nbx,1);
for jj=2:(Nbx-1)
    if by(jj)>by(jj-1) && by(jj)>by(jj+1) && by(jj)>peak_threshold && bx(jj)>short_distance && bx(jj)<max_fitting_distance
        I_localmaxima(jj)=1;
    end
end
I_localmaxima=logical(I_localmaxima);
bxpeak=bx(I_localmaxima);
bypeak=by(I_localmaxima);
Npeak_pre=size(bxpeak,1);
%% Fluctuation killer once
Npeak_after=Npeak_pre-100;
while(Npeak_after~=Npeak_pre)
    Npeak_pre=size(bxpeak,1);
    for jj=2:Npeak_pre
        if (bxpeak(jj)-bxpeak(jj-1))<surrounding
            if bypeak(jj)>bypeak(jj-1)
                bxpeak(jj-1)=bxpeak(jj);
                bypeak(jj-1)=bypeak(jj);
            else
                if jj<Npeak_pre
                    bxpeak(jj)=bxpeak(jj+1);
                    bypeak(jj)=bypeak(jj+1);
                end
            end
        end
    end
    bxpeak=unique(bxpeak,'stable');
    bypeak=unique(bypeak,'stable');
    Npeak_after=size(bxpeak,1);
end

%% Exponential fit
weight=diag(ones(size(bxpeak)));
Outliers=logical(zeros(size(bxpeak)));
Init=[1 xi_guess];
[fit_values,confidence_intervals]=exponential_decay_fit(bxpeak,bypeak,weight,Init,Outliers);
a=fit_values(1);
b=fit_values(2);
c=1;
y_fit=a*exp(-(x_interest/b).^c);
xi_fit=b;
xi_error=(confidence_intervals(2,2)-confidence_intervals(1,2))/2;
end