%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Script to Calculate CH Inequality
%(actual script used to obtain value in write up)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% set number display format
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];

%% reshaping data into columns ++; -+; +-;--
% rows are a sequence of a0b0, a0b1, a1b0, a1b1 times the number of
% experimental runs

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

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

% single counts for A and B
singles_A = reshape(data(:, 2), 4, 4*NO_OF_SET);
singles_B = reshape(data(:, 3), 4, 4*NO_OF_SET);

% average single rates
ra = mean(data(:, 2)/60);
rb = mean(data(:, 3)/60);
ra_err = sqrt(sum(data(:, 2)/60^2))/length(data(:, 2));
rb_err = sqrt(sum(data(:, 3)/60^2))/length(data(:, 3));

ra_err = std(data(:, 2)/60);
rb_err = std(data(:, 3)/60);

acc = ra*rb*1.2e-9;
acc_err = 1.2e-9*sqrt((ra*rb_err)^2+(rb*ra_err)^2);

fprintf('Calculated accidental rate: %f +/- %f s-1\n', acc, acc_err)

% coincidence counts
coinc = reshape(data(:, 4), 4, 4*NO_OF_SET);
rc_min = min(min(coinc))/60;
rc_max = max(max(coinc))/60;

%% calculating S
% getting normalization for all elements
coinc_total = [sum(coinc, 1);...
                    sum(coinc, 1);...
                    sum(coinc, 1);...
                    sum(coinc, 1)];

% normalization of the counts
coinc_normalised = coinc./coinc_total;

% correlations for every setting and every experimental run
correlation = sum(coinc_normalised([1, 4], :)) - sum(coinc_normalised([2, 3], :));

% calculating S for every experimental run
S = abs(sum([correlation(1:4:end); correlation(3:4:end); correlation(4:4:end)], 1) - correlation(2:4:end));

% final value of S
S_value = mean(S);

% distance of the measured S from the tsirelson's bound
Delta_T = 2*sqrt(2)- S_value;

% histogram of the distribution of the distance of S from the tsirelson' bound
% for experimental run
histfit(S-2*sqrt(2), 30);

%% error propagation and normalisation

% Error coming from the poissonian statistics of the process
poissonian_error = sqrt(coinc);

% Error due to the jitter of the total acquisition time of the counter
counter_err = coinc * 1e-7/60;

% Error due to the dead time of the detectors
dead_err = coinc .* sqrt(singles_A+singles_B) * 1.6e-6/60;

% Error due to the precision in the angular position of the polarizers
a = [cot(pi/8)
    cot(3*pi/8)
    cot(3*pi/8)
    cot(pi/8)];
b = [cot(3*pi/8)
    cot(pi/8)
    cot(pi/8)
    cot(3*pi/8)];
ab16 = cat(2, a,b, a, a);
ab312 = ab16;
for i = 2:312
    ab312 = cat(2, ab312, ab16);
end
pol_err = coinc * sqrt(2) .* ab312 * .1/180*pi; 

% Combination of the systematic and statistical errors in the counts
coinc_err = sqrt(poissonian_error.^2 + counter_err.^2 + dead_err.^2 + pol_err.^2);

%% Error in the final value of S
coinc_total_err = [sqrt(sum(coinc_err.^2, 1));...
                           sqrt(sum(coinc_err.^2, 1));...
                           sqrt(sum(coinc_err.^2, 1));...
                           sqrt(sum(coinc_err.^2, 1))];

coinc_normalised_err = coinc_normalised...
                            .*sqrt((coinc_err./coinc).^2 ...
                            + (coinc_total_err./coinc_total).^2);
                        

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

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

fprintf('Experimental value of S: %f +/- %f\n',S_value,mean_S_std)
fprintf('Distance from T bound: %f +/- %f\n', Delta_T, mean_S_std)
%% do not run

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 = 0
% 
%     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