function [pTime] = pickerRatioSING(x, S, L, rCut, stage)
% USAGE:
%       [pTime] = pickerRatio(D, sWindow, lWindow, rCutoff, nInterval, test)
%
% pickerRatio.m: function to use an sta/lta ratio on a hilbert transformed
%                waveform to determine picks
%
% INPUTS:
%   D: coral structure for picking (one channel at a time)
%   sWindow: short term window for picker in seconds (0.8 s is suggested)
%   lWindow: long term window for picker in seconds (7 s is suggested)
%   rCutoff: ratio cutoff for declaring event. Positive numbers are cutoffs
%         with respect to the actual ratio.  Negative numbers are cutoffs
%         with respect to standard deviations from the mean of the noise.
%         (2.5 or -2 is suggested)
%         ex.
%           2 means that the ratio lWindow/sWindow is two
%          -2 fits a normal distribution to the ratios of the time series,
%             (assumes that noise is normally distributed) and picks the
%             cutoff to be where the ratios are 2 standard deviations above
%             the mean
%   nInterval: Number of seconds between events (10 s is suggested)
%   test: Boolean, plot out the triggers? 1 = yes, 0 = no
%
% OUTPUT:
%   pTime: Cell array of pick times in coral format
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%test = 1;
%XX = x;  % Your crappy ass data
%S = 50;     % Short Window in samples
%L = 100;    % Long Window in samples
%rCut = 2.5;  % ratio cutoff



%x = abs(real(hilbert(D.data)));
%S = sWindow/D.recSampInt;
%L = lWindow/D.recSampInt;
% Calculate ratio (from Creager)
sL = zeros(size(x));  % initialize output vector
 if length(x) > S+L;
   X       = cumsum(x);       % cumulative sum 
   k       = L+2:length(X)-S+1; % index vector
   sL(k)   = abs([X(k+S-1)-X(k-1)] ./ [X(k-1)-X(k-L-1)] * L/S);
   sL(L+1) = (X(L+S)-X(L))/X(L) *L/S; % special case for first point in time series
 end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate stats on ratio
%[muhat, sigmahat] = normfit(sL);
% Find ratio
%if rCutoff < 1
%    rCut = muhat + abs(rCutoff)*sigmahat;
%else
%    rCut = rCutoff;
%end

% Find events above ratio
r = find(sL >= rCut);

% Get triggers
eTrace = zeros(length(sL),1);

eTrace(r) = 1;
for i = 1 : length(r)
        % Cut out appropriate window
 %       if (r(i)+nInterval*1/D.recSampInt < length(eTrace))
 %           tSamp = eTrace(r(i):r(i)+nInterval*1/D.recSampInt);
 %           tSampV = sL(r(i):r(i)+nInterval*1/D.recSampInt);
 %       else
            tSamp = eTrace(r(i):end);
            tSampV = sL(r(i):end);
 %       end
        k = find(tSamp == 1);
        % if multiple triggers
 %       if (length(k) > 1)
 %           kk = find(tSampV~=max(tSampV));
 %           tSamp(kk) = 0;
 %           ind = r(i)+kk-1;
 %           eTrace(ind) = 0;
 %       end
 %       if (r(i)+nInterval*1/D.recSampInt < length(eTrace))
 %           eTrace(r(i):r(i)+nInterval*1/D.recSampInt) = tSamp;
 %       else
 %           eTrace(r(i):end) = tSamp;
 %       end
end

%disp(sprintf('Number of boundaries found: %d', length(find(eTrace ==1))));
% Get absolute picktimes
rr = find(eTrace == 1);
%for i = 1 : length(rr)
%    pTime{i} = timeadd(D.recStartTime, rr(i)*D.recSampInt);
    pTime = rr;
%end
%figure;plot( eTrace );
if length(rr) < 1
%    disp('No amygdales found on a cross-section');
%     pTime = {}; %Wes original
    pTime = 0; %Sanjoy modified
    return;
end

% Testing plot
if stage == 1
    figure; subplot(2,1,1);
%    coralPlot(D);
    %plot(x); hold on;
    xlabel('pixel across');
    ylabel('ratio (STA/LTA)');
    %title('Original data (demeaned)')
    
%    tvec = [0:length(eTrace)-1]*D.recSampInt;
    tvec = [0:length(eTrace)-1];
%    plot((x/max(x) * 0.5), '-b');
%    for i = 1 : length(rr)
%        plot(tvec(rr(i)), 0, 'or');
%    end
%    xlim([0 tvec(end)]);
%    ylabel('Counts');
    %ylim([-0.5 0.5]);
    %legend('Waveform', 'Hilbert', 'Trigger', 'Location', 'NorthEast');
    
    subplot(2,1,2);
    plot(tvec, sL, '-k');
    hold on;
    for i = 1 : length(rr)
        plot(tvec(rr(i)), rCut, 'or');
    end
    xlim([0 tvec(end)]);
    xlabel('pixel across');
    ylabel('ratio (STA/LTA)');
    %title('Original data (demeaned) with detected boundaries')
    axis tight
%    legend('STA/LTA', 'Trigger', 'Location', 'NorthEast');
%     subplot(3,1,3);
%     title('Histogram of sta/lta values');
%     h = hist(sL, [0:0.1:max(sL)]);
%     hist(sL, [0:0.1:max(sL)]);
%     hold on;
% %     nn = (1/(sigmahat*sqrt(2*pi)))*exp((-0.5*((muhat+sigmahat*randn(length(sL),1))-muhat).^2)./(muhat.^2));
%     plot([rCut rCut], [0 max(h)], '-r', 'LineWidth', 3);
%     xlabel('STA/LTA Ratio');
%     ylabel('Occurences');
%     text(rCut+0.2, mean([0 max(h)]), 'cutoff', 'HorizontalAlignment', 'center', 'Rotation', 90, 'Color', 'r');
%     %text(rCut+0.2, mean([0 max(h)]));
%     axis tight;
end
    pTime=pTime-1; % time = pixel - 1; when pixel = 1, time = 0. remember 
                  % that this is from a seismic picker...

