function stats = lpt_error_stats(ref, meas, opts)
	% calculate statistics of errors in particle tracking
	%
	% Input
	%	ref		structure containing reference data
	%		x	3xN		position
	%		u	3xN		velocity
	%	meas	structure containing measured data
	%		x	3xN		position
	%		u	3xN		velocity
	%	opts	options structure with parameters for statistics
	%		R_ghost_tol distance for identifying ghosts
	
	n_ref			= size(ref.x , 2);
	n_meas			= size(meas.x, 2);
	
	%% count ghosts
	% ghost particles are defined as particles in reconstruction which do
	% not have a nearest neighbour within opts.R_ghost_tol
	
	% find nearest neighbour
	[i_nn,d_nn]		= knnsearch(ref.x.', meas.x.', 'K', 1);
	b_ghost			= d_nn > opts.R_ghost_tol;
	b_real			= ~b_ghost;
	n_ghost			= nnz(b_ghost);
	n_real			= nnz(b_real);
	
	% define a subset of the reference
	i_nn_real		= i_nn(b_real);
	x_ref			= ref.x(:, i_nn_real);
	u_ref			= ref.u(:, i_nn_real);
	x_real			= meas.x(:, b_real);
	u_real			= meas.u(:, b_real);
	
	%% measure RMS position and velocity error
	u_res			= u_real - u_ref;
	x_res			= x_real - x_ref;
	x_res_rms		= rms(x_res, 2);
	u_res_rms		= rms(u_res, 2);
	x_res_med		= median(abs(x_res),2);
	u_res_med		= median(abs(u_res),2);
	
	%% save statistics in structure
	stats	= struct('n_ghost', n_ghost, ...
				'n_ref', n_ref, 'n_meas', n_meas, 'n_real', n_real, ...
				'b_ghost', b_ghost, ...
				'x_res_rms', x_res_rms, 'x_res_med', x_res_med, ...
				'u_res_rms', u_res_rms, 'u_res_med', u_res_med);
end