%pool spike triggered averages of all units from one cell together
 close all
clearvars
load('20160825_Stage12_output_v6.mat','goodcells')

load('20160825_Stage8_V6_output.mat')

w1=-4;% w1 second delay before each spike and w2 s delay after
w2=4;
res2=0.04;
interval=w1:res2:w2;

for i=1:length(spktrig_EODr)
    %calculate the mean spktrigEOD
    temp=nanmean(spktrig_EODr{i});
    temp_SD=nanstd(spktrig_EODr{i});
    temp_shuff=nanmean(spktrig_EODr_shuff{i});
    temp_SD_shuff=nanstd(spktrig_EODr_shuff{i});
    
    
    temp_speed=nanmean(spktrig_speed{i});
    temp_speed_SD=nanstd(spktrig_speed{i});
    temp_speed_shuff=nanmean(spktrig_speed_shuff{i});
    temp_speed_SD_shuff=nanstd(spktrig_speed_shuff{i});
    
    
    temp_sampdens=nanmean(spktrig_sampdens{i});
    temp_sampdens_SD=nanstd(spktrig_sampdens{i});
    temp_sampdens_shuff=nanmean(spktrig_sampdens_shuff{i});
    temp_sampdens_SD_shuff=nanstd(spktrig_sampdens_shuff{i});
    %find the largest peak in the mean
    [Pks,Ind]=findpeaks(temp);
    AvePkInd=Ind(Pks==max(Pks));
    %is this peak significantly larger than the mean
    [a,b]=size(spktrig_EODr{i});
    for j=1:a
        amp(j)=spktrig_EODr{i}(j,AvePkInd);
    end
    allmean=nanmedian(nanmedian(spktrig_EODr{i}));
    [p,h]=signrank(amp,allmean)
    clear amp
    EOD_P(i,2)=p;
    EOD_P(i,1)=interval(AvePkInd);
    %compare the median to the median of the shuffled spikes
    trial_mean=nanmedian(spktrig_EODr{i},2);
    allmean_shuff=nanmedian(nanmedian(spktrig_EODr_shuff{i}))
    [p,h]=signrank(trial_mean,allmean_shuff)
    EOD_Shuff_P(i)=p;
    %do the same for speed peak
    
    [Pks,Ind]=findpeaks(temp_speed);
    AvePkInd=Ind(Pks==max(Pks));
    %is this peak significantly larger than the mean
    [a,b]=size(spktrig_speed{i});
    for j=1:a
        amp(j)=spktrig_speed{i}(j,AvePkInd);
    end
    allmean=nanmedian(nanmedian(spktrig_speed{i}));
    [p,h]=signrank(amp,allmean)
    clear amp
    speed_P(i,2)=p;
    speed_P(i,1)=interval(AvePkInd);
    
    if EOD_P(i,2)<0.05
        figure(i)
        subplot(1,3,1)
        plot(interval,temp,'r')
        hold on
        plot(interval, temp+(temp_SD/sqrt(a)),'r:')
        plot(interval, temp-(temp_SD/sqrt(a)),'r:')
        plot(interval,temp_shuff,'k')
        hold on
        plot(interval, temp_shuff+temp_SD_shuff/sqrt(a*100),'k:')
        plot(interval, temp_shuff-temp_SD_shuff/sqrt(a*100),'k:')
        xlim([-4,4])
        
        %plot speed     
        subplot(1,3,2)
        plot(interval,temp_speed,'r')
        hold on
        plot(interval, temp_speed+temp_speed_SD/sqrt(a),'r:')
        plot(interval, temp_speed-temp_speed_SD/sqrt(a),'r:')
        plot(interval,temp_speed_shuff,'k')
        hold on
        plot(interval, temp_speed_shuff+temp_speed_SD_shuff/sqrt(a*100),'k:')
        plot(interval, temp_speed_shuff-temp_speed_SD_shuff/sqrt(a*100),'k:')
        xlim([-4,4])

        %plot sampdens
        subplot(1,3,3)
        plot(interval,temp_sampdens,'r')
        hold on
        plot(interval, temp_sampdens+temp_sampdens_SD/sqrt(a),'r:')
        plot(interval, temp_sampdens-temp_sampdens_SD/sqrt(a),'r:')
        plot(interval,temp_sampdens_shuff,'k')
        hold on
        plot(interval, temp_sampdens_shuff+temp_sampdens_SD_shuff/sqrt(a*100),'k:')
        plot(interval, temp_sampdens_shuff-temp_sampdens_SD_shuff/sqrt(a*100),'k:')
        xlim([-4,4])
        %title(strcat('Spk Trig Aves U',goodcells(i),'p=',num2str(EOD_P(i,2))))
        %saveas(gcf,strcat('Spk Trig Aves U',goodcells(i),'.pdf'))
    end

end

%save('20160825_Stage8_V6_output.mat','EOD_P','speed_P','EOD_Shuff_P','goodcells','-append')