%% Get postions from optitrack OR and calculate PCs

clearvars;
dirData = 'A:\ORproject\';
    % only considering ses with clustering and good SWR
animals = {'Wmaze2\OR15\','Wmaze2\OR18\','Wmaze3\OR22\','Wmaze3\OR21\','Wmaze3\OR23\'};% 'OR15\',  
sessions = {{'day1','day2','day3','day4','day10'},{'day1','day2','day3'},{'day1','day4','day3','day5'},...
            {'day2','day4'},{'day1','day5'}};
% DONE: OR22-1,3,4,5; OR21-4,2; OR15-1,2,3,4,10; OR18-2,3; OR23-1,5
% PROB: OR22-3    ;OR18-1        
%% Manually get tracking for indiv sessions
for rat = 1:length(animals)
    for ses = 1:length(sessions{rat})
        cd([dirData animals{rat} sessions{rat}{ses}]);
        
        basename = bz_BasenameFromBasepath(pwd); 
        load([basename '.sessionInfo.mat']);sr = sessionInfo.rates.wideband;
        
        % manual (go to subfolder...)
        bz_processConvertOptitrack2Behav('OR18day1_2','syncSampFq',sr);
        % 
    end
end

%% Get tracking positions and trials
for rat = 2%:length(animals)
    for ses = 1:length(sessions{rat})
        cd([dirData animals{rat} sessions{rat}{ses}]);
        basename = bz_BasenameFromBasepath(pwd); 

        if ~exist([basename '.MergePoints.events.mat'])
           bz_ConcatenatedTimes(pwd);
        end
        load([basename '.MergePoints.events.mat']);

        % get tracking from individual sessions
        withTracking = [];
        for i = 1:length(MergePoints.foldernames)
            cd(MergePoints.foldernames{i});
            if ~isempty(dir('*.tracking.behavior.mat'))
               fileT = dir('*.tracking.behavior.mat');
               trackingSes{i} = load(fileT.name);
               withTracking = cat(1,withTracking,i);
            end
            clear tracking
            cd ..
        end
        
        % cat tracking and correct offset
        if exist('trackingSes','var')
           tracking.timestamps = [];
           tracking.orientation.rx=[];tracking.orientation.ry=[];tracking.orientation.rz=[];tracking.orientation.rw=[]; 
           tracking.position.x=[];tracking.position.y=[];tracking.position.z=[];
           for i = 1:length(withTracking)
               timestampsOffset = trackingSes{withTracking(i)}.tracking.timestamps + MergePoints.timestamps(withTracking(i),1);
               tracking.timestamps = cat(1,tracking.timestamps,timestampsOffset); clear timestampsOffset;
               
               tracking.orientation.rx = cat(1,tracking.orientation.rx,trackingSes{withTracking(i)}.tracking.orientation.rx);
               tracking.orientation.ry = cat(1,tracking.orientation.ry,trackingSes{withTracking(i)}.tracking.orientation.ry);
               tracking.orientation.rz = cat(1,tracking.orientation.rz,trackingSes{withTracking(i)}.tracking.orientation.rz);
               tracking.position.x = cat(1,tracking.position.x,trackingSes{withTracking(i)}.tracking.position.x);
               tracking.position.y = cat(1,tracking.position.y,trackingSes{withTracking(i)}.tracking.position.y);
               tracking.position.z = cat(1,tracking.position.z,trackingSes{withTracking(i)}.tracking.position.z);
           end
           save([basename '.tracking.behavior.mat'],'tracking'); 
           %figure;plot(tracking.position.x(1:end-10000),tracking.position.z(1:end-10000));
           
        % get trials   
          load([basename '.sessionInfo.mat']);sr = sessionInfo.rates.wideband;
          [SCORE, trials] = wMazeScoreTrials(pwd, sr, 0); 
          save([basename '.WmazeTrials.mat'],'trials');
          clear SCORE trials sr sessionInfo tracking;
          
        else
           warning([animals{rat} sessions{rat}{ses} ': NO tracking found']); 
        end
        
    end
end

% separte postions for trial type 
for rat = 2:length(animals)
    for ses = 1:length(sessions{rat})
        cd([dirData animals{rat} sessions{rat}{ses}]);
        basename = bz_BasenameFromBasepath(pwd); 

        if exist([basename '.tracking.behavior.mat'])
           load([basename '.tracking.behavior.mat']);
           load([basename '.WmazeTrials.mat']);
           %save([basename '.WmazeTrialsCorrected.mat'],'trials');
           
           % trial types
           iL=0;    iR=0;   oL=0;   oR=0; clear posTrials;
           for i = 1:length(trials.int)
               if trials.outbound(i) == 1 && trials.left(i) == 1
                   oL = oL +1;
                  posTrials{3}.int(oL,:) =  trials.int(i,:);
                  posTrials{3}.correct(oL,1) =  trials.correct(i,1);      
                  posTrials{3}.type = 'outboundLeft';
                  % remove trials from side to side arms
                  if i > 1 && trials.correct(i) == 0 && trials.outbound(i-1) == 1
                     posTrials{3}.int(oL,:) = nan(1,2); 
                  end
                elseif trials.outbound(i) == 1 && trials.right(i) == 1
                   oR = oR +1;
                  posTrials{4}.int(oR,:) =  trials.int(i,:);
                  posTrials{4}.correct(oR,1) =  trials.correct(i,1);
                  posTrials{4}.type = 'outboundRight';
                  if i > 1 && trials.correct(i) == 0 && trials.outbound(i-1) == 1
                     posTrials{4}.int(oR,:) = nan(1,2); 
                  elseif i == 1 && posTrials{t}.correct(i) == 0  
                     posTrials{4}.int(oR,:) = nan(1,2);  
                  end               
               elseif trials.inbound(i) == 1 && trials.left(i) == 1
                  iL = iL+1;
                  posTrials{1}.int(iL,:) =  trials.int(i,:);
                  posTrials{1}.correct(iL,1) =  trials.correct(i,1);
                  posTrials{1}.type = 'inboundLeft';
               elseif trials.inbound(i) == 1 && trials.right(i) == 1
                  iR = iR+1;
                  posTrials{2}.int(iR,:) =  trials.int(i,:);
                  posTrials{2}.correct(iR,1) =  trials.correct(i,1); 
                  posTrials{2}.type = 'inboundRight';
               end
           end
           
           posXY = [tracking.timestamps tracking.position.x tracking.position.z];
           % remove outlayers
           for i = 1:size(posXY,1) % Wmaze3
               if posXY(i,2) > 750 || posXY(i,2) < -750 || ...
                  posXY(i,3) > 1600 || posXY(i,3) < -200
                  posXY(i,2:3) = nan(1,2);
               end
           end
           for i = 1:size(posXY,1) % Wmaze2
               if posXY(i,2) > 0.8 || posXY(i,2) < -0.8 || ...
                  posXY(i,3) > 1.6 || posXY(i,3) < -0.2
                  posXY(i,2:3) = nan(1,2);
               end
           end           
           % normalize positons 
           temp2 = ZeroToOne(posXY(:,2));temp3 = ZeroToOne(posXY(:,3));
           posXY(:,2) = temp2; posXY(:,3) = temp3; clear temp2 temp3;
           figure; plot(posXY(:,2),posXY(:,3));
           
           % positions per trial type 
           for t = 1:4
               for i = 1:length(posTrials{t}.int)
                   posTrials{t}.pos = Restrict(posXY,posTrials{t}.int); 
               end
           end
           figure;
           for t = 1:4
               subplot(2,2,t); plot(posTrials{t}.pos(:,2),posTrials{t}.pos(:,3),'.');hold on;
               title(posTrials{t}.type);%xlim([-600 650]);ylim([-200 1500]);
           end
           suptitle([animals{rat} sessions{rat}{ses}]);
           saveas(gcf,'trialTypes.fig');
           
           % linearize positions- arctang
           centers = {[0.75 0.5],[0.3 0.5],[0.75 0.5],[0.3 0.5]};
           for t = 1:4
               posTrials{t}.linpos2(:,1) = posTrials{t}.pos(:,1);
               angleTmp = atan2(posTrials{t}.pos(:,3)-centers{t}(2), posTrials{t}.pos(:,2)-centers{t}(1));
               if t == 1 || t == 4
                  posTrials{t}.linpos2(:,2) = ZeroToOne(wrapTo2Pi(-angleTmp+pi/2));  clear angleTmp;
               elseif t == 2 || t == 3
                  posTrials{t}.linpos2(:,2) = ZeroToOne(wrapTo2Pi(angleTmp-pi/2));  clear angleTmp;
               end
           end
 
           % lineariz posic - add segments - Wmaze3 (0.18 hztal for OR21)
           for t = 1 % inb left 
               posTrials{t}.linpos(:,1) = posTrials{t}.pos(:,1);
               for i = 1:length(posTrials{t}.pos)
                   if posTrials{t}.pos(i,3) > 0.13 && posTrials{t}.pos(i,2) > 0.7 % first arm
                      if posTrials{t}.pos(i,3) > 0.9 
                          posTrials{t}.pos(i,3) = 0.9;
                      end
                      posTrials{t}.linpos(i,2) = 1-posTrials{t}.pos(i,3)-0.1; %[0 0.77] ;invert and rest 0.1 to remove 
                   elseif posTrials{t}.pos(i,3) < 0.13 % hztal arm
                       if posTrials{t}.pos(i,3) > 0.95
                          posTrials{t}.pos(i,3) = 0.95;
                       elseif posTrials{t}.pos(i,2) < 0.49
                          posTrials{t}.pos(i,2) = 0.49;
                       end
                      posTrials{t}.linpos(i,2) = 1-posTrials{t}.pos(i,2)+0.77-0.05;% [0.77 1.23]
                   elseif posTrials{t}.pos(i,3) > 0.13 && posTrials{t}.pos(i,2) < 0.7 % second arm
                       if posTrials{t}.pos(i,3) > 0.9
                          posTrials{t}.pos(i,3) = 0.9;
                       end
                       posTrials{t}.linpos(i,2) = posTrials{t}.pos(i,3)+1.23-0.13; %[1.23 2]
                   end
               end 
               posTrials{t}.rewZone = [0 0.3];
               posTrials{t}.waitZone = [NaN NaN];
           end
           for t = 2 % inb right 
               posTrials{t}.linpos(:,1) = posTrials{t}.pos(:,1);
               for i = 1:length(posTrials{t}.pos)           
                   if posTrials{t}.pos(i,3) > 0.13 && posTrials{t}.pos(i,2) < 0.3 % first arm
                      if posTrials{t}.pos(i,3) > 0.9 
                          posTrials{t}.pos(i,3) = 0.9;
                      end       
                      posTrials{t}.linpos(i,2) = 1-posTrials{t}.pos(i,3)-0.1; %[0 0.77] ;invert and rest 0.1 to remove 
                   elseif posTrials{t}.pos(i,3) < 0.13 % hztal arm   
                       if posTrials{t}.pos(i,3) > 0.55
                          posTrials{t}.pos(i,3) = 0.55;
                       elseif posTrials{t}.pos(i,2) < 0.09
                          posTrials{t}.pos(i,2) = 0.09;
                       end   
                     posTrials{t}.linpos(i,2) = posTrials{t}.pos(i,2)+0.77-0.09; %[0.77 1.23]  
                   elseif posTrials{t}.pos(i,3) > 0.13 && posTrials{t}.pos(i,2) > 0.3 % second arm
                       if posTrials{t}.pos(i,3) > 0.9
                          posTrials{t}.pos(i,3) = 0.9;
                       end
                       posTrials{t}.linpos(i,2) = posTrials{t}.pos(i,3)+1.23-0.13; %[1.23 2]
                   end                       
               end
               posTrials{t}.rewZone = [0 0.3];
               posTrials{t}.waitZone = [NaN NaN];               
           end
           for t = 3 % outb left 
               posTrials{t}.linpos(:,1) = posTrials{t}.pos(:,1);
               for i = 1:length(posTrials{t}.pos)       
                    if posTrials{t}.pos(i,3) > 0.13 && posTrials{t}.pos(i,2) < 0.7 % first arm
                      if posTrials{t}.pos(i,3) > 0.9 
                          posTrials{t}.pos(i,3) = 0.9;
                      end          
                      posTrials{t}.linpos(i,2) = 1-posTrials{t}.pos(i,3)-0.1; %[0 0.77] ;invert and rest 0.1 to remove 
                   elseif posTrials{t}.pos(i,3) < 0.13 % hztal arm   
                       if posTrials{t}.pos(i,3) > 0.96
                          posTrials{t}.pos(i,3) = 0.96;
                       elseif posTrials{t}.pos(i,2) < 0.5
                          posTrials{t}.pos(i,2) = 0.5;
                       end   
                     posTrials{t}.linpos(i,2) = posTrials{t}.pos(i,2)+0.77-0.5;% [0.77 1.23]          
                   elseif posTrials{t}.pos(i,3) > 0.13 && posTrials{t}.pos(i,2) > 0.7 % second arm
                       if posTrials{t}.pos(i,3) > 0.9
                          posTrials{t}.pos(i,3) = 0.9;
                       end
                       posTrials{t}.linpos(i,2) = posTrials{t}.pos(i,3)+1.23-0.13; %[1.23 2]
                    end    
               posTrials{t}.rewZone = [NaN NaN];
               posTrials{t}.waitZone = [0 0.3];               
               end   
           end
           for t = 4 % outb right
               posTrials{t}.linpos(:,1) = posTrials{t}.pos(:,1);
               for i = 1:length(posTrials{t}.pos)                     
                    if posTrials{t}.pos(i,3) > 0.13 && posTrials{t}.pos(i,2) > 0.3 % first arm
                      if posTrials{t}.pos(i,3) > 0.9 
                          posTrials{t}.pos(i,3) = 0.9;
                      end               
                      posTrials{t}.linpos(i,2) = 1-posTrials{t}.pos(i,3)-0.1; %[0 0.77] ;invert and rest 0.1 to remove 
                    elseif posTrials{t}.pos(i,3) < 0.13 % hztal arm 
                       if posTrials{t}.pos(i,3) > 0.55
                          posTrials{t}.pos(i,3) = 0.55;
                       elseif posTrials{t}.pos(i,2) < 0.09
                          posTrials{t}.pos(i,2) = 0.09;
                       end   
                     posTrials{t}.linpos(i,2) = 1-posTrials{t}.pos(i,2)+0.77-0.45; %[0.77 1.23]  
                   elseif posTrials{t}.pos(i,3) > 0.13 && posTrials{t}.pos(i,2) < 0.3 % second arm
                       if posTrials{t}.pos(i,3) > 0.9
                          posTrials{t}.pos(i,3) = 0.9;
                       end
                       posTrials{t}.linpos(i,2) = posTrials{t}.pos(i,3)+1.23-0.13; %[1.23 2]
                    end                                            
               posTrials{t}.rewZone = [NaN NaN];
               posTrials{t}.waitZone = [0 0.3];               
               end  
           end
           
           %%% lineariz posic Kathryn method
            % ...
            
           %%% FOR WMAZE2
           for t = 1:4
               posTrials{t}.linpos = posTrials{t}.linpos2;
           end
            
           figure;
           for t = 1:4
               subplot(4,1,t); plot(posTrials{t}.linpos(:,1),posTrials{t}.linpos(:,2),'.');hold on;
               title(posTrials{t}.type);%xlim([-600 650]);ylim([-200 1500]);
           end           
           saveas(gcf,'trialTypesLin.fig');
           save([basename '.posTrials.mat'],'posTrials');
        end
    end
end

%% Get place cells
for rat = 2%:length(animals)
    for ses = 3:length(sessions{rat})
        cd([dirData animals{rat} sessions{rat}{ses}]);
        basename = bz_BasenameFromBasepath(pwd); 

        if exist([basename 'posTrials.mat'])
           load([basename '.posTrials.mat']);
           if rat == 2
           load([basename '.spikesHP.cellinfo.mat']);
           else
           load([basename '.spikes.cellinfo.mat']); 
           end
         % Calculate PCs
        for t = 1:4 % with added segments linpos
            for i = 1:length(spikes.times)
                [curve{i}{t},stats{i}{t}] = FiringCurve(posTrials{t}.linpos,spikes.times{i},...
                                   'nBins',320,'threshold',0.15,'minSize',15,'minPeak',2,'smooth',5); 
            end
        end   
        for t = 1:4 % with arctan linpos
            for i = 1:length(spikes.times)
                [curve2{i}{t},stats2{i}{t}] = FiringCurve(posTrials{t}.linpos2,spikes.times{i},...
                                   'nBins',320,'threshold',0.15,'minSize',15,'minPeak',2,'smooth',5); 
            end
        end   
        
        secondaryPF = 0.6; 
        mkdir([pwd '\PCs\']);
        for i = 1:length(spikes.times)
            figure('visible','off');
            for t = 1:4
                subplot(2,2,t);plot(curve{i}{t}.x,curve{i}{t}.rate);
                title(posTrials{t}.type);ylabel('FR(Hz)');xlabel('track (cm)');hold on;
                if t == 1 % inbL
                   set(gca,'XTick',[0 0.6 1],'XTickLabel',{'leftR','T',''});
                elseif t == 2 % inbR
                   set(gca,'XTick',[0 0.6 1],'XTickLabel',{'rightR','T',''});
                else % outb
                   set(gca,'XTick',[0 0.4 1],'XTickLabel',{'centerR','T',''});            
                end
                if stats{i}{t}.peak ~= 0
                   for k=1:length(stats{i}{t}.peak) % for each PF of this cell
                       if stats{i}{t}.peak(k)>stats{i}{t}.peak(1)*secondaryPF % destroy small secondary PFs
                           peakLoc(k)=curve{i}{t}.x(stats{i}{t}.x(k));
                           scatter(peakLoc(k),stats{i}{t}.peak(k),'MarkerFaceColor','r','LineWidth',1);
                           fieldIni(k)=curve{i}{t}.x(stats{i}{t}.fieldX(k,1));
                           fieldFin(k)=curve{i}{t}.x(stats{i}{t}.fieldX(k,2));
                           x1=fieldIni(k);y1=get(gca,'ylim');hold on; plot([x1 x1],y1,'-.k','LineWidth',1);
                           x2=fieldFin(k);y1=get(gca,'ylim');hold on; plot([x2 x2],y1,'-.k','LineWidth',1);
                       else
                           stats{i}{t}.x(k) = NaN;      stats{i}{t}.y(k) = NaN; 
                           stats{i}{t}.fieldX(k) = 0;    stats{i}{t}.size(k) = 0; 
                           stats{i}{t}.peak(k) = 0;     stats{i}{t}.mean(k) = 0;  
                           stats{i}{t}.fieldX(k,1:2) = [NaN,NaN]; stats{i}{t}.fieldY(k,1:2) = [NaN,NaN];                 
                       end
                   end
                   for k=1:length(stats{i}{t}.peak)% destroy too large PF (CHANGE for recalculate!!)
                       if stats{i}{t}.peak(k)~=0 && ... 
                               curve{i}{t}.x(stats{i}{t}.fieldX(k,2))-curve{i}{t}.x(stats{i}{t}.fieldX(k,1)) > 0.6
                           stats{i}{t}.x(k) = NaN;      stats{i}{t}.y(k) = NaN; 
                           stats{i}{t}.fieldX(k) = 0;    stats{i}{t}.size(k) = 0; 
                           stats{i}{t}.peak(k) = 0;     stats{i}{t}.mean(k) = 0;  
                           stats{i}{t}.fieldX(k,1:2) = [NaN,NaN]; stats{i}{t}.fieldY(k,1:2) = [NaN,NaN];   
                       end     
                   end
                   clear temp peakLoc fieldIni fieldFin x1 x2 y1 y2 k
                end
            end
             suptitle([basename ' cell' num2str(i)]);
             saveas(gcf,[pwd '\PCs\cell_' num2str(i) '.png'],'png'); close all;     
        end
        save([basename '.PCs.mat'],'curve','stats');
        for i = 1:length(spikes.times)
            figure('visible','off');
            for t = 1:4
                subplot(2,2,t);plot(curve2{i}{t}.x,curve2{i}{t}.rate);
                title(posTrials{t}.type);ylabel('FR(Hz)');xlabel('track (cm)');hold on;
                if t == 1 % inbL
                   set(gca,'XTick',[0 0.6 1],'XTickLabel',{'leftR','T',''});
                elseif t == 2 % inbR
                   set(gca,'XTick',[0 0.6 1],'XTickLabel',{'rightR','T',''});
                else % outb
                   set(gca,'XTick',[0 0.4 1],'XTickLabel',{'centerR','T',''});            
                end
                if stats{i}{t}.peak ~= 0
                   for k=1:length(stats2{i}{t}.peak) % for each PF of this cell
                       if stats2{i}{t}.peak(k)>stats2{i}{t}.peak(1)*secondaryPF % destroy small secondary PFs
                           peakLoc(k)=curve2{i}{t}.x(stats2{i}{t}.x(k));
                           scatter(peakLoc(k),stats2{i}{t}.peak(k),'MarkerFaceColor','r','LineWidth',1);
                           fieldIni(k)=curve2{i}{t}.x(stats2{i}{t}.fieldX(k,1));
                           fieldFin(k)=curve2{i}{t}.x(stats2{i}{t}.fieldX(k,2));
                           x1=fieldIni(k);y1=get(gca,'ylim');hold on; plot([x1 x1],y1,'-.k','LineWidth',1);
                           x2=fieldFin(k);y1=get(gca,'ylim');hold on; plot([x2 x2],y1,'-.k','LineWidth',1);
                       else
                           stats2{i}{t}.x(k) = NaN;      stats2{i}{t}.y(k) = NaN; 
                           stats2{i}{t}.fieldX(k) = 0;    stats2{i}{t}.size(k) = 0; 
                           stats2{i}{t}.peak(k) = 0;     stats2{i}{t}.mean(k) = 0;  
                           stats2{i}{t}.fieldX(k,1:2) = [NaN,NaN]; stats2{i}{t}.fieldY(k,1:2) = [NaN,NaN];                 
                       end
                   end
                   for k=1:length(stats2{i}{t}.peak)% destroy too large PF (CHANGE for recalculate!!)
                       if stats2{i}{t}.peak(k)~=0 && ... 
                               curve2{i}{t}.x(stats2{i}{t}.fieldX(k,2))-curve2{i}{t}.x(stats2{i}{t}.fieldX(k,1)) > 0.6
                           stats2{i}{t}.x(k) = NaN;      stats2{i}{t}.y(k) = NaN; 
                           stats2{i}{t}.fieldX(k) = 0;    stats2{i}{t}.size(k) = 0; 
                           stats2{i}{t}.peak(k) = 0;     stats2{i}{t}.mean(k) = 0;  
                           stats2{i}{t}.fieldX(k,1:2) = [NaN,NaN]; stats2{i}{t}.fieldY(k,1:2) = [NaN,NaN];   
                       end     
                   end
                   clear temp peakLoc fieldIni fieldFin x1 x2 y1 y2 k
                end
            end
             suptitle([basename ' cell' num2str(i)]);
             saveas(gcf,[pwd '\PCs\cell_' num2str(i) '_2.png'],'png'); close all;     
        end
        save([basename '.PCs2.mat'],'curve2','stats2');
        
        % Compile PC stats
        for i = 1:length(spikes.times)
            for t = 1:4
                if spikes.shankID(i) < 6
                if stats{i}{t}.peak(1) ~= 0 && sum(~isnan(stats{i}{t}.x)) < 3 % discard cells >2 PFs
                   sesStats.isPC(i,t) = 1;
                   sesStats.nPF(i,t) = sum(~isnan(stats{i}{t}.x));
                   sesStats.specif(i,t) = stats{i}{t}.specificity;
                   sesStats.size(i,t) = stats{i}{t}.size(1); % first PF
                   sesStats.peakFR(i,t) = stats{i}{t}.peak(1);
                   sesStats.meanFR(i,t) = stats{i}{t}.mean(1);
                   if t < 3
                      sesStats.dist2rew(i,t) = stats{i}{t}.x(1);
                   else
                      sesStats.dist2rew(i,t) = 320-stats{i}{t}.x(1);
                   end     
                   if sesStats.dist2rew(i,t) < 48 %50 before; PF close to rew (last 0.3 units)
                      sesStats.rewPF(i,t) = 1;
                      sesStats.waitPF(i,t) = NaN;
                      sesStats.tPF(i,t) = NaN;
                   elseif sesStats.dist2rew(i,t) > 272 %270 before; PF close to wait area
                      sesStats.waitPF(i,t) = 1;
                      sesStats.rewPF(i,t) = NaN;
                      sesStats.tPF(i,t) = NaN;
        %            elseif %t < 3 && sesStats.dist2rew(i,t) > 167 && sesStats.dist2rew(i,t) < 217 || ...
        %                   t > 2 && sesStats.dist2rew(i,t) > 103 && sesStats.dist2rew(i,t) < 153 % 50 around T
                   elseif sesStats.dist2rew(i,t) > 123-30 && sesStats.dist2rew(i,t) < 123+20 % 123 is end of arm
                      sesStats.tPF(i,t) = 1;
                      sesStats.waitPF(i,t) = NaN;
                      sesStats.rewPF(i,t) = NaN;                
                   else
                      sesStats.waitPF(i,t) = NaN;
                      sesStats.rewPF(i,t) = NaN;   
                      sesStats.tPF(i,t) = NaN;
                   end
                else
                   sesStats.isPC(i,t) = 0;
                   sesStats.nPF(i,t) = NaN;
                   sesStats.specif(i,t) = NaN;
                   sesStats.size(i,t) = NaN;
                   sesStats.peakFR(i,t) = NaN;
                   sesStats.meanFR(i,t) = NaN;
                   sesStats.dist2rew(i,t) = NaN;
                   sesStats.waitPF(i,t) = NaN;
                   sesStats.rewPF(i,t) = NaN; 
                   sesStats.tPF(i,t) = NaN;
                end
                end
            end
        end
        save([basename '.PCsStats.mat'],'sesStats');
        for i = 1:length(spikes.times)
            for t = 1:4
                if stats2{i}{t}.peak(1) ~= 0 && sum(~isnan(stats2{i}{t}.x)) < 3 % discard cells >2 PFs
                   sesStats2.isPC(i,t) = 1;
                   sesStats2.nPF(i,t) = sum(~isnan(stats2{i}{t}.x));
                   sesStats2.specif(i,t) = stats2{i}{t}.specificity;
                   sesStats2.size(i,t) = stats2{i}{t}.size(1); % first PF
                   sesStats2.peakFR(i,t) = stats2{i}{t}.peak(1);
                   sesStats2.meanFR(i,t) = stats2{i}{t}.mean(1);
                   if t < 3
                      sesStats2.dist2rew(i,t) = stats2{i}{t}.x(1);
                   else
                      sesStats2.dist2rew(i,t) = 320-stats2{i}{t}.x(1);
                   end     
                   if sesStats2.dist2rew(i,t) < 50 %50 before; PF close to rew (last 0.3 units)
                      sesStats2.rewPF(i,t) = 1;
                      sesStats2.waitPF(i,t) = NaN;
                      sesStats2.tPF(i,t) = NaN;
                   elseif sesStats2.dist2rew(i,t) > 270 %270 before; PF close to wait area
                      sesStats2.waitPF(i,t) = 1;
                      sesStats2.rewPF(i,t) = NaN;
                      sesStats2.tPF(i,t) = NaN;
                    elseif t < 3 && sesStats.dist2rew(i,t) > 167 && sesStats.dist2rew(i,t) < 217 || ...
                           t > 2 && sesStats.dist2rew(i,t) > 103 && sesStats.dist2rew(i,t) < 153 % 50 around T
                      sesStats2.tPF(i,t) = 1;
                      sesStats2.waitPF(i,t) = NaN;
                      sesStats2.rewPF(i,t) = NaN;                
                   else
                      sesStats2.waitPF(i,t) = NaN;
                      sesStats2.rewPF(i,t) = NaN;   
                      sesStats2.tPF(i,t) = NaN;
                   end
                else
                   sesStats2.isPC(i,t) = 0;
                   sesStats2.nPF(i,t) = NaN;
                   sesStats2.specif(i,t) = NaN;
                   sesStats2.size(i,t) = NaN;
                   sesStats2.peakFR(i,t) = NaN;
                   sesStats2.meanFR(i,t) = NaN;
                   sesStats2.dist2rew(i,t) = NaN;
                   sesStats2.waitPF(i,t) = NaN;
                   sesStats2.rewPF(i,t) = NaN; 
                   sesStats2.tPF(i,t) = NaN;
                end
            end
        end
        save([basename '.PCsStats2.mat'],'sesStats2');
        clear ripSpkT curve stats curve2 stats2 sesStats sesStats2 spikes 
       
        end
    end
end

%% Get PCs excluding ripSpk/ pulSpk
for rat = 1:length(animals)
    for ses = 1:length(sessions{rat})
        cd([dirData animals{rat} sessions{rat}{ses}]);
        basename = bz_BasenameFromBasepath(pwd); 

        if exist([basename '.posTrials.mat'])
           load([basename '.posTrials.mat']);
           
           load([basename '.spkSWRTask.mat'],'ripSpkTask');
           load([basename '.spkCLpulses.mat'],'pulSpkCL');
           if rat == 2
           load([basename '.spikesHP.cellinfo.mat']);
           else
           load([basename '.spikesHP.cellinfo.mat']);    
           end
         % cat all ripSpk & pulSpk  
           temp1 = cell2mat(reshape(ripSpkTask.timesAbs,[],1));
           temp2 = cell2mat(reshape(pulSpkCL.timesAbs,[],1));
           ripSpkT = cat(1,temp1,temp2);
           if exist([basename '.spkDpulses.mat'])
              load([basename '.spkDpulses.mat'],'pulSpkD');
              temp3 = cell2mat(reshape(pulSpkD.timesAbs,[],1));
              ripSpkT = cat(1,ripSpkT,temp3);
           end
           clear ripSpkTask pulSpkCL pulSpkD temp1 temp2 temp3
           
         % exclude ripSpk & pulSpk
         for i = 1:length(spikes.times)
             spikes.times{i} = setdiff(spikes.times{i},ripSpkT);
         end
           
         % Calculate PCs
        for t = 1:4 % with added segments linpos
            for i = 1:length(spikes.times)
                [curve{i}{t},stats{i}{t}] = FiringCurve(posTrials{t}.linpos,spikes.times{i},...
                                   'nBins',320,'threshold',0.15,'minSize',15,'minPeak',2,'smooth',5); 
            end
        end  
        try
        for t = 1:4 % with arctan linpos
            for i = 1:length(spikes.times)
                [curve2{i}{t},stats2{i}{t}] = FiringCurve(posTrials{t}.linpos2,spikes.times{i},...
                                   'nBins',320,'threshold',0.15,'minSize',15,'minPeak',2,'smooth',5); 
            end
        end
        catch
        end
        
        secondaryPF = 0.6; 
        mkdir([pwd '\PCsNoRip\']);
        for i = 1:length(spikes.times)
            figure('visible','off');
            for t = 1:4
                subplot(2,2,t);plot(curve{i}{t}.x,curve{i}{t}.rate);
                title(posTrials{t}.type);ylabel('FR(Hz)');xlabel('track (cm)');hold on;
                if t == 1 % inbL
                   set(gca,'XTick',[0 0.6 1],'XTickLabel',{'leftR','T',''});
                elseif t == 2 % inbR
                   set(gca,'XTick',[0 0.6 1],'XTickLabel',{'rightR','T',''});
                else % outb
                   set(gca,'XTick',[0 0.4 1],'XTickLabel',{'centerR','T',''});            
                end
                if stats{i}{t}.peak ~= 0
                   for k=1:length(stats{i}{t}.peak) % for each PF of this cell
                       if stats{i}{t}.peak(k)>stats{i}{t}.peak(1)*secondaryPF % destroy small secondary PFs
                           peakLoc(k)=curve{i}{t}.x(stats{i}{t}.x(k));
                           scatter(peakLoc(k),stats{i}{t}.peak(k),'MarkerFaceColor','r','LineWidth',1);
                           fieldIni(k)=curve{i}{t}.x(stats{i}{t}.fieldX(k,1));
                           fieldFin(k)=curve{i}{t}.x(stats{i}{t}.fieldX(k,2));
                           x1=fieldIni(k);y1=get(gca,'ylim');hold on; plot([x1 x1],y1,'-.k','LineWidth',1);
                           x2=fieldFin(k);y1=get(gca,'ylim');hold on; plot([x2 x2],y1,'-.k','LineWidth',1);
                       else
                           stats{i}{t}.x(k) = NaN;      stats{i}{t}.y(k) = NaN; 
                           stats{i}{t}.fieldX(k) = 0;    stats{i}{t}.size(k) = 0; 
                           stats{i}{t}.peak(k) = 0;     stats{i}{t}.mean(k) = 0;  
                           stats{i}{t}.fieldX(k,1:2) = [NaN,NaN]; stats{i}{t}.fieldY(k,1:2) = [NaN,NaN];                 
                       end
                   end
                   for k=1:length(stats{i}{t}.peak)% destroy too large PF (CHANGE for recalculate!!)
                       if stats{i}{t}.peak(k)~=0 && ... 
                               curve{i}{t}.x(stats{i}{t}.fieldX(k,2))-curve{i}{t}.x(stats{i}{t}.fieldX(k,1)) > 0.6
                           stats{i}{t}.x(k) = NaN;      stats{i}{t}.y(k) = NaN; 
                           stats{i}{t}.fieldX(k) = 0;    stats{i}{t}.size(k) = 0; 
                           stats{i}{t}.peak(k) = 0;     stats{i}{t}.mean(k) = 0;  
                           stats{i}{t}.fieldX(k,1:2) = [NaN,NaN]; stats{i}{t}.fieldY(k,1:2) = [NaN,NaN];   
                       end     
                   end
                   clear temp peakLoc fieldIni fieldFin x1 x2 y1 y2 k
                end
            end
             suptitle([basename ' cell' num2str(i)]);
             saveas(gcf,[pwd '\PCs\cell_' num2str(i) '.png'],'png'); close all;     
        end
        save([basename '.PCsNoRip.mat'],'curve','stats');
        try
        for i = 1:length(spikes.times)
            figure('visible','off');
            for t = 1:4
                subplot(2,2,t);plot(curve2{i}{t}.x,curve2{i}{t}.rate);
                title(posTrials{t}.type);ylabel('FR(Hz)');xlabel('track (cm)');hold on;
                if t == 1 % inbL
                   set(gca,'XTick',[0 0.6 1],'XTickLabel',{'leftR','T',''});
                elseif t == 2 % inbR
                   set(gca,'XTick',[0 0.6 1],'XTickLabel',{'rightR','T',''});
                else % outb
                   set(gca,'XTick',[0 0.4 1],'XTickLabel',{'centerR','T',''});            
                end
                if stats{i}{t}.peak ~= 0
                   for k=1:length(stats2{i}{t}.peak) % for each PF of this cell
                       if stats2{i}{t}.peak(k)>stats2{i}{t}.peak(1)*secondaryPF % destroy small secondary PFs
                           peakLoc(k)=curve2{i}{t}.x(stats2{i}{t}.x(k));
                           scatter(peakLoc(k),stats2{i}{t}.peak(k),'MarkerFaceColor','r','LineWidth',1);
                           fieldIni(k)=curve2{i}{t}.x(stats2{i}{t}.fieldX(k,1));
                           fieldFin(k)=curve2{i}{t}.x(stats2{i}{t}.fieldX(k,2));
                           x1=fieldIni(k);y1=get(gca,'ylim');hold on; plot([x1 x1],y1,'-.k','LineWidth',1);
                           x2=fieldFin(k);y1=get(gca,'ylim');hold on; plot([x2 x2],y1,'-.k','LineWidth',1);
                       else
                           stats2{i}{t}.x(k) = NaN;      stats2{i}{t}.y(k) = NaN; 
                           stats2{i}{t}.fieldX(k) = 0;    stats2{i}{t}.size(k) = 0; 
                           stats2{i}{t}.peak(k) = 0;     stats2{i}{t}.mean(k) = 0;  
                           stats2{i}{t}.fieldX(k,1:2) = [NaN,NaN]; stats2{i}{t}.fieldY(k,1:2) = [NaN,NaN];                 
                       end
                   end
                   for k=1:length(stats2{i}{t}.peak)% destroy too large PF (CHANGE for recalculate!!)
                       if stats2{i}{t}.peak(k)~=0 && ... 
                               curve2{i}{t}.x(stats2{i}{t}.fieldX(k,2))-curve2{i}{t}.x(stats2{i}{t}.fieldX(k,1)) > 0.6
                           stats2{i}{t}.x(k) = NaN;      stats2{i}{t}.y(k) = NaN; 
                           stats2{i}{t}.fieldX(k) = 0;    stats2{i}{t}.size(k) = 0; 
                           stats2{i}{t}.peak(k) = 0;     stats2{i}{t}.mean(k) = 0;  
                           stats2{i}{t}.fieldX(k,1:2) = [NaN,NaN]; stats2{i}{t}.fieldY(k,1:2) = [NaN,NaN];   
                       end     
                   end
                   clear temp peakLoc fieldIni fieldFin x1 x2 y1 y2 k
                end
            end
             suptitle([basename ' cell' num2str(i)]);
             saveas(gcf,[pwd '\PCs\cell_' num2str(i) '_2.png'],'png'); close all;     
        end
        save([basename '.PCsNoRip2.mat'],'curve2','stats2');
        catch
        end
        
        % Compile PC stats
        for i = 1:length(spikes.times)
            for t = 1:4
                if spikes.shankID(i) < 6
                if stats{i}{t}.peak(1) ~= 0 && sum(~isnan(stats{i}{t}.x)) < 3 % discard cells >2 PFs
                   sesStats.isPC(i,t) = 1;
                   sesStats.nPF(i,t) = sum(~isnan(stats{i}{t}.x));
                   sesStats.specif(i,t) = stats{i}{t}.specificity;
                   sesStats.size(i,t) = stats{i}{t}.size(1); % first PF
                   sesStats.peakFR(i,t) = stats{i}{t}.peak(1);
                   sesStats.meanFR(i,t) = stats{i}{t}.mean(1);
                   if t < 3
                      sesStats.dist2rew(i,t) = stats{i}{t}.x(1);
                   else
                      sesStats.dist2rew(i,t) = 320-stats{i}{t}.x(1);
                   end     
                   if sesStats.dist2rew(i,t) < 48 %50 before; PF close to rew (last 0.3 units)
                      sesStats.rewPF(i,t) = 1;
                      sesStats.waitPF(i,t) = NaN;
                      sesStats.tPF(i,t) = NaN;
                   elseif sesStats.dist2rew(i,t) > 272 %270 before; PF close to wait area
                      sesStats.waitPF(i,t) = 1;
                      sesStats.rewPF(i,t) = NaN;
                      sesStats.tPF(i,t) = NaN;
        %            elseif %t < 3 && sesStats.dist2rew(i,t) > 167 && sesStats.dist2rew(i,t) < 217 || ...
        %                   t > 2 && sesStats.dist2rew(i,t) > 103 && sesStats.dist2rew(i,t) < 153 % 50 around T
                   elseif sesStats.dist2rew(i,t) > 123-30 && sesStats.dist2rew(i,t) < 123+20 % 123 is end of arm
                      sesStats.tPF(i,t) = 1;
                      sesStats.waitPF(i,t) = NaN;
                      sesStats.rewPF(i,t) = NaN;                
                   else
                      sesStats.waitPF(i,t) = NaN;
                      sesStats.rewPF(i,t) = NaN;   
                      sesStats.tPF(i,t) = NaN;
                   end
                else
                   sesStats.isPC(i,t) = 0;
                   sesStats.nPF(i,t) = NaN;
                   sesStats.specif(i,t) = NaN;
                   sesStats.size(i,t) = NaN;
                   sesStats.peakFR(i,t) = NaN;
                   sesStats.meanFR(i,t) = NaN;
                   sesStats.dist2rew(i,t) = NaN;
                   sesStats.waitPF(i,t) = NaN;
                   sesStats.rewPF(i,t) = NaN; 
                   sesStats.tPF(i,t) = NaN;
                end
                end
            end
        end
        save([basename '.PCsStatsNoRip.mat'],'sesStats');
        try
        for i = 1:length(spikes.times)
            for t = 1:4
                if stats2{i}{t}.peak(1) ~= 0 && sum(~isnan(stats2{i}{t}.x)) < 3 % discard cells >2 PFs
                   sesStats2.isPC(i,t) = 1;
                   sesStats2.nPF(i,t) = sum(~isnan(stats2{i}{t}.x));
                   sesStats2.specif(i,t) = stats2{i}{t}.specificity;
                   sesStats2.size(i,t) = stats2{i}{t}.size(1); % first PF
                   sesStats2.peakFR(i,t) = stats2{i}{t}.peak(1);
                   sesStats2.meanFR(i,t) = stats2{i}{t}.mean(1);
                   if t < 3
                      sesStats2.dist2rew(i,t) = stats2{i}{t}.x(1);
                   else
                      sesStats2.dist2rew(i,t) = 320-stats2{i}{t}.x(1);
                   end     
                   if sesStats2.dist2rew(i,t) < 50 %50 before; PF close to rew (last 0.3 units)
                      sesStats2.rewPF(i,t) = 1;
                      sesStats2.waitPF(i,t) = NaN;
                      sesStats2.tPF(i,t) = NaN;
                   elseif sesStats2.dist2rew(i,t) > 270 %270 before; PF close to wait area
                      sesStats2.waitPF(i,t) = 1;
                      sesStats2.rewPF(i,t) = NaN;
                      sesStats2.tPF(i,t) = NaN;
                    elseif t < 3 && sesStats.dist2rew(i,t) > 167 && sesStats.dist2rew(i,t) < 217 || ...
                           t > 2 && sesStats.dist2rew(i,t) > 103 && sesStats.dist2rew(i,t) < 153 % 50 around T
                      sesStats2.tPF(i,t) = 1;
                      sesStats2.waitPF(i,t) = NaN;
                      sesStats2.rewPF(i,t) = NaN;                
                   else
                      sesStats2.waitPF(i,t) = NaN;
                      sesStats2.rewPF(i,t) = NaN;   
                      sesStats2.tPF(i,t) = NaN;
                   end
                else
                   sesStats2.isPC(i,t) = 0;
                   sesStats2.nPF(i,t) = NaN;
                   sesStats2.specif(i,t) = NaN;
                   sesStats2.size(i,t) = NaN;
                   sesStats2.peakFR(i,t) = NaN;
                   sesStats2.meanFR(i,t) = NaN;
                   sesStats2.dist2rew(i,t) = NaN;
                   sesStats2.waitPF(i,t) = NaN;
                   sesStats2.rewPF(i,t) = NaN; 
                   sesStats2.tPF(i,t) = NaN;
                end
            end
        end
        save([basename '.PCsStatsNoRip.mat'],'sesStats2');
        catch
        end
        
        clear ripSpkT curve stats curve2 stats2 sesStats sesStats2 spikes 
        disp([animals{rat} sessions{rat}{ses} ' DONE']); 
        else
        warning([animals{rat} sessions{rat}{ses} ' PROBLEM']);
        end
    end
end


%% test CL particip vs PCs
load([basename '.spkSWRTask.mat']);
load([basename '.spkCLpulses.mat']);

% clasify cells according to particip prob in pul 
frI = 5;
count1 = 0; count2 = 0;
for unit = 1:length(pul2SpkCL.FRavg)
    if pul2SpkCL.FRavg(unit) < frI && pul2SpkCL.FRall(unit) < 10 % 
        if pul2SpkCL.particip2(unit) > ripSpkTask.particip(unit)
           count1 = count1 +1;
           prefPul(count1) = unit;
        else %pul2SpkCL.particip2(unit) < ripSpkTask.particip(unit)
           count2 = count2 +1;
           prefNoPul(count2) = unit;
        end
    end
end
count1 = 0; count2 = 0;
for unit = 1:length(pul2SpkCL.FRavg)
    if pul2SpkCL.FRavg(unit) < frI && pul2SpkCL.FRall(unit) < 10
        if pul2SpkCL.particip2(unit) > pul2SpkCL.particip1(unit)
           count1 = count1 +1;
           prefPul2(count1) = unit;
        else  %pul2SpkCL.particip2(unit) < pul2SpkCL.particip1(unit)
           count2 = count2 +1;
           prefPul1(count2) = unit;
        end
    end
end
unitPref{1}=prefNoPul;unitPref{2}=prefPul;
unitPref{3}=prefPul1;unitPref{4}=prefPul2;

% % PCs
for i = 1:length(unitPref)
    propPC(i,:) = sum(sesStats2.isPC(unitPref{i},:))/length(unitPref{i});
    nPF(i,:) = nansum(sesStats2.nPF(unitPref{i},:))/length(unitPref{i});
    rewPF(i,:) = nansum(sesStats2.rewPF(unitPref{i},:))/length(unitPref{i});
    waitPF(i,:) = nansum(sesStats2.waitPF(unitPref{i},:))/length(unitPref{i});
    tPF(i,:) = nansum(sesStats2.tPF(unitPref{i},:))/length(unitPref{i});
    specif(i,:) = nanmean(sesStats2.specif(unitPref{i},:));
    dist2rew(i,:) = nanmean(sesStats2.dist2rew(unitPref{i},:));
end
figure;
subplot(1,4,1);bar(propPC);title('% PCs');set(gca,'XTickLabel',{'prefNoPul','prefPul','prefPul1','prefPul2'});
subplot(1,4,2);bar(rewPF);title('n PFs rewards');set(gca,'XTickLabel',{'prefNoPul','prefPul','prefPul1','prefPul2'});
subplot(1,4,3);bar(waitPF);title('n PFs wait area');set(gca,'XTickLabel',{'prefNoPul','prefPul','prefPul1','prefPul2'});
subplot(1,4,4);bar(tPF);title('n PFs T');set(gca,'XTickLabel',{'prefNoPul','prefPul','prefPul1','prefPul2'});
saveas(gcf,'PCsPulNoPul.fig');

% figure;
% bar(specif);title('specificity');set(gca,'XTickLabel',{'prefNoPul','prefPul','prefPul1','prefPul2'});
% figure;
% bar(dist2rew);title('dist to rew');set(gca,'XTickLabel',{'prefNoPul','prefPul','prefPul1','prefPul2'});

%% test D particip vs PCs
frI = 5;
count1 = 0; count2 = 0;
for unit = 1:length(pul2SpkD.FRavg)
    if pul2SpkCL.FRavg(unit) < frI && pul2SpkCL.FRall(unit) < 10 % 
        if pul2SpkD.particip2(unit) > 2*ripSpkTask.particip(unit)
           count1 = count1 +1;
           prefPulD(count1) = unit;
        else %pul2SpkCL.particip2(unit) < ripSpkTask.particip(unit)
           count2 = count2 +1;
           prefNoPulD(count2) = unit;
        end
    end
end
count1 = 0; count2 = 0;
for unit = 1:length(pul2SpkD.FRavg)
    if pul2SpkCL.FRavg(unit) < frI && pul2SpkCL.FRall(unit) < 10
        if pul2SpkD.particip2(unit) > 2*pul2SpkD.particip1(unit)
           count1 = count1 +1;
           prefPul2D(count1) = unit;
        else  %pul2SpkCL.particip2(unit) < pul2SpkCL.particip1(unit)
           count2 = count2 +1;
           prefPul1D(count2) = unit;
        end
    end
end
unitPrefD{1}=prefNoPulD;unitPrefD{2}=prefPulD;
unitPrefD{3}=prefPul1D;unitPrefD{4}=prefPul2D;

% % PCs
for i = 1:length(unitPrefD)
    propPCD(i,:) = sum(sesStats.isPC(unitPrefD{i},:))/length(unitPrefD{i});
    nPFD(i,:) = nansum(sesStats.nPF(unitPrefD{i},:))/length(unitPrefD{i});
    rewPFD(i,:) = nansum(sesStats.rewPF(unitPrefD{i},:))/length(unitPrefD{i});
    waitPFD(i,:) = nansum(sesStats.waitPF(unitPrefD{i},:))/length(unitPrefD{i});
    tPFD(i,:) = nansum(sesStats.tPF(unitPrefD{i},:))/length(unitPrefD{i});
    specifD(i,:) = nanmean(sesStats.specif(unitPrefD{i},:));
    dist2rewD(i,:) = nanmean(sesStats.dist2rew(unitPrefD{i},:));
end
figure;
subplot(1,4,1);bar(propPCD);title('% PCs');set(gca,'XTickLabel',{'prefNoPul','prefPul','prefPul1','prefPul2'});
subplot(1,4,2);bar(rewPFD);title('n PFs rewards');set(gca,'XTickLabel',{'prefNoPul','prefPul','prefPul1','prefPul2'});
subplot(1,4,3);bar(waitPFD);title('n PFs wait area');set(gca,'XTickLabel',{'prefNoPul','prefPul','prefPul1','prefPul2'});
subplot(1,4,4);bar(tPFD);title('n PFs T');set(gca,'XTickLabel',{'prefNoPul','prefPul','prefPul1','prefPul2'});

figure;
bar(specifD);title('specificity');set(gca,'XTickLabel',{'prefNoPul','prefPul','prefPul1','prefPul2'});
figure;
bar(dist2rewD);title('dist to rew');set(gca,'XTickLabel',{'prefNoPul','prefPul','prefPul1','prefPul2'});

