function [avg_psth, bins] = compute_count_time(threshold_amp, direction, f_sample, binsize, rec_length, file_struct) 
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function computes the average PSTH from a set of preprocessed datasets.
% 
% threshold_amp [mV] = the amplitude threshold that separates high
% amplitude from low amplitude events;
% direction = selects which events to count: 'low' / 'high' / 'all'
% f_sample [Hz] = sampling frequency of event timestamps
% binsize [s] = length of PSTH bin
% rec_length [s] = how long to compute the PSTH  
% files = the file structure of preprocessed datasets: 
% file_struct(file_number, feature_number) = [session_number recording_number
% baseline/response control/target dpf errors]
% 
% Note that:
% session_number = number of the session - should be present in the
% folder name;
% recording_number = number of the recording - should be present in the
% last digits of the recording name;
% baseline/response = 0 if the trace is recorded only during spontaneous
% activity; 1 if trace contains response to PTZ 
% control/target = 0 if the trace is recorded from a Mismatch control fish;
% 1 if trace is recorded from a DEPDC5 animal
%
% file_struct(1,:) = [11 2 1 0 4 0];
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% define bin edges for PSTH
bins = 0:binsize:rec_length;

% Initialize structures: 
psth(length(file_struct)) = struct;           % will hold the psth for each dataset
filenames(length(file_struct(:,1))) = struct; % will hold the filenames for each dataset
sum_psth = zeros(1,length(bins));       % will hold the sum of the PSTHs across all datasets

% find number of files
no_sets = size(file_struct, 1);

% iterate through datasets
for i_dataset=1:no_sets
    
    % find dataset in file structure
    session = ['Session ' num2str(file_struct(i_dataset,1))];
    list = dir(session);
    
    if length(list)>2
        filename_base = list(3,1).name;
        clear list;
    end
    if file_struct(i_dataset,2)<10
        filename = [filename_base(1:7) num2str(file_struct(i_dataset,2))];
    else
        filename = [filename_base(1:6) num2str(file_struct(i_dataset,2))];
    end
    filenames(i_dataset).name = filename;
    
    % load preprocessed peak matrix 
    load([session '\' filename '_peaks' ]);

    % select desired events based on amplitude, as defined by 'direction'
    % and 'threshold amp'
    if strcmp(direction,'low')
        amp_int_peaks = peak_matrix(peak_matrix(:,3)< threshold_amp,:);
    elseif strcmp(direction,'high')
        amp_int_peaks = peak_matrix(peak_matrix(:,3)> threshold_amp,:);
    else
        amp_int_peaks = peak_matrix;
    end
    
    % compute PSTH
    if (amp_int_peaks)
        times = (amp_int_peaks(:,1))/f_sample;
        psth(i_dataset).psth = histc(times', bins);
    else
        psth(i_dataset).psth = zeros(1,length(bins));
    end
    
    % add PSTH to the sum
    sum_psth = sum_psth+psth(i_dataset).psth;

    clear peak_matrix trig times amp_int_peaks;
end

% compute the average PSTH
avg_psth = sum_psth/length(file_struct(:,1));
%--------------------------------------------------------------------------
