%this program analyzes behavior and puts it in the Lick struct
%can change what it's saved as here to later compare Int and Blocks
clear all;
load ('RAW.mat');
load ('R_2R.mat');

%which lick intervals to look at?
startinterval=1;
endinterval=30;

%time to look at difference in licking activity
earlylicks=[1 4.5];

%% conduct lick analysis
global Dura Tm BSIZE Tbin
tic
path='C:\Users\dottenh2\Documents\MATLAB\David\2Rewards Nex Files\2RLick_paper.xls';
if ~exist('RAW'), load ('RAW.mat'); end

%Main settings
BSIZE=0.01; %Do not change
Dura=[-22 20]; Tm=Dura(1):BSIZE:Dura(2);
Tbin=-0.5:0.005:0.5; %window used to determine the optimal binsize
MinNumTrials=5;
Lick=[];Lick.Ninfo={};LL=0;NSess=0;

%Smoothing
Smoothing=1; %0 for raw and 1 for smoothing
SmoothTYPE='lowess'; %can change this between lowess and rlowess (more robust, ignores outliers more)
SmoothSPAN=50; %percentage of total data points
if Smoothing~=1, SmoothTYPE='NoSmoothing';SmoothSPAN=NaN; end

% List of events to analyze and analysis windows EXTRACTED from excel file
[~,Erefnames]=xlsread(path,'Windows','a3:a10'); % cell that contains the event names

%Finds the total number of suc vs mal sessions
for i=1:length(RAW)
    if strcmp('NA',RAW(i).Type(1:2)) | strcmp('VP',RAW(i).Type(1:2))
        NSess=NSess+1;
        Lick.Linfo(NSess,1)=RAW(i).Ninfo(1,1);
    end
end
Lick.Erefnames= Erefnames;


%preallocating for behavior stats
Lick.RewDurR1=NaN(NSess,35);
Lick.RewDurR2=NaN(NSess,35);
Lick.NumLicksR1=NaN(NSess,35);
Lick.NumLicksR2=NaN(NSess,35);
Lick.LicksEarlyR1=NaN(NSess,35);
Lick.LicksEarlyR2=NaN(NSess,35);
Lick.LicksR1P1=NaN(NSess,35);
Lick.LicksR1P2=NaN(NSess,35);
Lick.LicksR2P1=NaN(NSess,35);
Lick.LicksR2P2=NaN(NSess,35);
for i=1:endinterval
    Lick.ILIR1{i,1}=NaN(NSess,35);
    Lick.ILIR2{i,1}=NaN(NSess,35);
end

%preallocating the result matrix
for k=1:length(Erefnames)
    Lick.Ev(k).PSTHraw(1:NSess,1:length(Tm))=NaN(NSess,length(Tm));
    Lick.Ev(k).BW(1:NSess,1)=NaN;
    Lick.Ev(k).NumberTrials(1:NSess,1)=NaN;
end

for i=1:length(RAW) %loops through sessions
    if strcmp('NA',RAW(i).Type(1:2)) | strcmp('VP',RAW(i).Type(1:2))
        LL=LL+1; %lick session counter
        for k=1:length(Erefnames) %loops thorough the events
            EvInd=strcmp(Erefnames(k),RAW(i).Einfo(:,2)); %find the event id number from RAW
            LickInd=strcmp('Licks',RAW(i).Einfo(:,2)); %find the event id number from RAW
            if sum(EvInd)==0
                fprintf('HOWDY, CANT FIND EVENTS FOR ''%s''\n',Erefnames{k});
            end

            Lick.Ev(k).NumberTrials(LL,1)=length(RAW(i).Erast{EvInd});
            if  ~isempty(EvInd) && Lick.Ev(k).NumberTrials(LL,1)>MinNumTrials %avoid analyzing sessions where that do not have enough trials
                [PSR1,N1]=MakePSR04(RAW(i).Erast(LickInd),RAW(i).Erast{EvInd},Dura,{1});% makes collpased rasters. PSR1 is a cell(neurons)
                if ~isempty(PSR1{1}) %to avoid errors, added on 12/28 2011
                    %forcing 100ms bin size to keep it consistent across
                    %sessions (no reason is should be different for licks)
                    [PTH1,BW1,~]=MakePTH07(PSR1,repmat(N1, size(RAW(i).Erast{LickInd},1),1),{2,0,0.1});%these values force bin size to be 100ms
                    PTH1=smooth(PTH1,SmoothSPAN,SmoothTYPE)';

                    %------------- Fills the R.Ev(k) fields --------------
                    Lick.Ev(k).BW(LL,1)=BW1;
                    Lick.Ev(k).PSTHraw(LL,1:length(Tm))=PTH1;
                end
            end
        end
    




%% this section finds the ILIs, number of licks, and lick duration


        RD1=strmatch('RD1',RAW(i).Einfo(:,2),'exact'); %finds the index of the event
        RD2=strmatch('RD2',RAW(i).Einfo(:,2),'exact'); %finds the index of the event
        RD1P1=strmatch('RD1P1',RAW(i).Einfo(:,2),'exact'); %finds the index of the event
        RD1P2=strmatch('RD1P2',RAW(i).Einfo(:,2),'exact'); %finds the index of the event
        RD2P1=strmatch('RD2P1',RAW(i).Einfo(:,2),'exact'); %finds the index of the event
        RD2P2=strmatch('RD2P2',RAW(i).Einfo(:,2),'exact'); %finds the index of the event
        Rind=strmatch('BoutEndRLast',RAW(i).Einfo(:,2),'exact');  %finds the index of the event
        Lind=strmatch('Licks',RAW(i).Einfo(:,2),'exact');  %finds the index of the event


        %trying to make the ISI calculator Peter recommended
        %unfortunately, there's one trial of maltodextrin with only one
        %lick
        %so I need to put this lame code in to avoid an error
    %         NL2{1,4}=NaN(40,1);
        %except it doesn't work!

        %find durations of reward consumption for all, R1, and R2 trls
        RDu1=MakePSR04(RAW(i).Erast(Rind),RAW(i).Erast{RD1},[0 15],{2,'first'});
        RDu2=MakePSR04(RAW(i).Erast(Rind),RAW(i).Erast{RD2},[0 15],{2,'first'});
        Lick.RewDurR1(LL,1:length(RDu1))=cell2mat(RDu1);
        Lick.RewDurR2(LL,1:length(RDu2))=cell2mat(RDu2);
        %find number of licks for all, R1, and R2 trls
        NL1=MakePSR04(RAW(i).Erast(Lind),RAW(i).Erast{RD1},[0 15],{2});
        NL2=MakePSR04(RAW(i).Erast(Lind),RAW(i).Erast{RD2},[0 15],{2});
        %find number of licks in 1-3 seconds post-reward
        EL1=MakePSR04(RAW(i).Erast(Lind),RAW(i).Erast{RD1},[earlylicks(1) earlylicks(2)],{2});
        EL2=MakePSR04(RAW(i).Erast(Lind),RAW(i).Erast{RD2},[earlylicks(1) earlylicks(2)],{2});
        %break that down by previous reward
        NL1P1=MakePSR04(RAW(i).Erast(Lind),RAW(i).Erast{RD1P1},[earlylicks(1) earlylicks(2)],{2});
        NL1P2=MakePSR04(RAW(i).Erast(Lind),RAW(i).Erast{RD1P2},[earlylicks(1) earlylicks(2)],{2});
        NL2P1=MakePSR04(RAW(i).Erast(Lind),RAW(i).Erast{RD2P1},[earlylicks(1) earlylicks(2)],{2});
        NL2P2=MakePSR04(RAW(i).Erast(Lind),RAW(i).Erast{RD2P2},[earlylicks(1) earlylicks(2)],{2});

        %make lick matrices
        for j=1:length(NL1)
            Lick.NumLicksR1(LL,j)=sum(NL1{1,j}>0);
            Lick.LicksEarlyR1(LL,j)=sum(EL1{1,j}>0);

            for k=startinterval:endinterval
                if length(NL1{1,j})>k
                    Lick.ILIR1{k,1}(LL,j)= NL1{1,j}(k+1,1)-NL1{1,j}(k,1);
                else
                    Lick.ILIR1{k,1}(LL,j)= NaN;
                end
            end
        end
        for j=1:length(NL2)
            Lick.NumLicksR2(LL,j)=sum(NL2{1,j}>0);
            Lick.LicksEarlyR2(LL,j)=sum(EL2{1,j}>0);

            for k=startinterval:endinterval
                if length(NL2{1,j})>k
                    Lick.ILIR2{k,1}(LL,j)= NL2{1,j}(k+1,1)-NL2{1,j}(k,1);
                else
                    Lick.ILIR2{k,1}(LL,j)= NaN;
                end
            end
        end

        %now for the previous reward situation
        for j=1:length(NL1P1)
            Lick.LicksR1P1(LL,j)=sum(NL1P1{1,j}>0);
        end
        for j=1:length(NL1P2)
            Lick.LicksR1P2(LL,j)=sum(NL1P2{1,j}>0);
        end
        for j=1:length(NL2P1)
            Lick.LicksR2P1(LL,j)=sum(NL2P1{1,j}>0);
        end
        for j=1:length(NL2P2)
            Lick.LicksR2P2(LL,j)=sum(NL2P2{1,j}>0);
        end

        %get session and rat #
        Lick.Sess(LL,1:35)=i;
        A=char(RAW(i).Ninfo(1,1));
        Lick.Rat(LL,1:35)=cellstr(A(1:3));

        fprintf('Session # %d\n',LL);
    end
end

%% computing stats

%setting up result stat matrix
R_2R.BehAvs{2,1}='RewDur';
R_2R.BehAvs{3,1}='NumLicks';
R_2R.BehAvs{4,1}='NumLicksEarly';
R_2R.BehAvs{1,2}='R1 Avg';
R_2R.BehAvs{1,3}='R2 Avg';

%putting all values into one row array and removing NaN, doing 3-way ANOVA
%with fixed effect of reward and random effects of session nested in rat

%reward duration
Lick.RDR1=reshape(Lick.RewDurR1,1,numel(Lick.RewDurR1));
Lick.RDR1Sess=reshape(Lick.Sess,1,numel(Lick.Sess));
Lick.RDR1Rat=reshape(Lick.Rat,1,numel(Lick.Rat));
Lick.RDR1Sess(isnan(Lick.RDR1))=[];
Lick.RDR1Rat(isnan(Lick.RDR1))=[];
Lick.RDR1(isnan(Lick.RDR1))=[];
Lick.RDR2=reshape(Lick.RewDurR2,1,numel(Lick.RewDurR2));
Lick.RDR2Sess=reshape(Lick.Sess,1,numel(Lick.Sess));
Lick.RDR2Rat=reshape(Lick.Rat,1,numel(Lick.Rat));
Lick.RDR2Sess(isnan(Lick.RDR2))=[];
Lick.RDR2Rat(isnan(Lick.RDR2))=[];
Lick.RDR2(isnan(Lick.RDR2))=[];
[~,R_2R.RDstats,~]=anovan(cat(2,Lick.RDR1,Lick.RDR2),{(cat(2,zeros(1,length(Lick.RDR1)),ones(1,length(Lick.RDR2)))),(cat(2,Lick.RDR1Sess,Lick.RDR2Sess)),(cat(2,Lick.RDR1Rat,Lick.RDR2Rat))},'varnames',{'reward','session','rat'},'nested',[0 0 0;0 0 1;0 0 0],'random',[2 3],'display','off','model','full'); %nested ANOVA


%number of licks
Lick.NLR1=reshape(Lick.NumLicksR1,1,numel(Lick.NumLicksR1));
Lick.NLR1Sess=reshape(Lick.Sess,1,numel(Lick.Sess));
Lick.NLR1Rat=reshape(Lick.Rat,1,numel(Lick.Rat));
Lick.NLR1Sess(isnan(Lick.NLR1))=[];
Lick.NLR1Rat(isnan(Lick.NLR1))=[];
Lick.NLR1(isnan(Lick.NLR1))=[];
Lick.NLR2=reshape(Lick.NumLicksR2,1,numel(Lick.NumLicksR2));
Lick.NLR2Sess=reshape(Lick.Sess,1,numel(Lick.Sess));
Lick.NLR2Rat=reshape(Lick.Rat,1,numel(Lick.Rat));
Lick.NLR2Sess(isnan(Lick.NLR2))=[];
Lick.NLR2Rat(isnan(Lick.NLR2))=[];
Lick.NLR2(isnan(Lick.NLR2))=[];
[~,R_2R.LicksStats,~]=anovan(cat(2,Lick.NLR1,Lick.NLR2),{(cat(2,zeros(1,length(Lick.NLR1)),ones(1,length(Lick.NLR2)))),(cat(2,Lick.NLR1Sess,Lick.NLR2Sess)),(cat(2,Lick.NLR1Rat,Lick.NLR2Rat))},'varnames',{'reward','session','rat'},'nested',[0 0 0;0 0 1;0 0 0],'random',[2 3],'display','off','model','full'); %nested ANOVA

%number of licks 1-3sec
Lick.ELR1=reshape(Lick.LicksEarlyR1,1,numel(Lick.LicksEarlyR1));
Lick.ELR1Sess=reshape(Lick.Sess,1,numel(Lick.Sess));
Lick.ELR1Rat=reshape(Lick.Rat,1,numel(Lick.Rat));
Lick.ELR1Sess(isnan(Lick.ELR1))=[];
Lick.ELR1Rat(isnan(Lick.ELR1))=[];
Lick.ELR1(isnan(Lick.ELR1))=[];
Lick.ELR2=reshape(Lick.LicksEarlyR2,1,numel(Lick.LicksEarlyR2));
Lick.ELR2Sess=reshape(Lick.Sess,1,numel(Lick.Sess));
Lick.ELR2Rat=reshape(Lick.Rat,1,numel(Lick.Rat));
Lick.ELR2Sess(isnan(Lick.ELR2))=[];
Lick.ELR2Rat(isnan(Lick.ELR2))=[];
Lick.ELR2(isnan(Lick.ELR2))=[];
[~,R_2R.LicksEarlyStats,~]=anovan(cat(2,Lick.ELR1,Lick.ELR2),{(cat(2,zeros(1,length(Lick.ELR1)),ones(1,length(Lick.ELR2)))),(cat(2,Lick.ELR1Sess,Lick.ELR2Sess)),(cat(2,Lick.ELR1Rat,Lick.ELR2Rat))},'varnames',{'reward','session','rat'},'nested',[0 0 0;0 0 1;0 0 0],'random',[2 3],'display','off','model','full'); %nested ANOVA

%number of licks 1-3sec
Lick.NLR1P1=reshape(Lick.LicksR1P1,1,numel(Lick.LicksR1P1));
Lick.NLR1P1Sess=reshape(Lick.Sess,1,numel(Lick.Sess));
Lick.NLR1P1Rat=reshape(Lick.Rat,1,numel(Lick.Rat));
Lick.NLR1P1Sess(isnan(Lick.NLR1P1))=[];
Lick.NLR1P1Rat(isnan(Lick.NLR1P1))=[];
Lick.NLR1P1(isnan(Lick.NLR1P1))=[];
Lick.NLR1P2=reshape(Lick.LicksR1P2,1,numel(Lick.LicksR1P2));
Lick.NLR1P2Sess=reshape(Lick.Sess,1,numel(Lick.Sess));
Lick.NLR1P2Rat=reshape(Lick.Rat,1,numel(Lick.Rat));
Lick.NLR1P2Sess(isnan(Lick.NLR1P2))=[];
Lick.NLR1P2Rat(isnan(Lick.NLR1P2))=[];
Lick.NLR1P2(isnan(Lick.NLR1P2))=[];
Lick.NLR2P1=reshape(Lick.LicksR2P1,1,numel(Lick.LicksR2P1));
Lick.NLR2P1Sess=reshape(Lick.Sess,1,numel(Lick.Sess));
Lick.NLR2P1Rat=reshape(Lick.Rat,1,numel(Lick.Rat));
Lick.NLR2P1Sess(isnan(Lick.NLR2P1))=[];
Lick.NLR2P1Rat(isnan(Lick.NLR2P1))=[];
Lick.NLR2P1(isnan(Lick.NLR2P1))=[];
Lick.NLR2P2=reshape(Lick.LicksR2P2,1,numel(Lick.LicksR2P2));
Lick.NLR2P2Sess=reshape(Lick.Sess,1,numel(Lick.Sess));
Lick.NLR2P2Rat=reshape(Lick.Rat,1,numel(Lick.Rat));
Lick.NLR2P2Sess(isnan(Lick.NLR2P2))=[];
Lick.NLR2P2Rat(isnan(Lick.NLR2P2))=[];
Lick.NLR2P2(isnan(Lick.NLR2P2))=[];
[~,R_2R.LicksPrevRewStats,stats]=anovan(cat(2,Lick.NLR1P1,Lick.NLR1P2,Lick.NLR2P1,Lick.NLR2P2),{cat(2,zeros(1,length(cat(2,Lick.NLR1P1,Lick.NLR1P2))),ones(1,length(cat(2,Lick.NLR2P1,Lick.NLR2P2)))),...
    cat(2,zeros(1,length(Lick.NLR1P1)),ones(1,length(Lick.NLR1P2)),zeros(1,length(Lick.NLR2P1)),ones(1,length(Lick.NLR2P2))),cat(2,Lick.NLR1P1Sess,Lick.NLR1P2Sess,Lick.NLR2P1Sess,Lick.NLR2P2Sess),...
    cat(2,Lick.NLR1P1Rat,Lick.NLR1P2Rat,Lick.NLR2P1Rat,Lick.NLR2P2Rat)},'varnames',{'current reward','prev reward','session','rat'},'nested',[0 0 0 0;0 0 0 0;0 0 0 1;0 0 0 0],'random',[3 4],'display','off','model','full'); %nested ANOVA

% %without random factor of session
% [~,R_2R.LicksPrevRewStats,stats]=anovan(cat(2,Lick.NLR1P1,Lick.NLR1P2,Lick.NLR2P1,Lick.NLR2P2),{cat(2,zeros(1,length(cat(2,Lick.NLR1P1,Lick.NLR1P2))),ones(1,length(cat(2,Lick.NLR2P1,Lick.NLR2P2)))),...
%     cat(2,zeros(1,length(Lick.NLR1P1)),ones(1,length(Lick.NLR1P2)),zeros(1,length(Lick.NLR2P1)),ones(1,length(Lick.NLR2P2))),...
%     cat(2,Lick.NLR1P1Rat,Lick.NLR1P2Rat,Lick.NLR2P1Rat,Lick.NLR2P2Rat)},'varnames',{'current reward','prev reward','rat'},'random',[3],'display','off','model','full');


%ILI
for m=startinterval:endinterval
    %first get session averages so I can do an ANOVA across all ILI
    Lick.ILIR1sessav{m,1}=nanmean(Lick.ILIR1{m,1},2);
    Lick.ILIR2sessav{m,1}=nanmean(Lick.ILIR2{m,1},2);
    Lick.ILIsesssess{m,1}=Lick.Sess(:,1);
    Lick.ILIsessrat{m,1}=Lick.Rat(:,1);
    
    
    %now do it with all trials
    Lick.ILIR1{m,1}=reshape(Lick.ILIR1{m,1},1,numel(Lick.ILIR1{m,1}));
    Lick.ILIR1Sess{m,1}=reshape(Lick.Sess,1,numel(Lick.ILIR1{m,1}));
    Lick.ILIR1Rat{m,1}=reshape(Lick.Rat,1,numel(Lick.ILIR1{m,1}));
    Lick.ILIR1Sess{m,1}(isnan(Lick.ILIR1{m,1}))=[];
    Lick.ILIR1Rat{m,1}(isnan(Lick.ILIR1{m,1}))=[];
    Lick.ILIR1{m,1}(isnan(Lick.ILIR1{m,1}))=[];
    Lick.ILIR2{m,1}=reshape(Lick.ILIR2{m,1},1,numel(Lick.ILIR2{m,1}));
    Lick.ILIR2Sess{m,1}=reshape(Lick.Sess,1,numel(Lick.ILIR2{m,1}));
    Lick.ILIR2Rat{m,1}=reshape(Lick.Rat,1,numel(Lick.ILIR2{m,1}));
    Lick.ILIR2Sess{m,1}(isnan(Lick.ILIR2{m,1}))=[];
    Lick.ILIR2Rat{m,1}(isnan(Lick.ILIR2{m,1}))=[];
    Lick.ILIR2{m,1}(isnan(Lick.ILIR2{m,1}))=[];
end

%getting effects of interval #, reward, rat, and session on ILI
Lick.ILIR1all=[];
Lick.ILIR1all=[];
Lick.ILIR1Ratall=[];
Lick.ILIR1Sessall=[];
Lick.ILIR1Rewall=[];
Lick.ILIR1Intall=[];
Lick.ILIR2all=[];
Lick.ILIR2all=[];
Lick.ILIR2Ratall=[];
Lick.ILIR2Sessall=[];
Lick.ILIR2Rewall=[];
Lick.ILIR2Intall=[];

%to do session averages for comparing across all intervals
for m=startinterval:endinterval
    Lick.ILIR1all=cat(1,Lick.ILIR1all,Lick.ILIR1sessav{m,1});
    Lick.ILIR1Ratall=cat(1,Lick.ILIR1Ratall,Lick.ILIsessrat{m,1});
    Lick.ILIR1Sessall=cat(1,Lick.ILIR1Sessall,Lick.ILIsesssess{m,1});
    Lick.ILIR1Rewall=cat(1,Lick.ILIR1Rewall,zeros(length(Lick.ILIR1sessav{m,1}),1));
    Lick.ILIR1Intall=cat(1,Lick.ILIR1Intall,(m*ones(length(Lick.ILIR1sessav{m,1}),1)));
    Lick.ILIR2all=cat(1,Lick.ILIR2all,Lick.ILIR2sessav{m,1});
    Lick.ILIR2Ratall=cat(1,Lick.ILIR2Ratall,Lick.ILIsessrat{m,1});
    Lick.ILIR2Sessall=cat(1,Lick.ILIR2Sessall,Lick.ILIsesssess{m,1});
    Lick.ILIR2Rewall=cat(1,Lick.ILIR2Rewall,ones(length(Lick.ILIR2sessav{m,1}),1));
    Lick.ILIR2Intall=cat(1,Lick.ILIR2Intall,(m*ones(length(Lick.ILIR2sessav{m,1}),1)));    
end

%no factor of session
[~,R_2R.ILItrialstats,~]=anovan(cat(1,Lick.ILIR1all,Lick.ILIR2all),{(cat(1,Lick.ILIR1Rewall,Lick.ILIR2Rewall)),(cat(1,Lick.ILIR1Ratall,Lick.ILIR2Ratall)),(cat(1,Lick.ILIR1Intall,Lick.ILIR2Intall))},'varnames',{'reward','rat','interval'},'random',[2],'display','off','model','full');

%finding values for R_2R.BehAvs
R_2R.BehAvs{2,2}=mean(Lick.RDR1);
R_2R.BehAvs{2,3}=mean(Lick.RDR2);
R_2R.BehAvs{3,2}=mean(Lick.NLR1);
R_2R.BehAvs{3,3}=mean(Lick.NLR2);
R_2R.BehAvs{4,2}=mean(Lick.ELR1);
R_2R.BehAvs{4,3}=mean(Lick.ELR2);

%Finding ILI
for m=startinterval:endinterval
    Lick.ILI1(m,1)=nanmean(Lick.ILIR1{m,1}(:));
    Lick.ILI1(m,2)=nanstd(Lick.ILIR1{m,1}(:))/sqrt(sum(~isnan(Lick.ILIR1{m,1}(:))));
    Lick.ILI2(m,1)=nanmean(Lick.ILIR2{m,1}(:));
    Lick.ILI2(m,2)=nanstd(Lick.ILIR2{m,1}(:))/sqrt(sum(~isnan(Lick.ILIR2{m,1}(:))));
end

%% plotting licks in each reward condition
Xaxis=[-2 12];
Ishow=find(Tm>=Xaxis(1) & Tm<=Xaxis(2));
time1=Tm(Ishow);
c=[-1000 7500];ClimE=sign(c).*abs(c).^(1/4);%ColorMapExc
colormap(jet);
sucrose=[1  0.6  0.1];
maltodextrin=[.9  0.3  .9];
water=[0.00  0.75  0.75];
total=[0.3  0.1  0.8];
exc=[0 113/255 188/255];
inh=[240/255 0 50/255];

Ev1=strcmp('RD1', Lick.Erefnames);
Ev2=strcmp('RD2', Lick.Erefnames);

psth1=nanmean(Lick.Ev(Ev1).PSTHraw(:,Ishow),1);
sem1=nanste(Lick.Ev(Ev1).PSTHraw(:,Ishow),1); %calculate standard error of the mean
up1=psth1+sem1;
down1=psth1-sem1;

psthE=nanmean(Lick.Ev(Ev2).PSTHraw(:,Ishow),1);
semE=nanste(Lick.Ev(Ev2).PSTHraw(:,Ishow),1); %calculate standard error of the mean
upE=psthE+semE;
downE=psthE-semE;

%plotting
subplot(2,2,1);
hold on;
plot(time1,psth1,'Color',sucrose,'linewidth',1);
plot(time1,psthE,'Color',maltodextrin,'linewidth',1);
patch([time1,time1(end:-1:1)],[up1,down1(end:-1:1)],sucrose,'EdgeColor','none');alpha(0.5);
patch([time1,time1(end:-1:1)],[upE,downE(end:-1:1)],maltodextrin,'EdgeColor','none');alpha(0.5);
title('Mean lick rate');
%plot([-2 10],[0 0],':','color','k');
plot([0 0],[-2 8],':','color','k','linewidth',0.75);
plot([earlylicks(1) earlylicks(2)],[7.3 7.3],'color','k','linewidth',0.75);
axis([-2 12 0 7.5]);
xlabel('Seconds from reward delivery');
ylabel('Licks/s');
legend('Sucrose trials','Maltodextrin trials');


%% plotting ILI
sucrose=[1  0.6  0.1];
maltodextrin=[.9  0.3  .9];
water=[0.00  0.75  0.75];
total=[0.3  0.1  0.8];
exc=[0 113/255 188/255];
inh=[240/255 0 50/255];
xaxis=(startinterval:endinterval);

subplot(2,2,3);
hold on;
errorbar(xaxis,Lick.ILI1(startinterval:endinterval,1),Lick.ILI1(startinterval:endinterval,2),'*','Color',sucrose,'linewidth',1.5);
errorbar(xaxis,Lick.ILI2(startinterval:endinterval,1),Lick.ILI2(startinterval:endinterval,2),'*','Color',maltodextrin,'linewidth',1.5);
xlabel('Lick Interval #');
ylabel('Duration (s)');
title(['Interlick intervals'])
axis([1 30 0.13 0.19]);
legend('Sucrose trials','Maltodextrin trials','location','northwest');

subplot(4,4,15);
hold on;
errorbar(2,nanmean(Lick.ILI1(startinterval:endinterval,1)),nanste(Lick.ILI1(startinterval:endinterval,1),1),'*','Color',sucrose,'linewidth',1.5);
errorbar(2,nanmean(Lick.ILI2(startinterval:endinterval,1)),nanste(Lick.ILI2(startinterval:endinterval,1),1),'*','Color',maltodextrin,'linewidth',1.5);
if R_2R.LicksEarlyStats{2,7}<0.05
    plot(2,nanmean(Lick.ILI2(startinterval:endinterval,1))+0.005,'*','Color','k','MarkerSize',5);
end
set(gca,'xtick',[]);
axis([1.5 2.5 0.14 0.17]);


%% plotting sucrose preferences

%Plot sucrose preference

%colors
NAc=[0.5  0.1  0.8];
VP=[0.3  0.7  0.1];

%data
x=[1 2];
NAcPref{1}=[81.7 81.3];
NAcPref{2}=[77 73.9];
NAcPref{3}=[37.4 74.6];
NAcPref{4}=[71 95.6];
NAcPref{5}=[89 83.6];
NAcPref{6}=[94 93];
VPPref{1}=[59.6 78.8];
VPPref{2}=[88.8 85.8];
VPPref{3}=[71.2 89];
VPPref{4}=[87 76.6];
VPPref{5}=[85.2 84.1];

subplot(2,4,3);
hold on;
%do the first 2 first in order to get the legend right
plot(x,NAcPref{1},'Marker','o','MarkerFaceColor',NAc,'LineWidth',1,'Color',NAc);
plot(x,VPPref{1},'Marker','o','MarkerFaceColor',VP,'LineWidth',1,'Color',VP);

for i=2:length(NAcPref)
    plot(x,NAcPref{i},'Marker','o','MarkerFaceColor',NAc,'LineWidth',1,'Color',NAc);
end

for i=2:length(VPPref)
    plot(x,VPPref{i},'Marker','o','MarkerFaceColor',VP,'LineWidth',1,'Color',VP);
end

axis([0.7 2.2 0 100]);
title('Two-bottle choice');
xlabel('Initial       Final');
ylabel('% sucrose consumption');
legend('NAc rat','VP rat','location','southeast');
plot([0 3],[50 50],':','color','k','linewidth',0.75);
set(gca,'xtick',[])


%% licks for each rat
figure;

Xaxis=[-2 12];
Ishow=find(Tm>=Xaxis(1) & Tm<=Xaxis(2));
time1=Tm(Ishow);

%colors
sucrose=[1  0.6  0.1];
maltodextrin=[.9  0.3  .9];
water=[0.00  0.75  0.75];
total=[0.3  0.1  0.8];
inh=[0.1 0.021154 0.6];
exc=[0.9 0.75 0.205816];
NAc=[0.5  0.1  0.8];
VP=[0.3  0.7  0.1];

%split data up into individual rats
for i=1:length(Lick.Linfo)
    A=char(Lick.Linfo(i,1));
    Rats(i,:)=cellstr(A(1:3));
end

C=unique(Rats);

for i=1:length(C)
    Sel=strcmp(C(i),Rats);
    Ev1=strcmp('RD1', Lick.Erefnames);
    Ev2=strcmp('RD2', Lick.Erefnames);

    psth1=nanmean(Lick.Ev(Ev1).PSTHraw(Sel,Ishow),1);
    sem1=nanste(Lick.Ev(Ev1).PSTHraw(Sel,Ishow),1); %calculate standard error of the mean
    up1=psth1+sem1;
    down1=psth1-sem1;

    psthE=nanmean(Lick.Ev(Ev2).PSTHraw(Sel,Ishow),1);
    semE=nanste(Lick.Ev(Ev2).PSTHraw(Sel,Ishow),1); %calculate standard error of the mean
    upE=psthE+semE;
    downE=psthE-semE;

    %plotting
    subplot(4,3,i);
    hold on;
    plot(time1,psth1,'Color',sucrose,'linewidth',1);
    plot(time1,psthE,'Color',maltodextrin,'linewidth',1);
    patch([time1,time1(end:-1:1)],[up1,down1(end:-1:1)],sucrose,'EdgeColor','none');alpha(0.5);
    patch([time1,time1(end:-1:1)],[upE,downE(end:-1:1)],maltodextrin,'EdgeColor','none');alpha(0.5);
    title([cell2mat(C(i)) ' (n=' num2str(sum(Sel)) ')' ]);
    %plot([-2 10],[0 0],':','color','k');
    plot([0 0],[-2 8],':','color','k','linewidth',0.75);
    patch([earlylicks(1) earlylicks(1) earlylicks(2) earlylicks(2)],[5 8 8 5],'k','facecolor','none','edgecolor',[0.7 0.7 0.7]);alpha(0.5);
    axis([-2 12 0 8]);
    xlabel('Seconds from reward delivery');
    ylabel('Licks/s');
    if i==1 legend('Sucrose trials','Maltodextrin trials'); end

end

save('R_2R.mat','R_2R');