%% 

clear all; clc; 
global Dura Baseline Tm Tbase BSIZE Tbin DuraITI TmITI
tic
path='E:\Youna\Post-doc\MATLAB\Youna Matlab\DT Nex Files\RESULTfs5.xls'; %excel file with windows etc


%Main settings
BSIZE=0.01; %Do not change
Dura=[-25 20]; Tm=Dura(1):BSIZE:Dura(2);
DuraITI=[-50 20]; TmITI=DuraITI(1):BSIZE:DuraITI(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.01; %for comparing pre versus post windows, or event A versus event B
MinNumTrials=10;
R=[];R.Ninfo={};NN=0;
Nneurons=[];
Neur=0;
BinSize=0.05;

load('RAW_FS5.mat')
RAW=RAW;

%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:a8'); % for FS5
prewin =  [-1 0]; %BL window relative to LP1, for all events
EventWin =  xlsread(path,'Windows','b3:c8'); %for FR5
PreEventWin = xlsread(path,'Windows','d3:e8');
PostEventWin = xlsread(path,'Windows','f3:g8');
[~,Session] = xlsread(path,'Windows','c23:c29');

%% 

R.Erefnames= Erefnames;
R.Session= Session;
SessionName = {RAW.SessionName};
Ncum=0;

%Finds the total number of neurons accross all sessions
for i=1:length(Session) %loops through sessions
    Ses=strcmp(Session(i),SessionName);
    SesIndex=find(Ses==1);    
    NN=0; Neur=0;
    R.Ses(i).Ninfo=[];
    
    for y=1:length(SesIndex)
        R.Ses(i).Ninfo=cat(1,R.Ses(i).Ninfo,RAW(SesIndex(y)).Ninfo);
        Neur=Neur+size(RAW(SesIndex(y)).Nrast,1);
        Nneurons(i)=Neur;
    end
    
        load('Coord_FS5.mat'); 
        R.Ses(i).Coord=part(i+5).Coord;
        R.Ses(i).Structure=part(i+5).Coord(:,4);
        R.Ses(i).Ninfo=cat(2, R.Ses(i).Ninfo, mat2cell([1:Nneurons(i)]',ones(1,Nneurons(i)),1));
        
    for k=1:length(Erefnames)
        R.Ses(i).Ev(k).PSTHraw(1:Nneurons(i),1:length(Tm))=NaN(Nneurons(i),length(Tm));
        R.Ses(i).Ev(k).PSTHz(1:Nneurons(i),1:length(Tm))=NaN(Nneurons(i),length(Tm));
        R.Ses(i).Ev(k).Meanraw(1:Nneurons(i),1)=NaN;
        
        R.Ses(i).Ev(k).PSTHrawBL(1:Nneurons(i),1:length(Tbase))=NaN(Nneurons(i),length(Tbase));
        R.Ses(i).Ev(k).Meanz(1:Nneurons(i),1)=NaN;
        R.Ses(i).Ev(k).rawMeanz(1:Nneurons(i),1)=NaN;
        R.Ses(i).Ev(k).rawMeanzPre(1:Nneurons(i),1)=NaN;
        R.Ses(i).Ev(k).ttestPreEvent(1:Nneurons(i),1)=NaN;
        R.Ses(i).Ev(k).ttestPostEvent(1:Nneurons(i),1)=NaN;
        R.Ses(i).Ev(k).RespDirPre(1:Nneurons(i),1)=NaN;
        R.Ses(i).Ev(k).RespDirPost(1:Nneurons(i),1)=NaN;
        R.Ses(i).Ev(k).NumberTrials(1:Nneurons(i),1)=NaN;
    end   
end


% preallocating
R.Param.Tm=Tm;
R.Param.Tbin=Tbin;
R.Param.Dura=Dura;
R.Param.Baseline=Baseline;
R.Param.PStat=PStat;
R.Param.MinNumTrials=MinNumTrials;
R.Param.path=path;
R.Param.prewin=prewin;
R.Param.postwinPreEvent=PreEventWin;
R.Param.postwinPostEvent=PostEventWin;
R.Param.SmoothTYPE=SmoothTYPE;
R.Param.SmoothSPAN=SmoothSPAN;


%% runs the main routine

NS=0;
TNN=0;             
for i=1:length(Session) %loops through sessions
    Ses=strcmp(Session(i),SessionName);
    SesIndex=find(Ses==1);    
    NN=0;
    NS=NS+1;
    
    for y=1:length(SesIndex)
        for j= 1:size(RAW(SesIndex(y)).Nrast,1) %Number of neurons per session
            NN=NN+1; %neuron counter
            TNN=TNN+1; %neuron counter across all session
            
            if R.Ses(i).Structure(NN)~=0 %to avoid analyzing excluded neurons              
                ev2=length(RAW(SesIndex(y)).Nrast{j})/RAW(SesIndex(y)).Nrast{j}(end);% FR computed across whole session.
                R.WF(TNN,:)=RAW(SesIndex(y)).waveforms(j,:);
                R.WFanalysis(TNN,1)=ev2; 
                proportionISI=0;
                for k=3:1:size(RAW(SesIndex(y)).Nrast{j,1},1) % loop thru timestamps waveform
                    ISIn=RAW(SesIndex(y)).Nrast{j,1}(k)-RAW(SesIndex(y)).Nrast{j,1}(k-1); 
                    ISIo=RAW(SesIndex(y)).Nrast{j,1}(k-1)-RAW(SesIndex(y)).Nrast{j,1}(k-2);
                    calcul(k-2)=2*abs(ISIn-ISIo)/(ISIn+ISIo);
                    if ISIo>0.5
                        proportionISI=proportionISI+1;
                    end
                end
                R.CV2(TNN,1)=mean(calcul);%CV2
                R.CV2(TNN,2)=100*proportionISI/(size(RAW(SesIndex(y)).Nrast{j,1},1)-1);% of ISI>0.5s
                R.SESSION(TNN,1)=NS;
                
                LP1Ind=strcmp('First_LPcomp',RAW(SesIndex(y)).Einfo(:,2)); %find the event id number from RAW
                for k=1:length(Erefnames) %loops thorough the events
                    EvInd=strcmp(Erefnames(k),RAW(SesIndex(y)).Einfo(:,2)); %find the event id number from RAW
                    if sum(EvInd)==0
                        fprintf('HOWDY, CANT FIND EVENTS FOR ''%s''\n',Erefnames{k});
                    end
                    if k==1
                        FirstPressInd=EvInd;
                    end
                    if ~isempty(RAW(SesIndex(y)).Erast{EvInd})
                        
                        R.Ses(i).Ev(k).NumberTrials(NN,1)=length(RAW(SesIndex(y)).Erast{EvInd});
                        Tbase=prewin(1,1):BSIZE:prewin(1,2); % for FR5 and FS, same prewin for all events
                        [PSR0,N0]=MakePSR04(RAW(SesIndex(y)).Nrast(j),RAW(SesIndex(y)).Erast{LP1Ind},prewin(1,:),{1});
                        [PSR1,N1]=MakePSR04(RAW(SesIndex(y)).Nrast(j),RAW(SesIndex(y)).Erast{EvInd},Dura,{1});% makes collpased rasters. PSR1 is a cell(neurons)
                        [PSR2,N2]=MakePSR04(RAW(SesIndex(y)).Nrast(j),RAW(SesIndex(y)).Erast{EvInd},Dura,{2});% makes trial by trail rasters. PSR1 is a cell(neurons, trials)
                        BLcell=MakePSR04(RAW(SesIndex(y)).Erast(FirstPressInd),RAW(SesIndex(y)).Erast{EvInd},[-1800 1],{2,'Last'});
                        
                        if ~isempty(PSR0{1}) || ~isempty(PSR1{1}) %to avoid errors, added on 12/28 2011
                            %Fixed bin size
                            [PTH1,BW1,~]=MakePTH07(PSR1,repmat(N1, size(RAW(SesIndex(y)).Nrast{j},1),1),{2,0,BinSize});%-----DP used here
                            [PTH0,~,~]=MakePTH07(PSR0,repmat(N0, size(RAW(SesIndex(y)).Nrast{j},1),1),{1,0,BinSize});%BW1 reinjected here to make sure PTH0 & PTH1 have the same BW
                            % calculate MeanZ from PSTH before smoothing
                            
                            if sum(PTH0,2)~=0
                                unsmoothPSTHz(1,1:length(Tm))=normalize(PTH1,PTH0,0);
                                R.Ses(i).Ev(k).rawMeanz(NN,1)=nanmean(unsmoothPSTHz(1,Tm>EventWin(k,1) & Tm<EventWin(k,2)),2);
                                R.Ses(i).Ev(k).rawMeanzPre(NN,1)=nanmean(unsmoothPSTHz(1,Tm>PreEventWin(k,1) & Tm<PreEventWin(k,2)),2);
                            else
                                R.Ses(i).Ev(k).rawMeanz(NN,1)=NaN;
                                R.Ses(i).Ev(k).rawMeanzPre(NN,1)=NaN;
                            end
                            
                            PTH1=smooth(PTH1,SmoothSPAN,SmoothTYPE)';
                            PTH0=smooth(PTH0,SmoothSPAN,SmoothTYPE)';
                            
                            %------------- Fills the R.Ses(i).Ev(k) fields --------------
                            R.Ses(i).Ev(k).PSTHraw(NN,1:length(Tm))=PTH1;
                            R.Ses(i).Ev(k).PSTHrawBL(NN,1:length(Tbase))=PTH0;
                            R.Ses(i).Ev(k).Meanraw(NN,1)=nanmean(R.Ses(i).Ev(k).PSTHraw(NN,Tm>PostEventWin(k,1) & Tm<PostEventWin(k,2)),2);
                            if sum(PTH0,2)~=0
                                R.Ses(i).Ev(k).PSTHz(NN,1:length(Tm))=normalize(PTH1,PTH0,0);
                                R.Ses(i).Ev(k).Meanz(NN,1)=nanmean(R.Ses(i).Ev(k).PSTHz(NN,Tm>EventWin(k,1) & Tm<EventWin(k,2)),2);
                            else
                                R.Ses(i).Ev(k).PSTHz(NN,1:length(Tm))=NaN(1,length(Tm));
                                R.Ses(i).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(SesIndex(y)).Erast{EvInd},1),1);
                            ev(k).PreEvent=NaN(size(RAW(SesIndex(y)).Erast{EvInd},1),1);
                            ev(k).PostEvent=NaN(size(RAW(SesIndex(y)).Erast{EvInd},1),1);
                            for m=1:size(RAW(SesIndex(y)).Erast{EvInd},1) %loops through trials
%                                 ev(k).pre(m)=sum(PSR2{m}<prewin(1,2) & PSR2{m}>prewin(1,1))/(prewin(1,2)-prewin(1,1));
                                ev(k).pre(m)=sum(PSR2{m}<(BLcell{1,m}+prewin(1,2)) & PSR2{m}>(BLcell{1,m}+prewin(1,1)))/(prewin(1,2)-prewin(1,1));
                                ev(k).PreEvent(m)=sum(PSR2{m}<PreEventWin(k,2) & PSR2{m}>PreEventWin(k,1))/(PreEventWin(k,2)-PreEventWin(k,1));
                                ev(k).PostEvent(m)=sum(PSR2{m}<PostEventWin(k,2) & PSR2{m}>PostEventWin(k,1))/(PostEventWin(k,2)-PostEventWin(k,1));
                            end
                            %---------------------------Response detection w/ SignRank on pre/post windows---------------------------
                            
                            if  ~isempty(EvInd) && R.Ses(i).Ev(k).NumberTrials(NN,1)>= MinNumTrials  && nanmean(ev(k).pre(:,1),1)>0.5 %avoid analyzing sessions where that do not have enough trials %avoid analyzing sessions where that do not have enough trials
                                
                                [~,R.Ses(i).Ev(k).ttestPreEvent(NN,1)]=signrank(ev(k).pre, ev(k).PreEvent); %Signrank used here because it is a dependant sample test
                                if R.Ses(i).Ev(k).ttestPreEvent(NN,1)<PStat
                                    R.Ses(i).Ev(k).RespDirPre(NN,1)=sign(mean(ev(k).PreEvent)-mean(ev(k).pre));
                                else R.Ses(i).Ev(k).RespDirPre(NN,1)=0;
                                end
                                
                                [~,R.Ses(i).Ev(k).ttestPostEvent(NN,1)]=signrank(ev(k).pre, ev(k).PostEvent); %Signrank used here because it is a dependant sample test
                                if R.Ses(i).Ev(k).ttestPostEvent(NN,1)<PStat
                                    R.Ses(i).Ev(k).RespDirPost(NN,1)=sign(mean(ev(k).PostEvent)-mean(ev(k).pre));
                                else R.Ses(i).Ev(k).RespDirPost(NN,1)=0;
                                end
                                
                                if R.Ses(i).Ev(k).RespDirPre(NN,1)>0 || R.Ses(i).Ev(k).RespDirPost(NN,1) >0
                                    Search=find(Tm>=PreEventWin(k,1) & Tm<=PostEventWin(k,2));
                                    [R.Ses(i).Ev(k).MaxVal(NN,1),MaxInd]=max(R.Ses(i).Ev(k).PSTHraw(NN,Search));
                                    R.Ses(i).Ev(k).MaxTime(NN,1)=Tm(Search(1)+MaxInd-1);
                                else
                                    Search=find(Tm>=PreEventWin(k,1) & Tm<PostEventWin(k,2));
                                    [R.Ses(i).Ev(k).MaxVal(NN,1),MaxInd]=min(R.Ses(i).Ev(k).PSTHraw(NN,Search));
                                    R.Ses(i).Ev(k).MaxTime(NN,1)=Tm(Search(1)+MaxInd-1);
                                end

                            end %if EvInd=0 OR n(trials) < MinNumTrials fills with NaN
                        end %if ~isempty(PSR0{1}) || ~isempty(PSR1{1})
                    end
                end
                
                fprintf('Neuron ID # %d\n',NN);
               elseif R.Ses(i).Structure(NN)~=0
            end %exclusion: IF R.Structure(NN)~=0 to avoid analyzing excluded neurons
        end %neurons: FOR j= 1:size(RAW(SesIndex(y)).Nrast,1)
    end %SesIndex
end %sessions: FOR i=1:length(RAW)

for k=1:size(R.Ses,2)
    tableTRN=[];
    for i=1:length(Erefnames)
        tableTRN(:,i)=R.Ses(k).Ev(i).RespDirPre(:,1);
        tableTRN(:,i+length(Erefnames))=R.Ses(k).Ev(i).RespDirPost(:,1);
    end
    R.Ses(k).TRN(:,1)=nansum(abs(tableTRN),2);
end
save('R_FS5','R')
toc