%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Script to Calculate CH Inequality
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

format long g;

% define data processing constants
SET_SIZE = 16;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% load data file

data_1 = load('chsh_h_optimise_20141129.dat');
data_2 = load('chsh_h_optimise_20141201.dat');
data_3 = load('chsh_h_optimise_20141202.dat');
data = [data_1; data_2; data_3];

% evaluate the number of complete set
NO_OF_SET = floor(size(data, 1)/SET_SIZE);

% truncating incomplete set
data(NO_OF_SET*SET_SIZE + 1 :end, :) = [];

%data_store = data;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Splitting data set

%index = 1:26*16:length(data)

%for i = 1:12
    
%data = data_store(i:i + 26*16 - 1, :);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%
data_singles = data(:, [2, 3]);
data_2_fold = reshape(data(:, 4), 4, 4*NO_OF_SET);

index = 4:4:length(data_2_fold);
i = 312;
data_2_fold_store = data_2_fold;
data_2_fold = [sum(data_2_fold_store(:, 1:4:index(i)), 2), sum(data_2_fold_store(:, 2:4:index(i)), 2), sum(data_2_fold_store(:, 3:4:index(i)), 2), sum(data_2_fold_store(:, 4:4:index(i)), 2)];

NO_OF_SET = 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% correcting for setting dependent detection efficiencies

correction_factor = load('correction_factor.dat');
delta_correction_factor = load('correction_factor_error.dat');

correction_factor_matrix = kron(reshape(correction_factor, 4, 4), ones(1, NO_OF_SET));
correction_factor_matrix_error = kron(reshape(delta_correction_factor, 4, 4), ones(1, NO_OF_SET));

%data_2_fold = correction_factor_matrix.*data_2_fold;
data_2_fold_error = sqrt(data_2_fold.^2.*correction_factor_matrix_error.^2 + correction_factor_matrix.^2.*sqrt(data_2_fold));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

data_2_fold_total = [sum(data_2_fold, 1); sum(data_2_fold, 1);...
    sum(data_2_fold, 1); sum(data_2_fold, 1)];

data_2_fold_normalised = data_2_fold./data_2_fold_total;

correlation = sum(data_2_fold_normalised([1, 4], :)) - sum(data_2_fold_normalised([2, 3], :));

S = sum([correlation(1:4:end); correlation(3:4:end); correlation(4:4:end)], 1) - correlation(2:4:end);

%end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% error propagation for normalisation
%data_2_fold_error = sqrt(data_2_fold);
data_2_fold_total_error = [sqrt(sum(data_2_fold_error.^2, 1));...
    sqrt(sum(data_2_fold_error.^2, 1)); sqrt(sum(data_2_fold_error.^2, 1));...
    sqrt(sum(data_2_fold_error.^2, 1))];

data_2_fold_normalised_error = data_2_fold_normalised.*sqrt((data_2_fold_error./data_2_fold).^2 + (data_2_fold_total_error./data_2_fold_total).^2);

S_error = sqrt(sum(reshape(data_2_fold_normalised_error, SET_SIZE, NO_OF_SET).^2, 1));

S = abs(S);

mean(S)
mean_S_std = sqrt(sum(S_error.^2))/NO_OF_SET;
sum(sum(data_2_fold))
NO_OF_SET

display_index = 0:NO_OF_SET + 3;
p1 = errorbar(1:NO_OF_SET, S, S_error, 'bx');
hold on;

p2 = plot(display_index, mean(S)*ones(size(display_index)), 'b');
plot(display_index, (mean(S) + mean_S_std)*ones(size(display_index)), 'b');
plot(display_index, (mean(S) - mean_S_std)*ones(size(display_index)), 'b');

p3 = plot(display_index, 2.83*ones(size(display_index)), 'm');
plot(display_index, (2.83 + 0.02)*ones(size(display_index)), 'm');
plot(display_index, (2.83 - 0.02)*ones(size(display_index)), 'm');

p4 = plot(display_index, 2.7252*ones(size(display_index)), 'k');
plot(display_index, (2.7252 + 0.0006)*ones(size(display_index)), 'k');
plot(display_index, (2.7252 - 0.0006)*ones(size(display_index)), 'k');

p5 = plot(display_index, 2*sqrt(2)*ones(size(display_index)), 'r');
grid;
hold off;
ylabel('S Value');
xlabel('No of Measurements');
title('S Value versus No. of Measurements');
legend([p1, p2, p3, p4, p5], {'Measurements', 'Measurement Average', 'Kwait, PRL 111, 130406 (2013)', 'Kwait, Opt. Exp. Vol. 12, No. 22, 8951','Max Theoretical Violation'});

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

no_points = cumsum(sum(reshape(data(:, 4), 16, NO_OF_SET)));

plot(no_points, sqrt(cumsum(S_error.^2))./(1:NO_OF_SET)),...
    xlabel('No. of Events'),...
    ylabel('Standard Error'),...
    title('Standard Error versus No. of Events'),...
    grid;

plot(no_points, abs(cumsum(S)./(1:NO_OF_SET) - mean(S))),...
    xlabel('No. of Events'),...
    ylabel('Deviation From Sample Mean'),...
    title('Deviation From Sample Mean versus No. of Events'),...
    grid;

plot((1:NO_OF_SET)*16, sqrt(cumsum(S_error.^2))./(1:NO_OF_SET)),...
    xlabel('Time Taken (min)'),...
    ylabel('Standard Error'),...
    title('Standard Error versus Time Taken (min)'),...
    grid;

plot((1:NO_OF_SET)*16, abs(cumsum(S)./(1:NO_OF_SET) - mean(S))),...
    xlabel('Time Taken (min)'),...
    ylabel('Absolute Deviation From Sample Mean'),...
    title('Absolute Deviation From Sample Mean versus Time Taken (min)'),...
    grid;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

plot(1:NO_OF_SET, S - mean(S), 'bx'), grid;

%A0 = 0
%A1 = 

B0H = 22.9227 - 22.5;
B0V = 113.0425 - 112.5;
B1H = 67.5669 - 67.5;
B1V = 157.2468 - 157.5;

%2.82762 pm 0.00092
%2.82767 pm 0.00092
%2.82775 pm 0.00091
%2.82768 pm 0.00092
%2.82780 pm 0.00090
%2.82776 pm 0.00086
%2.82763 pm 0.00084
%2.82760 pm 0.00083

div = 2:311;
remainder = rem(312*ones(1, 310), 2:311);

div(find(remainder ~=0)) = [];
div(end) = [];

no_set = div;
S_error = zeros(size(no_set));

for i = 1:length(no_set)
    
  test = reshape(S, no_set(i), 312/no_set(i));
  S_error(i) = std(sum(test')/(312/no_set(i)));
      length(sum(test')/(312/no_set(i)))
end;

plot(no_set, S_error, 'x-'),...
    xlabel('No. of Sets (no of sets 312 measurments are divided into)'),...
    ylabel('delta S'),...
    title('delta S versus No. of Sets'),...    
grid;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% plotting standard error versus number of events

N_tot = 1:100:33184329; % total number of events
N = N_tot/4; % total number of events for evaluation of each correlation function
Delta_S = sqrt(6./N);

semilogy(N_tot, Delta_S),...
    xlabel('No. of Events'),...
    ylabel('\DeltaS'),...
    title('\DeltaS versus No. of Events (Poissonian)'),...
    grid;
axis([N_tot(1) N_tot(end) Delta_S(1) Delta_S(end)]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

cum_total = cumsum(sum(reshape(sum(data_2_fold, 1), 4, 312)));

index = 1:312;

S_error = zeros(length(index), 1);
S_mean = zeros(length(index), 1);
S_mean_man = zeros(length(index), 1);

for i = 1:length(index)
    
    S_error(i) = std(S(1:index(i)));
    S_mean(i) = mean(S(1:index(i)));
    S_error_man(i) = sqrt(sum((S(1:index(i)) - S_mean(i)).^2)/index(i));
    
end;

plot(cum_total, S_error, N_tot, Delta_S),...
    grid;

start = 1;
plot(cum_total(start:end), S_error(start:end)),...
    xlabel('No of Events'),...
    ylabel('Delta S'),...
    title('Delta S versus No of Events'),...
    grid;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

index = 1:4:length(data_2_fold);
S = zeros(size(index));

for i = 1:length(index)
    
    data_1 = sum(data_2_fold(:, (1:4:index(i)) + 0), 2)/(sum(sum(data_2_fold(:, (1:4:index(i)) + 0), 2)));
    data_2 = sum(data_2_fold(:, (1:4:index(i)) + 1), 2)/(sum(sum(data_2_fold(:, (1:4:index(i)) + 1), 2)));
    data_3 = sum(data_2_fold(:, (1:4:index(i)) + 2), 2)/(sum(sum(data_2_fold(:, (1:4:index(i)) + 2), 2)));
    data_4 = sum(data_2_fold(:, (1:4:index(i)) + 3), 2)/(sum(sum(data_2_fold(:, (1:4:index(i)) + 3), 2)));

    e1 = sum(data_1([1, 4])) - sum(data_1([2, 3]));
    e2 = sum(data_2([1, 4])) - sum(data_2([2, 3]));
    e3 = sum(data_3([1, 4])) - sum(data_3([2, 3]));
    e4 = sum(data_4([1, 4])) - sum(data_4([2, 3]));
 
    S(i) = abs(e1 - e2 + e3 + e4);
    
end;

%abs(S + 2.82758700716579)

plot(index, S, 'x-'),...
    xlabel(''),...
    ylabel(''),...
    title(''),...
    grid;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
