function [ThresVel]=getFiringThesholdNewVel(peak,output1Vel,output2Vel,FigHandlVel)
%%
%peak = on which side the activity is 1 is right 0 is left         
% output1Vel= 1st output from averageEye_V3vel                                               
% output2Vel= 2nd output from averageEye_V3vel                                              
% ThresVel= calculated firing threshold                               
% ThresIndexVel= index of calculated threshold                        
% FigHandlVel = figure handle to plot firing thres in                 
%% calculate statistical firing threshold 

ThresVel=nan;

if peak==0 % max left
    iterbasesta=isnan(output1Vel(:,1));
    iterbasesta2=find(iterbasesta==0,1,'last');
    pvaluesIter=nan(size(output1Vel,1),1); 
    hvaluesIter=nan(size(output1Vel,1),1);
    for g=1:(iterbasesta2-2-(find(iterbasesta==0,1,'first'))) %highest x value with big enough sample size (iterbasesta2))  - 2 (from initial Baseline) - lowest xvalue with enough sample size = number of times I can run ranksum
    iterbase=reshape(output2Vel(:, iterbasesta2-1-g:iterbasesta2), [], 1); % highest x value with big enough sample size -1 -g (starts with 1)= (for g=1) 3 samples as baseline + increasing with loop 
        try
        [pvaluesIter(iterbasesta2-2-g,1) hvaluesIter(iterbasesta2-2-g,1)]=...
            ranksum(output2Vel(:,iterbasesta2-2-g), iterbase, 'tail', 'right','alpha', 0.05/g); % always one more right than baseline
        catch
            pvaluesIter(iterbasesta2-2-g,1)=NaN;
            hvaluesIter(iterbasesta2-2-g,1)=NaN;
        end
    end
elseif peak==1 % max right
    iterbasesta=isnan(output1Vel(:,1));
    iterbasesta2=find(iterbasesta==0,1,'first');
    pvaluesIter=nan(size(output1Vel,1),1); 
    hvaluesIter=nan(size(output1Vel,1),1);
    for g=1:(find(iterbasesta==0,1,'last')-iterbasesta2-2) %highest x value with big enough sample size  - 2 (from initial Baseline) - lowest xvalue with enough sample size = number of times I can run ranksum
    iterbase=reshape(output2Vel(:,iterbasesta2:iterbasesta2+1+g),[],1);
        try
        [pvaluesIter(iterbasesta2+2+g,1) hvaluesIter(iterbasesta2+2+g,1)]=...
            ranksum(output2Vel(:,iterbasesta2+2+g),iterbase, 'tail', 'right','alpha', 0.05/g); % always one more right than baseline
        catch
            pvaluesIter(iterbasesta2+2+g,1)=NaN;
            hvaluesIter(iterbasesta2+2+g,1)=NaN;
        end 
    end
end

if ~isnan(peak)
for k=1:size(hvaluesIter,1) % set 0 as NaN's so that only significant values := 1
    if hvaluesIter(k,1)==0
    hvaluesIter(k,1)=NaN;
    end
    if hvaluesIter(k,1)==1
       hvaluesIter(k,1)=output1Vel(k,1);
    end
end
end

%get x/y values to plot firing thres/ significant bins in figure
if peak==1
indexvector=find(~isnan(hvaluesIter));
threshrunning=0;
 for ivk=1:size(indexvector,1)
     if ivk~=size(indexvector,1)
        if (ismember(indexvector(ivk)+1,indexvector) ||...
            ismember(indexvector(ivk)+2,indexvector)) &&...
            threshrunning==0
        indexvectorEyeVel=-4.75:0.5:4.75;
        ThresVel=indexvectorEyeVel(indexvector(ivk));
        threshrunning=threshrunning+1;
        ThresIndexVel=indexvector(ivk);
        end
     elseif ivk==size(indexvector,1) 
         if threshrunning==0 && length(indexvector)>=3
            indexvectorEyeVel=-4.75:0.5:4.75;
            ThresVel=indexvectorEyeVel(indexvector(ivk));
            threshrunning=threshrunning+1;
            ThresIndexVel=indexvector(ivk);
         end
     else
         FORCEDERROR
     end
 end
elseif peak==0  
indexvector=find(~isnan(hvaluesIter));
threshrunning=0;
 for ivk=size(indexvector,1):-1:1
     if ivk~=1
        if (ismember(indexvector(ivk)-1,indexvector) ||...
            ismember(indexvector(ivk)-2,indexvector)) &&...
            threshrunning==0
            indexvectorEyeVel=-4.75:0.5:4.75;
            ThresVel=indexvectorEyeVel(indexvector(ivk));
            threshrunning=threshrunning+1;
            ThresIndexVel=indexvector(ivk);
        end
     elseif ivk==1
        if threshrunning==0 && length(indexvector)>=3
            indexvectorEyeVel=-4.75:0.5:4.75;
            ThresVel=indexvectorEyeVel(indexvector(ivk));
            threshrunning=threshrunning+1;
            ThresIndexVel=indexvector(ivk);
        end
     else 
         FORCEDERROR
     end
 end
end
    %plot x for significant bins
    if ~isnan(peak) && exist('ThresVel','var') && exist('ThresIndexVel','var')
    subplot(FigHandlVel)
    scatter(-4.75:0.5:4.75,hvaluesIter+0.75,'Marker','x',...
    'MarkerEdgeColor','k')
    %plot circle for firing threshold
    scatter(ThresVel,hvaluesIter(ThresIndexVel)+0.75)
    end

end % function end
