% Analysis of Rowan's temporal multiplexing data with 2x2 loop

clear all

load('LoopNoLoopCAR.mat');										% data file

numpower = 10;													% number of pump power levels
indexp = (1:1:10);

ts = 120;														% total time counted for (s)

reprate = 10e6;													% pump laser rep rate
aomratio = 2;													% divide down ratio of AOM
mpdepth = 4;													% Multiplexing depth
binrate = reprate/(aomratio*mpdepth);							% rep rate of final bins into which photons shifted

rawcountn = {P0N P1N P2N P3N P4N P5N P6N P7N P8N P9N};			% raw counts, each row per second, without loop
rawcountl = {P0L P1L P2L P3L P4L P5L P6L P7L P8L P9L};			% raw counts, each row per second, with loop

scountn = zeros(numpower,24);
scountl = zeros(numpower,24);
mcountn = zeros(numpower,24);
mcountl = zeros(numpower,24);

for n = 1:numpower
	scountn(n,:) = sum(cell2mat(rawcountn(n)))';				% total counts in 120s without loop
	scountl(n,:) = sum(cell2mat(rawcountl(n)))';				% total counts in 120s with loop
	mcountn(n,:) = mean(cell2mat(rawcountn(n)))';				% mean count rates per second without loop
	mcountl(n,:) = mean(cell2mat(rawcountl(n)))';				% mean count rates per second with loop
end

% col 1: raw herald from silicon
% col 2: "computed herald"
% col 3: idler singles from InGaAs
% col 4: coincidences between computed herald and InGaAs
scountn = [scountn(:,1) scountn(:,2) scountn(:,6) scountn(:,12)];	
scountl = [scountl(:,1) scountl(:,2) scountl(:,6) scountl(:,12)];
mcountn = [mcountn(:,1) mcountn(:,2) mcountn(:,6) mcountn(:,12)];
mcountl = [mcountl(:,1) mcountl(:,2) mcountl(:,6) mcountl(:,12)];

clear -regexp ^P ^C ^raw;										% clear unwanted variables

% detection probability in target bin
probn = mcountn(:,4)/binrate;
proberrn = sqrt(scountn(:,4))/(binrate*ts);					% error

probl = mcountl(:,4)/binrate;
proberrl = sqrt(scountl(:,4))/(binrate*ts);					% error

% maximum improvement
floop = probl(numpower,1)/probn(numpower,1);
flooperr = floop*sqrt((proberrn(numpower,1)/probn(numpower,1))^2 +...
						(proberrl(numpower,1)/probl(numpower,1))^2);

% heralding efficiency
heffn = mcountn(:,4)./mcountn(:,1);
hefferrn = heffn.*sqrt(((sqrt(scountn(:,4))/ts)./mcountn(:,4)).^2 +...
				((sqrt(scountn(:,1))/ts)./mcountn(:,1)).^2);

heffl = mcountl(:,4)./mcountl(:,1);
hefferrl = heffl.*sqrt(((sqrt(scountl(:,4))/ts)./mcountl(:,4)).^2 +...
				((sqrt(scountl(:,1))/ts)./mcountl(:,1)).^2);

% Herald-Idler cross-correlation (CAR?) using computed herald counts
xcorrn = binrate*mcountn(:,4)./(mcountn(:,2).*mcountn(:,3));

xcorrl = binrate*mcountl(:,4)./(mcountl(:,2).*mcountl(:,3));

%%

fs = 20;
lw = 2;
ms = 20;

% plot count rates against power index
% figure(1)
% clf
% subplot(2,2,1); plot(indexp,mcountn(:,1),indexp,mcountl(:,1))
% subplot(2,2,2); plot(indexp,mcountn(:,2),indexp,mcountl(:,2))
% subplot(2,2,3); plot(indexp,mcountn(:,3),indexp,mcountl(:,3))
% subplot(2,2,4); plot(indexp,mcountn(:,4),indexp,mcountl(:,4))

% plot count rates against herald rate
figure(2)
clf
subplot(2,2,1); plot(mcountn(:,1),mcountn(:,1),mcountl(:,1),mcountl(:,1))
subplot(2,2,2); plot(mcountn(:,1),mcountn(:,2),mcountl(:,1),mcountl(:,2))
subplot(2,2,3); plot(mcountn(:,1),mcountn(:,3),mcountl(:,1),mcountl(:,3))
subplot(2,2,4); plot(mcountn(:,1),mcountn(:,4),mcountl(:,1),mcountl(:,4))

figure(3)
clf
set(gcf,'position',[100 200 1000 350])
% Single photon delivery probability per bin
subplot(1,2,1)
hold on
errorbar(mcountn(:,1),probn,proberrn,...
	'k','linestyle','-','marker','.','linewidth',lw,'markersize',ms)
errorbar(mcountl(:,1),probl,proberrl,...
	'b','linestyle','-','marker','.','linewidth',lw,'markersize',ms)
set(gca, 'FontName', 'Helvetica');
set(gca, 'FontSize', fs);
xlabel('Herald count rate (s^{-1})')
ylabel('Probability per bin')
% Heralding efficiency
subplot(1,2,2)
hold on
errorbar(mcountn(:,1),heffn,hefferrn,...
	'k','linestyle','-','marker','.','linewidth',lw,'markersize',ms)
errorbar(mcountl(:,1),heffl,hefferrl,...
	'b','linestyle','-','marker','.','linewidth',lw,'markersize',ms)
set(gca, 'FontName', 'Helvetica');
set(gca, 'FontSize', fs);
xlabel('Herald count rate (s^{-1})')
ylabel('Heralding efficiency')

% CAR plot
figure(4)
clf
set(gcf,'position',[100 200 600 400])

hold on
plot(mcountn(:,4),xcorrn,...
	'k','linestyle','-','marker','.','linewidth',lw,'markersize',ms)
plot(mcountl(:,4),xcorrl,...
	'b','linestyle','-','marker','.','linewidth',lw,'markersize',ms)
set(gca, 'FontName', 'Helvetica');
set(gca, 'FontSize', fs);
xlabel('Mean coincidence count rate (s^{-1})')
ylabel('g^{h,i}(0)')


%% Improvement factors with loss backed out of simplex data

p1n = (1/10^-0.1)*mcountn(numpower,4)/binrate;
p1l = mcountl(numpower,4)/binrate;

flooptrue = p1l/p1n;

N = (1:1:20);

countNn = reprate*p1n.^N;
countNl = binrate*p1l.^N;

fN = countNl./countNn;

figure(5)
clf
set(gcf,'position',[100 200 400 300])

hold on
plot(N,fN,...
	'k','linestyle','-','marker','.','linewidth',lw,'markersize',ms)
%set(gca,'xlim',[min(N) max(N)]);
set(gca, 'FontName', 'Helvetica');
set(gca, 'FontSize', fs);
xlabel('N')
ylabel('f^{(m)}_N')
line([min(N) max(N)],[1 1],'color','r','linewidth',lw);

%% Master plot


figure(6)
clf
set(gcf,'position',[100 200 1000 700])
% Single photon delivery probability per bin
subplot(2,2,1)
hold on
errorbar(mcountn(:,1),probn,proberrn,...
	'k','linestyle','-','marker','.','linewidth',lw,'markersize',ms)
errorbar(mcountl(:,1),probl,proberrl,...
	'b','linestyle','-','marker','.','linewidth',lw,'markersize',ms)
set(gca, 'FontName', 'Helvetica');
set(gca, 'FontSize', fs);
xlabel('Herald count rate (s^{-1})')
ylabel('Probability per bin')
text(-0.55e5,4.3e-4,'(a)','FontName', 'Helvetica','FontSize', fs)

% Heralding efficiency
subplot(2,2,2)
hold on
errorbar(mcountn(:,1),heffn,hefferrn,...
	'k','linestyle','-','marker','.','linewidth',lw,'markersize',ms)
errorbar(mcountl(:,1),heffl,hefferrl,...
	'b','linestyle','-','marker','.','linewidth',lw,'markersize',ms)
set(gca, 'FontName', 'Helvetica');
set(gca, 'FontSize', fs);
xlabel('Herald count rate (s^{-1})')
ylabel('Heralding efficiency')
text(-0.6e5,1.85e-3,'(b)','FontName', 'Helvetica','FontSize', fs)

% CAR plot
subplot(2,2,3)
hold on
plot(mcountn(:,4),xcorrn,...
	'k','linestyle','-','marker','.','linewidth',lw,'markersize',ms)
plot(mcountl(:,4),xcorrl,...
	'b','linestyle','-','marker','.','linewidth',lw,'markersize',ms)
set(gca, 'FontName', 'Helvetica');
set(gca, 'FontSize', fs);
xlabel('Mean coincidence count rate (s^{-1})')
ylabel('g^{h,i}(0)')
text(-70,55,'(c)','FontName', 'Helvetica','FontSize', fs)

% improvement factor
subplot(2,2,4)
hold on
plot(N,fN,...
	'b','linestyle','-','marker','.','linewidth',lw,'markersize',ms)
%set(gca,'xlim',[min(N) max(N)]);
set(gca, 'FontName', 'Helvetica');
set(gca, 'FontSize', fs);
xlabel('N')
ylabel('f^{(m)}_N')
line([min(N) max(N)],[1 1],'color','r','linewidth',lw);
text(-3.5,8.5,'(d)','FontName', 'Helvetica','FontSize', fs)


%% Master plot version 2

figure(7)
clf
set(gcf,'position',[100 200 1000 700])
% Single photon delivery probability per bin
subplot(2,2,1)
hold on
errorbar(mcountn(:,1),probn,proberrn,...
	'k','linestyle','--','marker','.','linewidth',lw,'markersize',ms)
errorbar(mcountl(:,1),probl,proberrl,...
	'b','linestyle','--','marker','.','linewidth',lw,'markersize',ms)
set(gca, 'FontName', 'Helvetica');
set(gca, 'FontSize', fs);
axis([0 2.2e5 0 3.2e-4]);
xlabel('Herald count rate (s^{-1})')
ylabel('Probability per bin')
text(-0.3e5,3.5e-4,'(a)','FontName', 'Helvetica','FontSize', fs)

%fit(mcountn(:,4),xcorrn,'a1./(a2.*x)');

% CAR plot
subplot(2,2,2)
hold on
plot(mcountn(:,4),xcorrn,...
	'k','linestyle','--','marker','.','linewidth',lw,'markersize',ms)
plot(mcountl(:,4),xcorrl,...
	'b','linestyle','--','marker','.','linewidth',lw,'markersize',ms)
set(gca, 'FontName', 'Helvetica');
set(gca, 'FontSize', fs);
xlabel('Mean coincidence count rate (s^{-1})')
ylabel('g^{h,i}(0)')
text(-70,55,'(b)','FontName', 'Helvetica','FontSize', fs)

% improvement factor
subplot(2,2,3.5)
hold on
plot(N,fN,...
	'b','linestyle','none','marker','.','linewidth',lw,'markersize',ms)
%set(gca,'xlim',[min(N) max(N)]);
set(gca, 'FontName', 'Helvetica');
set(gca, 'FontSize', fs);
xlabel('N')
ylabel('f^{(m)}_N')
line([min(N) max(N)],[1 1],'color','r','linewidth',lw);
text(-3.5,8.5,'(c)','FontName', 'Helvetica','FontSize', fs)