%% input parameters
src_dir_A		= 'L:\array\tcf\30C_E\tracks\';
src_pat_A		= '%s/z001_t001.mat';
src_dir_B		= 'L:\array\tcf\30C_F\tracks\';
src_pat_B		= '%s/z001_t001.mat';

ref_dir			= 'L:\array\tcf\rescaled\';
ref_pat			= '%s/z001_t001.mat';
font_size		= 24;

%% load reference
ref_file		= sprintf(ref_pat, ref_dir);
fprintf('reference data %s\n', ref_file);
ref				= load(ref_file);

%% load data
src_file_A		= sprintf(src_pat_A, src_dir_A);
fprintf('source data %s\n', src_file_A);
src_A			= load(src_file_A);
s_A				= struct('x',[src_A.tracks.x],'u',[src_A.tracks.u]);

src_file_B		= sprintf(src_pat_B, src_dir_B);
fprintf('source data %s\n', src_file_B);
src_B			= load(src_file_B);
s_B				= struct('x',[src_B.tracks.x],'u',[src_B.tracks.u]);

%% merge datasets to find particles in both A and B
re_tol			= src_A.opts.matching.R_tolerance_obj;
fprintf('tolerance for matching = %0.3f mm\n', re_tol);

[iAnearB, dist]= knnsearch(s_A.x', s_B.x.', 'K', 1');
b_close			= dist < re_tol;
iB				= find(b_close);
matches			= [iAnearB(b_close,1), iB(:)];
n_AB			= size(matches, 1);
n_A				= size(s_A.x,2);
n_B				= size(s_B.x,2);
n_ref			= size(ref.x,2);

fprintf('A   has %d particles\n', n_A);
fprintf('B   has %d particles\n', n_B);
fprintf('A&B has %d particles in common\n', n_AB);
% compare 
x_pair_A		= s_A.x(:,matches(:,1));
u_pair_A		= s_A.u(:,matches(:,1));
x_pair_B		= s_B.x(:,matches(:,2));
u_pair_B		= s_B.u(:,matches(:,2));
s_AB			= struct();
s_AB.x			= 1/2 * (x_pair_A + x_pair_B);
s_AB.u			= 1/2 * (u_pair_A + u_pair_B);

%% evaluate statistics
stat_opts		= struct('R_ghost_tol', re_tol);
% error statistics for conventional ptv
statsA			= lpt_error_stats(ref, s_A, stat_opts);
statsB			= lpt_error_stats(ref, s_B, stat_opts);
statsAB			= lpt_error_stats(ref, s_AB, stat_opts);

%% report statistics
fprintf('particles: %6d / %6d / %6d\n', statsA.n_meas, statsB.n_meas, statsAB.n_meas);
fprintf('ghosts:    %6d / %6d / %6d\n', statsA.n_ghost, statsB.n_ghost, statsAB.n_ghost);
fprintf('reals:     %6d / %6d / %6d\n', statsA.n_real, statsB.n_real, statsAB.n_real);
fprintf('A pos err: %6.3f %6.3f %6.3f mm rms\n', statsA.x_res_rms);
fprintf('           %6.3f %6.3f %6.3f mm med\n', statsA.x_res_med);
fprintf('B pos err: %6.3f %6.3f %6.3f mm rms\n', statsB.x_res_rms);
fprintf('           %6.3f %6.3f %6.3f mm med\n', statsB.x_res_med);
fprintf('AB pos err:%6.3f %6.3f %6.3f mm rms\n', statsAB.x_res_rms);
fprintf('           %6.3f %6.3f %6.3f mm med\n', statsAB.x_res_med);

fprintf('A vel err: %6.2f %6.2f %6.2f mm/s rms\n', statsA.u_res_rms);
fprintf('           %6.2f %6.2f %6.2f mm/s med\n', statsA.u_res_med);
fprintf('B vel err: %6.2f %6.2f %6.2f mm/s rms\n', statsB.u_res_rms);
fprintf('           %6.2f %6.2f %6.2f mm/s med\n', statsB.u_res_med);
fprintf('AB vel err:%6.2f %6.2f %6.2f mm/s rms\n', statsAB.u_res_rms);
fprintf('           %6.2f %6.2f %6.2f mm/s med\n', statsAB.u_res_med);

%% bar chart comparison
figure(1);
layout=tiledlayout(1,3);

% number of particles reconstructed
nexttile;
cats=categorical({'Real','Ghosts'});
bar(categorical({'Real','Ghosts'}), ...
	[[statsA.n_real; statsA.n_ghost], [statsB.n_real; statsB.n_ghost], [statsAB.n_real; statsAB.n_ghost]])
hold on;
plot(cats, n_ref*[1 1 ], 'k--');
ylabel('\# particles', 'interpreter', 'latex');
set(gca,'fontsize',font_size,'ticklabelinterpreter','latex');
pbaspect([1 1 1]);
set(gca, 'yscale', 'log');
ylim([100 1e6]);

% velocity error
nexttile;
scale_px_to_mm=0.1;
bar(categorical({'u','v','w'}), [statsA.u_res_med, statsB.u_res_med, statsAB.u_res_med]*delta_t/scale_px_to_mm)
ylabel('Median velocity error [px]', 'interpreter', 'latex');
set(gca,'fontsize',font_size,'ticklabelinterpreter','latex');
pbaspect([1 1 1]);
legend('A','B','A\&B');
set(legend,'location','west','interpreter','latex');

% position error
nexttile;
bar(categorical({'x','y','z'}), [statsA.x_res_med, statsB.x_res_med, statsAB.x_res_med]/scale_px_to_mm)
ylabel('Median position error [px]', 'interpreter', 'latex');
set(gca,'fontsize',font_size,'ticklabelinterpreter','latex');
pbaspect([1 1 1]);
