%This program calculates spike triggered landmark plot
%for a position to be considered, I need to have at least 3 spikes at that
%position
clear all
close all
exp_date=input('please type in the experiment date : ','s');
load(strcat(exp_date,'_Stage1_output.mat'));
load(strcat(exp_date,'_Stage2_output.mat'));
%load(strcat(exp_date,'_Stage8_output.mat'));
ppc=7.5;%pixels per cm
for tn=1:length(hd_crct) %for all trials
    %% load the landmark locations data
    load(strcat('landmark_coordinates_',filename{tn}(1:8),'.mat'));
    tankfilt=filter;
    clear filter
    y=0:10:1600;
    x=0:10:1200;
    resy=1600/10;
    resx=1200/10;
    LMoccup=zeros(resx,resy);
    filtgrid=polygrid(1200-tankfilt(:,2),tankfilt(:,1),0.01);
    hcgrid=polygrid(1200-hc(:,2),hc(:,1),0.01);
    stpgrid=polygrid(1200-stp(:,2),stp(:,1),0.01);
    ltpgrid=polygrid(1200-ltp(:,2),ltp(:,1),0.01);
    sqpgrid=polygrid(1200-sqp(:,2),sqp(:,1),0.01);
%%
    %for each spike find the center of the fish and set the grid around
    %that
    ndata=length(nuraw_count{tn});        
    for i=1:(resx)
        for j=1:(resy)
            %set the value of the grid points that are occupied by the tank
            %walls to 1
            temp1=((find((1200-tank_circle(2,:))>=x(i)&(1200-tank_circle(2,:))<x(i+1)&(tank_circle(1,:))>=y(j)&(tank_circle(1,:)<y(j+1)))));
            if ~isempty(temp1)
                LMoccup(i,j)=1;
            end
            %generate points that are inside the polygons belonging to
            %landmarks
            temp1=((find((filtgrid(:,1))>=x(i)&(filtgrid(:,1))<x(i+1)&(filtgrid(:,2))>=y(j)&(filtgrid(:,2)<y(j+1)))));
            if ~isempty(temp1)
                LMoccup(i,j)=1;
            end
            
            %             %check for other landmarks
            if filename{tn}(11)=='2'
                temp1=((find((hcgrid(:,1))>=x(i)&(hcgrid(:,1))<x(i+1)&(hcgrid(:,2))>=y(j)&(hcgrid(:,2)<y(j+1)))));
                if ~isempty(temp1)
                    LMoccup(i,j)=1;
                end
            elseif filename{tn}(11)~='7'
                temp1=((find((hcgrid(:,1))>=x(i)&(hcgrid(:,1))<x(i+1)&(hcgrid(:,2))>=y(j)&(hcgrid(:,2)<y(j+1)))));
                if ~isempty(temp1)
                    LMoccup(i,j)=1;
                end
                temp1=((find((stpgrid(:,1))>=x(i)&(stpgrid(:,1))<x(i+1)&(stpgrid(:,2))>=y(j)&(stpgrid(:,2)<y(j+1)))));
                if ~isempty(temp1)
                    LMoccup(i,j)=1;
                end
                temp1=((find((ltpgrid(:,1))>=x(i)&(ltpgrid(:,1))<x(i+1)&(ltpgrid(:,2))>=y(j)&(ltpgrid(:,2)<y(j+1)))));
                if ~isempty(temp1)
                    LMoccup(i,j)=1;
                end
                temp1=((find((sqpgrid(:,1))>=x(i)&(sqpgrid(:,1))<x(i+1)&(sqpgrid(:,2))>=y(j)&(sqpgrid(:,2)<y(j+1)))));
                if ~isempty(temp1)
                    LMoccup(i,j)=1;
                end
                temp1=((find((1200-sc(2,:))>=x(i)&(1200-sc(2,:))<x(i+1)&(sc(1,:))>=y(j)&(sc(1,:)<y(j+1)))));
                if ~isempty(temp1)
                    LMoccup(i,j)=1;
                end
                temp1=((find((1200-lc(2,:))>=x(i)&(1200-lc(2,:))<x(i+1)&(lc(1,:))>=y(j)&(lc(1,:)<y(j+1)))));
                if ~isempty(temp1)
                    LMoccup(i,j)=1;
                end
            else
                temp1=((find((stpgrid(:,1))>=x(i)&(stpgrid(:,1))<x(i+1)&(stpgrid(:,2))>=y(j)&(stpgrid(:,2)<y(j+1)))));
                if ~isempty(temp1)
                    LMoccup(i,j)=1;
                end
                temp1=((find((ltpgrid(:,1))>=x(i)&(ltpgrid(:,1))<x(i+1)&(ltpgrid(:,2))>=y(j)&(ltpgrid(:,2)<y(j+1)))));
                if ~isempty(temp1)
                    LMoccup(i,j)=1;
                end
                temp1=((find((sqpgrid(:,1))>=x(i)&(sqpgrid(:,1))<x(i+1)&(sqpgrid(:,2))>=y(j)&(sqpgrid(:,2)<y(j+1)))));
                if ~isempty(temp1)
                    LMoccup(i,j)=1;
                end
                temp1=((find((1200-sc(2,:))>=x(i)&(1200-sc(2,:))<x(i+1)&(sc(1,:))>=y(j)&(sc(1,:)<y(j+1)))));
                if ~isempty(temp1)
                    LMoccup(i,j)=1;
                end
                temp1=((find((1200-lc(2,:))>=x(i)&(1200-lc(2,:))<x(i+1)&(lc(1,:))>=y(j)&(lc(1,:)<y(j+1)))));
                if ~isempty(temp1)
                    LMoccup(i,j)=1;
                end
            end
            
            
        end
    end
    %%calculate the position triggered landmark average. This will be
    %%calculated for each trial and used to normalize spike triggered LM to
    %%account for time spent at a given location'
    head=downsample(hd_crct{tn},20);
    tail=downsample(COM_crct{tn},20);%time between frames is 1.33 s
    clear dist
    for i=1:length(head)
            dist(i)=pdist([head(i,1) head(i,2);tail(i,1) tail(i,2)],'euclidean');
            mid(i,:)=[(head(i,1)+tail(i,1))/2,(head(i,2)+tail(i,2))/2];
            h2tvec=[(head(i,1)-tail(i,1)),(head(i,2)-tail(i,2)),0];
            %calculate the translation vector that brings the midpoint to
            %construct a grid around the center point
            %calculate the heading angle of the fish relative to the north
            %pointing arrow
            n0=[800,400];
            n1=[800,600];
            N=[n1(1)-n0(1),n1(2)-n0(2),0];
            crossvec=cross(h2tvec,N);
            ThetaInDegreesL(i) =atan2d(norm(crossvec),dot(N,h2tvec))*sign(crossvec(3));
            x_center=800;
            y_center=600;
            center = repmat([x_center; y_center], 1, 1);
            R = [cosd(ThetaInDegreesL(i)) -sind(ThetaInDegreesL(i)); sind(ThetaInDegreesL(i)) cosd(ThetaInDegreesL(i))];
            midshift=mid(i,:)'-center;
            s0=R*midshift;
            midshift=s0+center;
            transvect=[(800-midshift(1))/10,(600-midshift(2))/10,];
            %construct the landmark configuration around the fish
            RotateImage=imrotate(LMoccup, -1*ThetaInDegreesL(i),'bilinear','crop');
            TranslateImageL{i}=imtranslate(RotateImage, transvect);
    end
    range1=mean(dist)+3*std(dist);
    range2=mean(dist)-3*std(dist);
    index=find(dist(dist<range1 & dist>range2));
    M=zeros(length(index),resx,resy);
    for i=1:length(index)
            M(i,:,:)=cell2mat(TranslateImageL(index(i)));
    end
    PosTrigLM{tn}=(squeeze(sum(M)));
    figure
    imagesc(PosTrigLM{tn})
    colorbar
    set(gca,'YDir','normal')
    hold on
    %%saveas(gcf,strcat(filename{tn},'_Postriglm.pdf'))
    max_val=max(reshape(PosTrigLM{tn},[1,120*160]));
    Lindex(tn)=length(index);

    for i=1:120
        for j=1:160
            %if PosTrigLM{tn}(i,j)<(length(index)*0.05) || i<40 || i>70 || j<70 || j>90
            if i<40 || i>70 || j<70 || j>90

                PosTrigLM{tn}(i,j)=NaN;
            end
        end
    end   
        
    figure
    imagesc(PosTrigLM{tn})
    colorbar
    set(gca,'YDir','normal')
    hold on    

    
    for w=1:ndata
            clear dist    

        for k=1:nuraw_count{tn}(w)
            %figure
            %for each spike find the mid point between head and tail
            dist(k)=pdist([Place_int_unit{tn,w,1}(k) Place_int_unit{tn,w,2}(k);Place_int_unit_tail{tn,w,1}(k) Place_int_unit_tail{tn,w,2}(k)]);
            mid=[(Place_int_unit{tn,w,1}(k)+Place_int_unit_tail{tn,w,1}(k))/2,(Place_int_unit{tn,w,2}(k)+(Place_int_unit_tail{tn,w,2}(k)))/2];
            h2tvec=[Place_int_unit{tn,w,1}(k)-Place_int_unit_tail{tn,w,1}(k),Place_int_unit{tn,w,2}(k)-(Place_int_unit_tail{tn,w,2}(k)),0];
            %calculate the translation vector that brings the midpoint to
            %construct a grid around the center point
            %calculate the heading angle of the fish relative to the north
            %pointing arrow
            n0=[800,400];
            n1=[800,600];
            N=[n1(1)-n0(1),n1(2)-n0(2),0];
            crossvec=cross(h2tvec,N);
            ThetaInDegrees(k) =atan2d(norm(crossvec),dot(N,h2tvec))*sign(crossvec(3));
            x_center=800;
            y_center=600;
            center = repmat([x_center; y_center], 1, 1);
            R = [cosd(ThetaInDegrees(k)) -sind(ThetaInDegrees(k)); sind(ThetaInDegrees(k)) cosd(ThetaInDegrees(k))];
            midshift=mid'-center;
            s0=R*midshift;
            midshift=s0+center;
            transvect=[(800-midshift(1))/10,(600-midshift(2))/10,];
            %construct the landmark configuration around the fish
            RotateImage=imrotate(LMoccup, -1*ThetaInDegrees(k),'bilinear','crop');
            TranslateImage{k,w}=imtranslate(RotateImage, transvect);
        end
        range1=mean(dist)+3*std(dist);
        range2=mean(dist)-3*std(dist);
        index=find(dist(dist<range1 & dist>range2));

         M=zeros(length(index),resx,resy);
        for i=1:length(index)
            M(i,:,:)=cell2mat(TranslateImage(index(i),w));
        end
        SpkTrigLM{tn,w}=(squeeze(sum(M)));
        figure
        imagesc(SpkTrigLM{tn,w})
        colorbar
        set(gca,'YDir','normal')
        hold on
        line([80 80],[47.5 62.5])
%         axis equal

        title(strcat('trial number',num2str(tn),'cell',uname{w,tn}))
        %%saveas(gcf,strcat(filename{tn},'_',uname{w,tn},'_spktriglm.pdf'))
        %plot normalized spike triggered LM
        
        figure
        NormSpkTrigLM{tn,w}=(SpkTrigLM{tn,w}./(PosTrigLM{tn}));
        %plot a up to 5 cm distance fom the fish
        rad=3*ppc;
        imagesc(NormSpkTrigLM{tn,w})
        colorbar
        set(gca,'YDir','normal')
        hold on
        axis equal

        hold on

        title(strcat('trial number',num2str(tn),'cell',uname{w,tn},'Normalized'))
        %%saveas(gcf,strcat(filename{tn},'_',uname{w,tn},'_Normspktriglm.pdf'))

        
        %calculate the location of the peak in the lmtriggered average
        [maxval,idx]=(max(SpkTrigLM{tn,w}(:)));
        [row,col]=ind2sub(size(SpkTrigLM{tn,w}), idx);
        %calculate distance to the head of the fish in pixels
        %a point between fish's head and its center of mass is rowf=80,colf=60

        dist_cm{tn,w}=sqrt((row-60)^2+(col-80)^2)/ppc
        
        %calculate the location of the peak in the normalized lmtriggered average
        [maxval,idx]=(max(NormSpkTrigLM{tn,w}(:)));
        [row,col]=ind2sub(size(NormSpkTrigLM{tn,w}), idx);
        %calculate distance to the head of the fish in pixels
        %a point between fish's head and its center of mass is rowf=80,colf=60

        dist_cm_norm{tn,w}=sqrt((row-60)^2+(col-80)^2)/ppc

    


    end


end
close all
LtotalU1=0;
for i=1:6
 if i~=5   
    SpkTrigAllU1(i,:,:)=cell2mat(SpkTrigLM(i,1));
 else
    SpkTrigAllU1(i,:,:)=cell2mat(SpkTrigLM(i,5));
 end
    PosTrigAllU1(i,:,:)=cell2mat(PosTrigLM(i));
    LtotalU1=LtotalU1+Lindex(i);
end

sumSpkTrigLMU1=(squeeze(nansum(SpkTrigAllU1)));
figure
imagesc(sumSpkTrigLMU1)
set(gca,'YDir','normal')

sumPosTrigLMU1=squeeze(nansum(PosTrigAllU1));
figure
imagesc(sumPosTrigLMU1)
set(gca,'YDir','normal')
PosMean=prctile(reshape(sumPosTrigLMU1(40:70,70:90),1,31*21),10);
for i=1:120
    for j=1:160
        %if sumPosTrigLMU1(i,j)<(length(LtotalU1)*0.01) || i<40 || i>70 || j<70 || j>90
        if sumPosTrigLMU1(i,j)<PosMean || i<40 || i>70 || j<70 || j>90            
            sumPosTrigLMU1(i,j)=NaN;
        end
    end
end


NormSpkTrigLMU1=(sumSpkTrigLMU1./sumPosTrigLMU1);
% figure
% imagesc(NormSpkTrigLMU1)
% set(gca,'YDir','normal')

temp=imgaussfilt(NormSpkTrigLMU1);
[maxval,idx]=max(temp(:));
[row,col]=ind2sub(size(temp),idx);
MaxLoc(1,:)=[row-57,col-80];%subtract the fish com
stlmUname{1}='U1';
temp=temp/maxval;

reg_max(1,1)=max(reshape(temp(58:70,81:90),[1,numel(temp(58:70,81:90))]))
reg_max(1,2)=max(reshape(temp(58:70,70:80),[1,numel(temp(58:70,70:80))]));
reg_max(1,3)=max(reshape(temp(40:57,70:80),[1,numel(temp(40:57,70:80))]));
reg_max(1,4)=max(reshape(temp(40:57,81:90),[1,numel(temp(40:57,81:90))]));
%calculate an index head vs tail bias (:,1) up down preference, + up , -
%down, (:,2) left right preference + right - left
index_side_pref(1,1)=(max(reg_max(1,1:2))-max(reg_max(1,3:4)))/nansum([max(reg_max(1,1:2)),max(reg_max(1,3:4))]);
index_side_pref(1,2)=(max(reg_max(1,[1,4]))-max(reg_max(1,[2,3])))/nansum([max(reg_max(1,[1,4])),max(reg_max(1,[2,3]))]);

figure
imagesc(temp)
hold on
title('U1')
colorbar
set(gca,'YDir','normal')
hold on
line([80 80],[47.5 62.5])
xlim([70,90])
ylim([40,70])
axis equal
%%saveas(gcf,'NormU1spktriglm.pdf')

%for cell 2
LtotalU2=0;
for i=1:5
    if i~=5
        SpkTrigAllU2(i,:,:)=cell2mat(SpkTrigLM(i,2));
    else
        SpkTrigAllU2(i,:,:)=cell2mat(SpkTrigLM(i,1));
    end
    PosTrigAllU2(i,:,:)=cell2mat(PosTrigLM(i));
    LtotalU2=LtotalU2+Lindex(i);

end
sumSpkTrigLMU2=(squeeze(nansum(SpkTrigAllU2)));
% figure
% imagesc(sumSpkTrigLMU1)
% set(gca,'YDir','normal')

sumPosTrigLMU2=squeeze(nansum(PosTrigAllU2));
% figure
% imagesc(sumPosTrigLMU2)
% set(gca,'YDir','normal')
PosMean=prctile(reshape(sumPosTrigLMU2(40:70,70:90),1,31*21),10);

for i=1:120
    for j=1:160
       % if sumPosTrigLMU2(i,j)<(length(LtotalU2)*0.01) || i<40 || i>70 || j<70 || j>90
        if sumPosTrigLMU2(i,j)<PosMean || i<40 || i>70 || j<70 || j>90
            
            sumPosTrigLMU2(i,j)=NaN;
        end
    end
end


NormSpkTrigLMU2=(sumSpkTrigLMU2./sumPosTrigLMU2);
% figure
% imagesc(NormSpkTrigLMU2)
% set(gca,'YDir','normal')

temp=imgaussfilt(NormSpkTrigLMU2);
[maxval,idx]=max(temp(:))
[row,col]=ind2sub(size(temp), idx);
MaxLoc(2,:)=[row-57,col-80];
stlmUname{2}='U2';
temp=temp/maxval;

reg_max(2,1)=max(reshape(temp(58:70,81:90),[1,numel(temp(58:70,81:90))]));
reg_max(2,2)=max(reshape(temp(58:70,70:80),[1,numel(temp(58:70,70:80))]));
reg_max(2,3)=max(reshape(temp(40:57,70:80),[1,numel(temp(40:57,70:80))]));
reg_max(2,4)=max(reshape(temp(40:57,81:90),[1,numel(temp(40:57,81:90))]));
%calculate an index head vs tail bias (:,1) up down preference, + up , -
%down, (:,2) left right preference + right - left
index_side_pref(2,1)=(max(reg_max(2,1:2))-max(reg_max(2,3:4)))/nansum([max(reg_max(2,1:2)),max(reg_max(2,3:4))]);
index_side_pref(2,2)=(max(reg_max(2,[1,4]))-max(reg_max(2,[2,3])))/nansum([max(reg_max(2,[1,4])),max(reg_max(2,[2,3]))]);

figure
imagesc(temp)
hold on
title('U2')
colorbar
set(gca,'YDir','normal')
hold on
line([80 80],[47.5 62.5])
xlim([70,90])
ylim([40,70])
axis equal

%%saveas(gcf,'NormU2spktriglm.pdf')


sumSpkTrigLMU3=cell2mat(SpkTrigLM(2,3));
% figure
% imagesc(sumSpkTrigLMU3)
% colorbar
sumPosTrigLMU3=cell2mat(PosTrigLM(2));

% figure
% imagesc(sumPosTrigLMU3)
% colorbar
PosMean=prctile(reshape(sumPosTrigLMU3(40:70,70:90),1,31*21),10);

for i=1:120
    for j=1:160
       % if sumPosTrigLMU2(i,j)<(length(LtotalU2)*0.01) || i<40 || i>70 || j<70 || j>90
        if sumPosTrigLMU3(i,j)<PosMean || i<40 || i>70 || j<70 || j>90
            
            sumPosTrigLMU3(i,j)=NaN;
        end
    end
end
% figure
% imagesc(sumPosTrigAllU3)

NormSpkTrigLMU3=sumSpkTrigLMU3./sumPosTrigLMU3;
% figure
% imagesc(NormSpkTrigLMU3)

temp=imgaussfilt(NormSpkTrigLMU3);
[maxval,idx]=max(temp(:))
[row,col]=ind2sub(size(temp), idx);
MaxLoc(3,:)=[row-57,col-80];
stlmUname{3}='U3';
temp=temp/maxval;
reg_max(3,1)=max(reshape(temp(58:70,81:90),[1,numel(temp(58:70,81:90))]));
reg_max(3,2)=max(reshape(temp(58:70,70:80),[1,numel(temp(58:70,70:80))]));
reg_max(3,3)=max(reshape(temp(40:57,70:80),[1,numel(temp(40:57,70:80))]));
reg_max(3,4)=max(reshape(temp(40:57,81:90),[1,numel(temp(40:57,81:90))]));
%calculate an index head vs tail bias (:,1) up down preference, + up , -
%down, (:,2) left right preference + right - left
index_side_pref(3,1)=(max(reg_max(3,1:2))-max(reg_max(3,3:4)))/nansum([max(reg_max(3,1:2)),max(reg_max(3,3:4))]);
index_side_pref(3,2)=(max(reg_max(3,[1,4]))-max(reg_max(3,[2,3])))/nansum([max(reg_max(3,[1,4])),max(reg_max(3,[2,3]))]);

figure
imagesc(temp)
title('U3')
colorbar
set(gca,'YDir','normal')
hold on
line([80 80],[47.5 62.5])
xlim([70,90])
ylim([40,70])
axis equal
%%saveas(gcf,'NormU3spktriglm.pdf')




for i=1:4
    if i~=2
        SpkTrigAllU4(i,:,:)=cell2mat(SpkTrigLM(i+1,4));
    else
        SpkTrigAllU4(i,:,:)=cell2mat(SpkTrigLM(3,3));
    end
        PosTrigAllU4(i,:,:)=cell2mat(PosTrigLM(i+1));

end

sumSpkTrigLMU4=(squeeze(nansum(SpkTrigAllU4)));
sumPosTrigLMU4=squeeze(nansum(PosTrigAllU4));
PosMean=prctile(reshape(sumPosTrigLMU4(40:70,70:90),1,31*21),10);

for i=1:120
    for j=1:160
       % if sumPosTrigLMU2(i,j)<(length(LtotalU2)*0.01) || i<40 || i>70 || j<70 || j>90
        if sumPosTrigLMU4(i,j)<PosMean || i<40 || i>70 || j<70 || j>90
            
            sumPosTrigLMU4(i,j)=NaN;
        end
    end
end
NormSpkTrigLMU4=sumSpkTrigLMU4./sumPosTrigLMU4;

temp=imgaussfilt(NormSpkTrigLMU4);
[maxval,idx]=max(temp(:))
[row,col]=ind2sub(size(temp), idx);
MaxLoc(4,:)=[row-57,col-80];
stlmUname{4}='U4';
temp=temp/maxval;
reg_max(4,1)=max(reshape(temp(58:70,81:90),[1,numel(temp(58:70,81:90))]));
reg_max(4,2)=max(reshape(temp(58:70,70:80),[1,numel(temp(58:70,70:80))]));
reg_max(4,3)=max(reshape(temp(40:57,70:80),[1,numel(temp(40:57,70:80))]));
reg_max(4,4)=max(reshape(temp(40:57,81:90),[1,numel(temp(40:57,81:90))]));
%calculate an index head vs tail bias (:,1) up down preference, + up , -
%down, (:,2) left right preference + right - left
index_side_pref(4,1)=(max(reg_max(4,1:2))-max(reg_max(4,3:4)))/nansum([max(reg_max(4,1:2)),max(reg_max(4,3:4))]);
index_side_pref(4,2)=(max(reg_max(4,[1,4]))-max(reg_max(4,[2,3])))/nansum([max(reg_max(4,[1,4])),max(reg_max(4,[2,3]))]);


figure
imagesc(temp)
hold on
title('U4')
colorbar
set(gca,'YDir','normal')
line([80 80],[47.5 62.5])
xlim([70,90])
ylim([40,70])
axis equal
%%saveas(gcf,'NormU4spktriglm.pdf')


sumSpkTrigLMU5=(SpkTrigLM{3,4});
sumPosTrigLMU5=cell2mat(PosTrigLM(3));
PosMean=prctile(reshape(sumPosTrigLMU5(40:70,70:90),1,31*21),10);

for i=1:120
    for j=1:160
       % if sumPosTrigLMU2(i,j)<(length(LtotalU2)*0.01) || i<40 || i>70 || j<70 || j>90
        if sumPosTrigLMU5(i,j)<PosMean || i<40 || i>70 || j<70 || j>90
            
            sumPosTrigLMU5(i,j)=NaN;
        end
    end
end
NormSpkTrigLMU5=sumSpkTrigLMU5./sumPosTrigLMU5;

temp=imgaussfilt(NormSpkTrigLMU5);
[maxval,idx]=max(temp(:))
[row,col]=ind2sub(size(temp), idx);
MaxLoc(5,:)=[row-57,col-80];
stlmUname{5}='U5';
temp=temp/maxval;
reg_max(5,1)=max(reshape(temp(58:70,81:90),[1,numel(temp(58:70,81:90))]));
reg_max(5,2)=max(reshape(temp(58:70,70:80),[1,numel(temp(58:70,70:80))]));
reg_max(5,3)=max(reshape(temp(40:57,70:80),[1,numel(temp(40:57,70:80))]));
reg_max(5,4)=max(reshape(temp(40:57,81:90),[1,numel(temp(40:57,81:90))]));
%calculate an index head vs tail bias (:,1) up down preference, + up , -
%down, (:,2) left right preference + right - left
index_side_pref(5,1)=(max(reg_max(5,1:2))-max(reg_max(5,3:4)))/nansum([max(reg_max(5,1:2)),max(reg_max(5,3:4))]);
index_side_pref(5,2)=(max(reg_max(5,[1,4]))-max(reg_max(5,[2,3])))/nansum([max(reg_max(5,[1,4])),max(reg_max(5,[2,3]))]);


figure
imagesc(temp)
hold on
title('U5')
colorbar
set(gca,'YDir','normal')
hold on
line([80 80],[47.5 62.5])
xlim([70,90])
ylim([40,70])
axis equal

%%saveas(gcf,'NormU5spktriglm.pdf')


for i=1:5
    if i<3
        SpkTrigAllU6(i,:,:)=cell2mat(SpkTrigLM(i+1,5));
    elseif i==3
        SpkTrigAllU6(i,:,:)=cell2mat(SpkTrigLM(4,3));
    else
        SpkTrigAllU6(i,:,:)=cell2mat(SpkTrigLM(i+1,2));
    end
    PosTrigAllU6(i,:,:)=cell2mat(PosTrigLM(i+1));

end

sumSpkTrigLMU6=(squeeze(nansum(SpkTrigAllU6)));
figure
imagesc(sumSpkTrigLMU6)
set(gca,'YDir','normal')
sumPosTrigLMU6=squeeze(nansum(PosTrigAllU6));
ZZ=zscore(sumPosTrigLMU6);
figure
imagesc(sumPosTrigLMU6)
set(gca,'YDir','normal')
figure
imagesc(ZZ)
set(gca,'YDir','normal')
PosMean=prctile(reshape(sumPosTrigLMU6(40:70,70:90),1,31*21),10);

for i=1:120
    for j=1:160
       % if sumPosTrigLMU2(i,j)<(length(LtotalU2)*0.01) || i<40 || i>70 || j<70 || j>90
        if sumPosTrigLMU6(i,j)<PosMean || i<40 || i>70 || j<70 || j>90
            
            sumPosTrigLMU6(i,j)=NaN;
        end
    end
end
NormSpkTrigLMU6=sumSpkTrigLMU6./sumPosTrigLMU6;
figure
imagesc(NormSpkTrigLMU6)
set(gca,'YDir','normal')

temp=imgaussfilt(NormSpkTrigLMU6);
[maxval,idx]=max(temp(:))
[row,col]=ind2sub(size(temp), idx);
MaxLoc(6,:)=[row-57,col-80];
stlmUname{6}='U6';
temp=temp/maxval;
reg_max(6,1)=max(reshape(temp(58:70,81:90),[1,numel(temp(58:70,81:90))]));
reg_max(6,2)=max(reshape(temp(58:70,70:80),[1,numel(temp(58:70,70:80))]));
reg_max(6,3)=max(reshape(temp(40:57,70:80),[1,numel(temp(40:57,70:80))]));
reg_max(6,4)=max(reshape(temp(40:57,81:90),[1,numel(temp(40:57,81:90))]));
%calculate an index head vs tail bias (:,1) up down preference, + up , -
%down, (:,2) left right preference + right - left
index_side_pref(6,1)=(max(reg_max(6,1:2))-max(reg_max(6,3:4)))/nansum([max(reg_max(6,1:2)),max(reg_max(6,3:4))]);
index_side_pref(6,2)=(max(reg_max(6,[1,4]))-max(reg_max(6,[2,3])))/nansum([max(reg_max(6,[1,4])),max(reg_max(6,[2,3]))]);



figure
imagesc(temp)
hold on
title('U6')
colorbar
set(gca,'YDir','normal')
hold on
line([80 80],[47.5 62.5])
xlim([70,90])
ylim([40,70])
axis equal


%%saveas(gcf,'NormU6spktriglm.pdf')

%UA
sumSpkTrigLMUA=(cell2mat(SpkTrigLM(2,6)));
sumPosTrigLMUA=cell2mat(PosTrigLM(2));
PosMean=prctile(reshape(sumPosTrigLMUA(40:70,70:90),1,31*21),10);

for i=1:120
    for j=1:160
       % if sumPosTrigLMU2(i,j)<(length(LtotalU2)*0.01) || i<40 || i>70 || j<70 || j>90
        if sumPosTrigLMUA(i,j)<PosMean || i<40 || i>70 || j<70 || j>90
            
            sumPosTrigLMUA(i,j)=NaN;
        end
    end
end


NormSpkTrigLMUA=sumSpkTrigLMUA./sumPosTrigLMUA;

temp=imgaussfilt(NormSpkTrigLMUA);
[maxval,idx]=max(temp(:))
[row,col]=ind2sub(size(temp), idx);
MaxLoc(7,:)=[row-57,col-80];
stlmUname{7}='UA';
temp=temp/maxval;
reg_max(7,1)=max(reshape(temp(58:70,81:90),[1,numel(temp(58:70,81:90))]));
reg_max(7,2)=max(reshape(temp(58:70,70:80),[1,numel(temp(58:70,70:80))]));
reg_max(7,3)=max(reshape(temp(40:57,70:80),[1,numel(temp(40:57,70:80))]));
reg_max(7,4)=max(reshape(temp(40:57,81:90),[1,numel(temp(40:57,81:90))]));
%calculate an index head vs tail bias (:,1) up down preference, + up , -
%down, (:,2) left right preference + right - left
index_side_pref(7,1)=(max(reg_max(7,1:2))-max(reg_max(7,3:4)))/nansum([max(reg_max(7,1:2)),max(reg_max(7,3:4))]);
index_side_pref(7,2)=(max(reg_max(7,[1,4]))-max(reg_max(7,[2,3])))/nansum([max(reg_max(7,[1,4])),max(reg_max(7,[2,3]))]);



figure
imagesc(temp)
hold on
title('UA')
colorbar
set(gca,'YDir','normal')
hold on
line([80 80],[47.5 62.5])
xlim([70,90])
ylim([40,70])
axis equal
%%saveas(gcf,'NormUAspktriglm.pdf')

[t,r]=cart2pol(MaxLoc(:,2),MaxLoc(:,1));
ThetaLM=t*180/pi;

figure
histogram(cell2mat(dist_cm(:)),10, 'Normalization','Probability')
xlabel('Distance to LM peak (cm)')
ylabel('Probability')
figure
histogram(cell2mat(dist_cm_norm(:)),10, 'Normalization','Probability')
xlabel('Distance to LM peak of Normalized plot(cm)')
ylabel('Probability')

%%save(strcat(exp_date,'_Stage9_output_v6.mat'),'index_side_pref','reg_max','SpkTrigLM','dist_cm','dist_cm_norm','MaxLoc','ThetaLM','stlmUname');
