%% Data import

% set number display format
format long g;
% define data processing constants
SET_SIZE = 16;

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, :) = [];

%% reshaping data
singles_A = reshape(data(:, 2), 16, NO_OF_SET)';
singles_B = reshape(data(:, 3), 16, NO_OF_SET)';

coinc = reshape(data(:, 4), 16, NO_OF_SET)';

%% error in coincidences counts

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

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

% Error due to the on-board clock
clock_err = 0;%coinc * .1*1e-6;

% Error due to the dead time of the detectors
dead_err = 0;%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(1, a, b, a, a)';
ab312 = ab16;
for i = 2:312
    ab312 = cat(1, ab312, ab16);
end
pol_err = 0;%coinc * sqrt(2) .* ab312 * .1/180*pi;

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

%% combination of the 312 runs
ens = sum(coinc, 1);
ens_err = sqrt(sum(coinc_err.^2,1));

%% calculating E
E = zeros(4,1);
E_err = zeros(4,1);
for i=1:4
    N = ens(4*i-3)+ens(4*i)+ens(4*i-1)+ens(4*i-2);
    N2 = N^2;
    E(i) = (ens(4*i-3)+ens(i*4)-ens(4*i-2)-ens(4*i-1))/N;
    E_err(i) = sqrt(...
           (2*(ens(4*i-1)+ens(4*i-2))/N2)^2 * ens_err(4*i-3)^2+...
           (2*(ens(4*i-3)+ens(4*i))/N2)^2 * ens_err(4*i-2)^2+...
           (2*(ens(4*i-3)+ens(4*i))/N2)^2 * ens_err(4*i-1)^2+...
           (2*(ens(4*i-1)+ens(4*i-2))/N2)^2 * ens_err(4*i)^2 ...
                );
end

%% calculating S

S = sum(abs(E));
S_err = sqrt(sum(E_err.^2));
Delta_T = 2*sqrt(2)- S;
S_G = 2.825372;
Delta_G = abs((S_G-S)/S_err);

fprintf('Experimental value of S: %f +/- %f\n',S,S_err)
fprintf('Distance from T bound: %f +/- %f\n', Delta_T, S_err)
fprintf('Distance from G bound: %f sigmas\n', Delta_G)
