function spikeTrigAvg_postsyn_curveset(files, varargin)
%% Calculate spike triggered average of the postsynaptic response
% For each spike in the stimulated cell the membrane potential of the not
% stimulated cell is gathered and summed up from 5 ms before to 30 ms after
% the spike. Devided by the number of spikes the averaged postsynaptic
% membrane potential is obtained. If there is a clear peak (due to EPSPs 
% for example), the time of this peak is determined as the latency of the 
% reaction to the spike in the other cell. The latency and the number of 
% used spikes are returned. Alternatively the result can be plotted. The 
% function can handle multiple files at a time.
% -------------------------------------------------------------------------
% Input     [additional arguments are positional, 
% -----      optional arguments are name-value pairs]
%
%   file:           number(s) of the '.mat' file(s) in the directory,
%                       type: integer or vector
%   trials:         number(s) of the trials to be analysed, type: integer
%   (additional)        of vector, default: 'all'
%   Filter:         select between three filter options: 
%   (optional)          'hn'    applys a hum-noise-remove filter
%                       'hnhfn' additonally applys a low pass filter to
%                              remove high frequency noise
%                       If not set, no filter will be applyed (default).
%   MinPeakPromT:   set the value for 'MinPeakProminence' which is used in
%   (optional)          the 'findpeaks' function to find T-cell peaks. 
%                       Default: 17
%   StimCell:       set whether the data for a stimulated T cell or a 
%   (optional)          stimulated interneuron should be used.
%                       'INT' for the interneuron
%                       'T'   for the T-cell (default)
% ------
% Output
% ------
%   varargout{1}:   latency from spike of the stimulated cell to peak in 
%                   average response of the not stimullated cell (no plot)
%                       type: float (single file) or vector (multiple
%                             files)
%   varargout{2}:   number of spikes used from the stimulated cell
%                       type: float (single file) or vector (multiple
%                             files)
%
%   visually output of the plot, if no output variable set
% -------------------------------------------------------------------------
% Program by: Bjarne Schultze [last modified 09.03.2022]
% -------------------------------------------------------------------------

%% Preparations
% parse the user input
[file_nums, trials, filter_choice, ~, MinPeakPromT, ...
    stimCell, ~,~,~] = parse_input_peakfun(files, varargin{:});
% extract the requested data
recData = extract_data(file_nums, filter_choice);

% load a file with color values (shades)
load("../AdditionalFiles/colorShades.mat", "colorShades")

% assign variables to check the min and max values of the y-data
minY = 100; maxY = -100;

%% Main calculations
% iterate through the single files as requested by the input
for file = 1:length(file_nums)
    % get the number or trials for the current file
    trial_count = size(recData{file,1}.recordingT,2);
    % set the length of the time frame to gather the triggered data
    trig_len = 0.03;                     % in seconds
    trig_len_dp = trig_len * 10000;      % in datapoints
    
    % distinguish between a stimulated T cell and a stimulated interneuron
    if isequal(stimCell, "INT")
        % get the columns with data for a stimulated interneuron
        stimulated = setdiff(1:trial_count, recData{file,1}.stimulatedT(:), ...
            'stable');
        % select the requested trials, ensure only valid trials
        if ~isequal(trials, 'all')
            stimulated = intersect(stimulated,trials);
        end
        % set an appropriate plot title
        plot_title = sprintf("Postsynaptic Spike Response - File: " + ...
            "%02.f - Interneuron stimulated", file_nums(file));
        % assign the data to working variables
        recDataNSTIM = recData{file,1}.recordingT(:,stimulated);
        recDataSTIM = recData{file,1}.recordingINT(:,stimulated);
    else
        % get the columns with data for a stimulated T-cell
        stimulated = recData{file,1}.stimulatedT(:);
        % select the requested trials, ensure only valid trials
        if ~isequal(trials, 'all')
            stimulated = intersect(stimulated,trials);
        end
        % set an appropriate plot title
        plot_title = sprintf("Postsynaptic Spike Response - File: " + ...
            "%02.f - T-cell stimulated", file_nums(file));
        % assign the data to working variables
        recDataSTIM = recData{file,1}.recordingT(:,stimulated);
        recDataNSTIM = recData{file,1}.recordingINT(:,stimulated);
        MinPeakPromSTIM = MinPeakPromT;
    end
    
    %% Main calculations
    % pre-define a variable to sum up the triggered membrane potentials
    triggered_sum = 0;                                         
    spike_counter = 0;
    % set up a figure window for the curvesets of each file
    figure();
    % iterate through the trials
    for exp = 1:length(stimulated)
        % assign a variable for the curve set data
        triggered_data_trial = 0; spike_counter_trial = 0;
        % search for spikes in the stimulated cell
        if isequal(stimCell, 'INT')
            [~,spike_indexSTIM_all] = findIntSpikes(recDataSTIM(:,exp));
        else
            [~, spike_indexSTIM_all] = findpeaks(recDataSTIM(:,exp), ...
                'MinPeakProminence', MinPeakPromSTIM);
        end
        % exclude spikes at the very beginning and the end of one recording
        % to ensure a complete time frame for the triggered data
        spike_indexSTIM_sel = spike_indexSTIM_all > 50 & ...
            spike_indexSTIM_all < (size(recDataSTIM,1) - trig_len_dp);
        spike_indexSTIM = spike_indexSTIM_all(spike_indexSTIM_sel);
        
        % if no spikes are found the next trial is examined
        if isempty(spike_indexSTIM)
            continue
        else
            % select the membrane potentail of the not-stimulated cell in 
            % an index frame of the length 'trig_len_dp' datapoints for 
            % each spike of the stimulated cell
            for spike = 1:length(spike_indexSTIM)
                % select the data (triggered by a spike)
                triggered_data = ...
                    recDataNSTIM((spike_indexSTIM(spike)-50):...
                    (spike_indexSTIM(spike) + trig_len_dp), exp);
                % sum up the data
                triggered_sum = triggered_sum + triggered_data;

                % sum up the data for each trial, then restart
                triggered_data_trial = triggered_data_trial + ...
                    triggered_data;
            end
            % count the spikes
            spike_counter = spike_counter + spike;

            % count the spikes for one trial, then restart
            spike_counter_trial = spike_counter_trial + spike;
        end
        % create the curve set
        afterSpkAvg_trial = triggered_data_trial ./ spike_counter_trial;
        % create an appropriate time vector
        timeVector = [sort(-recData{file,1}.timeVector(1:50)); 0;...
                recData{file,1}.timeVector(1:trig_len_dp)];
        timeVector = timeVector * 1000;     % in ms
        % set the color shades, getting brighter for later trials
        if exp > length(colorShades.purple(:,1))
            plot_color = colorShades.purple(end,:);
        else
            plot_color = colorShades.blue(exp,:);
        end
        % plot the data
        hold on;
        plot(timeVector,afterSpkAvg_trial,"Color", plot_color, ...
            'LineWidth', 1)
        hold off; 
        % check for the min and max of the y-data
        if min(afterSpkAvg_trial) < minY
            minY = min(afterSpkAvg_trial);
        end
        if max(afterSpkAvg_trial) > maxY
            maxY = max(afterSpkAvg_trial);
        end
        ylim([minY,maxY])
        if isequal(exp,1)
            xline(0, 'Color', [0.6 0.6 0.6])
        end
    end
    % calculate the averaged postsynaptic membrane potential if more than 8
    % spikes were found in the T cell
    if spike_counter < 8
        fprintf("Not enough T cell spikes (n < 8)!\n")
        title(sprintf("File: %g - Not engough T cell spikes (n < 8)!",...
            file_nums(file)))
        xlabel("time after spike [ms]"); 
        ylabel("average membrane potentail [mV]");
    else
        afterSpkAvg = triggered_sum ./ spike_counter;
        hold on;
        plot(timeVector, afterSpkAvg, 'Color', colorShades.green(1,:), ...
            'LineWidth', 2);
        hold off;
        title(plot_title);
        xlabel("time after spike [ms]"); 
        ylabel("average membrane potentail [mV]");
        
    end
end

% end of function definition
end