% use this function to detect the PSPs in your data using a PSP finding algorithm
% params = detectPSPs(data, PSPsDown, windows, windowSize, parameter pairs)
% or params = detectPSPs(data, PSPsDown, windows, windowSize)
% or params = detectPSPs(data, PSPsDown)
% if either pairs of parameter names and values xor a parameters structure is passed they must be the last arguement

% The algorithm used is based on that of Cohen and Miles 2000
% Briefly, the data are filtered and the derivative is taken.
% A function is generated based on the rule that a point is equal to the
% sum of the derivative at that point and the previous point only if the
% derivative is greater than (for upward PSPs) zero otherwise the function
% is zero at that point.
% This is then thresholded and sent for fitting by alpha functions and one
% variable is returned:
% Params(1,:) = amplitude of psp
% Params(2,:) = rise time
% Params(3,:) = x-offset of alpha function
% Params(4,:) = decay time

% To untangle multiple stacked PSPs as after a stim the readjustment of
% minYOffset/maxYOffset and also an iterative mechanism to break down giant

function [params decayFits] = detectPSPs(data, varargin)

warning off MATLAB:divideByZero
warning off all % to avoid a lot of errors with the best fitting

%%%
% Parse through the inputs
%%%

% pick the orientation that we like
if nargin > 0
if size(data, 1) > size(data, 2)
data = data';
end
else
error('Must specify data on which to operate');
end

% if more than one data set then reject
if size(data, 1) > 1
error('Input must be a 1 x N vector')
end

% look for the structure of input parameters
hasParameters = 0;
for i = 1:nargin - 1
if isstruct(varargin{i})
hasParameters = i;
parameters = varargin{i};
break
end
end

% find the first pair of input properties
firstChar = nargin + 1;
if ~hasParameters
for i = 1:nargin - 1
if ischar(varargin{i})
firstChar = i + 1;
break
end
end
end

switch firstChar - 1 - (hasParameters > 0)
case 0
windows = 1; % so that we include the whole time
windowSize = length(data);
% determine whether events are up or down
if median(data) > mean(data)
% events are down
PSPsDown = 1;
else
% events are up
PSPsDown = 0;
end
case 1
windows = 1; % so that we include the whole time
windowSize = length(data);
% determine whether events are up or down
if median(data) > mean(data)
% events are down
PSPsDown = 1;
else
% events are up
PSPsDown = 0;
end
case 2
PSPsDown = varargin{1};
windows = 1; % so that we include the whole time
windowSize = length(data);
case 3
% determine whether events are up or down
if median(data) > mean(data)
% events are down
PSPsDown = 1;
else
% events are up
PSPsDown = 0;
end
windows = varargin{1};
windowSize = varargin{2};
case 4
if ismember(varargin{1}, [0 1])
PSPsDown = varargin{1};
windows = varargin{2};
windowSize = varargin{3};
else
% determine whether events are up or down
if median(data) > mean(data)
% events are down
PSPsDown = 1;
else
% events are up
PSPsDown = 0;
end
windows = varargin{1};
windowSize = varargin{2};
parameters = varargin{3};
end
case 5
PSPsDown = varargin{1};
windows = varargin{2};
windowSize = varargin{3};
parameters = varargin{4};
end

if ~hasParameters
if PSPsDown
parameters = struct('minAmp', (-5000),... % minimum allowable amplitude for alpha functions (in units of samples)
'maxAmp', (-5),... % maximum allowable amplitude
'minTau', (5),... % minimum allowable tau for alpha functions (in units of samples)
'maxTau', (100),... %maximum allowable tau
'minYOffset', (-100),... % minimum allowable yOffset for alpha functions (in units of mV)
'maxYOffset', (-30),... % maximum allowable yOffset
'minDecay', (5),... % minimum allowable decay tau
'maxDecay', (500),... % maximum allowable decay tau
'derThresh', (1),... % threshold used to determine if the change of derivative is great enough to separate out two alpha functions
'closestEPSPs', (5),... % second EPSP is not recognized if it is closer than this value to the previous one (in units of samples)
'errThresh', (0.08),... % threshold for standard error above which a fit with multiple alphas is attempted
'dataFilterType', 1,... % 0 = none 1 = windowFilter, 2 = medianFilter, 3 = savitsky-golay
'derFilterType', 3,... % 0 = none 1 = windowFilter, 2 = medianFilter, 3 = savitsky-golay
'dataFilterLength', 5,... % length of data filter
'derFilterLength', 7,... % length of derivative filter
'debugging', (0),... % if set to 1 then debugging figures appear
'outputAxis', 0,... % axis to which any output will be written
'dataStart', 1,... % index of first data point
'forceDisplay', 0,... % forces a graphical output even if other outputs are taken
'alphaFit', 0,... % alpha function fit to the first part of the response
'decayFit', 0,... % exponential fit to the decay
'riseFit', 0); % exponential fit to the rise
else
parameters = struct('minAmp', (5),... % minimum allowable amplitude for alpha functions (in units of samples)
'maxAmp', (5000),... % maximum allowable amplitude
'minTau', (5),... % minimum allowable tau for alpha functions (in units of samples)
'maxTau', (100),... %maximum allowable tau
'minYOffset', (-100),... % minimum allowable yOffset for alpha functions (in units of mV)
'maxYOffset', (-30),... % maximum allowable yOffset
'minDecay', (5),... % minimum allowable decay tau
'maxDecay', (500),... % maximum allowable decay tau
'derThresh', (1),... % threshold used to determine if the change of derivative is great enough to separate out two alpha functions
'closestEPSPs', (5),... % second EPSP is not recognized if it is closer than this value to the previous one (in units of samples)
'errThresh', (0.004),... % threshold for standard error above which a fit with multiple alphas is attempted
'dataFilterType', 1,... % 0 = none 1 = windowFilter, 2 = medianFilter, 3 = savitsky-golay
'derFilterType', 3,... % 0 = none 1 = windowFilter, 2 = medianFilter, 3 = savitsky-golay
'dataFilterLength', 5,... % length of data filter
'derFilterLength', 7,... % length of derivative filter
'debugging', (0),... % if set to 1 then debugging figures appear
'outputAxis', 0,... % axis to which any output will be written
'dataStart', 1,... % index of first data point
'forceDisplay', 0,... % forces a graphical output even if other outputs are taken
'alphaFit', 0,... % alpha function fit to the first part of the response
'decayFit', 0,... % exponential fit to the decay
'riseFit', 0); % exponential fit to the rise
end
end

if firstChar < nargin
if ~mod(nargin - firstChar, 2)
disp('Inappropriate input');
return
else
% specific parameters where passed so figure out which they were
strArray = fieldnames(parameters);
for paramIndex = firstChar - 1:2:numel(varargin)
whichParam = strfind(strArray, varargin{paramIndex});
if find(cell2mat(whichParam))
parameters.(varargin{paramIndex}) = varargin{paramIndex + 1};
else
error([varargin{paramIndex} ' is not a valid parameter'])
end
end
end
end

if nargin > 0 && isa(data, 'char') || nargin == 0
if nargin == 0
else
end

[selection,okTrue] = listdlg('PromptString','Select a channel:',...
'SelectionMode','single',...
'ListSize', [200 80],...
'ListString',zData.protocol.channelNames);
if okTrue
data = zData.traceData(:, selection)';
else
disp('Must select input data');
params = -1;
return
end
end

% make sure that the filter length is odd
if parameters.dataFilterLength / 2 == fix(parameters.dataFilterLength / 2)
parameters.dataFilterLength = parameters.dataFilterLength + 1;
end

if parameters.derFilterLength / 2 == fix(parameters.derFilterLength / 2)
parameters.derFilterLength = parameters.derFilterLength + 1;
end

% reset the error threshold to reflect the noise in the signal
sortedData = sort(data);
if PSPsDown
dataMean = mean(sortedData(size(data, 2) * .7:size(data, 2) * .9)); % mean of top 10-30% of data
dataDev = std(sortedData(size(data, 2) * .6:size(data, 2) * .8)); % standard deviation of top 20-40% of data
% or noise = data - dataFilt; dataDev = std(noise(1:1000));
% dataDev = std(highPass(data)) * 3;
else
dataMean = mean(sortedData(size(data, 2) * .1:size(data, 2) * .3)); % mean of bottom 10-30% of data
dataDev = std(sortedData(size(data, 2) * .2:size(data, 2) * .4)); % standard deviation of bottom 20-40% of data
end

% determine if data is current or voltage clamp
if nargin < 5 && ~hasParameters
if max(data) > 50 || min(data) < -100
% this must be voltage clamp
parameters.minYOffset = dataMean - 200 * dataDev;
parameters.maxYOffset = dataMean + 200 * dataDev;
else
% probably current clamp
parameters.derThresh = .3;
parameters.minAmp = 0.2;
parameters.maxAmp = 50;
parameters.minYOffset = dataMean - 3;
parameters.maxYOffset = dataMean + 3;
end
end

% if the person mixed up the max and min values then switch them
if (PSPsDown && (parameters.minAmp > 0 || parameters.maxAmp > 0)) || (~PSPsDown && (parameters.minAmp < 0 || parameters.maxAmp < 0))
tempAmp = min(sign(-PSPsDown + 0.5) * abs([parameters.minAmp parameters.maxAmp]));
parameters.maxAmp = max(sign(-PSPsDown + 0.5) * abs([parameters.minAmp parameters.maxAmp]));
parameters.minAmp = tempAmp;
end

% filter the raw data
switch parameters.dataFilterType
case 0 % no filtering
dataFilt = data;
case 1 % window filtering
dataFilt = movingAverage(data, parameters.dataFilterLength);
case 2 % median filtering
dataFilt = medfilt1(data, parameters.dataFilterLength);
case 3 % savitsky-golay
dataFilt = sgolayfilt(data, 2, parameters.dataFilterLength);
end

% filter the derivative
tempDiff = diff(dataFilt);
switch parameters.derFilterType
case 0 % no filtering
case 1 % window filtering
case 2 % median filtering
case 3 % savitsky-golay
end
clear tempDiff;

% filter as per Cohen and Miles 2000
outData = zeros(size(dataFilt));
if PSPsDown
for index = 2:length(dataFilt)
if dataDerFilt(index - 1) < 0
outData(index) = outData(index - 1) + dataDerFilt(index - 1);
end
end
else
for index = 2:length(dataFilt)
if dataDerFilt(index - 1) > 0
outData(index) = outData(index - 1) + dataDerFilt(index - 1);
end
end
end

if PSPsDown
% find where derivative of this function is changing from negative to positive
functionDer = diff(outData);
peaks = find((functionDer(2:length(functionDer)) ./ functionDer(1:length(functionDer) -1) < 0 | functionDer(2:length(functionDer)) == 0) & functionDer(1:length(functionDer) - 1) < 0);
else
% find where derivative of this function is changing from positive to negative
functionDer = diff(outData);
peaks = find((functionDer(2:length(functionDer)) ./ functionDer(1:length(functionDer) -1) < 0 | functionDer(2:length(functionDer)) == 0) & functionDer(1:length(functionDer) - 1) > 0);
end

% for each such value greater than derThresh find where the function last
% began to deviate from 0 and call that an event start
numStarts = 0;
whereStarts = ones(length(peaks), 1);
wherePeaks = whereStarts;
for index = 1:length(peaks)
if abs(outData(peaks(index))) > parameters.derThresh
numStarts = numStarts + 1;
whereStarts(numStarts) = peaks(index);
while outData(whereStarts(numStarts)) ~= 0
whereStarts(numStarts) = whereStarts(numStarts) - 1;
end
wherePeaks(numStarts) = peaks(index);
end
end
whereStarts(numStarts + 1) = length(outData);
whereStarts(numStarts + 2:end) = [];
wherePeaks(numStarts + 1:end) = [];

if parameters.debugging
figure
plot(dataFilt);
line(1:length(outData), outData + mean(dataFilt), 'Color', [1 0 0]);
line(1:length(functionDer), functionDer + mean(dataFilt), 'Color', [0 1 0]);
line(wherePeaks, dataFilt(wherePeaks), 'Color', [1 0 1], 'linestyle', 'none', 'marker', '+');
line(whereStarts, dataFilt(whereStarts), 'Color', [0 0 0], 'linestyle', 'none', 'marker', '+');
if PSPsDown
line([1 length(dataFilt)], -[parameters.derThresh parameters.derThresh] + mean(dataFilt), 'Color', [1 0 0], 'linestyle', ':');
else
line([1 length(dataFilt)], [parameters.derThresh parameters.derThresh] + mean(dataFilt), 'Color', [1 0 0], 'linestyle', ':');
end
fitLine = line([0 1], [dataMean dataMean], 'Color', [0 0 0]);
decayLine = line([0 1], [dataMean dataMean], 'Color', [1 0 1]);
riseStart = line([0 1], [dataMean dataMean], 'linestyle', 'none', 'marker', 'x', 'Color', [0.5430 0.2695 0.0742]);
fitProps = annotation('textbox',[.8 .1 .1 .1], 'backgroundcolor', [1 1 1]);
legend('Data', 'Up-Only Function', 'Derivative of U-O Fn', 'Event Peaks', 'Events Starts', 'Threshold for choosing', 'Fit line', 'Decay Fit', 'Rise Line');
end

% set a counter for events
numEvents = 0;
% create output array for alpha best fit parameters
params = zeros(1,4);

for indexEPSP = 1:length(whereStarts) - 1
if whereStarts(indexEPSP + 1) - whereStarts(indexEPSP) < parameters.closestEPSPs
continue % skip this iteration if the next fitting spot is really close
end

if any(~(whereStarts(indexEPSP)  < windows  - windowSize | whereStarts(indexEPSP)  > windows +  2 * windowSize)) % skip this iteration if the found event is outside of the given windows + a buffer of windowSize on each side in case the location is imprecise
% must get tau within +/- 50% of the true value or the fitting
% algorithm does not converge (norm of step goes below boundary
% condition)
tauGuess = max(wherePeaks(indexEPSP) - whereStarts(indexEPSP), parameters.minTau);
ampGuess = dataFilt(wherePeaks(indexEPSP)) - dataFilt(whereStarts(indexEPSP));
whereStop = whereStarts(indexEPSP + 1) - whereStarts(indexEPSP);
if whereStop > 4 * tauGuess
whereStop = 2 * tauGuess;
end

if parameters.alphaFit
[bestParams, residuals] = nlinfit(whereStarts(indexEPSP):whereStarts(indexEPSP) + whereStop, dataFilt(whereStarts(indexEPSP):whereStarts(indexEPSP) + whereStop), @alpha, [ampGuess, tauGuess, whereStarts(indexEPSP), dataFilt(whereStarts(indexEPSP))], statset('MaxIter', 1000, 'FunValCheck', 'off'));
%calculate standard error
bestStdErr = sqrt(residuals' * residuals) / (length(residuals) - 2) / abs(bestParams(1));

%show the fit line
if parameters.debugging
set(fitLine, 'Xdata', whereStarts(indexEPSP):whereStarts(indexEPSP) + whereStop, 'Ydata', alpha(bestParams, whereStarts(indexEPSP):whereStarts(indexEPSP) + whereStop));
xBounds = get(gca, 'Xlim');
if xBounds(2) < whereStarts(indexEPSP) + whereStop
set(gca, 'Xlim', [whereStarts(indexEPSP) - 50 whereStarts(indexEPSP) + 250]);
end
end
else
bestParams = [0 parameters.minTau + 1 whereStarts(indexEPSP) dataFilt(whereStarts(indexEPSP))];
if PSPsDown
bestParams(1) = min(dataFilt(whereStarts(indexEPSP):whereStarts(indexEPSP) + whereStop));
else
bestParams(1) = max(dataFilt(whereStarts(indexEPSP):whereStarts(indexEPSP) + whereStop));
end
bestParams(1:2) = bestParams(1:2) - [dataFilt(whereStarts(indexEPSP)) 0];
bestStdErr = 0;
end

if ~(parameters.minAmp < bestParams(1)  && parameters.maxAmp * 2.7183 > bestParams(1)  &&... % is the amplitude ok?
parameters.minTau < bestParams(2)  && parameters.maxTau > bestParams(2)  &&... % is the tau ok?
whereStarts(indexEPSP + 1) - bestParams(3) > bestParams(2) * 0.5 &&... % do we really have enough to fit a function to it?
parameters.minYOffset < bestParams(4) && parameters.maxYOffset > bestParams(4)) ||...% is it in the right place?
~any(~(bestParams(3)  < windows  | bestParams(3)  > windows + windowSize)) % is it in the window that we're looking in
bestStdErr = inf; % set high so that any other fit is accepted since this one has parameters out of range
end

% refit variables
if bestStdErr < parameters.errThresh

bestParams(1) = dataFilt(wherePeaks(indexEPSP)) - dataFilt(whereStarts(indexEPSP));
bestParams(4) = dataFilt(whereStarts(indexEPSP));
% find where the trace is at 20% of max rise
howLong = min([45 wherePeaks(indexEPSP) - whereStarts(indexEPSP)]);
[loc loc] = max(1 ./ (dataFilt(whereStarts(indexEPSP):whereStarts(indexEPSP) + howLong) - (bestParams(4) + .10 * bestParams(1))));
% make sure that this point is the last time we have this situation (is on the rising phase)

% linearly interpolate
pointDrop = dataFilt(loc + whereStarts(indexEPSP) - 1) - dataFilt(loc + whereStarts(indexEPSP));

% sometimes the pointDrop is zero and this creates an infinite rise time
if pointDrop == 0 || pointDrop * (PSPsDown - .5) < 0 % if this product is negative then we're dealing with a jagged spot
loc = loc + 1;
pointDrop = dataFilt(loc + whereStarts(indexEPSP) - 1) - dataFilt(loc + whereStarts(indexEPSP));
twentyRise = wherePeaks(indexEPSP) - whereStarts(indexEPSP) - (loc + 1 - ((bestParams(4) + .10 * bestParams(1) - dataFilt(whereStarts(indexEPSP) + loc)) / pointDrop));
else
twentyRise = wherePeaks(indexEPSP) - whereStarts(indexEPSP) - (loc + 1 - ((bestParams(4) + .10 * bestParams(1) - dataFilt(whereStarts(indexEPSP) + loc)) / pointDrop));
end

if parameters.riseFit
startPoint = [];
twentyRise = round(twentyRise);
if wherePeaks(indexEPSP) - twentyRise <= whereStarts(indexEPSP)
whereStarts(indexEPSP) = wherePeaks(indexEPSP) - twentyRise - 1;
end
%                 if twentyRise > 10 % if this is a slow one then fit it with a more advanced method
%                     % fit the portion of the rise that we think is definitley clean
%                     [estimates estimates estimates estimates] = fitDecayDouble(dataFilt(wherePeaks(indexEPSP) - twentyRise:wherePeaks(indexEPSP) - 2));
%                     if ~isnan(estimates)
%                         xdata = ((whereStarts(indexEPSP) - wherePeaks(indexEPSP)) + twentyRise:twentyRise - 2);
%                         errorVector = (estimates(1) + sign(bestParams(3)) * estimates(2) .* exp(xdata./estimates(3)) + sign(bestParams(3)) * estimates(4) .* exp(xdata./estimates(5)) - dataFilt(round(xdata + wherePeaks(indexEPSP) - twentyRise))) .^ 2;
%                         errBound = mean(errorVector(length(xdata) - twentyRise + 1:end)) + 3 * std(errorVector(length(xdata) - twentyRise + 1:end));
%
%                         % find the first two points where we deviate from the curve more than we do in the clean part
%                         startPoint = find(errorVector(1:length(xdata) - twentyRise) > errBound & errorVector(2:length(xdata) - twentyRise + 1) > errBound, 1, 'last') + 1;
%                         if ~isempty(startPoint) && parameters.debugging
%                             xdata = xdata(startPoint:end);
%                             FittedRise = estimates(1) + sign(bestParams(3)) * estimates(2) .* exp(xdata./estimates(3)) + sign(bestParams(3)) * estimates(4) .* exp(xdata./estimates(5));
%                         end
%                     end
%                 end

% if the above fit failed or we didn't have enough point to attempt it then
if isempty(startPoint) && twentyRise > 3
% fit the portion of the rise that we think is definitley clean
[estimates estimates estimates] = fitDecaySingle(dataFilt(wherePeaks(indexEPSP) - twentyRise:wherePeaks(indexEPSP) - 2));
if ~isnan(estimates)
xdata = ((whereStarts(indexEPSP) - wherePeaks(indexEPSP)) + twentyRise:twentyRise - 2);
errorVector = (estimates(1) + sign(bestParams(3)) * estimates(2) .* exp(xdata./estimates(3)) - dataFilt(round(xdata + wherePeaks(indexEPSP) - twentyRise))) .^ 2;
errBound = mean(errorVector(length(xdata) - twentyRise + 1:end)) + 3 * std(errorVector(length(xdata) - twentyRise + 1:end));

% find the first two points where we deviate from the curve more than we do in the clean part
startPoint = find(errorVector(1:length(xdata) - twentyRise) > errBound & errorVector(2:length(xdata) - twentyRise + 1) > errBound, 1, 'last') + 1;
if ~isempty(startPoint) && parameters.debugging
xdata = xdata(startPoint:end);
FittedRise = estimates(1) + sign(bestParams(3)) * estimates(2) .* exp(xdata./estimates(3));
end
end
end

if isempty(startPoint)
startPoint = loc;
FittedRise = ones(twentyRise, 1) * dataFilt(whereStarts(indexEPSP));
errBound = 0;
end
bestParams(3) = startPoint  + whereStarts(indexEPSP);
else
bestParams(3) = whereStarts(indexEPSP);
end
bestParams(3) = whereStarts(indexEPSP);
bestParams(2) = wherePeaks(indexEPSP) - bestParams(3); %whereStarts(indexEPSP);
if parameters.decayFit
% change bestParams(4) to the decay tau
whichData = round(bestParams(3) + bestParams(2) * 1.5:min([bestParams(3) + bestParams(2) * 1.5 + 100 whereStarts(indexEPSP + 1)]));
if length(whichData) > 5
[bestParams(4) FittedDecay] = fitDecaySingle(dataFilt(whichData));
residuals = FittedDecay' - dataFilt(whichData);
if sqrt(residuals' * residuals) / (length(whichData) - 2) / abs(bestParams(1)) > parameters.errThresh
bestParams(4) = NaN;
FittedDecay = NaN;
end
else
bestParams(4) = NaN;
FittedDecay = NaN;
end

if isnan(bestParams(4))
% make a manual guess
if PSPsDown == 0
tempDecay = find(dataFilt(wherePeaks(indexEPSP) + 2:min([wherePeaks(indexEPSP) + 90 whereStarts(indexEPSP + 1)])) < (dataFilt(wherePeaks(indexEPSP)) - log(2) * (dataFilt(wherePeaks(indexEPSP)) - min(dataFilt(wherePeaks(indexEPSP) + 2:min([wherePeaks(indexEPSP) + 90 whereStarts(indexEPSP + 1)]))))), 1, 'first');
else
tempDecay = find(dataFilt(wherePeaks(indexEPSP) + 2:min([wherePeaks(indexEPSP) + 90 whereStarts(indexEPSP + 1)])) > (dataFilt(wherePeaks(indexEPSP)) - log(2) * (dataFilt(wherePeaks(indexEPSP)) - max(dataFilt(wherePeaks(indexEPSP) + 2:min([wherePeaks(indexEPSP) + 90 whereStarts(indexEPSP + 1)]))))), 1, 'first');
end
if ~isempty(tempDecay)
bestParams(4) = tempDecay;
FittedDecay = ones([1 length(whichData)]) * dataFilt(wherePeaks(indexEPSP));
else
bestParams(4) = parameters.minDecay - 1;
end
end

if parameters.debugging
set(fitProps, 'String',  [{['\bf Amp: ' num2str(bestParams(1))]};{[' Rise: ' num2str(bestParams(2))]};{[' Start: ' num2str(bestParams(3))]};{[' Decay: ' num2str(bestParams(4))]}]);
set(decayLine, 'Xdata', whichData, 'Ydata', FittedDecay);
set(gca, 'ylimmode', 'auto');
set(gca, 'xlim', [max([1 whereStarts(indexEPSP) - 10]) min([max([1 whereStarts(indexEPSP) - 10]) + max([bestParams(2) * 5 200]) length(data)])]);
savedLims = get(gca, 'ylim');
set(riseStart, 'xdata', bestParams(3), 'ydata', dataFilt(round(bestParams(3))));
set(gca, 'ylim', savedLims);
pause
end
else
bestParams(4) = parameters.minDecay + 1;
end

if bestParams(4) >= parameters.minDecay && bestParams(4) <= parameters.maxDecay &&...
bestParams(1) >= parameters.minAmp && bestParams(1) <= parameters.maxAmp
%                 bestParams(2) = wherePeaks(indexEPSP) - whereStarts(indexEPSP);
params(numEvents + 1, :) = bestParams;
if parameters.decayFit
decayFits{numEvents + 1} = FittedDecay;
end
%dataFilt  = dataFilt - alpha([bestParams(:, 1:3) zeros(size(bestParams, 1), 1)], 1:length(dataFilt)); % subtract off found data (0 is to avoid subtracting off the yOffset)
numEvents = numEvents + 1;
end
end
end
end

if ~parameters.decayFit
decayFits = nan(size(params));
params(:,4) = nan;
end

% show how well our found trace fits the original data
if nargout == 0 || parameters.forceDisplay
if parameters.outputAxis == 0 || ~ishandle(parameters.outputAxis)
if max(data) > 50 || min(data) < -100
parentHandle = newScope(data, 1:length(data), 'Trace I');
else
parentHandle = newScope(data, 1:length(data), 'Trace V');
end
parentHandle = parentHandle.axes;
else
parentHandle = parameters.outputAxis;
end

set(0, 'currentfigure', get(parentHandle, 'parent'));
set(gcf, 'currentAxes', parentHandle);

% set the time variables to be in the time units of the current axis
timeDiff = get(gca, 'child');
timeDiff = timeDiff(end - 2);
timeDiff = get(timeDiff, 'xdata');
timeDiff = timeDiff(2) - timeDiff(1);
params(:, 2:4) = params(:, 2:4) * timeDiff;
parameters.dataStart = parameters.dataStart * timeDiff;

try
for i = 1:size(params, 1)
% plot starts
line(parameters.dataStart + params(i, 3), data(round(params(i, 3)/timeDiff)), 'Color', [0 0 1], 'linestyle', 'none', 'marker', '+', 'buttondownfcn', ['set(get(get(gcbo, ''parent''), ''userdata''), ''string'', ''Start time = ' sprintf('%6.1f', params(i,3)) ''')']);
% plot amplitudes
line(parameters.dataStart + [params(i, 3) + params(i, 2) params(i, 3) + params(i, 2)], [data(round(params(i, 3)/ timeDiff + params(i, 2) / timeDiff)) - params(i, 1) data(round(params(i, 3)/ timeDiff + params(i, 2)/ timeDiff))], 'Color', [0 1 0], 'buttondownfcn', ['set(get(get(gcbo, ''parent''), ''userdata''), ''string'', {''Rise time = ' sprintf('%4.1f', params(i,2)) '''; ''Amplitude = ' sprintf('%7.2f', params(i,1)) '''})']);
if parameters.decayFit
% plot decay tau
line(parameters.dataStart + ((params(i, 3)/timeDiff + params(i, 2)/timeDiff * 1.5:params(i, 3)/timeDiff + params(i, 2)/timeDiff * 1.5 + numel(decayFits{i}) - 1) * timeDiff), decayFits{i}, 'Color', [1 0 0], 'buttondownfcn', ['set(get(get(gcbo, ''parent''), ''userdata''), ''string'', ''Decay tau = ' sprintf('%6.2f', params(i,4)) ''')']);
end
end
catch
% at least return something

end
end