function rho = mf_snr(filename, Sh, scalefac, mkplot, T)
% function rho = mf_snr(filename, Sh, scalefac, mkplot, T)

% make a diagnostic plot?
try
  mkplot;
catch
  mkplot=true;
end

% scale injection?
try
  scalefac;
catch
  scalefac = 1;
end

% number of trials
try
  T;
catch
  T=100;
end

% load injection file
h = load(filename);

% parse injection file
t = h(:,1);
hp = scalefac*h(:,2);
hx = scalefac*h(:,3);
sampleRate = 1 / (t(2)-t(1));
duration = t(end);

% EHT: look at just hp for now
% EHT: antenna factors
s = hp;

% create signal
s = repmat(s, 1, T);
[sf f] = fft_eht(s, sampleRate);

% apply cut
cut = f>10 & f<1800;
sf = sf(cut, :);
f = f(cut);

% create noise
NN = length(t);
freqs = Sh(:,1);
P = Sh(:,2);
P2 = interp1(freqs, P, f);
Sh2 = [f P2];
n = gaussian_noise(Sh2, sampleRate, duration, T, NN);
[nf fn] = fft_eht(n, sampleRate);
nf = nf(cut, :);

% create MC data
h = s + n;
[hf fh] = fft_eht(h, sampleRate);
hf = hf(cut, :);

% calculate matched filtering SNR
% see Eq. 2.3 in Owen and Sathyaprakash paper
deltaF = f(2)-f(1);
S = real(inner_product(hf, sf, P2));
N = real(inner_product(nf, sf, P2));

% calculate SNR
rho = mean(S)/std(N);
rhoN = mean(N)/std(N);
rhoNmax = max(N)/std(N);

% print rho to screen
fprintf('rho = %1.3e\n', rho);
if rho>=4.5
  fprintf('d = %1.2e * d_ref for rho=5\n', sqrt(rho/5) );
end
%fprintf('rho_noise_avg = %1.3e\n', rhoN);
%fprintf('rho_noise_max = %1.3e\n', rhoNmax);

if mkplot
  % make a plot of signal, noise, and theoretical noise spectra
  figure;
  loglog(f, abs(nf(:,1)).^2, 'r');
  hold on;
  loglog(f, P2(:,1), 'k');
  hold on;
  loglog(f, abs(sf(:,1)).^2, 'b');
  xlabel('f (Hz)')
  ylabel('strain^2');
  pretty;
  print('-dpng', 'mf_snr');
  % histogram of SNR
  x = -10:0.5:10;
  nS = hist(S/std(N), x);
  nN = hist(N/std(N), x);
  figure;
  plot(x, nS, 'b');
  hold on;
  plot(x, nN, 'r');
  xlabel('\rho');
  pretty;
  print('-dpng', 'mf_histo');
end

return
