%% initial analysis of event-evoked activity

%need to save excel sheets in the specified folder
%the excel sheets say which events to make PSTHs for and what windows to
%use for statistical testing

%clear all; clc;
global Dura Baseline Tm Tbase BSIZE Tbin
tic

if ~exist('RAW'), load ('RAW.mat'); end

%Main settings
SAVE_FLAG=1;
BSIZE=0.01; %Do not change
Dura=[-5 12]; Tm=Dura(1):BSIZE:Dura(2);
%Baseline=[-22 0]; Tbase=Baseline(1):BSIZE:Baseline(2); %now defined line 98
Tbin=-0.5:0.005:0.5; %window used to determine the optimal binsize
PStat=0.05; %for comparing pre versus post windows, or event A versus event B
MinNumTrials=5; %how many trials of event there need to be to conduct analysis
BinSize=0.05;

%Smoothing
Smoothing=1; %0 for raw and 1 for smoothing
SmoothTYPE='lowess';
SmoothSPAN=5;
if Smoothing~=1, SmoothTYPE='NoSmoothing';SmoothSPAN=NaN; end

%for bins analysis
BinBase=[-22 -12]; %For normalizing activity
BinDura=[0 0.6]; %size of bin
bins = 66; %number of bins
binint = 0.1; %spacing of bins
binstart = -2; %start time of first bin relative to event

%%  Standard sucrose and maltodextrin sessions

%start fresh
R_2R=[];R_2R.Ninfo={};NN=0;Nneurons=0;

% List of events to analyze and analysis windows EXTRACTED from excel file
path='C:\Users\dottenh2\Documents\MATLAB\David\2Rewards Nex Files\2R_paper.xls';
[~,Erefnames]=xlsread(path,'Windows','a3:a15'); % cell that contains the event names
prewin =  xlsread(path,'Windows','b3:c15');
postwin = xlsread(path,'Windows','d3:e15');

%saves event names for reference later
R_2R.Erefnames=Erefnames;

%Finds the total number of neurons in 2R and marks them by region/session
for i=1:length(RAW)
    if strcmp('NA',RAW(i).Type(1:2)) | strcmp('VP',RAW(i).Type(1:2))
        R_2R.Ninfo=cat(1,R_2R.Ninfo,RAW(i).Ninfo);
        Nneurons=Nneurons+size(RAW(i).Nrast,1);
    end
end
for i=1:Nneurons
    Session=string(R_2R.Ninfo(i,1));
    Name=char(Session);
    R_2R.Ninfo(i,3)=cellstr(Name(1:2));
    R_2R.Ninfo(i,4)=cellstr(Name(1:3));
end

% preallocating
R_2R.Param.Tm=Tm;
R_2R.Param.Tbin=Tbin;
R_2R.Param.Dura=Dura;
R_2R.Param.Baseline=Baseline;
R_2R.Param.PStat=PStat;
R_2R.Param.MinNumTrials=MinNumTrials;
R_2R.Param.path=path;
R_2R.Param.prewin=prewin;
R_2R.Param.postwin=postwin;
R_2R.Param.SmoothTYPE=SmoothTYPE;
R_2R.Param.SmoothSPAN=SmoothSPAN;
R_2R.Param.BinBase=BinBase;
R_2R.Param.BinDura=BinDura;
R_2R.Param.bins = bins;
R_2R.Param.binint = binint;
R_2R.Param.binstart = binstart;

for k=1:length(Erefnames)
    R_2R.Ev(k).PSTHraw(1:Nneurons,1:length(Tm))=NaN(Nneurons,length(Tm));
    R_2R.Ev(k).PSTHz(1:Nneurons,1:length(Tm))=NaN(Nneurons,length(Tm));
    R_2R.Ev(k).Meanraw(1:Nneurons,1)=NaN;
    R_2R.Ev(k).Meanz(1:Nneurons,1)=NaN;
    R_2R.Ev(k).BW(1:Nneurons,1)=NaN;
    R_2R.Ev(k).signrank(1:Nneurons,1)=NaN;
    R_2R.Ev(k).RespDir(1:Nneurons,1)=NaN;
    R_2R.Ev(k).NumberTrials(1:Nneurons,1)=NaN;
end

% runs the main routine

for i=1:length(RAW) %loops through sessions
    if strcmp('NA',RAW(i).Type(1:2)) | strcmp('VP',RAW(i).Type(1:2))
        for j= 1:size(RAW(i).Nrast,1) %Number of neurons per session
            NN=NN+1; %neuron 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

                if sum(EvInd)==0
                    fprintf('HOWDY, CANT FIND EVENTS FOR ''%s''\n',Erefnames{k});
                end

                R_2R.Ev(k).NumberTrials(NN,1)=length(RAW(i).Erast{EvInd});
                Tbase=prewin(k,1):BSIZE:prewin(k,2);
                if  ~isempty(EvInd) && R_2R.Ev(k).NumberTrials(NN,1)>=MinNumTrials %avoid analyzing sessions where that do not have enough trials
                    [PSR0,N0]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},prewin(k,:),{1});% makes collapsed rasters for baseline for use in normalizing
                    [PSR1,N1]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},Dura,{1});% makes collpased rasters. PSR1 is a cell(neurons)
                    [PSR2,N2]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},Dura,{2});% makes trial by trial rasters. PSR1 is a cell(neurons, trials)
                    [PSR3,N0]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},prewin(k,:),{2});% makes trial by trial rasters for baseline for use in response detection
                    if ~isempty(PSR0{1}) || ~isempty(PSR1{1}) %to avoid errors, added on 12/28 2011
                        %Optimal bin size
                        %[PTH1,BW1,~]=MakePTH07(PSR1,repmat(N1, size(RAW(i).Nrast{j},1),1),{2,1});%-----DP used here
                        %[PTH0,~,~]=MakePTH07(PSR0,repmat(N0, size(RAW(i).Nrast{j},1),1),{1,0,BW1});%BW1 reinjected here to make sure PTH0 & PTH1 have the same BW
                        %Fixed bin size
                        [PTH1,BW1,~]=MakePTH07(PSR1,repmat(N1, size(RAW(i).Nrast{j},1),1),{2,0,BinSize});%-----Fixed bin used here
                        [PTH0,~,~]=MakePTH07(PSR0,repmat(N0, size(RAW(i).Nrast{j},1),1),{1,0,BinSize});%BW1 reinjected here to make sure PTH0 & PTH1 have the same BW
                        PTH1=smooth(PTH1,SmoothSPAN,SmoothTYPE)';
                        PTH0=smooth(PTH0,SmoothSPAN,SmoothTYPE)';

                        %------------- Fills the R_2R.Ev(k) fields --------------
                        R_2R.Ev(k).BW(NN,1)=BW1;
                        R_2R.Ev(k).PSTHraw(NN,1:length(Tm))=PTH1;
                        R_2R.Ev(k).Meanraw(NN,1)=nanmean(R_2R.Ev(k).PSTHraw(NN,Tm>postwin(k,1) & Tm<postwin(k,2)),2);
                        if sum(PTH0,2)~=0
                            R_2R.Ev(k).PSTHz(NN,1:length(Tm))=normalize(PTH1,PTH0,0);
                            R_2R.Ev(k).Meanz(NN,1)=nanmean(R_2R.Ev(k).PSTHz(NN,Tm>postwin(k,1) & Tm<postwin(k,2)),2);
                        else
                            R_2R.Ev(k).PSTHz(NN,1:length(Tm))=NaN(1,length(Tm));
                            R_2R.Ev(k).Meanz(NN,1)=NaN;
                        end

                        %------------------ firing (in Hz) per trial in pre/post windows ------------------
                        %used to make the between events comparisons and Response detection in a single window----
                        ev(k).pre=NaN(size(RAW(i).Erast{EvInd},1),1);
                        ev(k).post=NaN(size(RAW(i).Erast{EvInd},1),1);
                        for m=1:size(RAW(i).Erast{EvInd},1) %loops through trials
                            ev(k).pre(m)=sum(PSR3{m}<prewin(k,2) & PSR3{m}>prewin(k,1))/(prewin(k,2)-prewin(k,1)); %CHANGED FROM PSR2 to PSR0 here 10/24/17
                            ev(k).post(m)=sum(PSR2{m}<postwin(k,2) & PSR2{m}>postwin(k,1))/(postwin(k,2)-postwin(k,1));
                        end

                        %-------------------- signrank on event and direction----
                        [R_2R.Ev(k).signrank(NN,1),~]=signrank(ev(k).pre, ev(k).post); %Signrank used here because it is a dependant sample test
                        if R_2R.Ev(k).signrank(NN,1)<PStat
                            R_2R.Ev(k).RespDir(NN,1)=sign(mean(ev(k).post)-mean(ev(k).pre));
                        else R_2R.Ev(k).RespDir(NN,1)=0;
                        end

                    end %if ~isempty(PSR0{1}) || ~isempty(PSR1{1})
                end %if EvInd=0 OR n(trials) < MinNumTrials fills with NaN
            end %Events


            %----------------------- Check for difference in cue responding after sucrose and maltodextrin trials -----------------------         
            CueP1=strcmp('CueP1', Erefnames);
            CueP2=strcmp('CueP2', Erefnames);
            [R_2R.CueP12Stat(NN,1),~]=ranksum(ev(CueP1).post,ev(CueP2).post); %Ranksum test used becasue it is an independent sample test


            %----------------------- Check for difference in PE responding after sucrose and maltodextrin trials -----------------------         
            PEP1=strcmp('PEP1', Erefnames);
            PEP2=strcmp('PEP2', Erefnames);
            [R_2R.PEP12Stat(NN,1),~]=ranksum(ev(PEP1).post,ev(PEP1).post); %Ranksum test used becasue it is an independent sample test

            %----------------------- Check for reward modulation across all the bins surrounding reward delivery-----------------------
            BHz=0;
            Trial=0;
            Reward=0;

            EV1=strmatch('RD1',RAW(i).Einfo(:,2),'exact');
            EV2=strmatch('RD2',RAW(i).Einfo(:,2),'exact');

            %get mean baseline firing for all reward trials
            [Bcell1,B1n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV1},BinBase,{2});% makes trial by trial rasters for baseline
            for y= 1:B1n
                BHz(1,y)=sum(Bcell1{1,y}>BinBase(1))/(BinBase(1,2)-BinBase(1,1));
                Reward(1,y)=1;
            end
            [Bcell2,B2n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV2},BinBase,{2});% makes trial by trial rasters for baseline
            for z= 1:B2n
                BHz(1,z+B1n)=sum(Bcell2{1,z}>BinBase(1))/(BinBase(1,2)-BinBase(1,1));
                Reward(1,z+B1n)=2;
            end

            for k=1:bins %go through each bin for each neuron
                EVHz=0;
                %get trial by trial firing rate for sucrose trials
                [EV1cell,EV1n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV1},[(BinDura(1)+(binstart - binint)+k*binint) (BinDura(2)+(binstart - binint)+k*binint)],{2});% makes trial by trial rasters for event
                for y= 1:EV1n
                    EVHz(1,y)=sum(EV1cell{1,y}>binstart)/(BinDura(1,2)-BinDura(1,1));
                end       

                %get trial by trial firing rate for maltodextrin trials
                [EV2cell,EV2n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV2},[(BinDura(1)+(binstart - binint)+k*binint) (BinDura(2)+(binstart - binint)+k*binint)],{2});% makes trial by trial rasters for event
                for z= 1:EV2n
                    EVHz(1,z+EV1n)=sum(EV2cell{1,z}>binstart)/(BinDura(1,2)-BinDura(1,1));
                end       


                %run ANOVA on the effects of pre vs post, cued or not,
                %kind of reward, and the interactions
                [~,Bintbl,~]=anovan([BHz EVHz],{cat(2,zeros(1,length(Reward)),ones(1,length(Reward))),[Reward Reward]},'varnames',{'Pre vs Post','Reward'},'model','full','display','off');
                R_2R.BinStatRwrd{k,1}(NN).IntSig=cell2mat(Bintbl(4,7));

                %get mean firing rates
                R_2R.BinStatRwrd{k,1}(NN).R1Mean=(nanmean(EVHz(Reward==1))-nanmean(BHz(Reward==1)))/std(BHz(Reward==1));
                R_2R.BinStatRwrd{k,1}(NN).R2Mean=(nanmean(EVHz(Reward==2))-nanmean(BHz(Reward==2)))/std(BHz(Reward==2));

                %signranks for checking if inhibitory or excitatory

                %sucrose trials
                signranksig=signrank(BHz(Reward==1), EVHz(Reward==1)); %Signrank used here because it is a dependant sample test
                        if signranksig<PStat
                            R_2R.BinStatRwrd{k,1}(NN).SucRespDir=sign(nanmean(EVHz(Reward==1))-mean(BHz(Reward==1)));
                        else
                            R_2R.BinStatRwrd{k,1}(NN).SucRespDir=0;
                        end

                %maltodextrin trials
                signranksig=signrank(BHz(Reward==2), EVHz(Reward==2)); %Signrank used here because it is a dependant sample test
                        if signranksig<PStat
                            R_2R.BinStatRwrd{k,1}(NN).MalRespDir=sign(nanmean(EVHz(Reward==2))-mean(BHz(Reward==2)));
                        else
                            R_2R.BinStatRwrd{k,1}(NN).MalRespDir=0;
                        end

            end %bins for bins analysis


            fprintf('2R Neuron #%d\n',NN);
        end %neurons: FOR j= 1:size(RAW(i).Nrast,1)
    end %if it's the right kind of session
end %sessions: FOR i=1:length(RAW)

if SAVE_FLAG
    save('R_2R.mat','R_2R')
end

toc


%%  water and maltodextrin sessions

%start fresh
R_W=[];R_W.Ninfo={};NN=0;Nneurons=0;

% List of events to analyze and analysis windows EXTRACTED from excel file
path='C:\Users\dottenh2\Documents\MATLAB\David\2Rewards Nex Files\2R_paper.xls';
[~,Erefnames]=xlsread(path,'Windows','a3:a15'); % cell that contains the event names
prewin =  xlsread(path,'Windows','b3:c15');
postwin = xlsread(path,'Windows','d3:e15');

%saves event names for reference later
R_W.Erefnames=Erefnames;

%Finds the total number of neurons in 2R and marks them by region/session
for i=1:length(RAW)
    if strcmp('WA',RAW(i).Type(1:2))
        R_W.Ninfo=cat(1,R_W.Ninfo,RAW(i).Ninfo);
        Nneurons=Nneurons+size(RAW(i).Nrast,1);
    end
end
for i=1:Nneurons
    Session=string(R_W.Ninfo(i,1));
    Name=char(Session);
    R_W.Ninfo(i,3)=cellstr(Name(1:2));
    R_W.Ninfo(i,4)=cellstr(Name(1:3));
end

% preallocating
R_W.Param.Tm=Tm;
R_W.Param.Tbin=Tbin;
R_W.Param.Dura=Dura;
R_W.Param.Baseline=Baseline;
R_W.Param.PStat=PStat;
R_W.Param.MinNumTrials=MinNumTrials;
R_W.Param.path=path;
R_W.Param.prewin=prewin;
R_W.Param.postwin=postwin;
R_W.Param.SmoothTYPE=SmoothTYPE;
R_W.Param.SmoothSPAN=SmoothSPAN;
R_W.Param.BinBase=BinBase;
R_W.Param.BinDura=BinDura;
R_W.Param.bins = bins;
R_W.Param.binint = binint;
R_W.Param.binstart = binstart;

for k=1:length(Erefnames)
    R_W.Ev(k).PSTHraw(1:Nneurons,1:length(Tm))=NaN(Nneurons,length(Tm));
    R_W.Ev(k).PSTHz(1:Nneurons,1:length(Tm))=NaN(Nneurons,length(Tm));
    R_W.Ev(k).Meanraw(1:Nneurons,1)=NaN;
    R_W.Ev(k).Meanz(1:Nneurons,1)=NaN;
    R_W.Ev(k).BW(1:Nneurons,1)=NaN;
    R_W.Ev(k).signrank(1:Nneurons,1)=NaN;
    R_W.Ev(k).RespDir(1:Nneurons,1)=NaN;
    R_W.Ev(k).NumberTrials(1:Nneurons,1)=NaN;
end

% runs the main routine

for i=1:length(RAW) %loops through sessions
    if strcmp('WA',RAW(i).Type(1:2))
        for j= 1:size(RAW(i).Nrast,1) %Number of neurons per session
            NN=NN+1; %neuron 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

                if sum(EvInd)==0
                    fprintf('HOWDY, CANT FIND EVENTS FOR ''%s''\n',Erefnames{k});
                end

                R_W.Ev(k).NumberTrials(NN,1)=length(RAW(i).Erast{EvInd});
                Tbase=prewin(k,1):BSIZE:prewin(k,2);
                if  ~isempty(EvInd) && R_W.Ev(k).NumberTrials(NN,1)>=MinNumTrials %avoid analyzing sessions where that do not have enough trials
                    [PSR0,N0]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},prewin(k,:),{1});% makes collapsed rasters for baseline for use in normalizing
                    [PSR1,N1]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},Dura,{1});% makes collpased rasters. PSR1 is a cell(neurons)
                    [PSR2,N2]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},Dura,{2});% makes trial by trial rasters. PSR1 is a cell(neurons, trials)
                    [PSR3,N0]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},prewin(k,:),{2});% makes trial by trial rasters for baseline for use in response detection
                    if ~isempty(PSR0{1}) || ~isempty(PSR1{1}) %to avoid errors, added on 12/28 2011
                        %optimal bin size
                        %[PTH1,BW1,~]=MakePTH07(PSR1,repmat(N1, size(RAW(i).Nrast{j},1),1),{2,1});%-----DP used here
                        %[PTH0,~,~]=MakePTH07(PSR0,repmat(N0, size(RAW(i).Nrast{j},1),1),{1,0,BW1});%BW1 reinjected here to make sure PTH0 & PTH1 have the same BW
                        %Fixed bin size
                        [PTH1,BW1,~]=MakePTH07(PSR1,repmat(N1, size(RAW(i).Nrast{j},1),1),{2,0,BinSize});%-----fixed bin size
                        [PTH0,~,~]=MakePTH07(PSR0,repmat(N0, size(RAW(i).Nrast{j},1),1),{1,0,BinSize});%BW1 reinjected here to make sure PTH0 & PTH1 have the same BW
                        PTH1=smooth(PTH1,SmoothSPAN,SmoothTYPE)';
                        PTH0=smooth(PTH0,SmoothSPAN,SmoothTYPE)';

                        %------------- Fills the R_W.Ev(k) fields --------------
                        R_W.Ev(k).BW(NN,1)=BW1;
                        R_W.Ev(k).PSTHraw(NN,1:length(Tm))=PTH1;
                        R_W.Ev(k).Meanraw(NN,1)=nanmean(R_W.Ev(k).PSTHraw(NN,Tm>postwin(k,1) & Tm<postwin(k,2)),2);
                        if sum(PTH0,2)~=0
                            R_W.Ev(k).PSTHz(NN,1:length(Tm))=normalize(PTH1,PTH0,0);
                            R_W.Ev(k).Meanz(NN,1)=nanmean(R_W.Ev(k).PSTHz(NN,Tm>postwin(k,1) & Tm<postwin(k,2)),2);
                        else
                            R_W.Ev(k).PSTHz(NN,1:length(Tm))=NaN(1,length(Tm));
                            R_W.Ev(k).Meanz(NN,1)=NaN;
                        end

                        %------------------ firing (in Hz) per trial in pre/post windows ------------------
                        %used to make the between events comparisons and Response detection in a single window----
                        ev(k).pre=NaN(size(RAW(i).Erast{EvInd},1),1);
                        ev(k).post=NaN(size(RAW(i).Erast{EvInd},1),1);
                        for m=1:size(RAW(i).Erast{EvInd},1) %loops through trials
                            ev(k).pre(m)=sum(PSR3{m}<prewin(k,2) & PSR3{m}>prewin(k,1))/(prewin(k,2)-prewin(k,1)); %CHANGED FROM PSR2 to PSR0 here 10/24/17
                            ev(k).post(m)=sum(PSR2{m}<postwin(k,2) & PSR2{m}>postwin(k,1))/(postwin(k,2)-postwin(k,1));
                        end

                        %-------------------- signrank on event and direction----
                        [R_W.Ev(k).signrank(NN,1),~]=signrank(ev(k).pre, ev(k).post); %Signrank used here because it is a dependant sample test
                        if R_W.Ev(k).signrank(NN,1)<PStat
                            R_W.Ev(k).RespDir(NN,1)=sign(mean(ev(k).post)-mean(ev(k).pre));
                        else R_W.Ev(k).RespDir(NN,1)=0;
                        end

                    end %if ~isempty(PSR0{1}) || ~isempty(PSR1{1})
                end %if EvInd=0 OR n(trials) < MinNumTrials fills with NaN
            end %Events


            %----------------------- Check for reward modulation across all the bins surrounding reward delivery-----------------------
            BHz=0;
            Trial=0;
            Reward=0;

            EV1=strmatch('RD1',RAW(i).Einfo(:,2),'exact');
            EV2=strmatch('RD2',RAW(i).Einfo(:,2),'exact');

            %get mean baseline firing for all reward trials
            [Bcell1,B1n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV1},BinBase,{2});% makes trial by trial rasters for baseline
            for y= 1:B1n
                BHz(1,y)=sum(Bcell1{1,y}>BinBase(1))/(BinBase(1,2)-BinBase(1,1));
                Reward(1,y)=1;
            end
            [Bcell2,B2n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV2},BinBase,{2});% makes trial by trial rasters for baseline
            for z= 1:B2n
                BHz(1,z+B1n)=sum(Bcell2{1,z}>BinBase(1))/(BinBase(1,2)-BinBase(1,1));
                Reward(1,z+B1n)=2;
            end

            for k=1:bins %go through each bin for each neuron
                EVHz=0;
                %get trial by trial firing rate for water trials
                [EV1cell,EV1n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV1},[(BinDura(1)+(binstart - binint)+k*binint) (BinDura(2)+(binstart - binint)+k*binint)],{2});% makes trial by trial rasters for event
                for y= 1:EV1n
                    EVHz(1,y)=sum(EV1cell{1,y}>binstart)/(BinDura(1,2)-BinDura(1,1));
                end       

                %get trial by trial firing rate for maltodextrin trials
                [EV2cell,EV2n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV2},[(BinDura(1)+(binstart - binint)+k*binint) (BinDura(2)+(binstart - binint)+k*binint)],{2});% makes trial by trial rasters for event
                for z= 1:EV2n
                    EVHz(1,z+EV1n)=sum(EV2cell{1,z}>binstart)/(BinDura(1,2)-BinDura(1,1));
                end       


                %run ANOVA on the effects of pre vs post, cued or not,
                %kind of reward, and the interactions
                [~,Bintbl,~]=anovan([BHz EVHz],{cat(2,zeros(1,length(Reward)),ones(1,length(Reward))),[Reward Reward]},'varnames',{'Pre vs Post','Reward'},'model','full','display','off');
                R_W.BinStatRwrd{k,1}(NN).IntSig=cell2mat(Bintbl(4,7));

                %get mean firing rates
                R_W.BinStatRwrd{k,1}(NN).R1Mean=(nanmean(EVHz(Reward==1))-nanmean(BHz(Reward==1)))/std(BHz(Reward==1));
                R_W.BinStatRwrd{k,1}(NN).R2Mean=(nanmean(EVHz(Reward==2))-nanmean(BHz(Reward==2)))/std(BHz(Reward==2));

                %signranks for checking if inhibitory or excitatory

                %water trials
                signranksig=signrank(BHz(Reward==1), EVHz(Reward==1)); %Signrank used here because it is a dependant sample test
                if signranksig<PStat
                    R_W.BinStatRwrd{k,1}(NN).WatRespDir=sign(nanmean(EVHz(Reward==1))-mean(BHz(Reward==1)));
                else
                    R_W.BinStatRwrd{k,1}(NN).WatRespDir=0;
                end

                %maltodextrin trials
                signranksig=signrank(BHz(Reward==2), EVHz(Reward==2)); %Signrank used here because it is a dependant sample test
                if signranksig<PStat
                    R_W.BinStatRwrd{k,1}(NN).MalRespDir=sign(nanmean(EVHz(Reward==2))-mean(BHz(Reward==2)));
                else
                    R_W.BinStatRwrd{k,1}(NN).MalRespDir=0;
                end

            end %bins for bins analysis


            fprintf('Water neuron #%d\n',NN);
        end %neurons: FOR j= 1:size(RAW(i).Nrast,1)
    end %if right kind of session
end %sessions: FOR i=1:length(RAW)

if SAVE_FLAG
    save('R_W.mat','R_W')
end

toc

%%  All three rewards sessions

%make this lower because there are only 3 of some current/previous reward combinations
MinNumTrials=3; %how many trials of event there need to be to conduct analysis


%start fresh
R_3R=[];R_3R.Ninfo={};NN=0;Nneurons=0;

% List of events to analyze and analysis windows EXTRACTED from excel file
path='C:\Users\dottenh2\Documents\MATLAB\David\2Rewards Nex Files\3R_paper.xls';
[~,Erefnames]=xlsread(path,'Windows','a3:a23'); % cell that contains the event names
prewin =  xlsread(path,'Windows','b3:c23');
postwin = xlsread(path,'Windows','d3:e23');

%saves event names for reference later
R_3R.Erefnames=Erefnames;

%Finds the total number of neurons in 2R and marks them by region/session
for i=1:length(RAW)
    if strcmp('TH',RAW(i).Type(1:2))
        R_3R.Ninfo=cat(1,R_3R.Ninfo,RAW(i).Ninfo);
        Nneurons=Nneurons+size(RAW(i).Nrast,1);
    end
end
for i=1:Nneurons
    Session=string(R_3R.Ninfo(i,1));
    Name=char(Session);
    R_3R.Ninfo(i,3)=cellstr(Name(1:2));
    R_3R.Ninfo(i,4)=cellstr(Name(1:3));
end

% preallocating
R_3R.Param.Tm=Tm;
R_3R.Param.Tbin=Tbin;
R_3R.Param.Dura=Dura;
R_3R.Param.Baseline=Baseline;
R_3R.Param.PStat=PStat;
R_3R.Param.MinNumTrials=MinNumTrials;
R_3R.Param.path=path;
R_3R.Param.prewin=prewin;
R_3R.Param.postwin=postwin;
R_3R.Param.SmoothTYPE=SmoothTYPE;
R_3R.Param.SmoothSPAN=SmoothSPAN;
R_3R.Param.BinBase=BinBase;
R_3R.Param.BinDura=BinDura;
R_3R.Param.bins = bins;
R_3R.Param.binint = binint;
R_3R.Param.binstart = binstart;

for k=1:length(Erefnames)
    R_3R.Ev(k).PSTHraw(1:Nneurons,1:length(Tm))=NaN(Nneurons,length(Tm));
    R_3R.Ev(k).PSTHz(1:Nneurons,1:length(Tm))=NaN(Nneurons,length(Tm));
    R_3R.Ev(k).Meanraw(1:Nneurons,1)=NaN;
    R_3R.Ev(k).Meanz(1:Nneurons,1)=NaN;
    R_3R.Ev(k).BW(1:Nneurons,1)=NaN;
    R_3R.Ev(k).signrank(1:Nneurons,1)=NaN;
    R_3R.Ev(k).RespDir(1:Nneurons,1)=NaN;
    R_3R.Ev(k).NumberTrials(1:Nneurons,1)=NaN;
end

% runs the main routine

for i=1:length(RAW) %loops through sessions
    if strcmp('TH',RAW(i).Type(1:2))
        for j= 1:size(RAW(i).Nrast,1) %Number of neurons per session
            NN=NN+1; %neuron 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

                if sum(EvInd)==0
                    fprintf('HOWDY, CANT FIND EVENTS FOR ''%s''\n',Erefnames{k});
                end

                R_3R.Ev(k).NumberTrials(NN,1)=length(RAW(i).Erast{EvInd});
                Tbase=prewin(k,1):BSIZE:prewin(k,2);
                if  ~isempty(EvInd) && R_3R.Ev(k).NumberTrials(NN,1)>=MinNumTrials %avoid analyzing sessions where that do not have enough trials
                    [PSR0,N0]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},prewin(k,:),{1});% makes collapsed rasters for baseline for use in normalizing
                    [PSR1,N1]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},Dura,{1});% makes collpased rasters. PSR1 is a cell(neurons)
                    [PSR2,N2]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},Dura,{2});% makes trial by trial rasters. PSR1 is a cell(neurons, trials)
                    [PSR3,N0]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EvInd},prewin(k,:),{2});% makes trial by trial rasters for baseline for use in response detection
                    if ~isempty(PSR0{1}) || ~isempty(PSR1{1}) %to avoid errors, added on 12/28 2011
                        %optimal bin size
                        %[PTH1,BW1,~]=MakePTH07(PSR1,repmat(N1, size(RAW(i).Nrast{j},1),1),{2,1});%-----DP used here
                        %[PTH0,~,~]=MakePTH07(PSR0,repmat(N0, size(RAW(i).Nrast{j},1),1),{1,0,BW1});%BW1 reinjected here to make sure PTH0 & PTH1 have the same BW
                        %Fixed bin size
                        [PTH1,BW1,~]=MakePTH07(PSR1,repmat(N1, size(RAW(i).Nrast{j},1),1),{2,0,BinSize});%-----Fixed bin used here
                        [PTH0,~,~]=MakePTH07(PSR0,repmat(N0, size(RAW(i).Nrast{j},1),1),{1,0,BinSize});%BW1 reinjected here to make sure PTH0 & PTH1 have the same BW
                        PTH1=smooth(PTH1,SmoothSPAN,SmoothTYPE)';
                        PTH0=smooth(PTH0,SmoothSPAN,SmoothTYPE)';

                        %------------- Fills the R_3R.Ev(k) fields --------------
                        R_3R.Ev(k).BW(NN,1)=BW1;
                        R_3R.Ev(k).PSTHraw(NN,1:length(Tm))=PTH1;
                        R_3R.Ev(k).Meanraw(NN,1)=nanmean(R_3R.Ev(k).PSTHraw(NN,Tm>postwin(k,1) & Tm<postwin(k,2)),2);
                        if sum(PTH0,2)~=0
                            R_3R.Ev(k).PSTHz(NN,1:length(Tm))=normalize(PTH1,PTH0,0);
                            R_3R.Ev(k).Meanz(NN,1)=nanmean(R_3R.Ev(k).PSTHz(NN,Tm>postwin(k,1) & Tm<postwin(k,2)),2);
                        else
                            R_3R.Ev(k).PSTHz(NN,1:length(Tm))=NaN(1,length(Tm));
                            R_3R.Ev(k).Meanz(NN,1)=NaN;
                        end

                        %------------------ firing (in Hz) per trial in pre/post windows ------------------
                        %used to make the between events comparisons and Response detection in a single window----
                        ev(k).pre=NaN(size(RAW(i).Erast{EvInd},1),1);
                        ev(k).post=NaN(size(RAW(i).Erast{EvInd},1),1);
                        for m=1:size(RAW(i).Erast{EvInd},1) %loops through trials
                            ev(k).pre(m)=sum(PSR3{m}<prewin(k,2) & PSR3{m}>prewin(k,1))/(prewin(k,2)-prewin(k,1)); %CHANGED FROM PSR2 to PSR0 here 10/24/17
                            ev(k).post(m)=sum(PSR2{m}<postwin(k,2) & PSR2{m}>postwin(k,1))/(postwin(k,2)-postwin(k,1));
                        end

                        %-------------------- signrank on event and direction----
                        [R_3R.Ev(k).signrank(NN,1),~]=signrank(ev(k).pre, ev(k).post); %Signrank used here because it is a dependant sample test
                        if R_3R.Ev(k).signrank(NN,1)<PStat
                            R_3R.Ev(k).RespDir(NN,1)=sign(mean(ev(k).post)-mean(ev(k).pre));
                        else R_3R.Ev(k).RespDir(NN,1)=0;
                        end

                    end %if ~isempty(PSR0{1}) || ~isempty(PSR1{1})
                end %if EvInd=0 OR n(trials) < MinNumTrials fills with NaN
            end %Events

            %----------------------- Check for reward modulation across all the bins surrounding reward delivery-----------------------

            EV1=strmatch('RD1',RAW(i).Einfo(:,2),'exact');
            EV2=strmatch('RD2',RAW(i).Einfo(:,2),'exact');
            EV3=strmatch('RD3',RAW(i).Einfo(:,2),'exact');

            Bcell1=0;
            Bcell2=0;
            Bcell3=0;
            BHz=0;
            Trial=0;
            Reward=0;

            %get mean baseline firing for all reward trials
            [Bcell1,B1n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV1},BinBase,{2});% makes trial by trial rasters for baseline
            for y= 1:B1n
                BHz(1,y)=sum(Bcell1{1,y}>BinBase(1))/(BinBase(1,2)-BinBase(1,1));
                Trial(1,y)=y;
                Reward(1,y)=1;
            end
            [Bcell2,B2n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV2},BinBase,{2});% makes trial by trial rasters for baseline
            for z= 1:B2n
                BHz(1,z+B1n)=sum(Bcell2{1,z}>BinBase(1))/(BinBase(1,2)-BinBase(1,1));
                Trial(1,z+B1n)=z+B1n;
                Reward(1,z+B1n)=2;
            end
            [Bcell3,B3n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV3},BinBase,{2});% makes trial by trial rasters for baseline
            for x= 1:B3n
                BHz(1,x+B2n+B1n)=sum(Bcell3{1,x}>BinBase(1))/(BinBase(1,2)-BinBase(1,1));
                Trial(1,x+B2n+B1n)=x+B2n+B1n;
                Reward(1,x+B2n+B1n)=3;
            end

            for k=1:bins %go through each bin for each neuron
                EVHz=0;
                EV1cell=0;
                EV2cell=0;
                EV3Cell=0;

                %get trial by trial firing rate for water trials
                [EV1cell,EV1n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV1},[(BinDura(1)+(binstart - binint)+k*binint) (BinDura(2)+(binstart - binint)+k*binint)],{2});% makes trial by trial rasters for event
                for y= 1:EV1n
                    EVHz(1,y)=sum(EV1cell{1,y}>binstart)/(BinDura(1,2)-BinDura(1,1));
                end       

                %get trial by trial firing rate for maltodextrin trials
                [EV2cell,EV2n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV2},[(BinDura(1)+(binstart - binint)+k*binint) (BinDura(2)+(binstart - binint)+k*binint)],{2});% makes trial by trial rasters for event
                for z= 1:EV2n
                    EVHz(1,z+EV1n)=sum(EV2cell{1,z}>binstart)/(BinDura(1,2)-BinDura(1,1));
                end

                %get trial by trial firing rate for water trials
                [EV3cell,EV3n]=MakePSR04(RAW(i).Nrast(j),RAW(i).Erast{EV3},[(BinDura(1)+(binstart - binint)+k*binint) (BinDura(2)+(binstart - binint)+k*binint)],{2});% makes trial by trial rasters for event
                for x= 1:EV3n
                    EVHz(1,x+EV1n+EV2n)=sum(EV3cell{1,x}>binstart)/(BinDura(1,2)-BinDura(1,1));
                end    

                [~,Bintbl,~]=anovan([BHz EVHz],{[Reward Reward],cat(2,zeros(1,length(Trial)),ones(1,length(Trial)))},'varnames',{'Reward','Pre vs Post'},'model','full','display','off');
                R_3R.BinStatRwrd{k,1}(NN).IntSig=cell2mat(Bintbl(4,7));
                
                %get mean firing rates
                R_3R.BinStatRwrd{k,1}(NN).R1Mean=(nanmean(EVHz(Reward==1))-nanmean(BHz(Reward==1)))/std(BHz(Reward==1));
                R_3R.BinStatRwrd{k,1}(NN).R2Mean=(nanmean(EVHz(Reward==2))-nanmean(BHz(Reward==2)))/std(BHz(Reward==2));
                R_3R.BinStatRwrd{k,1}(NN).R3Mean=(nanmean(EVHz(Reward==3))-nanmean(BHz(Reward==3)))/std(BHz(Reward==3));

                %signranks for checking if inhibitory or excitatory

                %sucrose trials
                signranksig=signrank(BHz(Reward==1), EVHz(Reward==1)); %Signrank used here because it is a dependant sample test
                if signranksig<PStat
                    R_3R.BinStatRwrd{k,1}(NN).SucRespDir=sign(nanmean(EVHz(Reward==1))-mean(BHz(Reward==1)));
                else
                    R_3R.BinStatRwrd{k,1}(NN).SucRespDir=0;
                end

                %maltodextrin trials
                signranksig=signrank(BHz(Reward==2), EVHz(Reward==2)); %Signrank used here because it is a dependant sample test
                if signranksig<PStat
                    R_3R.BinStatRwrd{k,1}(NN).MalRespDir=sign(nanmean(EVHz(Reward==2))-mean(BHz(Reward==2)));
                else
                    R_3R.BinStatRwrd{k,1}(NN).MalRespDir=0;
                end

                %water trials
                signranksig=signrank(BHz(Reward==3), EVHz(Reward==3)); %Signrank used here because it is a dependant sample test
                if signranksig<PStat
                    R_3R.BinStatRwrd{k,1}(NN).WatRespDir=sign(nanmean(EVHz(Reward==3))-mean(BHz(Reward==3)));
                else
                    R_3R.BinStatRwrd{k,1}(NN).WatRespDir=0;
                end

            end %bins for bins analysis


            
            fprintf('3R neuron #%d\n',NN);
        end %neurons: FOR j= 1:size(RAW(i).Nrast,1)
    end %if right kind of session
end %sessions: FOR i=1:length(RAW)

if SAVE_FLAG
    save('R_3R.mat','R_3R')
end

toc