%%

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


%Main settings
SAVE_FLAG=1;
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
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_DT5.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:a9'); % cell that contains the event names
prewin = [-1 0]; 
EventWin =  xlsread(path,'Windows','l3:m9'); %to compute MeanZ in relevant window
PreEventWin = xlsread(path,'Windows','d3:e9');
PostEventWin = xlsread(path,'Windows','f3:g9');
[~,Session] = xlsread(path,'Windows','c21:c30'); 

%%

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_DT5.mat');
    R.Ses(i).Coord=part(i+3).Coord;
    R.Ses(i).Structure=part(i+3).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).PSTHrawBL(1:Nneurons(i),1:501)=NaN(Nneurons(i),501);
        R.Ses(i).Ev(k).PSTHz(1:Nneurons(i),1:length(Tm))=NaN(Nneurons(i),length(Tm));
        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).Meanraw(1:Nneurons(i),1)=NaN;
        R.Ses(i).Ev(k).Meanz(1:Nneurons(i),1)=NaN;
        R.Ses(i).Ev(k).MeanzPRE(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).lowBL_marker(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;

NS=0;
TNN=0;

%% runs the main routine

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; % mean FR computed across whole session.
                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;
                
                LeverInsertionInd=strcmp('LeverInsertion',RAW(SesIndex(y)).Einfo(:,2)); %Ref event for BL is lever insertion
                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
                                        
                    clear PreEventLeverIns
                    EventTimestamp=RAW(SesIndex(y)).Erast{EvInd};
                    LI=RAW(SesIndex(y)).Erast{LeverInsertionInd};
                    
                    if k<7 && ~isempty(EventTimestamp) %lever press or PE event
                        for l=1:length(EventTimestamp) %loop thru trial
                            PreEventLeverIns(l,1)=LI(find(LI(:,1)<EventTimestamp(l),1,'last'));
                        end
                        CellLevIns{1,1}=PreEventLeverIns;
                    else
                        CellLevIns{1,1}=RAW(SesIndex(y)).Erast{LeverInsertionInd};
                    end
                    
                    if sum(EvInd)==0
                        fprintf('HOWDY, CANT FIND EVENTS FOR ''%s''\n',Erefnames{k});
                    end
                    
                    R.Ses(i).Ev(k).NumberTrials(NN,1)=length(RAW(SesIndex(y)).Erast{EvInd});
                    if  ~isempty(EvInd) && R.Ses(i).Ev(k).NumberTrials(NN,1)>= MinNumTrials
                        
                        Tbase=prewin(1,1):BSIZE:prewin(1,2);
                        [PSR0,N0]=MakePSR04(RAW(SesIndex(y)).Nrast(j),RAW(SesIndex(y)).Erast{LeverInsertionInd},prewin(1,:),{1});% makes collapsed rasters for baseline. BL=[-1 0] compared to LI
                        [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(CellLevIns,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(PTH0))=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);
                                R.Ses(i).Ev(k).MeanzPRE(NN,1)=nanmean(R.Ses(i).Ev(k).PSTHz(NN,Tm>PreEventWin(k,1) & Tm<PreEventWin(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;
                                R.Ses(i).Ev(k).MeanzPRE(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}<(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
                                
                                [~,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).PSTHz(NN,Search));
                                    R.Ses(i).Ev(k).MaxTime(NN,1)=Tm(Search(1)+MaxInd-1);
                                end
                            else R.Ses(i).Ev(k).lowBL_marker(NN,1)=1;
                            end
                            
                            
                        end %if ~isempty(PSR0{1}) || ~isempty(PSR1{1})
                    end %if EvInd=0 OR n(trials) < MinNumTrials fills with NaN
                end %Events
                
                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)-1
        tableTRN(:,i)=R.Ses(k).Ev(i).RespDirPre(:,1);
        tableTRN(:,i+length(Erefnames)-1)=R.Ses(k).Ev(i).RespDirPost(:,1);
    end
    R.Ses(k).TRN(:,1)=nansum(abs(tableTRN),2);
end
save('R_DT5.mat','R');
toc