%%
% Sept 26th, 2012: Added line 278. Now excludes Baseline calculations on excluded neurons
% May 21st, 2012: Added the time of the max(PSTH) and the max for signigicant responses
% May 15th, 2012: this version includes waveform analysis
% May 9th, 2012: The baseline used to compute z-scores is matched with the prewindow matrix
% May 9th, 2012: R.Erefnames added to avoid having to fetch it from the excel file
% May 2nd, 2012: this version uses parametric tests (ttest/ttest2 instead of SignRank and Ranksum)
% April 27th, 2012:
%   -handles both GoNoGo and 2Go exps
%   -uses the new cleaned-up ResDetectSignRank02 code
% April 13th, 2012:
%   -uses a different way of excluding neurons (it skips the excluded neurons instead of filtering then afterwards
%   -R.Structure added that stores the brain region (DMS, CORE, SHELL)

clear all; clc;
global Dura Baseline Tm Tbase BSIZE Tbin
tic
path='RESULTpas.xls';
load('RAWvppasALL.mat')
RAW=RAWvppasALL;

%Main settings
SAVE_FLAG=1;
BSIZE=0.01; %Do not change
Dura=[-22 20]; 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=1;
RvppasALL=[];RvppasALL.Ninfo={};NN=0;Nneurons=0;
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


% List of events to analyze and analysis windows EXTRACTED from excel file
[~,Erefnames]=xlsread(path,'Windows','a3:a12'); % cell that contains the event names
prewin =  xlsread(path,'Windows','b3:c12');
postwin = xlsread(path,'Windows','d3:e12');
BLWin=    xlsread(path,'Windows','f3:g12');
RespWin=  xlsread(path,'Windows','h3:i12');

%Settings for Response detection w/ signrank
WINin=[0.03,0.3]; %Onset and offset requirements in sec.
%WINin=[-2 -4]; %negative values define onset requirement by the number of consecutive bins and not by duration
CI=[.99,.99];

%%
%Finds the total number of neurons accross all sessions
for i=1:length(RAW)
    RvppasALL.Ninfo=cat(1,RvppasALL.Ninfo,RAW(i).Ninfo);
    Nneurons=Nneurons+size(RAW(i).Nrast,1);
end

% load('CoordPPSall.mat');
% RPPSall.Coord=CoordPPSall;
RvppasALL.Coord=10*ones(Nneurons,4); %comment this line and uncomment the 2 lines above when you have real coordinates.
RvppasALL.Structure=RvppasALL.Coord(:,4);
RvppasALL.Ninfo=cat(2, RvppasALL.Ninfo, mat2cell([1:Nneurons]',ones(1,Nneurons),1));
RvppasALL.Erefnames= Erefnames;

RvppasALL.Rat(1:Nneurons,1)=NaN;
RvppasALL.Session(1:Nneurons,1)=NaN;

for i=1:length(RvppasALL.Ninfo)
    RvppasALL.Rat(i,1)=str2num(RvppasALL.Ninfo{i,1}(4:5));
    RvppasALL.Session(i,1)=str2num(RvppasALL.Ninfo{i,1}(10:11));
end

% preallocating
RvppasALL.Param.Tm=Tm;
RvppasALL.Param.Tbin=Tbin;
RvppasALL.Param.Dura=Dura;
RvppasALL.Param.Baseline=Baseline;
RvppasALL.Param.PStat=PStat;
RvppasALL.Param.MinNumTrials=MinNumTrials;
RvppasALL.Param.path=path;
RvppasALL.Param.prewin=prewin;
RvppasALL.Param.postwin=postwin;
RvppasALL.Param.BLwin=BLWin;
RvppasALL.Param.RespWin=RespWin;
RvppasALL.Param.ResponseReq=WINin;
RvppasALL.Param.SmoothTYPE=SmoothTYPE;
RvppasALL.Param.SmoothSPAN=SmoothSPAN;
RvppasALL.CSPlusResponseStat(1:Nneurons,1)=NaN;
RvppasALL.CSStat(1:Nneurons,1)=NaN;


for k=1:length(Erefnames)
    RvppasALL.Ev(k).PSTHraw(1:Nneurons,1:length(Tm))=NaN(Nneurons,length(Tm));
    RvppasALL.Ev(k).PSTHz(1:Nneurons,1:length(Tm))=NaN(Nneurons,length(Tm));
    RvppasALL.Ev(k).Meanraw(1:Nneurons,1)=NaN;
    RvppasALL.Ev(k).Meanz(1:Nneurons,1)=NaN;
    RvppasALL.Ev(k).BW(1:Nneurons,1)=NaN;
    RvppasALL.Ev(k).ttest(1:Nneurons,1)=NaN;
    RvppasALL.Ev(k).RespDir(1:Nneurons,1)=NaN;
    RvppasALL.Ev(k).RFLAG(1:Nneurons,1)=NaN;
    RvppasALL.Ev(k).CFLAG(1:Nneurons,1)=NaN;
    RvppasALL.Ev(k).Onsets(1:Nneurons,:)=NaN(Nneurons,2);
    RvppasALL.Ev(k).Offsets(1:Nneurons,:)=NaN(Nneurons,2);
    RvppasALL.Ev(k).NumberTrials(1:Nneurons,1)=NaN;
    RvppasALL.Ev(k).MaxVal(1:Nneurons,1)=NaN;
    RvppasALL.Ev(k).MaxTime(1:Nneurons,1)=NaN;
end

%% runs the main routine

for i=1:length(RAW) %loops through sessions
    for j= 1:size(RAW(i).Nrast,1) %Number of neurons per session
        NN=NN+1; %neuron counter
        if RvppasALL.Structure(NN)~=0 %to avoid analyzing excluded neurons
            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
                
                RvppasALL.Ev(k).NumberTrials(NN,1)=length(RAW(i).Erast{EvInd});
                Tbase=prewin(k,1):BSIZE:prewin(k,2);
                if  ~isempty(EvInd) && RvppasALL.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.
                    [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 trail rasters. PSR1 is a cell(neurons, trials)
                    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.Ev(k) fields --------------
                        RvppasALL.Ev(k).BW(NN,1)=BW1;
                        RvppasALL.Ev(k).PSTHraw(NN,1:length(Tm))=PTH1;
                        RvppasALL.Ev(k).Meanraw(NN,1)=nanmean(RvppasALL.Ev(k).PSTHraw(NN,Tm>postwin(k,1) & Tm<postwin(k,2)),2);
                        if sum(PTH0,2)~=0
                            RvppasALL.Ev(k).PSTHz(NN,1:length(Tm))=normalize(PTH1,PTH0,0);
                            RvppasALL.Ev(k).Meanz(NN,1)=nanmean(RvppasALL.Ev(k).PSTHz(NN,Tm>postwin(k,1) & Tm<postwin(k,2)),2);
                        else
                            RvppasALL.Ev(k).PSTHz(NN,1:length(Tm))=NaN(1,length(Tm));
                            RvppasALL.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(PSR2{m}<prewin(k,2) & PSR2{m}>prewin(k,1))/(prewin(k,2)-prewin(k,1));
                            ev(k).post(m)=sum(PSR2{m}<postwin(k,2) & PSR2{m}>postwin(k,1))/(postwin(k,2)-postwin(k,1));
                        end
                        
                        %---------------------------Response detection w/ SignRank on pre/post windows---------------------------
                        [~,RvppasALL.Ev(k).ttest(NN,1)]=ttest(ev(k).pre, ev(k).post); %Signrank used here because it is a dependant sample test
                        if RvppasALL.Ev(k).ttest(NN,1)<PStat
                            RvppasALL.Ev(k).RespDir(NN,1)=sign(mean(ev(k).post)-mean(ev(k).pre));
                            if RvppasALL.Ev(k).RespDir(NN,1)>0
                                Search=find(Tm>=postwin(k,1) & Tm<=postwin(k,2));
                                [RvppasALL.Ev(k).MaxVal(NN,1),MaxInd]=max(RvppasALL.Ev(k).PSTHraw(NN,Search));
                                RvppasALL.Ev(k).MaxTime(NN,1)=Tm(Search(1)+MaxInd-1);
                            else
                                Search=find(Tm>=postwin(k,1) & Tm<=postwin(k,2));
                                [RvppasALL.Ev(k).MaxVal(NN,1),MaxInd]=min(RvppasALL.Ev(k).PSTHraw(NN,Search));
                                RvppasALL.Ev(k).MaxTime(NN,1)=Tm(Search(1)+MaxInd-1);
                            end
                        else RvppasALL.Ev(k).RespDir(NN,1)=0;
                        end
                        
                        %---------------------------Response detection w/ RESDETECT08 on running windows---------------------------
                        RD=ResDetect08(PTH1,PTH0,RespWin(k,:),0,WINin, CI);
                        
                        RvppasALL.Ev(k).RFLAG(NN,1)=RD.RFLAG;
                        RvppasALL.Ev(k).CFLAG(NN,1)=RD.CFLAG;
                        RvppasALL.Ev(k).Onsets(NN,:)=RD.Onsets';
                        RvppasALL.Ev(k).Offsets(NN,:)=RD.Offsets';
                        RvppasALL.Param.Paramnames=RD.Paramnames;
                        RvppasALL.Param.RDParam=RD.Params;
                        
                    end %if ~isempty(PSR0{1}) || ~isempty(PSR1{1})
                end %if EvInd=0 OR n(trials) < MinNumTrials fills with NaN
            end %Events
            
            
            %----------------------- CUES -----------------------
            CSPlus=strcmp('CSPlus', Erefnames);
            CSMinus=strcmp('CSMinus', Erefnames);
            if sum(CSPlus)~=0 %To avoid this analysis is Go and NoGo cues are not included in the analysis
                if sum(CSMinus)~=0
                    RvppasALL.CScriteria=linspace(0, max([ev(CSPlus).post(:);ev(CSMinus).post(:)]),200);%Determine the criteria range
                    RvppasALL.CSprecriteria=linspace(0, max([ev(CSPlus).pre(:);ev(CSMinus).pre(:)]),200);%Determine the criteria range
                    for s=1:length(RvppasALL.CScriteria);
                        RvppasALL.CSfalsepos(NN,s)=sum(ev(CSMinus).post>=RvppasALL.CScriteria(s))/length(ev(CSMinus).post); %Determine the probability that the baseline is greater than each criteria (x)
                        RvppasALL.CStruepos(NN,s)=sum(ev(CSPlus).post>=RvppasALL.CScriteria(s))/length(ev(CSPlus).post);%Determine the probability that the post-window firing firing is greater than critera(y)
                        RvppasALL.CSprefalsepos(NN,s)=sum(ev(CSMinus).pre>=RvppasALL.CScriteria(s))/length(ev(CSMinus).pre); %Determine the probability that the baseline is greater than each criteria (x)
                        RvppasALL.CSpretruepos(NN,s)=sum(ev(CSPlus).pre>=RvppasALL.CScriteria(s))/length(ev(CSPlus).pre);%Determine the probability that the pre-window firing firing is greater than critera(y)
                    end
                    RvppasALL.CSauROC(NN,1)=trapz(RvppasALL.CSfalsepos(NN,:),RvppasALL.CStruepos(NN,:));%Calculate area under the curve for this neuron for this event
                    RvppasALL.CSpreauROC(NN,1)=trapz(RvppasALL.CSprefalsepos(NN,:),RvppasALL.CSpretruepos(NN,:));%Calculate area under the curve for this neuron for this event
                    [~,RvppasALL.CSStat(NN,1)]=ttest2(ev(CSPlus).post,ev(CSMinus).post); %Ranksum test used becasue it is an independant sample test
                else RvppasALL.CSStat(NN,1)=NaN;
                     RvppasALL.CSauROC(NN,1)=NaN;
                     RvppasALL.CSpreauROC(NN,1)=NaN;
                end
            end
 
                                
            
            PECSPlus=strcmp('PECSPlus', Erefnames);
            PECSMinus=strcmp('PECSMinus', Erefnames);
            if sum(PECSPlus)~=0 %To avoid this analysis is Go and NoGo cues are not included in the analysis
                if sum(PECSMinus)~=0
                    [RvppasALL.PECueStat(NN,1),~]=ranksum(ev(PECSPlus).post,ev(PECSMinus).post); %Ranksum test used becasue it is an independant sample test
                else RvppasALL.PECueStat(NN,1)=NaN;
                end
            end
            
            PECSPlus=strcmp('PECSPlus', Erefnames);
            ITIPE=strcmp('ITIPE', Erefnames);
            if sum(PECSPlus)~=0 %To avoid this analysis is Go and NoGo cues are not included in the analysis
                if RvppasALL.Ev(PECSPlus).RespDir(NN,1)~=0 || RvppasALL.Ev(ITIPE).RespDir(NN,1)~=0
                    [RvppasALL.PERewStat(NN,1),~]=ranksum(ev(PECSPlus).post,ev(ITIPE).post); %Ranksum test used becasue it is an independant sample test
                else RvppasALL.PERewStat(NN,1)=NaN;
                end
            end
            
            %             PECSMinus=strcmp('PECSMinus', Erefnames);
            %             ITIPE=strcmp('ITIPE', Erefnames);
            %             if sum(PECSMinus)~=0 %To avoid this analysis is Go and NoGo cues are not included in the analysis
            %                 if R.Ev(PECSMinus).RespDir(NN,1)~=0 || R.Ev(ITIPE).RespDir(NN,1)~=0
            %                     [R.PEUnrewCueStat(NN,1),~]=ranksum(ev(PECSMinus).post,ev(ITIPE).post); %Ranksum test used becasue it is an independant sample test
            %                 else R.PEUnrewCueStat(NN,1)=NaN;
            %                 end
            %             end
            
            CSPluswPE=strcmp('CSPluswPE', Erefnames);
            CSPlusnoPE=strcmp('CSPlusnoPE', Erefnames);
            if sum(CSPlusnoPE)~=0 %To avoid this analysis is Go and NoGo cues are not included in the analysis
                if sum(CSPluswPE)~=0
                    RvppasALL.CSrespcriteria=linspace(0, max([ev(CSPlusnoPE).post(:);ev(CSPluswPE).post(:)]),200);%Determine the criteria range
                    RvppasALL.CSprerespcriteria=linspace(0, max([ev(CSPlusnoPE).pre(:);ev(CSPluswPE).pre(:)]),200);%Determine the criteria range
                    for s=1:length(RvppasALL.CSrespcriteria);
                        RvppasALL.CSrespfalsepos(NN,s)=sum(ev(CSPlusnoPE).post>=RvppasALL.CSrespcriteria(s))/length(ev(CSPlusnoPE).post); %Determine the probability that the baseline is greater than each criteria (x)
                        RvppasALL.CSresptruepos(NN,s)=sum(ev(CSPluswPE).post>=RvppasALL.CSrespcriteria(s))/length(ev(CSPluswPE).post);%Determine the probability that the post-window firing firing is greater than critera(y)
                        RvppasALL.CSprerespfalsepos(NN,s)=sum(ev(CSPlusnoPE).pre>=RvppasALL.CSrespcriteria(s))/length(ev(CSPlusnoPE).pre); %Determine the probability that the baseline is greater than each criteria (x)
                        RvppasALL.CSpreresptruepos(NN,s)=sum(ev(CSPluswPE).pre>=RvppasALL.CSrespcriteria(s))/length(ev(CSPluswPE).pre);%Determine the probability that the post-window firing firing is greater than critera(y)
                    end
                    RvppasALL.CSrespauROC(NN,1)=trapz(RvppasALL.CSrespfalsepos(NN,:),RvppasALL.CSresptruepos(NN,:));%Calculate area under the curve for this neuron for this event
                    RvppasALL.CSprerespauROC(NN,1)=trapz(RvppasALL.CSprerespfalsepos(NN,:),RvppasALL.CSpreresptruepos(NN,:));
                    [~,RvppasALL.CSPlusResponseStat(NN,1)]=ttest2(ev(CSPluswPE).post,ev(CSPlusnoPE).post); %Ranksum test used becasue it is an independant sample test
                else RvppasALL.CSPlusResponseStat(NN,1)=NaN;
                    RvppasALL.CSrespauROC(NN,1)=NaN;
                end
            else
                RvppasALL.CSrespauROC(NN,1)=NaN;
                RvppasALL.CSPlusResponseStat(NN,1)=NaN;
                RvppasALL.CSprerespauROC(NN,1)=NaN;
            end
            
            
            
            %
            CSMinuswPE=strcmp('CSMinuswPE', Erefnames);
            CSMinusnoPE=strcmp('CSMinusnoPE', Erefnames);
            if sum(CSMinuswPE)~=0 %To avoid this analysis is Go and NoGo cues are not included in the analysis
                if RvppasALL.Ev(CSMinuswPE).RespDir(NN,1)~=0 || RvppasALL.Ev(CSMinusnoPE).RespDir(NN,1)~=0
                    [RvppasALL.CSMinusResponseStat(NN,1),~]=ranksum(ev(CSMinuswPE).post,ev(CSMinusnoPE).post); %Ranksum test used becasue it is an independant sample test
                else RvppasALL.CSMinusResponseStat(NN,1)=NaN;
                end
            end
            %
    RvppasALL.CSPlusRatio(NN,1)=length(ev(CSPluswPE).post)/length(ev(CSPlus).post);
    RvppasALL.CSMinusRatio(NN,1)=length(ev(CSMinuswPE).post)/length(ev(CSMinus).post);
            fprintf('Neuron ID # %d\n',NN);
        elseif RvppasALL.Structure(NN)~=0
        end %exclusion: IF R.Structure(NN)~=0 to avoid analyzing excluded neurons
    end %neurons: FOR j= 1:size(RAW(i).Nrast,1)
end %sessions: FOR i=1:length(RAW)

RvppasALL.CSfalseposMEAN(1,:)=nanmean(RvppasALL.CSfalsepos(RvppasALL.Structure==10,:));
RvppasALL.CSfalseposSEM(1,:)=nanste(RvppasALL.CSfalsepos(RvppasALL.Structure==10,:),1);

RvppasALL.CStrueposMEAN(1,:)=nanmean(RvppasALL.CStruepos(RvppasALL.Structure==10,:));
RvppasALL.CStrueposSEM(1,:)=nanste(RvppasALL.CStruepos(RvppasALL.Structure==10,:),1);

RvppasALL.CSrespfalseposMEAN(1,:)=nanmean(RvppasALL.CSrespfalsepos(RvppasALL.Structure==10,:));
RvppasALL.CSrespfalseposSEM(1,:)=nanste(RvppasALL.CSrespfalsepos(RvppasALL.Structure==10,:),1);

RvppasALL.CSresptrueposMEAN(1,:)=nanmean(RvppasALL.CSresptruepos(RvppasALL.Structure==10,:));
RvppasALL.CSresptrueposSEM(1,:)=nanste(RvppasALL.CSresptruepos(RvppasALL.Structure==10,:),1);

% RPPSall.CSprefalseposMEAN(1,:)=nanmean(RPPSall.CSprefalsepos(RPPSall.Structure==10,:));
% RPPSall.CSprefalseposSEM(1,:)=nanste(RPPSall.CSprefalsepos(RPPSall.Structure==10,:),1);
% 
% RPPSall.CSpretrueposMEAN(1,:)=nanmean(RPPSall.CSpretruepos(RPPSall.Structure==10,:));
% RPPSall.CSpretrueposSEM(1,:)=nanste(RPPSall.CSpretruepos(RPPSall.Structure==10,:),1);


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

toc


%% Create a Summury excel spreadsheet

% RESULT=[];
% %RESULT=cat(2,RESULT, R.CueCorStat, R.CueCorRExtStat,R.CueIncorStat,R.CueIncorRExtStat,R.BaselineTrend, R.LExtTrend, R.BaselineLExtOne, R.BaselineLExtNo, R.DeltaLP, R.DeltaNo);
% 
% for k=1:length(Erefnames)
%     RESULT=cat(2,RESULT,RPPSall.Ev(k).RespDir, RPPSall.Ev(k).RFLAG, RPPSall.Ev(k).CFLAG,RPPSall.Ev(k).Onsets, RPPSall.Ev(k).Offsets);
% end
% 
% xlswrite(path,RPPSall.Ninfo,'RESULTsheet', 'A3'); % Writes Session - Neuron - ID#
% xlswrite(path,RPPSall.Structure,'RESULTsheet', 'D3');  % Writes the brain region
% xlswrite(path,RESULT,'RESULTsheet', 'F3'); % Writes RESULT matrix
% toc