Delete 'MatlabScripts/b_EventAnalysis.m'

MatlabScripts/b_EventAnalysis.m

@@ -1,681 +0,0 @@
-%% 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.02;
-
-%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
-                        [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});%-----DP 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});%-----DP 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_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});%-----DP 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
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-