Looking-to-cue-sides-and-modulation-of-value-in-OFC/code/neuralAnalysesDelayShiftAligned.m

addpath(genpath('../../../my-tools'));
addpath(genpath('../../../external-tools'));

clearvars -except subjsData poolstack* pooltw*;
%close all; clc;

reload=0; recomp_shLRalgn=1;
recomputelm=1;
Nlag=10; Nspan=200; %Nspan=200;
Nshf=10;


selflag=1; %LookL during offer1 is most common
%selflag=4; %LookR during offer2 is most common
%selflag=2; %LookR during offer1 is less common
%selflag=3; %LookL during offer2 is less common

%caseflag=1; % perform analysis saccade-aligned forcing stare after first shift
%caseflag=2; % perform analysis saccade-aligned ignoring gaze after first shift (removes data if second shifts are detected)
%caseflag=3; % perform analysis for unilateral staring only trials
%caseflag=4; % same as caseflag=1 but forcing offerLL/LR
%caseflag=5; % same as caseflag=2 but forcing offerLL/LR
%caseflag=6; % same as caseflag=3 but forcing offerLL/LR
for selflag=4%[1 4]
for caseflag=2
clearvars -except subjsData reload recomp_shLRalgn recomputelm Nlag Nspan Nshf selflag caseflag; close all; clc;

stareflag=(caseflag==1 | caseflag==4); %forces stare after saccade
casesstrs={'1_shift_NOoffLR','2_shift_NOstare_NOoffLR','3_NOshift_NOoffLR','4_shift','5_shift_NOstare','6_NOshift'};

disp(['running selflag=' num2str(selflag) '; caseflag=' num2str(caseflag) '; (' casesstrs{caseflag} ')']);


switch selflag
    case 1 
        %most_recent_folder='../../Figures/Temp/NoBaselineNorm/evLevR_off1L_del1LR_sacc/';
        most_recent_folder=['../../Figures/Temp/NoBaselineNorm/evLevR_off_del_focus_minus200ms_apr24/evLevR_off1L_del1LR_' casesstrs{caseflag} '/'];
        if ~exist(most_recent_folder,'dir'); mkdir(most_recent_folder); mkdir([most_recent_folder 'span' num2str(Nspan) 'ms_lag' num2str(Nlag) 'ms']); end
        %most_recent_folder='../../Figures/Temp/NoBaselineNorm/evLevR_off1L_del1LR_NOshift_NOoff1L/';
        selectedTimes='((ntce_ds(3,1)-46):ntce_ds(3,2))';
        %selectedCases{1}='find((cciso1LL) & (ccepx(:,tt)<0))';
        %selectedCases{2}='find((cciso1LL) & (ccepx(:,tt)>0))';
    case 2
        %most_recent_folder='../../Figures/Temp/NoBaselineNorm/evLevR_off1R_del1LR_sacc/';
        most_recent_folder='../../Figures/Temp/NoBaselineNorm/evLevR_off1R_del1LR_sacc/';
        selectedTimes='(ntce_ds(3,1):ntce_ds(3,2))';
        %selectedCases{1}='find((~cciso1LL) & (ccepx(:,tt)<0))';
        %selectedCases{2}='find((~cciso1LL) & (ccepx(:,tt)>0))';
    case 3
        %most_recent_folder='../../Figures/Temp/NoBaselineNorm/evLevR_off2L_del2LR_sacc/';
        most_recent_folder='../../Figures/Temp/NoBaselineNorm/evLevR_off2L_del2LR_shift/';
        %selectedTimes='(ntce_ds(5,1):ntce_ds(5,2))';
        selectedTimes='(ntce_ds(5,1):ntce_ds(6,1)-2)'; % FULL DELAY2
        %selectedCases{1}='find((cciso2LL) & (ccepx(:,tt)<0))';
        %selectedCases{2}='find((cciso2LL) & (ccepx(:,tt)>0))';
    case 4
        %most_recent_folder='../../Figures/Temp/NoBaselineNorm/evLevR_off2R_del2LR_sacc/';
        most_recent_folder=['../../Figures/Temp/NoBaselineNorm/evLevR_off_del_focus_minus200ms_apr24/evLevR_off2R_del2LR_' casesstrs{caseflag} '/'];
        if ~exist(most_recent_folder,'dir'); mkdir(most_recent_folder); mkdir([most_recent_folder 'span' num2str(Nspan) 'ms_lag' num2str(Nlag) 'ms']); end
        %most_recent_folder='../../Figures/Temp/NoBaselineNorm/evLevR_off2R_del2LR_NOshift_NOoff2R/';
        selectedTimes='((ntce_ds(5,1)-46):ntce_ds(5,2))';
        %selectedTimes='(ntce_ds(5,1):ntce_ds(6,1)-2)'; % FULL DELAY2
        %selectedCases{1}='find((~cciso2LL) & (ccepx(:,tt)<0))';
        %selectedCases{2}='find((~cciso2LL) & (ccepx(:,tt)>0))';
end

% for off1L del1L/R
% for off1R del1L/R
% for off2L del2L/R
% for off2R del2L/R


if reload
% if ~exist('subjsData','var')
%     subjsData(1).sbdata=getSubjectData(1);
%     subjsData(2).sbdata=getSubjectData(2);
% end

if ~exist('poolstackvars','var')
    load latestpoolstackvars_14feb22.mat
end
if ~exist('pooltweyepos','var')
    load latestpooltweyepos_14feb22.mat
end
if ~exist('pooltwspksq','var')
    load latestpooltwspks_14feb22.mat
end

twLabels={'preoffer1','offer1','delay1','offer2','delay2','refixate','choice-go','ch-hold'};

nccs=[ 51 59 24 29 18 32  9 26];
ntrs=arrayfun(@(xx) size(poolstackvars(xx).offer1sd,1),1:8);

Nccs=sum(nccs);
Ntrs12=max(ntrs);
Ntps12=arrayfun(@(tw) size(pooltweyepos(tw).posX,3), 1:8);
Ntps12(1)=400;

pooloffer1sd=[]; poolofferLev=[]; poolofferRev=[]; pooloffer1ev=[]; pooloffer2ev=[]; pooloffer1rw=[]; pooloffer2rw=[]; poolofferchR=[];
for sn=1:8
    pooloffer1sd=cat(2, pooloffer1sd, repmat(cat(1,poolstackvars(sn).offer1sd,nan(Ntrs12-ntrs(sn),1)),[1 nccs(sn)]));
    poolofferLev=cat(2, poolofferLev, repmat(cat(1,poolstackvars(sn).offerLev,nan(Ntrs12-ntrs(sn),1)),[1 nccs(sn)]));
    poolofferRev=cat(2, poolofferRev, repmat(cat(1,poolstackvars(sn).offerRev,nan(Ntrs12-ntrs(sn),1)),[1 nccs(sn)]));
    pooloffer1ev=cat(2, pooloffer1ev, repmat(cat(1,poolstackvars(sn).offer1ev,nan(Ntrs12-ntrs(sn),1)),[1 nccs(sn)]));
    pooloffer2ev=cat(2, pooloffer2ev, repmat(cat(1,poolstackvars(sn).offer2ev,nan(Ntrs12-ntrs(sn),1)),[1 nccs(sn)]));
    pooloffer1rw=cat(2, pooloffer1rw, repmat(cat(1,poolstackvars(sn).offer1rw,nan(Ntrs12-ntrs(sn),1)),[1 nccs(sn)]));
    pooloffer2rw=cat(2, pooloffer2rw, repmat(cat(1,poolstackvars(sn).offer2rw,nan(Ntrs12-ntrs(sn),1)),[1 nccs(sn)]));
    poolofferchR=cat(2, poolofferchR, repmat(cat(1,poolstackvars(sn).chooseR,nan(Ntrs12-ntrs(sn),1)),[1 nccs(sn)]));
end
pooleposxmat=cat(3,pooltweyepos(1).posX,pooltweyepos(2).posX,pooltweyepos(3).posX,pooltweyepos(4).posX,...
    pooltweyepos(5).posX,pooltweyepos(6).posX,pooltweyepos(7).posX,pooltweyepos(8).posX);
poolspksqmat=cat(3,pooltwspksq(1).spkseqs,pooltwspksq(2).spkseqs,pooltwspksq(3).spkseqs,pooltwspksq(4).spkseqs,...
    pooltwspksq(5).spkseqs,pooltwspksq(6).spkseqs,pooltwspksq(7).spkseqs,pooltwspksq(8).spkseqs);
poolspksqmat=cat(3,nan(size(pooltwspksq(1).spkseqs)),nan(size(pooltwspksq(2).spkseqs)),pooltwspksq(3).spkseqs,...
    nan(size(pooltwspksq(4).spkseqs)),pooltwspksq(5).spkseqs,nan(size(pooltwspksq(6).spkseqs)),nan(size(pooltwspksq(7).spkseqs)),nan(size(pooltwspksq(8).spkseqs)));

poolcovhist=cat(2,pooltwspksq(1).covhist,pooltwspksq(2).covhist,pooltwspksq(3).covhist,pooltwspksq(4).covhist,...
    pooltwspksq(5).covhist,pooltwspksq(6).covhist,pooltwspksq(7).covhist,pooltwspksq(8).covhist);
%clearvars pooltwspksq pooltweyepos

%fchg=find(poolcovhist>=.8); fchl=find(poolcovhist<.8);
fchg=find(poolcovhist>=.999); fchl=find(poolcovhist<.999);
ntcv=[fchg(diff([fchg sum(Ntps12)])>1)' fchl([diff(fchl) sum(Ntps12)]>1)'];
ntce=[[1; 401; ntcv(1:end-1,2)] [400; ntcv(:,1)]];
vtps=ones(1,sum(Ntps12)); for jj=1:7; vtps(ntcv(jj,1):ntcv(jj,2))=nan; end

%if sparsity==1 it means 1spike/ms that is 1000 spikes/s
%nanmean(vectorise((nansum(poolspksqmat,3))./(nansum(poolspksqmat>=0,3))))*1000
%nanstd(vectorise((nansum(poolspksqmat,3))./(nansum(poolspksqmat>=0,3))))*1000

%% smooth by movmean filtering, downsample to 20ms, 50Hz resolution
Ntps_ds=ceil(sum(Ntps12)/Nlag);
tmax_ds=nanmovmean(1:size(pooleposxmat,3),Nlag,Nlag);
poolspksqmat_ds=nan(Ntrs12,Nccs,length(tmax_ds));
pooleposxmat_ds=nan(Ntrs12,Nccs,length(tmax_ds));

for cc=1:Nccs
    poolspksqmat_ds(:,cc,:)=nanmovmean(squeeze(poolspksqmat(:,cc,:))',Nlag,[0 Nspan])';
    pooleposxmat_ds(:,cc,:)=nanmovmean(squeeze(pooleposxmat(:,cc,:))',Nlag,Nlag)';
end

clearvars poolspksqmat pooleposxmat


%%


poolcovhist_ds=nanmovmean(poolcovhist,Nlag,Nlag);
%fcgh_ds=find(poolcovhist_ds>=.8); fchl_ds=find(poolcovhist_ds<.8);
fcgh_ds=find(poolcovhist_ds>=.999); fchl_ds=find(poolcovhist_ds<.999);
ntcv_ds=[fcgh_ds(diff([fcgh_ds Ntps_ds])>1)' fchl_ds([diff(fchl_ds) Ntps_ds]>1)'];
vtps_ds=ones(1,Ntps_ds);  for jj=1:7; vtps_ds(ntcv_ds(jj,1):ntcv_ds(jj,2))=nan; end
ntce_ds=[[1 ceil(400/Nlag)+1 find(vtps_ds(2:end)+isnan(vtps_ds(1:end-1))==2)+1]' [ceil(400/Nlag) find(vtps_ds(1:end-1)+isnan(vtps_ds(2:end))==2)]'];
Nvtps_ds=sum(vtps_ds==1);

% Apply movmean independently across time windows
% checktps_ds=zeros(1,length(poolcovhist_ds));
% for tw=1:8
%     tpvec=zeros(1,length(poolcovhist));
%     tpvec(ntce(tw,1):ntce(tw,2))=1;
%     tpvec_ds=zeros(1,length(poolcovhist_ds));
%     tpvec_ds(ntce_ds(tw,1):ntce_ds(tw,2))=1;
%     for cc=1:Nccs
%         nnspktemp=nanmovmean(squeeze(poolspksqmat(:,cc,tpvec==1))',Nlag,[0 Nspan])';
%         poolspksqmat_ds(:,cc,tpvec_ds==1)=nnspktemp(:,1:sum(tpvec_ds==1));
%         nnepxtemp=nanmovmean(squeeze(pooleposxmat(:,cc,tpvec==1))',Nlag,Nlag)';
%         pooleposxmat_ds(:,cc,tpvec_ds==1)=nnepxtemp(:,1:sum(tpvec_ds==1));
%     end
%     checktps_ds(tpvec_ds==1)=checktps_ds(tpvec_ds==1)+tw;
% end

vtmax_ds=rmnans(tmax_ds.*vtps_ds);%find(rmnans(vtps_ds));
rtmax_ds=1:length(vtmax_ds);
rtmes_ds=arrayfun(@(jj) length(ntce_ds(jj,1):ntce_ds(jj,2)), 1:8);

% baseline normalization
%poolspksqmat_ds_bln=poolspksqmat_ds-repmat(nanmean(poolspksqmat_ds(:,:,1:floor(400/Nlag)),3),[1 1 Ntps_ds]);
%clearvars poolspksqmat_ds

% baseline normalization
%poolspksqmat_ds_bln=poolspksqmat_ds./repmat(nanmean(poolspksqmat_ds(:,:,1:floor(400/Nlag)),3)+eps,[1 1 Ntps_ds]);
%poolspksqmat_ds_bln=poolspksqmat_ds-repmat(nanmean(nanmean(poolspksqmat_ds(:,:,1:floor(400/Nlag)),3),1),[Ntrs12 1 Ntps_ds]);

% No baseline normalization
poolspksqmat_ds_bln=poolspksqmat_ds;
else
    load pool_all_light.mat
end

%%
% fullfig();
% subplot(1,2,1);
% spy(vectorisen(~isnan(pooleposxmat_ds_bln(1,cumsum(nccs),:)),1:2));
% subplot(1,2,2);
% spy(vectorisen(~isnan(poolspksqmat_ds_bln(1,cumsum(nccs),:)),1:2));
% 
% fulldel2time=ntce_ds(5,1):ntce_ds(6,1)-20-2;
% 
% fullfig();
% subplot(1,2,1);
% spy(vectorisen(~isnan(pooleposxmat_ds_bln(1,cumsum(nccs),fulldel2time)),1:2));
% subplot(1,2,2);
% spy(vectorisen(~isnan(poolspksqmat_ds_bln(1,cumsum(nccs),fulldel2time)),1:2));


%%
off1time=ntce_ds(2,2)-(19:-1:0);
del1time=ntce_ds(3,2)-(19:-1:0);
off2time=ntce_ds(4,2)-(19:-1:0);
del2time=ntce_ds(5,2)-(19:-1:0);

%fulldel1time=(ntce_ds(3,1)):ntce_ds(3,2);
fulldel1time=[(ntce_ds(3,1)-46):(ntce_ds(3,1)-7) ntce_ds(3,1):ntce_ds(3,2)];
%fulldel2time=(ntce_ds(5,1)):ntce_ds(6,1)-20-2; %due to spikecount/subsampling setup
fulldel2time=[(ntce_ds(5,1)-46):(ntce_ds(5,1)-7) ntce_ds(5,1):ntce_ds(6,1)-20-2];

num1stshiftsL=zeros(248,1);
num2ndshiftsR=zeros(248,1);
num1stshiftsR=zeros(248,1);
num2ndshiftsL=zeros(248,1);


if recomp_shLRalgn
   if selflag<=2
   %%%%%%%%%%%%%%
    pooleposxmat_shRalgn=nan(size(pooleposxmat_ds_bln));
    pooleposxmat_shLalgn=nan(size(pooleposxmat_ds_bln));
    poolspksqmat_shRalgn=nan(size(poolspksqmat_ds_bln));
    poolspksqmat_shLalgn=nan(size(poolspksqmat_ds_bln));
    pooleposxmat_shRonly=nan(size(pooleposxmat_ds_bln));
    pooleposxmat_shLonly=nan(size(pooleposxmat_ds_bln));
    poolspksqmat_shRonly=nan(size(poolspksqmat_ds_bln));
    poolspksqmat_shLonly=nan(size(poolspksqmat_ds_bln));
    isdelSR=nan(size(poolspksqmat_ds_bln,1:2));
    isdelSL=nan(size(poolspksqmat_ds_bln,1:2));
    for cc=1:248
        % firstL offer1 - thenL delay2
        % Detect saccades
        ccepx=squeeze(pooleposxmat_ds_bln(:,cc,fulldel1time));
        ccspk=squeeze(poolspksqmat_ds_bln(:,cc,fulldel1time));
        [nuR,~,~,ishiftsR]=nanmonotonicdrifts_forward_simplified(ccepx(:,40:end),2,0,+1,1,stareflag);
        % delay2 shiftR trials have sum(xdriftsR,2)>0
        isdelSR(:,cc)=double(~isnan(ishiftsR));
        ind_delSR=find(~isnan(ishiftsR));
        % saccade R times
        %stR=arrayfun(@(xx) find(xdriftsR(xx,:)==1,1,'first'),isSR);
        stR=ishiftsR(ind_delSR);
        num1stshiftsR(cc)=sum(nuR>0);
        for jj=1:length(ind_delSR)
            shRbins=(ishiftsR(ind_delSR(jj)):length(fulldel1time));
            %removes post-shift returns
            [n2L,~,~,ishifts2L]=nanmonotonicdrifts_forward_simplified(ccepx(ind_delSR(jj),shRbins(41:end)),2,0,-1,+1,stareflag);
            if n2L>0
                shRbins((41+ishifts2L):end)=[];
                num2ndshiftsL(cc)=num2ndshiftsL(cc)+1;
            end
            pooleposxmat_shRalgn(ind_delSR(jj),cc,fulldel1time)=[ccepx(ind_delSR(jj),shRbins) nan(1,length(fulldel1time)-sum(shRbins<=length(fulldel1time)))];
            poolspksqmat_shRalgn(ind_delSR(jj),cc,fulldel1time)=[ccspk(ind_delSR(jj),shRbins) nan(1,length(fulldel1time)-sum(shRbins<=length(fulldel1time)))];
            pooleposxmat_shRonly(ind_delSR(jj),cc,fulldel1time(1)-1+shRbins)=ccepx(ind_delSR(jj),shRbins);
            poolspksqmat_shRonly(ind_delSR(jj),cc,fulldel1time(1)-1+shRbins)=ccepx(ind_delSR(jj),shRbins);
        end
        %figure(); plot(1:76,squeeze(pooleposxmat_shRalgn(:,1,fulldel1time)))
        %%
        [nuL,~,~,ishiftsL]=nanmonotonicdrifts_forward_simplified(ccepx(:,40:end),2,0,-1,1,stareflag);
        % delay2 shiftL trials
        isdelSL(:,cc)=double(~isnan(ishiftsL));
        ind_delSL=find(~isnan(ishiftsL));
        % saccade L times
        %stL=arrayfun(@(xx) find(xdriftsL(xx,:)==1,1,'first'),isSL);
        stL=ishiftsL(ind_delSL);
        num1stshiftsL(cc)=sum(nuL>0);
        for jj=1:length(ind_delSL)
            shLbins=(ishiftsL(ind_delSL(jj)):length(fulldel1time));
            %removes post-shift returns
            [n2R,~,~,ishifts2R]=nanmonotonicdrifts_forward_simplified(ccepx(ind_delSL(jj),shLbins(41:end)),2,0,+1,+1,stareflag);
            if n2R>0
                shLbins((41+ishifts2R):end)=[];
                num2ndshiftsR(cc)=num2ndshiftsR(cc)+1;
            end
            pooleposxmat_shLalgn(ind_delSL(jj),cc,fulldel1time)=[ccepx(ind_delSL(jj),shLbins) nan(1,length(fulldel1time)-sum(shLbins<=length(fulldel1time)))];
            poolspksqmat_shLalgn(ind_delSL(jj),cc,fulldel1time)=[ccspk(ind_delSL(jj),shLbins) nan(1,length(fulldel1time)-sum(shLbins<=length(fulldel1time)))];
            pooleposxmat_shLonly(ind_delSL(jj),cc,fulldel1time(1)-1+shLbins)=ccepx(ind_delSL(jj),shLbins);
            poolspksqmat_shLonly(ind_delSL(jj),cc,fulldel1time(1)-1+shLbins)=ccepx(ind_delSL(jj),shLbins);
        end
        %figure(); plot(1:76,squeeze(pooleposxmat_shLalgn(:,1,fulldel1time)))
    end
    
    %save('pool_del1_shLRalgn.mat','pooleposxmat_shRalgn','poolspksqmat_shRalgn','pooleposxmat_shLalgn','poolspksqmat_shLalgn','-v7.3')
       
   elseif selflag>2    
    %%%%%%%%%%%%%%%%%%%%%%%
    pooleposxmat_shRalgn=nan(size(pooleposxmat_ds_bln));
    pooleposxmat_shLalgn=nan(size(pooleposxmat_ds_bln));
    poolspksqmat_shRalgn=nan(size(poolspksqmat_ds_bln));
    poolspksqmat_shLalgn=nan(size(poolspksqmat_ds_bln));
    pooleposxmat_shRonly=nan(size(pooleposxmat_ds_bln));
    pooleposxmat_shLonly=nan(size(pooleposxmat_ds_bln));
    poolspksqmat_shRonly=nan(size(poolspksqmat_ds_bln));
    poolspksqmat_shLonly=nan(size(poolspksqmat_ds_bln));
    isdelSR=nan(size(poolspksqmat_ds_bln,1:2));
    isdelSL=nan(size(poolspksqmat_ds_bln,1:2));
    for cc=1:248
        % firstL offer1 - thenL delay2
        % Detect saccades
        ccepx=squeeze(pooleposxmat_ds_bln(:,cc,fulldel2time));
        ccspk=squeeze(poolspksqmat_ds_bln(:,cc,fulldel2time));
        [nuR,~,~,ishiftsR]=nanmonotonicdrifts_forward_simplified(ccepx(:,40:end),2,0,+1,1,stareflag);
        % delay2 shiftR trials have sum(xdriftsR,2)>0
        isdelSR(:,cc)=double(~isnan(ishiftsR));
        ind_delSR=find(~isnan(ishiftsR));
        % saccade R times
        %stR=arrayfun(@(xx) find(xdriftsR(xx,:)==1,1,'first'),isSR);
        stR=ishiftsR(ind_delSR);
        num1stshiftsR(cc)=sum(nuR>0);
        for jj=1:length(ind_delSR)
            shRbins=(ishiftsR(ind_delSR(jj)):length(fulldel2time));
            %removes post-shift returns
            [n2L,~,~,ishifts2L]=nanmonotonicdrifts_forward_simplified(ccepx(ind_delSR(jj),shRbins(41:end)),2,0,-1,+1,stareflag);
            if n2L>0
                shRbins((41+ishifts2L):end)=[];
                num2ndshiftsL(cc)=num2ndshiftsL(cc)+1;
            end
            pooleposxmat_shRalgn(ind_delSR(jj),cc,fulldel2time)=[ccepx(ind_delSR(jj),shRbins) nan(1,length(fulldel2time)-sum(shRbins<=length(fulldel2time)))];
            poolspksqmat_shRalgn(ind_delSR(jj),cc,fulldel2time)=[ccspk(ind_delSR(jj),shRbins) nan(1,length(fulldel2time)-sum(shRbins<=length(fulldel2time)))];
            pooleposxmat_shRonly(ind_delSR(jj),cc,fulldel2time(1)-1+shRbins)=ccepx(ind_delSR(jj),shRbins);
            poolspksqmat_shRonly(ind_delSR(jj),cc,fulldel2time(1)-1+shRbins)=ccepx(ind_delSR(jj),shRbins);
        end
        
        [nuL,~,~,ishiftsL]=nanmonotonicdrifts_forward_simplified(ccepx(:,40:end),2,0,-1,1,stareflag);
        % delay2 shiftL trials
        isdelSL(:,cc)=double(~isnan(ishiftsL));
        ind_delSL=find(~isnan(ishiftsL));
        % saccade L times
        %stL=arrayfun(@(xx) find(xdriftsL(xx,:)==1,1,'first'),isSL);
        stL=ishiftsL(ind_delSL);
        num1stshiftsL(cc)=sum(nuL>0);
        for jj=1:length(ind_delSL)
            shLbins=(ishiftsL(ind_delSL(jj)):length(fulldel2time));
            %removes post-shift returns
            [n2R,~,~,ishifts2R]=nanmonotonicdrifts_forward_simplified(ccepx(ind_delSL(jj),shLbins(41:end)),2,0,+1,+1,stareflag);
            if n2R>0
                shLbins((41+ishifts2R):end)=[];
                num2ndshiftsR(cc)=num2ndshiftsR(cc)+1;
            end
            pooleposxmat_shLalgn(ind_delSL(jj),cc,fulldel2time)=[ccepx(ind_delSL(jj),shLbins) nan(1,length(fulldel2time)-sum(shLbins<=length(fulldel2time)))];
            poolspksqmat_shLalgn(ind_delSL(jj),cc,fulldel2time)=[ccspk(ind_delSL(jj),shLbins) nan(1,length(fulldel2time)-sum(shLbins<=length(fulldel2time)))];
            pooleposxmat_shLonly(ind_delSL(jj),cc,fulldel2time(1)-1+shLbins)=ccepx(ind_delSL(jj),shLbins);
            poolspksqmat_shLonly(ind_delSL(jj),cc,fulldel2time(1)-1+shLbins)=ccepx(ind_delSL(jj),shLbins);
        end
        %disp(cc)
    end
    
    %save('pool_del2_shLRalgn.mat','pooleposxmat_shRalgn','poolspksqmat_shRalgn','pooleposxmat_shLalgn','poolspksqmat_shLalgn','-v7.3')
   end
   
else
    if selflag<=2
        load pool_del1_shLRalgn.mat
    elseif selflag>2
        load pool_del2_shLRalgn.mat
    end
end

sum(num1stshiftsL)/sum(ntrs.*nccs)
sum(num1stshiftsR)/sum(ntrs.*nccs)
sum(num2ndshiftsR)/sum(num1stshiftsL)
sum(num2ndshiftsL)/sum(num1stshiftsR)
sum(num2ndshiftsR)/sum(ntrs.*nccs)
sum(num2ndshiftsL)/sum(ntrs.*nccs)




%%
pooleposxmat_L=nan(Ntrs12,Nccs,Ntps_ds); pooleposxmat_R=nan(Ntrs12,Nccs,Ntps_ds);
if any(caseflag==[1 2 4 5])
    pooleposxmat_L=pooleposxmat_shLalgn(:,:,:);
    pooleposxmat_R=pooleposxmat_shRalgn(:,:,:);
else
    dtime=eval(selectedTimes);
    inddelstL=arrayfun(@(xx) all((rmnans(squeeze(pooleposxmat_ds_bln(xx,1,dtime))))<=0), 1:size(pooleposxmat_ds_bln,1));
    inddelstR=arrayfun(@(xx) all((rmnans(squeeze(pooleposxmat_ds_bln(xx,1,dtime))))>=0), 1:size(pooleposxmat_ds_bln,1));
    pooleposxmat_L(1:sum(inddelstL),:,dtime)=pooleposxmat_ds_bln(inddelstL,:,dtime);
    pooleposxmat_R(1:sum(inddelstR),:,dtime)=pooleposxmat_ds_bln(inddelstR,:,dtime);    
end
if selflag <=2 
    hfig=fullfig();
    subplot(1,2,1); title('Shifts to Left screen side');
    plot(rtmax_ds,(squeeze(pooleposxmat_L(:,1,vtps_ds==1)))./32558); ylim([-1 1]); plothline(0);
    xlabel('time after delay1 aligned to gaze shift initiation (ms)'); ylabel('horizontal eye position (normalized)'); box off
    subplot(1,2,2); title('Shifts to Right screen side');
    plot(rtmax_ds,(squeeze(pooleposxmat_R(:,1,vtps_ds==1)))./32558); ylim([-1 1]); plothline(0);
    xlabel('time after delay1 aligned to gaze shift initiation (ms)'); ylabel('horizontal eye position (normalized)'); box off
    saveas(hfig,[most_recent_folder 'span' num2str(Nspan) 'ms_lag' num2str(Nlag) 'ms/sample_trials_off_del_shifts_monitor.png']);
    saveas(hfig,[most_recent_folder 'span' num2str(Nspan) 'ms_lag' num2str(Nlag) 'ms/sample_trials_off_del_shifts_monitor.fig']);
    saveas(hfig,[most_recent_folder 'span' num2str(Nspan) 'ms_lag' num2str(Nlag) 'ms/sample_trials_off_del_shifts_monitor.svg']);
    
    fullfig();
    subplot(1,2,1);
    plotmsem(fulldel1time,nanmedian(squeeze(pooleposxmat_L(:,1,fulldel1time))),'k');
    plot(fulldel1time,arrayfun(@(xx) prctile(rmnans(squeeze(pooleposxmat_L(:,1,xx))),99), fulldel1time),'color',[1 1 1].*.85);
    plot(fulldel1time,arrayfun(@(xx) prctile(rmnans(squeeze(pooleposxmat_L(:,1,xx))),01), fulldel1time),'color',[1 1 1].*.85);
    plot(fulldel1time,arrayfun(@(xx) prctile(rmnans(squeeze(pooleposxmat_L(:,1,xx))),95), fulldel1time),'color',[1 1 1].*.60);
    plot(fulldel1time,arrayfun(@(xx) prctile(rmnans(squeeze(pooleposxmat_L(:,1,xx))),05), fulldel1time),'color',[1 1 1].*.60);
    plot(fulldel1time,arrayfun(@(xx) prctile(rmnans(squeeze(pooleposxmat_L(:,1,xx))),90), fulldel1time),'color',[1 1 1].*.45);
    plot(fulldel1time,arrayfun(@(xx) prctile(rmnans(squeeze(pooleposxmat_L(:,1,xx))),10), fulldel1time),'color',[1 1 1].*.45);
    plot(fulldel1time,arrayfun(@(xx) prctile(rmnans(squeeze(pooleposxmat_L(:,1,xx))),75), fulldel1time),'color',[1 1 1].*.25);
    plot(fulldel1time,arrayfun(@(xx) prctile(rmnans(squeeze(pooleposxmat_L(:,1,xx))),25), fulldel1time),'color',[1 1 1].*.25);
    plot(fulldel1time,arrayfun(@(xx) prctile(rmnans(squeeze(pooleposxmat_L(:,1,xx))),65), fulldel1time),'color',[0 0 1].*.1);
    plot(fulldel1time,arrayfun(@(xx) prctile(rmnans(squeeze(pooleposxmat_L(:,1,xx))),35), fulldel1time),'color',[1 1 1].*.1); %ylim([-1 1]);
    plotvline([min(eval(selectedTimes)) max(eval(selectedTimes))]); plothline(0); xlim([min(eval(selectedTimes)) max(eval(selectedTimes))]);
    
    subplot(1,2,2);
    plotmsem(fulldel1time,nanmedian(squeeze(pooleposxmat_R(:,1,fulldel1time))),'k');
    plot(fulldel1time,arrayfun(@(xx) prctile(rmnans(squeeze(pooleposxmat_R(:,1,xx))),99), fulldel1time),'color',[1 1 1].*.85);
    plot(fulldel1time,arrayfun(@(xx) prctile(rmnans(squeeze(pooleposxmat_R(:,1,xx))),01), fulldel1time),'color',[1 1 1].*.85);
    plot(fulldel1time,arrayfun(@(xx) prctile(rmnans(squeeze(pooleposxmat_R(:,1,xx))),95), fulldel1time),'color',[1 1 1].*.60);
    plot(fulldel1time,arrayfun(@(xx) prctile(rmnans(squeeze(pooleposxmat_R(:,1,xx))),05), fulldel1time),'color',[1 1 1].*.60);
    plot(fulldel1time,arrayfun(@(xx) prctile(rmnans(squeeze(pooleposxmat_R(:,1,xx))),90), fulldel1time),'color',[1 1 1].*.45);
    plot(fulldel1time,arrayfun(@(xx) prctile(rmnans(squeeze(pooleposxmat_R(:,1,xx))),10), fulldel1time),'color',[1 1 1].*.45);
    plot(fulldel1time,arrayfun(@(xx) prctile(rmnans(squeeze(pooleposxmat_R(:,1,xx))),75), fulldel1time),'color',[1 1 1].*.25);
    plot(fulldel1time,arrayfun(@(xx) prctile(rmnans(squeeze(pooleposxmat_R(:,1,xx))),25), fulldel1time),'color',[1 1 1].*.25);
    plot(fulldel1time,arrayfun(@(xx) prctile(rmnans(squeeze(pooleposxmat_R(:,1,xx))),65), fulldel1time),'color',[0 0 1].*.1);
    plot(fulldel1time,arrayfun(@(xx) prctile(rmnans(squeeze(pooleposxmat_R(:,1,xx))),35), fulldel1time),'color',[1 1 1].*.1); %ylim([-1 1]);
    plotvline([min(eval(selectedTimes)) max(eval(selectedTimes))]); plothline(0); xlim([min(eval(selectedTimes)) max(eval(selectedTimes))]);
else
    hfig=fullfig();
    subplot(1,2,1); title('Shifts to Left screen side');
    plot(rtmax_ds,(squeeze(pooleposxmat_L(:,1,vtps_ds==1)))./32558); ylim([-1 1]); plothline(0);
    xlabel('time after delay2 aligned to gaze shift initiation (ms)'); ylabel('horizontal eye position (normalized)'); box off
    subplot(1,2,2); title('Shifts to Right screen side');
    plot(rtmax_ds,(squeeze(pooleposxmat_R(:,1,vtps_ds==1)))./32558); ylim([-1 1]); plothline(0);
    xlabel('time after delay2 aligned to gaze shift initiation (ms)'); ylabel('horizontal eye position (normalized)'); box off
    saveas(hfig,[most_recent_folder 'span' num2str(Nspan) 'ms_lag' num2str(Nlag) 'ms/sample_trials_off_del_shifts_monitor.png']);
    saveas(hfig,[most_recent_folder 'span' num2str(Nspan) 'ms_lag' num2str(Nlag) 'ms/sample_trials_off_del_shifts_monitor.fig']);
    saveas(hfig,[most_recent_folder 'span' num2str(Nspan) 'ms_lag' num2str(Nlag) 'ms/sample_trials_off_del_shifts_monitor.svg']);
    %%
    fullfig();
    subplot(1,2,1);
    plotmsem(fulldel2time,nanmedian(squeeze(pooleposxmat_L(:,1,fulldel2time))),'k');
    plot(fulldel2time,arrayfun(@(xx) prctile(rmnans(squeeze(pooleposxmat_L(:,1,xx))),99), fulldel2time),'color',[1 1 1].*.85);
    plot(fulldel2time,arrayfun(@(xx) prctile(rmnans(squeeze(pooleposxmat_L(:,1,xx))),01), fulldel2time),'color',[1 1 1].*.85);
    plot(fulldel2time,arrayfun(@(xx) prctile(rmnans(squeeze(pooleposxmat_L(:,1,xx))),95), fulldel2time),'color',[1 1 1].*.60);
    plot(fulldel2time,arrayfun(@(xx) prctile(rmnans(squeeze(pooleposxmat_L(:,1,xx))),05), fulldel2time),'color',[1 1 1].*.60);
    plot(fulldel2time,arrayfun(@(xx) prctile(rmnans(squeeze(pooleposxmat_L(:,1,xx))),90), fulldel2time),'color',[1 1 1].*.45);
    plot(fulldel2time,arrayfun(@(xx) prctile(rmnans(squeeze(pooleposxmat_L(:,1,xx))),10), fulldel2time),'color',[1 1 1].*.45);
    plot(fulldel2time,arrayfun(@(xx) prctile(rmnans(squeeze(pooleposxmat_L(:,1,xx))),75), fulldel2time),'color',[1 1 1].*.25);
    plot(fulldel2time,arrayfun(@(xx) prctile(rmnans(squeeze(pooleposxmat_L(:,1,xx))),25), fulldel2time),'color',[1 1 1].*.25);
    plot(fulldel2time,arrayfun(@(xx) prctile(rmnans(squeeze(pooleposxmat_L(:,1,xx))),65), fulldel2time),'color',[0 0 1].*.1);
    plot(fulldel2time,arrayfun(@(xx) prctile(rmnans(squeeze(pooleposxmat_L(:,1,xx))),35), fulldel2time),'color',[1 1 1].*.1);%ylim([-1 1]);
    plotvline([min(eval(selectedTimes)) max(eval(selectedTimes))]); plothline(0); xlim([min(eval(selectedTimes)) max(eval(selectedTimes))]);
    
    subplot(1,2,2);
    plotmsem(fulldel2time,nanmedian(squeeze(pooleposxmat_R(:,1,fulldel2time))),'k');
    plot(fulldel2time,arrayfun(@(xx) prctile(rmnans(squeeze(pooleposxmat_R(:,1,xx))),99), fulldel2time),'color',[1 1 1].*.85);
    plot(fulldel2time,arrayfun(@(xx) prctile(rmnans(squeeze(pooleposxmat_R(:,1,xx))),01), fulldel2time),'color',[1 1 1].*.85);
    plot(fulldel2time,arrayfun(@(xx) prctile(rmnans(squeeze(pooleposxmat_R(:,1,xx))),95), fulldel2time),'color',[1 1 1].*.60);
    plot(fulldel2time,arrayfun(@(xx) prctile(rmnans(squeeze(pooleposxmat_R(:,1,xx))),05), fulldel2time),'color',[1 1 1].*.60);
    plot(fulldel2time,arrayfun(@(xx) prctile(rmnans(squeeze(pooleposxmat_R(:,1,xx))),90), fulldel2time),'color',[1 1 1].*.45);
    plot(fulldel2time,arrayfun(@(xx) prctile(rmnans(squeeze(pooleposxmat_R(:,1,xx))),10), fulldel2time),'color',[1 1 1].*.45);
    plot(fulldel2time,arrayfun(@(xx) prctile(rmnans(squeeze(pooleposxmat_R(:,1,xx))),75), fulldel2time),'color',[1 1 1].*.25);
    plot(fulldel2time,arrayfun(@(xx) prctile(rmnans(squeeze(pooleposxmat_R(:,1,xx))),25), fulldel2time),'color',[1 1 1].*.25);
    plot(fulldel2time,arrayfun(@(xx) prctile(rmnans(squeeze(pooleposxmat_R(:,1,xx))),65), fulldel2time),'color',[0 0 1].*.1);
    plot(fulldel2time,arrayfun(@(xx) prctile(rmnans(squeeze(pooleposxmat_R(:,1,xx))),35), fulldel2time),'color',[1 1 1].*.1);%ylim([-1 1]);
    plotvline([min(eval(selectedTimes)) max(eval(selectedTimes))]); plothline(0); xlim([min(eval(selectedTimes)) max(eval(selectedTimes))]);
end

%%

% fullfig();
% subplot(1,4,1);
% spy(vectorisen(~isnan(poolspksqmat_shLalgn(:,cumsum(nccs),:)),1:2));
% subplot(1,4,2);
% spy(vectorisen(~isnan(poolspksqmat_shRalgn(:,cumsum(nccs),:)),1:2));
% subplot(1,4,3);
% spy(vectorisen(~isnan(pooleposxmat_shLalgn(:,cumsum(nccs),:)),1:2));
% subplot(1,4,4);
% spy(vectorisen(~isnan(pooleposxmat_shRalgn(:,cumsum(nccs),:)),1:2));
% title('Aligned to shifts');
%%
hfig=fullfig(); hfig.Position=[0 0 .5 1];
subplot(2,1,1); ylim([0 100]); hold on
plot(rtmax_ds,squeeze(sum(~isnan(pooleposxmat_L(:,cumsum(nccs),vtps_ds==1))))./repmat(ntrs',1,sum(vtps_ds>0))*100); hold on
for jj=1:8; plot([1 1].*sum(rtmes_ds(1:jj)),get(gca,'YLim'),'k:'); text(sum(rtmes_ds(1:jj-1)),max(get(gca,'YLim')),twLabels{jj}); end
subplot(2,1,2); ylim([0 100]); hold on
plot(rtmax_ds,squeeze(sum(~isnan(pooleposxmat_R(:,cumsum(nccs),vtps_ds==1))))./repmat(ntrs',1,sum(vtps_ds>0))*100); hold on
for jj=1:8; plot([1 1].*sum(rtmes_ds(1:jj)),get(gca,'YLim'),'k:'); text(sum(rtmes_ds(1:jj-1)),max(get(gca,'YLim')),twLabels{jj}); end
legend({'subj 1, sn 1','subj 1, sn 2','subj 1, sn 3','subj 1, sn 4','subj 2, sn 1','subj 2, sn 2','subj 2, sn 3','subj 2, sn 4'},'box','off');

switch selflag
    case 1
        subplot(2,1,1);
        title('Look Left Offer1 & Look Left Delay1'); ylim([0 100]); ylabel('% trials'); xlabel('Delay1 time bins (10 ms)');
        subplot(2,1,2);
        title('Look Left Offer1 & Look Right Delay1'); ylim([0 100]); ylabel('% trials'); xlabel('Delay1 time bins (10 ms)');
    case 2
        subplot(2,1,1);
        title('Look Right Offer1 & Look Left Delay1'); ylim([0 100]); ylabel('% trials'); xlabel('Delay1 time bins (10 ms)');
        subplot(2,1,2);
        title('Look Right Offer1 & Look Right Delay1'); ylim([0 100]); ylabel('% trials'); xlabel('Delay1 time bins (10 ms)');
    case 3
        subplot(2,1,1);
        title('Look Left Offer2 & Look Left Delay2'); ylim([0 100]); ylabel('% trials'); xlabel('Delay2 time bins (10 ms)');
        subplot(2,1,2);
        title('Look Left Offer2 & Look Right Delay2'); ylim([0 100]); ylabel('% trials'); xlabel('Delay2 time bins (10 ms)');
    case 4
        subplot(2,1,1);
        title('Look Right Offer2 & Look Left Delay2'); ylim([0 100]); ylabel('% trials'); xlabel('Delay2 time bins (10 ms)');
        subplot(2,1,2);
        title('Look Right Offer2 & Look Right Delay2'); ylim([0 100]); ylabel('% trials'); xlabel('Delay2 time bins (10 ms)');
end
saveas(hfig,[most_recent_folder 'span' num2str(Nspan) 'ms_lag' num2str(Nlag) 'ms/frac_trials_off_del.png']);
saveas(hfig,[most_recent_folder 'span' num2str(Nspan) 'ms_lag' num2str(Nlag) 'ms/frac_trials_off_del.fig']);
saveas(hfig,[most_recent_folder 'span' num2str(Nspan) 'ms_lag' num2str(Nlag) 'ms/frac_trials_off_del.svg']);

clearvars pooleposxmat_L pooleposxmat_R

%%
if any(caseflag==[1 2 4 5])

hfig=fullfig(); hfig.Position=[0 0 .5 1];
subplot(2,1,1); ylim([0 100]); hold on; gg=gca; 
%plot(1:Ntps_ds,squeeze(sum(~isnan(poolspksqmat_shLonly(:,cumsum(nccs),:))))./repmat(ntrs',1,Ntps_ds)*100); hold on
pp=plot(eval(selectedTimes),squeeze(sum(~isnan(poolspksqmat_shLonly(:,cumsum(nccs),eval(selectedTimes)))))./repmat(ntrs',1,length(eval(selectedTimes)))*100); hold on
for jj=1:8
plot(1:Ntps_ds,squeeze(sum(~isnan(poolspksqmat_shLonly(:,sum(nccs(1:jj)),:))))./repmat(ntrs(jj)',1,Ntps_ds)*100,':','color',pp(jj).Color); hold on
end
for jj=1:7; plotvline(ntce_ds(jj,1:2)); text(ntce_ds(jj,1),max(get(gca,'YLim')),twLabels{jj}); end
%xlim([189 438]);%xlim([1 Ntps_ds]);
subplot(2,1,2); ylim([0 100]); hold on
%plot(1:Ntps_ds,squeeze(sum(~isnan(poolspksqmat_shRonly(:,cumsum(nccs),:))))./repmat(ntrs',1,Ntps_ds)*100); hold on
pp=plot(eval(selectedTimes),squeeze(sum(~isnan(poolspksqmat_shRonly(:,cumsum(nccs),eval(selectedTimes)))))./repmat(ntrs',1,length(eval(selectedTimes)))*100); hold on
for jj=1:8
plot(1:Ntps_ds,squeeze(sum(~isnan(poolspksqmat_shRonly(:,sum(nccs(1:jj)),:))))./repmat(ntrs(jj)',1,Ntps_ds)*100,':','color',pp(jj).Color); hold on
end
    
for jj=1:7; plotvline(ntce_ds(jj,1:2)); tt=text(ntce_ds(jj,1),.85*max(get(gca,'YLim')),twLabels{jj},'Rotation',90); end
legend({'subj 1, sn 1','subj 1, sn 2','subj 1, sn 3','subj 1, sn 4','subj 2, sn 1','subj 2, sn 2','subj 2, sn 3','subj 2, sn 4'},'box','off');
%xlim([189 438]);%xlim([1 Ntps_ds]);
switch selflag
    case 1
        subplot(2,1,1);
        title('Look Left Offer1 & Look Left Delay1'); ylim([0 100]); ylabel('% trials'); xlabel('Delay1 time bins (10 ms)');
        subplot(2,1,2);
        title('Look Left Offer1 & Look Right Delay1'); ylim([0 100]); ylabel('% trials'); xlabel('Delay1 time bins (10 ms)');
    case 2
        subplot(2,1,1);
        title('Look Right Offer1 & Look Left Delay1'); ylim([0 100]); ylabel('% trials'); xlabel('Delay1 time bins (10 ms)');
        subplot(2,1,2);
        title('Look Right Offer1 & Look Right Delay1'); ylim([0 100]); ylabel('% trials'); xlabel('Delay1 time bins (10 ms)');
    case 3
        subplot(2,1,1);
        title('Look Left Offer2 & Look Left Delay2'); ylim([0 100]); ylabel('% trials'); xlabel('Delay2 time bins (10 ms)');
        subplot(2,1,2);
        title('Look Left Offer2 & Look Right Delay2'); ylim([0 100]); ylabel('% trials'); xlabel('Delay2 time bins (10 ms)');
    case 4
        subplot(2,1,1);
        title('Look Right Offer2 & Look Left Delay2'); ylim([0 100]); ylabel('% trials'); xlabel('Delay2 time bins (10 ms)');
        subplot(2,1,2);
        title('Look Right Offer2 & Look Right Delay2'); ylim([0 100]); ylabel('% trials'); xlabel('Delay2 time bins (10 ms)');
end
saveas(hfig,[most_recent_folder 'span' num2str(Nspan) 'ms_lag' num2str(Nlag) 'ms/frac_trials_off_del_shifts.png']);
saveas(hfig,[most_recent_folder 'span' num2str(Nspan) 'ms_lag' num2str(Nlag) 'ms/frac_trials_off_del_shifts.fig']);
saveas(hfig,[most_recent_folder 'span' num2str(Nspan) 'ms_lag' num2str(Nlag) 'ms/frac_trials_off_del_shifts.svg']);

end

% subplot(1,4,2);
% spy(vectorisen(~isnan(poolspksqmat_shRalgn(:,cumsum(nccs),:)),1:2));
% subplot(1,4,3);
% spy(vectorisen(~isnan(pooleposxmat_shLalgn(:,cumsum(nccs),:)),1:2));
% subplot(1,4,4);
% spy(vectorisen(~isnan(pooleposxmat_shRalgn(:,cumsum(nccs),:)),1:2));
% title('Aligned to shifts');


%  SET ALIGNMENT
%  figure('units','normalized','position',[0 0 1 1]);
%  subplot(2,1,1); shtr=1:length(isSR); plot(rerange(ccpex(isSR(shtr),:)',ccpex(:))); hold on
%  plot(stR(shtr),arrayfun(@(xx) rerange(ccpex(isSR(xx),stR(xx)),ccpex(:)), shtr),'rx');
%  plot(get(gca,'XLim'),[0 0],'k:'); ylim([-1 1]); title(sprintf('sample cell, shift R detection, n_s_h_R=%d', sum(nuR)));
%  
%  subplot(2,1,2); shtr=1:length(isSL); plot(rerange(ccpex(isSL(shtr),:)',ccpex(:))); hold on
%  plot(stL(shtr),arrayfun(@(xx) rerange(ccpex(isSL(xx),stL(xx)),ccpex(:)), shtr),'rx');
%  plot(get(gca,'XLim'),[0 0],'k:'); ylim([-1 1]); title(sprintf('sample cell, shift L detection, n_s_h_L=%d', sum(nuL)));
% 
%  %
%  figure('units','normalized','position',[0 0 1 1]); 
%  subplot(2,1,1); 
%  plot(vectorisen(pooleposxmat_shRalgn(:,end,:),1:2)'./max(abs(pooleposxmat_shRalgn(:)))); hold on;
%  plot(get(gca,'XLim'),[0 0],'k:'); ylim([-1 1]); title(sprintf('sample cell, shift R detection, n_s_h_R=%d', sum(nuR)));
%  
%  subplot(2,1,2); 
%  plot(vectorisen(pooleposxmat_shLalgn(:,end,:),1:2)'./max(abs(pooleposxmat_shLalgn(:)))); hold on;
%  plot(get(gca,'XLim'),[0 0],'k:'); ylim([-1 1]); title(sprintf('sample cell, shift L detection, n_s_h_L=%d', sum(nuL))); 
 
% firstL offer1 - thenR delay2
% firstR offer1 - thenL delay2
% firstR offer1 - thenR delay2

%%
if 0
    % re-defined on vtps_ds
    o1LLinterval=20; %x10ms=200ms
    d1LRinterval=20; %x10ms=200ms
    o2LRinterval=20; %x10ms=200ms
    d2LLinterval=20; %x10ms=200ms
    chLLinterval=44;
    Offer1LLtime=sum(rtmes_ds(1:2))-o1LLinterval+(1:o1LLinterval);
    Delay1LRtime=sum(rtmes_ds(1:3))-d1LRinterval+(1:d1LRinterval);
    Offer2LRtime=sum(rtmes_ds(1:4))-o2LRinterval+(1:o2LRinterval);
    Delay2LLtime=sum(rtmes_ds(1:5))-d2LLinterval+(1:d2LLinterval);
    ChoiceLLtime=sum(rtmes_ds(1:8))-chLLinterval+(1:chLLinterval);

    
    pooleposxmat_vds=pooleposxmat_ds(:,:,vtps_ds==1);
    pooleposxmat_vrds=rmnans(vectorisen(pooleposxmat_ds(:,cumsum(nccs),vtps_ds==1),1:2));
    xvro1LL=vectorise(nanmean(pooleposxmat_vrds(:,Offer1LLtime),2)<0);
    xvrd1LR=vectorise(nanmean(pooleposxmat_vrds(:,Delay1LRtime),2)>0);
    xvro2LR=vectorise(nanmean(pooleposxmat_vrds(:,Offer2LRtime),2)>0);
    xvrd2LL=vectorise(nanmean(pooleposxmat_vrds(:,Delay2LLtime),2)<0);
    xvrchLL=vectorise(nanmean(pooleposxmat_vrds(:,ChoiceLLtime),2)<0);

    pooleposxmat_vrdds=pooleposxmat_vrds;
    pooleposxmat_vrdds(rmnans(vectorise(pooloffer1sd(:,cumsum(nccs))))==0,:)=...
        -pooleposxmat_vrdds(rmnans(vectorise(pooloffer1sd(:,cumsum(nccs))))==0,:);
    
    xvrdo1LL=vectorise(nanmean(pooleposxmat_vrdds(:,Offer1LLtime),2)<0);
    xvrdd1LR=vectorise(nanmean(pooleposxmat_vrdds(:,Delay1LRtime),2)>0);
    xvrdo2LR=vectorise(nanmean(pooleposxmat_vrdds(:,Offer2LRtime),2)>0);
    xvrdd2LL=vectorise(nanmean(pooleposxmat_vrdds(:,Delay2LLtime),2)<0);
    xvrdchLL=vectorise(nanmean(pooleposxmat_vrdds(:,ChoiceLLtime),2)<0);
    
    xo1LL=vectorise(nanmean(pooleposxmat_vds(:,:,Offer1LLtime),3)<0);
    xo1LR=vectorise(nanmean(pooleposxmat_vds(:,:,Offer1LLtime),3)>0);
    xd1LR=vectorise(nanmean(pooleposxmat_vds(:,:,Delay1LRtime),3)>0);
    xo2LR=vectorise(nanmean(pooleposxmat_vds(:,:,Offer2LRtime),3)>0);
    xd2LL=vectorise(nanmean(pooleposxmat_vds(:,:,Delay2LLtime),3)<0);
    vvss=vectorisen(pooleposxmat_vds,1:2);
    vvrss=(pooleposxmat_vrds);
    iso1LL=nanmean(pooleposxmat_vds(:,:,Offer1LLtime),3)<0;
    isd1LR=nanmean(pooleposxmat_vds(:,:,Delay1LRtime),3)>0;
    iso2LR=nanmean(pooleposxmat_vds(:,:,Offer2LRtime),3)>0;
    isd2LL=nanmean(pooleposxmat_vds(:,:,Delay2LLtime),3)<0;
    iso1LLd2LL=  iso1LL &  isd2LL;
    iso1LLd1LR=  iso1LL &(~isd2LL);
    iso1LLd1LLd2LL=  iso1LL & (~isd1LR) &  isd2LL;
    iso1LLd1LLd2LR=  iso1LL & (~isd1LR) &(~isd2LL);
    
    iso2LRd2LL=  iso2LR &  isd2LL;
    iso2LRd2LR=  iso2LR &(~isd2LL);
    iso2LLd2LL=(~iso2LR)&  isd2LL;
    iso2LLd2LR=(~iso2LR)&(~isd2LL);
    xebins=linspace(-max(abs(vvss(:))),max(abs(vvss(:))),1e3);
    
    
    firstL  =logical(rmnans(vectorise(pooloffer1sd(:,cumsum(nccs))))==1);
    firstR  =logical(rmnans(vectorise(pooloffer1sd(:,cumsum(nccs))))==0);
    chooseR =logical(rmnans(vectorise(poolofferchR(:,cumsum(nccs))))==1);
    
    %%
    % offer1sd==1 if L; 0 if R.
    clc
    fprintf('LookL on offer1 %2.2f%% of the trials where first offer is Left.\n',   mean([double(xvrdo1LL(firstL))==1 ])*100);
    fprintf('LookR on offer1 %2.2f%% of the trials where first offer is Left.\n',   mean([double(xvrdo1LL(firstL))==0 ])*100);
    fprintf('LookL on offer1 %2.2f%% of the trials where first offer is Right.\n',  mean([double(xvrdo1LL(firstR))==1 ])*100);
    fprintf('LookR on offer1 %2.2f%% of the trials where first offer is Right.\n\n',mean([double(xvrdo1LL(firstR))==0 ])*100);
    
    fprintf('LookL on delay1 %2.2f%% of the trials where first offer is Left.\n',   mean([double(xvrdd1LR(firstL))==0 ])*100);
    fprintf('LookR on delay1 %2.2f%% of the trials where first offer is Left.\n',   mean([double(xvrdd1LR(firstL))==1 ])*100);
    fprintf('LookL on delay1 %2.2f%% of the trials where first offer is Right.\n',  mean([double(xvrdd1LR(firstR))==0 ])*100);
    fprintf('LookR on delay1 %2.2f%% of the trials where first offer is Right.\n\n',mean([double(xvrdd1LR(firstR))==1 ])*100);
    
    fprintf('LookL on offer2 %2.2f%% of the trials where first offer is Left.\n',   mean([double(xvrdo2LR(firstL))==0 ])*100);
    fprintf('LookR on offer2 %2.2f%% of the trials where first offer is Left.\n',   mean([double(xvrdo2LR(firstL))==1 ])*100);
    fprintf('LookL on offer2 %2.2f%% of the trials where first offer is Right.\n',  mean([double(xvrdo2LR(firstR))==0 ])*100);
    fprintf('LookR on offer2 %2.2f%% of the trials where first offer is Right.\n\n',mean([double(xvrdo2LR(firstR))==1 ])*100);
    
    fprintf('LookL on delay2 %2.2f%% of the trials where first offer is Left.\n',   mean([double(xvrdd2LL(firstL))==1 ])*100);
    fprintf('LookR on delay2 %2.2f%% of the trials where first offer is Left.\n',   mean([double(xvrdd2LL(firstL))==0 ])*100);
    fprintf('LookL on delay2 %2.2f%% of the trials where first offer is Right.\n',  mean([double(xvrdd2LL(firstR))==1 ])*100);
    fprintf('LookR on delay2 %2.2f%% of the trials where first offer is Right.\n\n',mean([double(xvrdd2LL(firstR))==0 ])*100);
    
    fprintf('LookL on choice %2.2f%% of the trials where first offer is Left.\n',   mean([double(xvrdchLL(firstL))==1 ])*100);
    fprintf('LookR on choice %2.2f%% of the trials where first offer is Left.\n',   mean([double(xvrdchLL(firstL))==0 ])*100);
    fprintf('LookL on choice %2.2f%% of the trials where first offer is Right.\n',  mean([double(xvrdchLL(firstR))==1 ])*100);
    fprintf('LookR on choice %2.2f%% of the trials where first offer is Right.\n\n',mean([double(xvrdchLL(firstR))==0 ])*100);
    
    %%
    Fisbest =logical(rmnans(vectorise(pooloffer1ev(:,cumsum(nccs)))) >rmnans(vectorise(pooloffer2ev(:,cumsum(nccs)))));
    Sisbest =logical(rmnans(vectorise(pooloffer1ev(:,cumsum(nccs)))) <rmnans(vectorise(pooloffer2ev(:,cumsum(nccs)))));
    FSssame =logical(rmnans(vectorise(pooloffer1ev(:,cumsum(nccs))))==rmnans(vectorise(pooloffer2ev(:,cumsum(nccs)))));
    
    Lisbest_reref= logical([Fisbest & firstL;  Fisbest & firstR]);
    Risbest_reref= logical([Sisbest & firstL;  Sisbest & firstR]);
    chR_reref=logical([chooseR==1 & firstL; chooseR==0 & firstR]);
    
    clc
    % LookL and Lis best means LookL(firstL&Fisbest) + LookR(firstR&Risbest) 
    % LL fL fLbest - LL fR fRbest - LL fL SRbest - LL fR SLbest  
    % LL fL fLbest - LL fL fLbest - LR fL fRbest - LL fL 
    fprintf('LookL on offer1 %2.2f%% of the trials where first, left offer is best.\n',  mean([double(xvrdo1LL(Fisbest & firstL)==1)  double(xvrdo1LL(Fisbest & firstR)==0)])*100);
    fprintf('LookL on offer1 %2.2f%% of the trials where first, left offer is worse.\n', mean([double(xvrdo1LL(Sisbest & firstL)==1)  double(xvrdo1LL(Sisbest & firstR)==0)])*100);
    fprintf('LookL on offer1 %2.2f%% of the trials where first, left offer is same.\n',  mean([double(xvrdo1LL(firstL  & FSssame)==1) double(xvrdo1LL(firstR  & FSssame)==0)])*100);
    fprintf('LookR on offer1 %2.2f%% of the trials where first, left offer is best.\n',  mean([double(xvrdo1LL(Fisbest & firstL)==0)  double(xvrdo1LL(Fisbest & firstR)==1)])*100);
    fprintf('LookR on offer1 %2.2f%% of the trials where first, left offer is worse.\n', mean([double(xvrdo1LL(Sisbest & firstL)==0)  double(xvrdo1LL(Sisbest & firstR)==1)])*100);
    fprintf('LookR on offer1 %2.2f%% of the trials where first, left offer is same.\n\n',mean([double(xvrdo1LL(firstL  & FSssame)==0) double(xvrdo1LL(firstR  & FSssame)==1)])*100);

    fprintf('LookL on delay1 %2.2f%% of the trials where first, left offer is best.\n',  mean([double(xvrdd1LR(Fisbest & firstL)==0)  double(xvrdd1LR(Fisbest & firstR)==1)])*100);
    fprintf('LookL on delay1 %2.2f%% of the trials where first, left offer is worse.\n', mean([double(xvrdd1LR(Sisbest & firstL)==0)  double(xvrdd1LR(Sisbest & firstR)==1)])*100);
    fprintf('LookL on delay1 %2.2f%% of the trials where first, left offer is same.\n',  mean([double(xvrdd1LR(firstL  & FSssame)==0) double(xvrdd1LR(firstR  & FSssame)==1)])*100);
    fprintf('LookR on delay1 %2.2f%% of the trials where first, left offer is best.\n',  mean([double(xvrdd1LR(Fisbest & firstL)==1)  double(xvrdd1LR(Fisbest & firstR)==0)])*100);
    fprintf('LookR on delay1 %2.2f%% of the trials where first, left offer is worse.\n', mean([double(xvrdd1LR(Sisbest & firstL)==1)  double(xvrdd1LR(Sisbest & firstR)==0)])*100);
    fprintf('LookR on delay1 %2.2f%% of the trials where first, left offer is same.\n\n',mean([double(xvrdd1LR(firstL  & FSssame)==1) double(xvrdd1LR(firstR  & FSssame)==0)])*100);
    
    fprintf('LookL on offer2 %2.2f%% of the trials where first, left offer is best.\n',  mean([double(xvrdo2LR(Fisbest & firstL)==0)  double(xvrdo2LR(Fisbest & firstR)==1)])*100);
    fprintf('LookL on offer2 %2.2f%% of the trials where first, left offer is worse.\n', mean([double(xvrdo2LR(Sisbest & firstL)==0)  double(xvrdo2LR(Sisbest & firstR)==1)])*100);
    fprintf('LookL on offer2 %2.2f%% of the trials where first, left offer is same.\n',  mean([double(xvrdo2LR(firstL  & FSssame)==0) double(xvrdo2LR(firstR  & FSssame)==1)])*100);
    fprintf('LookR on offer2 %2.2f%% of the trials where first, left offer is best.\n',  mean([double(xvrdo2LR(Fisbest & firstL)==1)  double(xvrdo2LR(Fisbest & firstR)==0)])*100);
    fprintf('LookR on offer2 %2.2f%% of the trials where first, left offer is worse.\n', mean([double(xvrdo2LR(Sisbest & firstL)==1)  double(xvrdo2LR(Sisbest & firstR)==0)])*100);
    fprintf('LookR on offer2 %2.2f%% of the trials where first, left offer is same.\n\n',mean([double(xvrdo2LR(firstL  & FSssame)==1) double(xvrdo2LR(firstR  & FSssame)==0)])*100);
    
    fprintf('LookL on delay2 %2.2f%% of the trials where first, left offer is best.\n',  mean([double(xvrdd2LL(Fisbest & firstL)==1)  double(xvrdd2LL(Fisbest & firstR)==0)])*100);
    fprintf('LookL on delay2 %2.2f%% of the trials where first, left offer is worse.\n', mean([double(xvrdd2LL(Sisbest & firstL)==1)  double(xvrdd2LL(Sisbest & firstR)==0)])*100);
    fprintf('LookL on delay2 %2.2f%% of the trials where first, left offer is same.\n',  mean([double(xvrdd2LL(firstL  & FSssame)==1) double(xvrdd2LL(firstR  & FSssame)==0)])*100);
    fprintf('LookR on delay2 %2.2f%% of the trials where first, left offer is best.\n',  mean([double(xvrdd2LL(Fisbest & firstL)==0)  double(xvrdd2LL(Fisbest & firstR)==1)])*100);
    fprintf('LookR on delay2 %2.2f%% of the trials where first, left offer is worse.\n', mean([double(xvrdd2LL(Sisbest & firstL)==0)  double(xvrdd2LL(Sisbest & firstR)==1)])*100);
    fprintf('LookR on delay2 %2.2f%% of the trials where first, left offer is same.\n\n',mean([double(xvrdd2LL(firstL  & FSssame)==0) double(xvrdd2LL(firstR  & FSssame)==1)])*100);
    
    fprintf('LookL on choice %2.2f%% of the trials where first, left offer is best.\n',  mean([double(xvrdchLL(Fisbest & firstL)==1)  double(xvrdchLL(Fisbest & firstR)==0)])*100);
    fprintf('LookL on choice %2.2f%% of the trials where first, left offer is worse.\n', mean([double(xvrdchLL(Sisbest & firstL)==1)  double(xvrdchLL(Sisbest & firstR)==0)])*100);
    fprintf('LookL on choice %2.2f%% of the trials where first, left offer is same.\n',  mean([double(xvrdchLL(firstL  & FSssame)==1) double(xvrdchLL(firstR  & FSssame)==0)])*100);
    fprintf('LookR on choice %2.2f%% of the trials where first, left offer is best.\n',  mean([double(xvrdchLL(Fisbest & firstL)==0)  double(xvrdchLL(Fisbest & firstR)==1)])*100);
    fprintf('LookR on choice %2.2f%% of the trials where first, left offer is worse.\n', mean([double(xvrdchLL(Sisbest & firstL)==0)  double(xvrdchLL(Sisbest & firstR)==1)])*100);
    fprintf('LookR on choice %2.2f%% of the trials where first, left offer is same.\n\n',mean([double(xvrdchLL(firstL  & FSssame)==0) double(xvrdchLL(firstR  & FSssame)==1)])*100);
    
   
    %%
    
    clc
    fprintf('LookL on offer1 %2.2f%% of the trials where chosen offer is Left.\n',   mean([double(xvrdo1LL(firstL & chooseR==0))==1 double(xvrdo1LL(firstR & chooseR==1))==0])*100);
    fprintf('LookR on offer1 %2.2f%% of the trials where chosen offer is Left.\n',   mean([double(xvrdo1LL(firstL & chooseR==0))==0 double(xvrdo1LL(firstR & chooseR==1))==1])*100);
    fprintf('LookL on offer1 %2.2f%% of the trials where chosen offer is Right.\n',  mean([double(xvrdo1LL(firstL & chooseR==1))==1 double(xvrdo1LL(firstR & chooseR==0))==0])*100);
    fprintf('LookR on offer1 %2.2f%% of the trials where chosen offer is Right.\n\n',mean([double(xvrdo1LL(firstL & chooseR==1))==0 double(xvrdo1LL(firstR & chooseR==0))==1])*100);
    
    fprintf('LookL on delay1 %2.2f%% of the trials where chosen offer is Left.\n',   mean([double(xvrdd1LR(firstL & chooseR==0))==0 double(xvrdd1LR(firstR & chooseR==1))==1])*100);
    fprintf('LookR on delay1 %2.2f%% of the trials where chosen offer is Left.\n',   mean([double(xvrdd1LR(firstL & chooseR==0))==1 double(xvrdd1LR(firstR & chooseR==1))==0])*100);
    fprintf('LookL on delay1 %2.2f%% of the trials where chosen offer is Right.\n',  mean([double(xvrdd1LR(firstL & chooseR==1))==0 double(xvrdd1LR(firstR & chooseR==0))==1])*100);
    fprintf('LookR on delay1 %2.2f%% of the trials where chosen offer is Right.\n\n',mean([double(xvrdd1LR(firstL & chooseR==1))==1 double(xvrdd1LR(firstR & chooseR==0))==0])*100);
    
    fprintf('LookL on offer2 %2.2f%% of the trials where chosen offer is Left.\n',   mean([double(xvrdo2LR(firstL & chooseR==0))==0 double(xvrdo2LR(firstR & chooseR==1))==1])*100);
    fprintf('LookR on offer2 %2.2f%% of the trials where chosen offer is Left.\n',   mean([double(xvrdo2LR(firstL & chooseR==0))==1 double(xvrdo2LR(firstR & chooseR==1))==0])*100);
    fprintf('LookL on offer2 %2.2f%% of the trials where chosen offer is Right.\n',  mean([double(xvrdo2LR(firstL & chooseR==1))==0 double(xvrdo2LR(firstR & chooseR==0))==1])*100);
    fprintf('LookR on offer2 %2.2f%% of the trials where chosen offer is Right.\n\n',mean([double(xvrdo2LR(firstL & chooseR==1))==1 double(xvrdo2LR(firstR & chooseR==0))==0])*100);

    fprintf('LookL on delay2 %2.2f%% of the trials where chosen offer is Left.\n',   mean([double(xvrdd2LL(firstL & chooseR==0))==1 double(xvrdd2LL(firstR & chooseR==1))==0])*100);
    fprintf('LookR on delay2 %2.2f%% of the trials where chosen offer is Left.\n',   mean([double(xvrdd2LL(firstL & chooseR==0))==0 double(xvrdd2LL(firstR & chooseR==1))==1])*100);
    fprintf('LookL on delay2 %2.2f%% of the trials where chosen offer is Right.\n',  mean([double(xvrdd2LL(firstL & chooseR==1))==1 double(xvrdd2LL(firstR & chooseR==0))==0])*100);
    fprintf('LookR on delay2 %2.2f%% of the trials where chosen offer is Right.\n\n',mean([double(xvrdd2LL(firstL & chooseR==1))==0 double(xvrdd2LL(firstR & chooseR==0))==1])*100);

    fprintf('LookL on choice %2.2f%% of the trials where chosen offer is Left.\n',   mean([double(xvrdchLL(firstL & chooseR==0))==1 double(xvrdchLL(firstR & chooseR==1))==0])*100);
    fprintf('LookR on choice %2.2f%% of the trials where chosen offer is Left.\n',   mean([double(xvrdchLL(firstL & chooseR==0))==0 double(xvrdchLL(firstR & chooseR==1))==1])*100);
    fprintf('LookL on choice %2.2f%% of the trials where chosen offer is Right.\n',  mean([double(xvrdchLL(firstL & chooseR==1))==1 double(xvrdchLL(firstR & chooseR==0))==0])*100);
    fprintf('LookR on choice %2.2f%% of the trials where chosen offer is Right.\n\n',mean([double(xvrdchLL(firstL & chooseR==1))==0 double(xvrdchLL(firstR & chooseR==0))==1])*100);

    
    %%
    fullfig();
    subplot(3,1,1);
    %SUBJECTS1&2
    for tt=1:Nvtps_ds; [hc(:,tt),he]=histcounts(vvrss(:,tt),xebins,'Normalization', 'probability'); end
    imagesc(rtmax_ds,he,hc); colormap(getdivergentcolmap(repmat(1:-.25:0,3,1)',[8 16 32 128])); hold on; view(2)
    for jj=1:8; plot([1 1].*sum(rtmes_ds(1:jj)),get(gca,'YLim'),'k:'); text(sum(rtmes_ds(1:jj-1)),max(get(gca,'YLim')),twLabels{jj}); end
    plot([0 sum(rtmes_ds)],[0 0],'k:');        
    xlabel('time bins (10 ms)'); ylabel('eye pos. (<0 is L; >0 is R)');
        
    subplot(3,1,2);
    %SUBJECT1
    for tt=1:Nvtps_ds; [hc(:,tt),he]=histcounts(vectorise(pooleposxmat_vds(:,1:sum(nccs(1:4)),tt)),xebins,'Normalization', 'probability'); end
    imagesc(rtmax_ds,he,hc); colormap(getdivergentcolmap(repmat(1:-.25:0,3,1)',[8 16 32 128])); hold on; view(2)
    for jj=1:8; plot([1 1].*sum(rtmes_ds(1:jj)),get(gca,'YLim'),'k:'); text(sum(rtmes_ds(1:jj-1)),max(get(gca,'YLim')),twLabels{jj}); end
    plot([0 sum(rtmes_ds)],[0 0],'k:');        
    xlabel('time bins (10 ms)'); ylabel('eye pos. (<0 is L; >0 is R)');
    
    subplot(3,1,3);
    %SUBJECT2
    for tt=1:Nvtps_ds; [hc(:,tt),he]=histcounts(vectorise(pooleposxmat_vds(:,sum(nccs(1:4))+1:end,tt)),xebins,'Normalization', 'probability'); end
    imagesc(rtmax_ds,he,hc); colormap(getdivergentcolmap(repmat(1:-.25:0,3,1)',[8 16 32 128])); hold on; view(2)
    for jj=1:8; plot([1 1].*sum(rtmes_ds(1:jj)),get(gca,'YLim'),'k:'); text(sum(rtmes_ds(1:jj-1)),max(get(gca,'YLim')),twLabels{jj}); end
    plot([0 sum(rtmes_ds)],[0 0],'k:');        
    xlabel('time bins (10 ms)'); ylabel('eye pos. (<0 is L; >0 is R)');
    
    
        %%
    fullfig();
    subplot(3,1,1);
    %SUBJECTS1&2
    for tt=1:Nvtps_ds; [hc(:,tt),he]=histcounts(vectorise(pooleposxmat_vrdds(:,tt)),xebins,'Normalization', 'probability'); end
    imagesc(rtmax_ds,he,hc); colormap(getdivergentcolmap(repmat(1:-.25:0,3,1)',[8 16 32 128])); hold on; view(2)
    for jj=1:8; plot([1 1].*sum(rtmes_ds(1:jj)),get(gca,'YLim'),'k:'); text(sum(rtmes_ds(1:jj-1)),max(get(gca,'YLim')),twLabels{jj}); end
    plot([0 sum(rtmes_ds)],[0 0],'k:');        
    xlabel('time bins (10 ms)'); ylabel('eye pos. (<0 is L; >0 is R)');
        
    subplot(3,1,2);
    %SUBJECT1
    for tt=1:Nvtps_ds; [hc(:,tt),he]=histcounts(vectorise(pooleposxmat_vrdds(1:sum(ntrs(1:4)),tt)),xebins,'Normalization', 'probability'); end
    imagesc(rtmax_ds,he,hc); colormap(getdivergentcolmap(repmat(1:-.25:0,3,1)',[8 16 32 128])); hold on; view(2)
    for jj=1:8; plot([1 1].*sum(rtmes_ds(1:jj)),get(gca,'YLim'),'k:'); text(sum(rtmes_ds(1:jj-1)),max(get(gca,'YLim')),twLabels{jj}); end
    plot([0 sum(rtmes_ds)],[0 0],'k:');        
    xlabel('time bins (10 ms)'); ylabel('eye pos. (<0 is L; >0 is R)');
    
    subplot(3,1,3);
    %SUBJECT2
    for tt=1:Nvtps_ds; [hc(:,tt),he]=histcounts(vectorise(pooleposxmat_vrdds(sum(ntrs(1:4))+1:end,tt)),xebins,'Normalization', 'probability'); end
    imagesc(rtmax_ds,he,hc); colormap(getdivergentcolmap(repmat(1:-.25:0,3,1)',[8 16 32 128])); hold on; view(2)
    for jj=1:8; plot([1 1].*sum(rtmes_ds(1:jj)),get(gca,'YLim'),'k:'); text(sum(rtmes_ds(1:jj-1)),max(get(gca,'YLim')),twLabels{jj}); end
    plot([0 sum(rtmes_ds)],[0 0],'k:');        
    xlabel('time bins (10 ms)'); ylabel('eye pos. (<0 is L; >0 is R)');
    
    
    %%
    fullfig();
    subplot(2,1,1);
    % LOOK LEFT DURING OFFER1
    for tt=1:Nvtps_ds; [hc(:,tt),he]=histcounts(vvrss(xvro1LL,tt),xebins,'Normalization', 'probability'); end
    imagesc(rtmax_ds,he,hc); colormap(getdivergentcolmap(repmat(1:-.25:0,3,1)',[8 16 32 128])); hold on; view(2)
    for jj=1:8; plot([1 1].*sum(rtmes_ds(1:jj)),get(gca,'YLim'),'k:'); text(sum(rtmes_ds(1:jj-1)),max(get(gca,'YLim')),twLabels{jj}); end
    plot([0 sum(rtmes_ds)],[0 0],'k:');        text(Offer1LLtime(1),-10000,sprintf('LookL:%2.2f%%',mean(xvro1LL)*100),'color','r');
    plot([1 1].*Offer1LLtime(1),get(gca,'YLim'),'r:');
    xlabel('time bins (10 ms)'); ylabel('eye pos. (<0 is L; >0 is R)');
    
    subplot(2,1,2);
    %LOOK RIGHT DURING OFFER1
    for tt=1:Nvtps_ds; [hc(:,tt),he]=histcounts(vvrss(~xvro1LL,tt),xebins,'Normalization', 'probability'); end
    imagesc(rtmax_ds,he,hc); colormap(getdivergentcolmap(repmat(1:-.25:0,3,1)',[8 16 32 128])); hold on; view(2)
    for jj=1:8; plot([1 1].*sum(rtmes_ds(1:jj)),get(gca,'YLim'),'k:'); text(sum(rtmes_ds(1:jj-1)),max(get(gca,'YLim')),twLabels{jj}); end
    plot([0 sum(rtmes_ds)],[0 0],'k:');         text(Offer1LLtime(1),+10000,sprintf('LookR:%2.2f%%',mean(~xvro1LL)*100),'color','r');
    plot([1 1].*Offer1LLtime(1),get(gca,'YLim'),'r:');
    xlabel('time bins (10 ms)'); ylabel('eye pos. (<0 is L; >0 is R)');
    supertitle('Compare offer1 side: LookL (top), LookR (bottom)');
    %%
    
    fullfig();
    subplot(2,1,1);
    % LOOK RIGHT DURING DELAY1
    for tt=1:Nvtps_ds; [hc(:,tt),he]=histcounts(vvrss(xvrd1LR,tt),xebins,'Normalization', 'probability'); end
    imagesc(rtmax_ds,he,hc); colormap(getdivergentcolmap(repmat(1:-.25:0,3,1)',[8 16 32 128])); hold on; view(2)
    for jj=1:8; plot([1 1].*sum(rtmes_ds(1:jj)),get(gca,'YLim'),'k:'); text(sum(rtmes_ds(1:jj-1)),max(get(gca,'YLim')),twLabels{jj}); end
    plot([0 sum(rtmes_ds)],[0 0],'k:');        text(Delay1LRtime(1),+10000,sprintf('LookR:%2.2f%%',mean(xvrd1LR)*100),'color','r');
    plot([1 1].*Delay1LRtime(1),get(gca,'YLim'),'r:');
    xlabel('time bins (10 ms)'); ylabel('eye pos. (<0 is L; >0 is R)');
    
    subplot(2,1,2);
    %LOOK LEFT DURING DELAY1
    for tt=1:Nvtps_ds; [hc(:,tt),he]=histcounts(vvrss(~xvrd1LR,tt),xebins,'Normalization', 'probability'); end
    imagesc(rtmax_ds,he,hc); colormap(getdivergentcolmap(repmat(1:-.25:0,3,1)',[8 16 32 128])); hold on; view(2)
    for jj=1:8; plot([1 1].*sum(rtmes_ds(1:jj)),get(gca,'YLim'),'k:'); text(sum(rtmes_ds(1:jj-1)),max(get(gca,'YLim')),twLabels{jj}); end
    plot([0 sum(rtmes_ds)],[0 0],'k:');         text(Delay1LRtime(1),-10000,sprintf('LookL:%2.2f%%',mean(~xvrd1LR)*100),'color','r');
    plot([1 1].*Delay1LRtime(1),get(gca,'YLim'),'r:');
    xlabel('time bins (10 ms)'); ylabel('eye pos. (<0 is L; >0 is R)');
    supertitle('Compare delay1 side: LookR (top), LookL (bottom)');
    
    %%
    fullfig();
    subplot(2,1,1);
    % LOOK LEFT DURING OFFER1 AND LEFT DURING DELAY2
    for tt=1:Nvtps_ds; [hc(:,tt),he]=histcounts(vvrss(xvro1LL&(~xvrd1LR)&xvrd2LL,tt),xebins,'Normalization', 'probability'); end
    imagesc(rtmax_ds,he,hc); colormap(getdivergentcolmap(repmat(1:-.25:0,3,1)',[8 16 32 128])); hold on; view(2)
    for jj=1:8; plot([1 1].*sum(rtmes_ds(1:jj)),get(gca,'YLim'),'k:'); text(sum(rtmes_ds(1:jj-1)),max(get(gca,'YLim')),twLabels{jj}); end
    plot([0 sum(rtmes_ds)],[0 0],'k:');        text(Delay2LLtime(1),-10000,sprintf('LookL:%2.2f%%',mean(xvro1LL&(~xvrd1LR)&xvrd2LL)*100),'color','r');
    plot([1 1].*Offer1LLtime(1),get(gca,'YLim'),'r:');
    plot([1 1].*Delay1LRtime(1),get(gca,'YLim'),'r:');
    plot([1 1].*Delay2LLtime(1),get(gca,'YLim'),'r:');
    xlabel('time bins (10 ms)'); ylabel('eye pos. (<0 is L; >0 is R)');
    
    subplot(2,1,2);
    %LOOK LEFT DURING OFFER1 AND RIGHT DURING DELAY2
    for tt=1:Nvtps_ds; [hc(:,tt),he]=histcounts(vvss(xvro1LL&(~xvrd1LR)&(~xvrd2LL),tt),xebins,'Normalization', 'probability'); end
    imagesc(rtmax_ds,he,hc); colormap(getdivergentcolmap(repmat(1:-.25:0,3,1)',[8 16 32 128])); hold on; view(2)
    for jj=1:8; plot([1 1].*sum(rtmes_ds(1:jj)),get(gca,'YLim'),'k:'); text(sum(rtmes_ds(1:jj-1)),max(get(gca,'YLim')),twLabels{jj}); end
    plot([0 sum(rtmes_ds)],[0 0],'k:');         text(Delay2LLtime(1),+10000,sprintf('LookR:%2.2f%%',mean(xvro1LL&(~xvrd1LR)&(~xvrd2LL))*100),'color','r');
    plot([1 1].*Offer1LLtime(1),get(gca,'YLim'),'r:');
    plot([1 1].*Delay1LRtime(1),get(gca,'YLim'),'r:');
    plot([1 1].*Delay2LLtime(1),get(gca,'YLim'),'r:');
    xlabel('time bins (10 ms)'); ylabel('eye pos. (<0 is L; >0 is R)');
    supertitle('offer1 LookL, delay1 LookL, compare delay2 side: LookL (top), LookR (bottom)');
    
    %%
    Vpl(1,:)=( xvro1LL)&(~xvrd1LR);% offer1=L  delay1=L
    Vpl(2,:)=( xvro1LL)&( xvrd1LR);% offer1=L  delay1=R
    Vpl(3,:)=(~xvro1LL)&(~xvrd1LR);% offer1=R  delay1=L
    Vpl(4,:)=(~xvro1LL)&( xvrd1LR);% offer1=R  delay1=R
    Vpl(5,:)=(~xvro2LR)&( xvrd2LL);% offer2=L  delay1=L
    Vpl(6,:)=(~xvro2LR)&(~xvrd2LL);% offer2=L  delay1=R
    Vpl(7,:)=( xvro2LR)&( xvrd2LL);% offer2=R  delay1=L
    Vpl(8,:)=( xvro2LR)&(~xvrd2LL);% offer2=R  delay1=R
    
    Spl={'Off1L,Del1L=','Off1L,Del1R=','Off1R,Del1L=','Off1R,Del1R=',...
         'Off2L,Del2L=','off2L,Del2R=','Off2R,Del2L=','Off2R,Del2R='};
    Ppl=[0 0; 0 1; 1 0; 1 1];
     
    for jj=[1 3 5 7]
    fullfig([Ppl((jj+1)/2,:)*.5 .5 .5]);
    subplot(2,1,1);
    % LOOK RIGHT DURING OFFER2 AND LEFT DURING DELAY2
    for tt=1:Nvtps_ds; [hc(:,tt),he]=histcounts(vvrss(Vpl(jj,:),tt),xebins,'Normalization', 'probability'); end
    imagesc(rtmax_ds,he,hc); colormap(getdivergentcolmap(repmat(1:-.25:0,3,1)',[8 16 32 128])); hold on; view(2)
    for tt=1:8; plot([1 1].*sum(rtmes_ds(1:tt)),get(gca,'YLim'),'k:'); text(sum(rtmes_ds(1:tt-1)),max(get(gca,'YLim')),twLabels{tt}); end
    plot([0 sum(rtmes_ds)],[0 0],'k:');        text(Delay2LLtime(1),-10000,sprintf('%s:%2.2f%%',Spl{jj},mean(Vpl(jj,:))*100),'color','r');
    plot([1 1].*Offer2LRtime(1),get(gca,'YLim'),'r:');
    plot([1 1].*Delay2LLtime(1),get(gca,'YLim'),'r:');
    xlabel('time bins (10 ms)'); ylabel('eye pos. (<0 is L; >0 is R)');
    
    subplot(2,1,2);
    %LOOK RIGHT DURING OFFER2 AND RIGHT DURING DELAY2
    for tt=1:Nvtps_ds; [hc(:,tt),he]=histcounts(vvrss(Vpl(jj+1,:),tt),xebins,'Normalization', 'probability'); end
    imagesc(rtmax_ds,he,hc); colormap(getdivergentcolmap(repmat(1:-.25:0,3,1)',[8 16 32 128])); hold on; view(2)
    for tt=1:8; plot([1 1].*sum(rtmes_ds(1:tt)),get(gca,'YLim'),'k:'); text(sum(rtmes_ds(1:tt-1)),max(get(gca,'YLim')),twLabels{tt}); end
    plot([0 sum(rtmes_ds)],[0 0],'k:');         text(Delay2LLtime(1),+10000,sprintf('%s:%2.2f%%',Spl{jj+1},mean(Vpl(jj+1,:))*100),'color','r');
    plot([1 1].*Offer2LRtime(1),get(gca,'YLim'),'r:');
    plot([1 1].*Delay2LLtime(1),get(gca,'YLim'),'r:');
    xlabel('time bins (10 ms)'); ylabel('eye pos. (<0 is L; >0 is R)');

    end
    
end

%% Regress EV vs Spiking
close all;
if recomputelm
    lmELpv=nan(Nccs,Ntps_ds);   lmELpv_shf=nan(Nccs,Ntps_ds,Nshf);   lmERpv=nan(Nccs,Ntps_ds);   lmERpv_shf=nan(Nccs,Ntps_ds,Nshf);
    lmELpvLL=nan(Nccs,Ntps_ds); lmELpvLL_shf=nan(Nccs,Ntps_ds,Nshf); lmERpvLL=nan(Nccs,Ntps_ds); lmERpvLL_shf=nan(Nccs,Ntps_ds,Nshf);
    lmELpvLR=nan(Nccs,Ntps_ds); lmELpvLR_shf=nan(Nccs,Ntps_ds,Nshf); lmERpvLR=nan(Nccs,Ntps_ds); lmERpvLR_shf=nan(Nccs,Ntps_ds,Nshf);
    lmELR2=nan(Nccs,Ntps_ds);   lmELR2_shf=nan(Nccs,Ntps_ds,Nshf);   lmERR2=nan(Nccs,Ntps_ds);   lmERR2_shf=nan(Nccs,Ntps_ds,Nshf);
    lmELR2LL=nan(Nccs,Ntps_ds); lmELR2LL_shf=nan(Nccs,Ntps_ds,Nshf); lmERR2LL=nan(Nccs,Ntps_ds); lmERR2LL_shf=nan(Nccs,Ntps_ds,Nshf);
    lmELR2LR=nan(Nccs,Ntps_ds); lmELR2LR_shf=nan(Nccs,Ntps_ds,Nshf); lmERR2LR=nan(Nccs,Ntps_ds); lmERR2LR_shf=nan(Nccs,Ntps_ds,Nshf);
    lmELb0=nan(Nccs,Ntps_ds);   lmELb0_shf=nan(Nccs,Ntps_ds,Nshf);   lmERb0=nan(Nccs,Ntps_ds);   lmERb0_shf=nan(Nccs,Ntps_ds,Nshf);
    lmELb0LL=nan(Nccs,Ntps_ds); lmELb0LL_shf=nan(Nccs,Ntps_ds,Nshf); lmERb0LL=nan(Nccs,Ntps_ds); lmERb0LL_shf=nan(Nccs,Ntps_ds,Nshf);
    lmELb0LR=nan(Nccs,Ntps_ds); lmELb0LR_shf=nan(Nccs,Ntps_ds,Nshf); lmERb0LR=nan(Nccs,Ntps_ds); lmERb0LR_shf=nan(Nccs,Ntps_ds,Nshf);
    lmELb1=nan(Nccs,Ntps_ds);   lmELb1_shf=nan(Nccs,Ntps_ds,Nshf);   lmERb1=nan(Nccs,Ntps_ds);   lmERb1_shf=nan(Nccs,Ntps_ds,Nshf);
    lmELb1LL=nan(Nccs,Ntps_ds); lmELb1LL_shf=nan(Nccs,Ntps_ds,Nshf); lmERb1LL=nan(Nccs,Ntps_ds); lmERb1LL_shf=nan(Nccs,Ntps_ds,Nshf);
    lmELb1LR=nan(Nccs,Ntps_ds); lmELb1LR_shf=nan(Nccs,Ntps_ds,Nshf); lmERb1LR=nan(Nccs,Ntps_ds); lmERb1LR_shf=nan(Nccs,Ntps_ds,Nshf);
    lmntrsLL=nan(Nccs,Ntps_ds); lmntrsLR=nan(Nccs,Ntps_ds); lmntrsLL_nn=nan(Nccs,Ntps_ds); lmntrsLR_nn=nan(Nccs,Ntps_ds); lmntrsLL_nn_shLalgn=nan(Nccs,Ntps_ds); lmntrsLR_nn_shRalgn=nan(Nccs,Ntps_ds);
    for cc=1:Nccs
        ccoff1rw=pooloffer1rw(:,cc);
        ccoff2rw=pooloffer2rw(:,cc);
        ccoff1ev=pooloffer1ev(:,cc);
        ccoff2ev=pooloffer2ev(:,cc);
        ccspk=squeeze(poolspksqmat_ds_bln(:,cc,:)); % For linear fit
        %%%%%    ccspk=squeeze(poolspksqmat_ds(:,cc,:));      % For Poisson fit
        ccepx=squeeze(pooleposxmat_ds_bln(:,cc,:));
        
        %%% shift aligned
        ccspk_shLalgn=squeeze(poolspksqmat_shLalgn(:,cc,:));
        ccspk_shRalgn=squeeze(poolspksqmat_shRalgn(:,cc,:));
        ccepx_shLalgn=squeeze(pooleposxmat_shLalgn(:,cc,:));
        ccepx_shRalgn=squeeze(pooleposxmat_shRalgn(:,cc,:));
        ccisdelSL=isdelSL(:,cc);
        ccisdelSR=isdelSR(:,cc);
        
        %     cciso1LLd2LL=iso1LLd2LL(:,cc);
        %     cciso1LLd2LR=iso1LLd2LR(:,cc);
        %cciso1LLd1LLd2LL=iso1LLd1LLd2LL(:,cc);
        %cciso1LLd1LLd2LR=iso1LLd1LLd2LR(:,cc);
        %     cciso2LRd2LL=iso2LRd2LL(:,cc);
        %     cciso2LLd2LL=iso2LLd2LL(:,cc);
        %     cciso2LRd2LR=iso2LRd2LR(:,cc);
        %     cciso2LLd2LR=iso2LLd2LR(:,cc);
        rmallnans=isnan(ccoff1ev);
        rmtrials=rmallnans;%(rmallnans | ccoff1rw~=1 | ccoff2rw~=1);
        %ccspk(rmallnans,:)=[];
        %ccepx(rmallnans,:)=[];
        %ccoff1ev(rmallnans)=[];
        %ccoff2ev(rmallnans)=[];
        ccspk(rmtrials,:)=[];
        ccepx(rmtrials,:)=[];
        ccoff1ev(rmtrials)=[];
        ccoff2ev(rmtrials)=[];
        ccoff1rw(rmtrials)=[];
        ccoff2rw(rmtrials)=[];
        %cciso2LRd2LL(rmtrials)=[];
        %cciso2LRd2LR(rmtrials)=[];
        %cciso2LLd2LL(rmtrials)=[];
        %cciso2LLd2LR(rmtrials)=[];
        %cciso1LLd1LLd2LL(rmtrials)=[];
        %cciso1LLd1LLd2LR(rmtrials)=[];
        ccisdelSL(rmtrials)=[];
        ccisdelSR(rmtrials)=[];
        ccspk_shLalgn(rmtrials,:)=[];
        ccspk_shRalgn(rmtrials,:)=[];
        ccepx_shLalgn(rmtrials,:)=[];
        ccepx_shRalgn(rmtrials,:)=[];
        ccntrs=size(ccspk,1);
        
        %    ccspk_shf=nan(Nshf,size(ccspk,1),size(ccspk,2));                        %time-shuffle
        %    for jj=1:Nshf                                                           %time-shuffle
        %        for ttr=1:ccntrs                                                    %time-shuffle
        %            isnn=find(~isnan(ccspk(ttr,:)));                                %time-shuffle
        %            ccspk_shf(jj,ttr,isnn)=ccspk(ttr,isnn(randperm(length(isnn)))); %time-shuffle
        %            %isvtps=find(vtps_ds==1);                                              %time-shuffle
        %            %ccspk_shf(jj,ttr,isvtps)=ccspk(ttr,isvtps(randperm(length(isvtps)))); %time-shuffle
        %        end                                                                 %time-shuffle
        %    end                                                                     %time-shuffle
        parfor tt=eval(selectedTimes) % delay1
        %for tt=eval(selectedTimes)
        %parfor tt=(ntce_ds(5,1):ntce_ds(5,2)) % delay2
            %for tt=1:Ntps_ds
            warning('off','all')
            if 1%vtps_ds(tt)==1 % discards low trials times
                ctspk=squeeze(ccspk(:,tt));
                ctepx=squeeze(ccepx(:,tt));
                ctspk_shLalgn=squeeze(ccspk_shLalgn(:,tt));
                ctspk_shRalgn=squeeze(ccspk_shRalgn(:,tt));
                
                ctlm1=fitlm(ccoff1ev,ctspk);
                ctlm2=fitlm(ccoff2ev,ctspk);
                %            ctlm1=fitglm(ccoff1ev,ctspk,'linear','Distribution','poisson');
                %            ctlm2=fitglm(ccoff2ev,ctspk,'linear','Distribution','poisson');
                lmELpv(cc,tt)=ctlm1.Coefficients.pValue(2);
                lmERpv(cc,tt)=ctlm2.Coefficients.pValue(2);
                lmELR2(cc,tt)=ctlm1.Rsquared.Ordinary;
                lmERR2(cc,tt)=ctlm2.Rsquared.Ordinary;
                lmELb0(cc,tt)=ctlm1.Coefficients.Estimate(1);
                lmERb0(cc,tt)=ctlm2.Coefficients.Estimate(1);
                lmELb1(cc,tt)=ctlm1.Coefficients.Estimate(2);
                lmERb1(cc,tt)=ctlm2.Coefficients.Estimate(2);
                
                cciso1LL=(nanmean(ccepx(:,off1time),2)<0);
                ccisd1LL=(nanmean(ccepx(:,del1time),2)<0);
                cciso2LL=(nanmean(ccepx(:,off2time),2)<0);
                ccisd2LL=(nanmean(ccepx(:,del2time),2)<0);
                
                %[nuL,~,~,ishiftsL]=nanmonotonicdrifts_forward_simplified(ccepx(:,fulldel1time),2,0,-1,1,1);
                %check1=find(~isnan(ishiftsL));
                %check0=find((ccisdelSL));
                                
%                ccisd1allLL=(arrayfun(@(xx) all((ccepx(xx,fulldel1time(1)-1+find(~isnan(ccepx(xx,fulldel1time)))))<=0),1:size(ccepx,1)));
%                ccisd1allLR=(arrayfun(@(xx) all((ccepx(xx,fulldel1time(1)-1+find(~isnan(ccepx(xx,fulldel1time)))))>=0),1:size(ccepx,1)));
                
                ccisd1allLL=(arrayfun(@(xx) all((ccepx(xx,fulldel1time(41)-1+find(~isnan(ccepx(xx,fulldel1time(41:end))))))<=0),1:size(ccepx,1)));
                ccisd1allLR=(arrayfun(@(xx) all((ccepx(xx,fulldel1time(41)-1+find(~isnan(ccepx(xx,fulldel1time(41:end))))))>=0),1:size(ccepx,1)));
                
                %fullfig(); subplot(1,2,1); plot(fulldel1time,ccepx(ccisd1allLL,fulldel1time));  
                %           subplot(1,2,2); plot(fulldel1time,ccepx(ccisd1allLR,fulldel1time)); 
                           
%                ccisd2allLL=(arrayfun(@(xx) all((ccepx(xx,fulldel2time(1)-1+find(~isnan(ccepx(xx,fulldel2time)))))<=0),1:size(ccepx,1)));
%                ccisd2allLR=(arrayfun(@(xx) all((ccepx(xx,fulldel2time(1)-1+find(~isnan(ccepx(xx,fulldel2time)))))>=0),1:size(ccepx,1)));

                ccisd2allLL=(arrayfun(@(xx) all((ccepx(xx,fulldel2time(41)-1+find(~isnan(ccepx(xx,fulldel2time(41:end))))))<=0),1:size(ccepx,1)));
                ccisd2allLR=(arrayfun(@(xx) all((ccepx(xx,fulldel2time(41)-1+find(~isnan(ccepx(xx,fulldel2time(41:end))))))>=0),1:size(ccepx,1)));
                %fullfig(); subplot(1,2,1); plot(fulldel2time,ccepx(ccisd2allLL,fulldel2time));  
                %           subplot(1,2,2); plot(fulldel2time,ccepx(ccisd2allLR,fulldel2time)); 
                
                switch selflag
                    case 1 % OFFER1 LookL and DELAY1 LookL/R
                        %isLL=find((cciso1LL) & (ccepx(:,tt)<0));
                        %isLR=find((cciso1LL) & (ccepx(:,tt)>0));
                        %isLL=find((cciso1LL) & (ccisdelSL));
                        %isLR=find((cciso1LL) & (ccisdelSR));
                        %isLL=find((ccisdelSL));
                        %isLR=find((ccisdelSR));
                        switch caseflag
                            case {1,2}
                                isLL=find((ccisdelSL));
                                isLR=find((ccisdelSR));
                            case 3
                                isLL=find((ccisd1allLL));
                                isLR=find((ccisd1allLR));
                            case {4,5}
                                isLL=find((cciso1LL) & (ccisdelSL));
                                isLR=find((cciso1LL) & (ccisdelSR));
                            case 6
                                isLL=find((cciso1LL) & (ccisd1allLL)');
                                isLR=find((cciso1LL) & (ccisd1allLR)');
                        end
                    case 2 % OFFER1 LookR and DELAY1 LookL/R
                        %isLL=find((~cciso1LL) & (ccepx(:,tt)<0));
                        %isLR=find((~cciso1LL) & (ccepx(:,tt)>0));
                        isLL=find((~cciso1LL) & (ccisdelSL));
                        isLR=find((~cciso1LL) & (ccisdelSR));
                    case 3 % OFFER2 LookL and DELAY2 LookL/R 
                        %isLL=find((cciso2LL) & (ccepx(:,tt)<0));
                        %isLR=find((cciso2LL) & (ccepx(:,tt)>0));
                        isLL=find((cciso2LL) & (ccisdelSL));
                        isLR=find((cciso2LL) & (ccisdelSR));
                    case 4 % OFFER2 LookR and DELAY2 LookL/R
                        %isLL=find((~cciso2LL) & (ccepx(:,tt)<0));
                        %isLR=find((~cciso2LL) & (ccepx(:,tt)>0));
                        %isLL=find((~cciso2LL) & (ccisdelSL));
                        %isLR=find((~cciso2LL) & (ccisdelSR));
                        %isLL=find((ccisdelSL));
                        %isLR=find((ccisdelSR));
                        switch caseflag
                            case {1,2}
                                isLL=find((ccisdelSL));
                                isLR=find((ccisdelSR));
                            case 3
                                isLL=find((ccisd2allLL));
                                isLR=find((ccisd2allLR));
                            case {4,5}
                                isLL=find((~cciso2LL) & (ccisdelSL));
                                isLR=find((~cciso2LL) & (ccisdelSR));
                            case 6
                                isLL=find((~cciso2LL) & (ccisd2allLL)');
                                isLR=find((~cciso2LL) & (ccisd2allLR)');
                        end
                        %isLL=find((~cciso2LL) & (ccisd2allLL)');
                        %isLR=find((~cciso2LL) & (ccisd2allLR)');
                        %isLL=find((ccisd2allLL));
                        %isLR=find((ccisd2allLR));
                end
                
                %Relevant cases:
                % Left in offer1 and Right in delay1
                % Right in offer2 and Left in delay2
                % control cases:
                % Right in offer1
                % Left in offer2
                niLL=length(isLL);
                if (caseflag==3 || caseflag==6)
                    ctlm1LL=fitlm(ccoff1ev(isLL),ctspk(isLL));
                    ctlm2LL=fitlm(ccoff2ev(isLL),ctspk(isLL));
                else
                    ctlm1LL=fitlm(ccoff1ev(isLL),ctspk_shLalgn(isLL));
                    ctlm2LL=fitlm(ccoff2ev(isLL),ctspk_shLalgn(isLL));
                end
                lmELpvLL(cc,tt)=ctlm1LL.Coefficients.pValue(2);
                lmERpvLL(cc,tt)=ctlm2LL.Coefficients.pValue(2);
                lmELR2LL(cc,tt)=ctlm1LL.Rsquared.Ordinary;
                lmERR2LL(cc,tt)=ctlm2LL.Rsquared.Ordinary;
                lmELb0LL(cc,tt)=ctlm1LL.Coefficients.Estimate(1);
                lmERb0LL(cc,tt)=ctlm2LL.Coefficients.Estimate(1);
                lmELb1LL(cc,tt)=ctlm1LL.Coefficients.Estimate(2);
                lmERb1LL(cc,tt)=ctlm2LL.Coefficients.Estimate(2);
                lmntrsLL(cc,tt)=niLL;
                lmntrsLL_nn(cc,tt)=sum(~isnan(ctspk(isLL)));
                lmntrsLL_nn_shLalgn(cc,tt)=sum(~isnan(ctspk_shLalgn(isLL)));
                
                niLR=length(isLR);
                if (caseflag==3 || caseflag==6)
                    ctlm1LR=fitlm(ccoff1ev(isLR),ctspk(isLR));
                    ctlm2LR=fitlm(ccoff2ev(isLR),ctspk(isLR));
                else
                    ctlm1LR=fitlm(ccoff1ev(isLR),ctspk_shRalgn(isLR));
                    ctlm2LR=fitlm(ccoff2ev(isLR),ctspk_shRalgn(isLR));
                end
                lmELpvLR(cc,tt)=ctlm1LR.Coefficients.pValue(2);
                lmERpvLR(cc,tt)=ctlm2LR.Coefficients.pValue(2);
                lmELR2LR(cc,tt)=ctlm1LR.Rsquared.Ordinary;
                lmERR2LR(cc,tt)=ctlm2LR.Rsquared.Ordinary;
                lmELb0LR(cc,tt)=ctlm1LR.Coefficients.Estimate(1);
                lmERb0LR(cc,tt)=ctlm2LR.Coefficients.Estimate(1);
                lmELb1LR(cc,tt)=ctlm1LR.Coefficients.Estimate(2);
                lmERb1LR(cc,tt)=ctlm2LR.Coefficients.Estimate(2);
                lmntrsLR(cc,tt)=niLR;
                lmntrsLR_nn(cc,tt)=sum(~isnan(ctspk(isLR)));
                lmntrsLR_nn_shRalgn(cc,tt)=sum(~isnan(ctspk_shRalgn(isLR)));
                
                for shf=1:Nshf
                    i1shf=randperm(ccntrs);                                   %trial-shuffle
                    i2shf=randperm(ccntrs);                                   %trial-shuffle
                    isshf=randperm(ccntrs);                                   %trial-shuffle
                    ctlm1_shf=fitlm(ccoff1ev(i1shf),ctspk(isshf));            %trial-shuffle
                    ctlm2_shf=fitlm(ccoff2ev(i2shf),ctspk(isshf));            %trial-shuffle
                    lmELpv_shf(cc,tt,shf)=ctlm1_shf.Coefficients.pValue(2);
                    lmERpv_shf(cc,tt,shf)=ctlm2_shf.Coefficients.pValue(2);
                    lmELR2_shf(cc,tt,shf)=ctlm1_shf.Rsquared.Ordinary;
                    lmERR2_shf(cc,tt,shf)=ctlm2_shf.Rsquared.Ordinary;
                    lmELb0_shf(cc,tt,shf)=ctlm1_shf.Coefficients.Estimate(1);
                    lmERb0_shf(cc,tt,shf)=ctlm2_shf.Coefficients.Estimate(1);
                    lmELb1_shf(cc,tt,shf)=ctlm1_shf.Coefficients.Estimate(2);
                    lmERb1_shf(cc,tt,shf)=ctlm2_shf.Coefficients.Estimate(2);
                    
                    i1LLshf=isLL(randperm(niLL));                                        %trial-shuffle
                    i2LLshf=isLL(randperm(niLL));                                        %trial-shuffle
                    isLLshf=isLL(randperm(niLL));                                        %trial-shuffle
                    if (caseflag==3 || caseflag==6)
                        ctlm1LL_shf=fitlm(ccoff1ev(i1LLshf),ctspk(isLLshf));         %trial-shuffle
                        ctlm2LL_shf=fitlm(ccoff2ev(i2LLshf),ctspk(isLLshf));         %trial-shuffle
                    else
                        ctlm1LL_shf=fitlm(ccoff1ev(i1LLshf),ctspk_shLalgn(isLLshf));         %trial-shuffle
                        ctlm2LL_shf=fitlm(ccoff2ev(i2LLshf),ctspk_shLalgn(isLLshf));         %trial-shuffle
                    end
                    lmELpvLL_shf(cc,tt,shf)=ctlm1LL_shf.Coefficients.pValue(2);
                    lmERpvLL_shf(cc,tt,shf)=ctlm2LL_shf.Coefficients.pValue(2);
                    lmELR2LL_shf(cc,tt,shf)=ctlm1LL_shf.Rsquared.Ordinary;
                    lmERR2LL_shf(cc,tt,shf)=ctlm2LL_shf.Rsquared.Ordinary;
                    lmELb0LL_shf(cc,tt,shf)=ctlm1LL_shf.Coefficients.Estimate(1);
                    lmERb0LL_shf(cc,tt,shf)=ctlm2LL_shf.Coefficients.Estimate(1);
                    lmELb1LL_shf(cc,tt,shf)=ctlm1LL_shf.Coefficients.Estimate(2);
                    lmERb1LL_shf(cc,tt,shf)=ctlm2LL_shf.Coefficients.Estimate(2);
                    
                    i1LRshf=isLR(randperm(niLR));                                        %trial-shuffle
                    i2LRshf=isLR(randperm(niLR));                                        %trial-shuffle
                    isLRshf=isLR(randperm(niLR));                                        %trial-shuffle
                    if (caseflag==3 || caseflag==6)
                        ctlm1LR_shf=fitlm(ccoff1ev(i1LRshf),ctspk(isLRshf));         %trial-shuffle
                        ctlm2LR_shf=fitlm(ccoff2ev(i2LRshf),ctspk(isLRshf));         %trial-shuffle
                    else
                        ctlm1LR_shf=fitlm(ccoff1ev(i1LRshf),ctspk_shRalgn(isLRshf));         %trial-shuffle
                        ctlm2LR_shf=fitlm(ccoff2ev(i2LRshf),ctspk_shRalgn(isLRshf));         %trial-shuffle
                    end
                    lmELpvLR_shf(cc,tt,shf)=ctlm1LR_shf.Coefficients.pValue(2);
                    lmERpvLR_shf(cc,tt,shf)=ctlm2LR_shf.Coefficients.pValue(2);
                    lmELR2LR_shf(cc,tt,shf)=ctlm1LR_shf.Rsquared.Ordinary;
                    lmERR2LR_shf(cc,tt,shf)=ctlm2LR_shf.Rsquared.Ordinary;
                    lmELb0LR_shf(cc,tt,shf)=ctlm1LR_shf.Coefficients.Estimate(1);
                    lmERb0LR_shf(cc,tt,shf)=ctlm2LR_shf.Coefficients.Estimate(1);
                    lmELb1LR_shf(cc,tt,shf)=ctlm1LR_shf.Coefficients.Estimate(2);
                    lmERb1LR_shf(cc,tt,shf)=ctlm2LR_shf.Coefficients.Estimate(2);
                end
            end
            warning('on','all')
        end
        disp([sprintf('%s',datetime('now','Format','HH:mm')) ' - done: ' num2str(cc) '/' num2str(Nccs) ', ' num2str(cc/Nccs.*100) '%']);
    end
    save([most_recent_folder 'fitlmstats_span_' num2str(Nspan) 'ms_lag_' num2str(Nlag) 'ms_' date '.mat'],'lmE*','lmntr*');
else
    %load(['fitlmstats_span_' num2str(Nspan) 'ms_lag_' num2str(Nlag) 'ms_13-May-2022.mat']);
end



%%
hfig=fullfig(); hfig.Position=[0 0 .5 1];
if any(caseflag==[1 2 4 5])
    subplot(2,1,1);
    plot(rtmax_ds,(lmntrsLL_nn_shLalgn(cumsum(nccs),vtps_ds==1))./repmat(ntrs',1,nansum(vtps_ds)).*100);
    subplot(2,1,2);
    plot(rtmax_ds,(lmntrsLR_nn_shRalgn(cumsum(nccs),vtps_ds==1))./repmat(ntrs',1,nansum(vtps_ds)).*100);
else
    subplot(2,1,1);
    plot(rtmax_ds,(lmntrsLL_nn(cumsum(nccs),vtps_ds==1))./repmat(ntrs',1,nansum(vtps_ds)).*100);
    subplot(2,1,2);
    plot(rtmax_ds,(lmntrsLR_nn(cumsum(nccs),vtps_ds==1))./repmat(ntrs',1,nansum(vtps_ds)).*100);
end
switch selflag
    case 1
        subplot(2,1,1);
        title('Look Left Offer1 & Look Left Delay1'); ylim([0 100]); ylabel('% trials'); xlabel('Delay1 time bins (10 ms)');
        subplot(2,1,2);
        title('Look Left Offer1 & Look Right Delay1'); ylim([0 100]); ylabel('% trials'); xlabel('Delay1 time bins (10 ms)');
    case 2
        subplot(2,1,1);
        title('Look Right Offer1 & Look Left Delay1'); ylim([0 100]); ylabel('% trials'); xlabel('Delay1 time bins (10 ms)');
        subplot(2,1,2);
        title('Look Right Offer1 & Look Right Delay1'); ylim([0 100]); ylabel('% trials'); xlabel('Delay1 time bins (10 ms)');
    case 3
        subplot(2,1,1);
        title('Look Left Offer2 & Look Left Delay2'); ylim([0 100]); ylabel('% trials'); xlabel('Delay2 time bins (10 ms)');
        subplot(2,1,2);
        title('Look Left Offer2 & Look Right Delay2'); ylim([0 100]); ylabel('% trials'); xlabel('Delay2 time bins (10 ms)');
    case 4
        subplot(2,1,1);
        title('Look Right Offer2 & Look Left Delay2'); ylim([0 100]); ylabel('% trials'); xlabel('Delay2 time bins (10 ms)');
        subplot(2,1,2);
        title('Look Right Offer2 & Look Right Delay2'); ylim([0 100]); ylabel('% trials'); xlabel('Delay2 time bins (10 ms)');
end
hl=legend(strsplit(sprintf('session %d-',1:8),'-'));
saveas(hfig,[most_recent_folder 'span' num2str(Nspan) 'ms_lag' num2str(Nlag) 'ms/frac_trials_off_del_selected.png']);
saveas(hfig,[most_recent_folder 'span' num2str(Nspan) 'ms_lag' num2str(Nlag) 'ms/frac_trials_off_del_selected.fig']);
saveas(hfig,[most_recent_folder 'span' num2str(Nspan) 'ms_lag' num2str(Nlag) 'ms/frac_trials_off_del_selected.svg']);




%% Plot colors
o1col = [ 55 188  52]/255;
o2col = [ 55 131 200]/255;
o3col = [121 121 121]/255;

o1LLcol=[ 50 148  48]/255;
o2LLcol=[ 55  97 135]/255;
o3LLcol=[ 43  43  43]/255;

o1LRcol=[135 215 133]/255;
o2LRcol=[144 187 226]/255;
o3LRcol=[168 168 168]/255;

o1LLcol=[212 87 28]/255;
o1LRcol=[251 179 144]/255;
o2LLcol=[119 30 152]/255;
o2LRcol=[185 118 210]/255;

%%

lmELpvd=double(lmELpv(:,vtps_ds==1)<0.05);
lmERpvd=double(lmERpv(:,vtps_ds==1)<0.05);
lmELpvdLL=double(lmELpvLL(:,vtps_ds==1)<0.05);
lmERpvdLL=double(lmERpvLL(:,vtps_ds==1)<0.05);
lmELpvdLR=double(lmELpvLR(:,vtps_ds==1)<0.05);
lmERpvdLR=double(lmERpvLR(:,vtps_ds==1)<0.05);
lmELpvd_shf=squeeze(double(lmELpv_shf(:,vtps_ds==1,:)<0.05));
lmERpvd_shf=squeeze(double(lmERpv_shf(:,vtps_ds==1,:)<0.05));
lmELpvdLL_shf=squeeze(double(lmELpvLL_shf(:,vtps_ds==1,:)<0.05));
lmERpvdLL_shf=squeeze(double(lmERpvLL_shf(:,vtps_ds==1,:)<0.05));
lmELpvdLR_shf=squeeze(double(lmELpvLR_shf(:,vtps_ds==1,:)<0.05));
lmERpvdLR_shf=squeeze(double(lmERpvLR_shf(:,vtps_ds==1,:)<0.05));

hfig=fullfig();
subplot(3,1,1);
plotcmapdots([o1col; o2col]); hold on
runlength(nanmean(lmELpvd),squeeze(nanmean(lmELpvd_shf)),1,rtmax_ds,o1col,0.005);
runlength(nanmean(lmERpvd),squeeze(nanmean(lmERpvd_shf)),1,rtmax_ds,o2col,0.002);
for jj=1:8; plot([1 1].*sum(rtmes_ds(1:jj)),[0 0.2],'k:'); text(sum(rtmes_ds(1:jj-1)),.21,twLabels{jj}); end
hl=legend({'E(L)','E(R)'}); ylim([0 0.2]); box(hl,'on');
xlabel('time bins (10 ms)'); ylabel('% signif. cells');

subplot(3,1,2);
plotcmapdots([o1LLcol; o1LRcol]); hold on
runlength(nanmean(lmELpvdLL),squeeze(nanmean(lmELpvdLL_shf)),1,rtmax_ds,o1LLcol,0.005);
runlength(nanmean(lmELpvdLR),squeeze(nanmean(lmELpvdLR_shf)),1,rtmax_ds,o1LRcol,0.002);
for jj=1:8; plot([1 1].*sum(rtmes_ds(1:jj)),[0 0.2],'k:'); text(sum(rtmes_ds(1:jj-1)),.21,twLabels{jj}); end
hl=legend({'E(L), LookL','E(L), LookR'}); ylim([0 0.2]); box(hl,'on');
xlabel('time bins (10 ms)'); ylabel('% signif. cells');

subplot(3,1,3);
plotcmapdots([o2LLcol; o2LRcol]);
runlength(nanmean(lmERpvdLL),squeeze(nanmean(lmERpvdLL_shf)),1,rtmax_ds,o2LLcol,0.005);
runlength(nanmean(lmERpvdLR),squeeze(nanmean(lmERpvdLR_shf)),1,rtmax_ds,o2LRcol,0.002);
for jj=1:8; plot([1 1].*sum(rtmes_ds(1:jj)),[0 0.2],'k:'); text(sum(rtmes_ds(1:jj-1)),.21,twLabels{jj}); end
hl=legend({'E(R), LookL','E(R), LookR'}); ylim([0 0.2]); box(hl,'on');
xlabel('time bins (10 ms)'); ylabel('% signif. cells');
saveas(hfig,[most_recent_folder 'span' num2str(Nspan) 'ms_lag' num2str(Nlag) 'ms/sig_in_time_pv.png']);
saveas(hfig,[most_recent_folder 'span' num2str(Nspan) 'ms_lag' num2str(Nlag) 'ms/sig_in_time_pv.fig']);
saveas(hfig,[most_recent_folder 'span' num2str(Nspan) 'ms_lag' num2str(Nlag) 'ms/sig_in_time_pv.svg']);
%%

lmELpvd=double(lmELpv(:,:)<0.05);
lmERpvd=double(lmERpv(:,:)<0.05);
lmELpvdLL=double(lmELpvLL(:,:)<0.05);
lmERpvdLL=double(lmERpvLL(:,:)<0.05);
lmELpvdLR=double(lmELpvLR(:,:)<0.05);
lmERpvdLR=double(lmERpvLR(:,:)<0.05);
lmELpvd_shf=squeeze(double(lmELpv_shf(:,:,:)<0.05));
lmERpvd_shf=squeeze(double(lmERpv_shf(:,:,:)<0.05));
lmELpvdLL_shf=squeeze(double(lmELpvLL_shf(:,:,:)<0.05));
lmERpvdLL_shf=squeeze(double(lmERpvLL_shf(:,:,:)<0.05));
lmELpvdLR_shf=squeeze(double(lmELpvLR_shf(:,:,:)<0.05));
lmERpvdLR_shf=squeeze(double(lmERpvLR_shf(:,:,:)<0.05));

hfig=fullfig();
subplot(3,1,1);
plotcmapdots([o1col; o2col]);
runlength(nanmean(lmELpvd),squeeze(nanmean(lmELpvd_shf)),1,1:Ntps_ds,o1col,0.005);
runlength(nanmean(lmERpvd),squeeze(nanmean(lmERpvd_shf)),1,1:Ntps_ds,o2col,0.002);
ylim([0 0.2]); for jj=1:7; plotvline(ntce_ds(jj,1:2)); text(ntce_ds(jj,1),max(get(gca,'YLim')),twLabels{jj}); end
hl=legend({'E(L)','E(R)'}); box(hl,'on'); box off; 
if selflag<=2; xlim(fulldel1time([1 end])); else; xlim(fulldel2time([1 end])); end
xlabel('time bins (10 ms)'); ylabel('% signif. cells');

subplot(3,1,2);
plotcmapdots([o1LLcol; o1LRcol]);
runlength(nanmean(lmELpvdLL),squeeze(nanmean(lmELpvdLL_shf)),1,1:Ntps_ds,o1LLcol,0.005);
runlength(nanmean(lmELpvdLR),squeeze(nanmean(lmELpvdLR_shf)),1,1:Ntps_ds,o1LRcol,0.002);
ylim([0 0.2]); for jj=1:7; plotvline(ntce_ds(jj,1:2)); text(ntce_ds(jj,1),max(get(gca,'YLim')),twLabels{jj}); end
hl=legend({'E(L), LookL','E(L), LookR'}); box(hl,'on'); box off;
if selflag<=2; xlim(fulldel1time([1 end])); else; xlim(fulldel2time([1 end])); end
xlabel('time bins (10 ms)'); ylabel('% signif. cells');

subplot(3,1,3);
plotcmapdots([o2LLcol; o2LRcol]);
runlength(nanmean(lmERpvdLL),squeeze(nanmean(lmERpvdLL_shf)),1,1:Ntps_ds,o2LLcol,0.005);
runlength(nanmean(lmERpvdLR),squeeze(nanmean(lmERpvdLR_shf)),1,1:Ntps_ds,o2LRcol,0.002);
ylim([0 0.2]); for jj=1:7; plotvline(ntce_ds(jj,1:2)); text(ntce_ds(jj,1),max(get(gca,'YLim')),twLabels{jj}); end
hl=legend({'E(R), LookL','E(R), LookR'}); box(hl,'on'); box off; 
if selflag<=2; xlim(fulldel1time([1 end])); else; xlim(fulldel2time([1 end])); end
xlabel('time bins (10 ms)'); ylabel('% signif. cells');
saveas(hfig,[most_recent_folder 'span' num2str(Nspan) 'ms_lag' num2str(Nlag) 'ms/sig_in_time_pv_deltime.png']);
saveas(hfig,[most_recent_folder 'span' num2str(Nspan) 'ms_lag' num2str(Nlag) 'ms/sig_in_time_pv_deltime.fig']);
saveas(hfig,[most_recent_folder 'span' num2str(Nspan) 'ms_lag' num2str(Nlag) 'ms/sig_in_time_pv_deltime.svg']);


end



end