function [PeakLocations]=FACSTC_PeakFind(filenames, col, Bins,Gates,Output)
%%
%%%
% SUMMARY:
% This function takes a filename of a .fcs file OR the data from the .fcs
% file (see below), lets the user select a gate from a plot of the
% side-scatter vs the forward scatter and returns the indices of the data
% points in the fcs data that are within this gate.
%
% INPUTS:
%   filenameORdata: This is either the name of a .fcs file to open, or the actual fcs FACS data 
%   (ie fcsdat from the call to  [fcsdat, fcshdr] = fca_readfcs(filename))
%   
%   col-column of data from the .fcs file that we are interested in (first
%   two columns for our machine are forward and side scatter, after that
%   the columns are the fluorescence intensity information for different
%   channels
%
%   Gates-indices of points to include (based on forward and side scatter)
%   returned by FindGate
%  

%   varargin:
%   varargin{1}: 0 if output is to be suppressed (default), 1 if output
%   (ie, plot of the gate and the stuff in the gate) IS to be plotted
%   varargin{2}: name for the title of the plot when data is sent in rather
%   than a filename
%   
%
% OUTPUTS:
%   PeakLocations: Locations of the largest two peaks in the histogram of the 
%   logarithm of the data (in column) for all files that filenameORdata
%   specifies.
%   
%    
% REQUIRES:
%   fca_readfcs
%   
%
% REFERENCES:
% Mathworks documentation on inpolygon:
%    http://www.mathworks.com/help/techdoc/ref/inpolygon.html
%
%    Written by Megan McClean, Ph.D.
%               Lewis-Sigler Institute for Integrative Genomics
%               Princeton University
%               120 Carl Icahn
%               Princeton, NJ 08544
%               mmcclean@princeton.edu
%
%    Last revised on February 22, 2012


%%
%%
if ~iscell(filenames)
    files=dir(filenames);
else
    files=struct('name',[]);    
    for i=1:length(filenames)
        files(i).name=filenames{i};
    end
end


PeakLocations=[];


for i=1:length(files)
  [fcsdat, fcshdr] = fca_readfcs(files(i).name);
  [n,x]=hist(log(fcsdat(Gates{i},col)),Bins);
  [P,L]=findpeaks(smooth(n,0.30,'loess'));
  Peaks=sortrows([P L]);
  if length(P)>1
    PeakLocations=[PeakLocations; x(Peaks(end,2)) x(Peaks(end-1,2))];
  else
    PeakLocations=[PeakLocations; x(Peaks(1,2)) x(Peaks(1,2))];
  end
  %Display purposes only
  if Output==1
      figure; hold;
      plot(x,n);
      plot(x,smooth(n,0.30,'loess'),'r');
      plot(x(L),P,'go')
      pause;
      close all;
  end
      



end