%% Recalcualte PC for all OR sessions

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'}};
sesD = [5:8 10:16];

%% Calculate PCs and stats       
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 '.spikes.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 with arctan linpos
        for t = 1:4 % 
            for i = 1:length(spikes.times)
                [curve{i}{t},stats{i}{t}] = FiringCurve(posTrials{t}.linpos2,spikes.times{i},...
                                   'nBins',320,'threshold',0.15,'minSize',15,'minPeak',3,'smooth',5); 
            end
        end   
        
        secondaryPF = 0.7; 
        mkdir([pwd '\newPCs\']);
        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
                           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.5
                           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)
                       if ~isnan(stats{i}{t}.x(k))
                       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);
                       end
                   end
                   clear temp peakLoc fieldIni fieldFin x1 x2 y1 y2 k
                end
            end
             suptitle([basename ' cell' num2str(i)]);
             saveas(gcf,[pwd '\newPCs\cell_' num2str(i) '.png'],'png'); close all;     
        end
        save([basename '.PCsNew.mat'],'curve','stats');

        % Calculate PC stats
        for i = 1:length(stats)
            % PF stats
            for t = 1:4
                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);
                   sesStats.x(i,t) = stats{i}{t}.x(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     
                   % PF rew/wait/T
                   if t<3 && stats{i}{t}.x(1) < 32 % 0.2 units, rew only in inb
                      sesStats.rewPF(i,t) = 1;
                      sesStats.waitPF(i,t) = NaN;
                      sesStats.tPF(i,t) = NaN;
                   elseif t>2 && stats{i}{t}.x(1) < 32 % 0.2 units, wait only in outb
                      sesStats.waitPF(i,t) = 1;
                      sesStats.rewPF(i,t) = NaN;
                      sesStats.tPF(i,t) = NaN;
                    elseif t < 3 && stats{i}{t}.x(1) > 192-30 && stats{i}{t}.x(1) < 192+20 || ...
                           t > 2 && stats{i}{t}.x(1) > 128-20 && stats{i}{t}.x(1) < 128+30 % +-25 around T
                      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;
                   sesStats.x(i,t) = NaN;
                end
            end
        
            % PF peak to 2D
            for t = 1 % inb left
                if stats{i}{t}.peak(1) ~= 0
                   if stats{i}{t}.x(1) < 128 % first arm
                      sesStats.x2D(i,t) = 0.9-0.05 + (0.1)*rand(1);
                      lin = linspace(0.9,0.1,128);
                      sesStats.y2D(i,t) = lin(stats{i}{t}.x(1));
                   elseif stats{i}{t}.x(1) >= 128 && stats{i}{t}.x(1) <= 192 % hztal arm
                      lin = linspace(0.5,0.9,64);
                      if stats{i}{t}.x(1) == 128;stats{i}{t}.x(1) = 129; end
                      sesStats.x2D(i,t) = lin(stats{i}{t}.x(1)-128);
                      sesStats.y2D(i,t) = 0.1-0.05 + (0.1)*rand(1);
                   elseif stats{i}{t}.x(1) > 192 % center arm
                      sesStats.x2D(i,t) = 0.5-0.05 + (0.1)*rand(1);
                      lin = linspace(0.1,0.9,128);
                      sesStats.y2D(i,t) = lin(stats{i}{t}.x(1)-192);
                   end
                else
                      sesStats.x2D(i,t) = NaN; sesStats.y2D(i,t) = NaN; 
                end
            end
            for t = 2 % inb right
                if stats{i}{t}.peak(1) ~= 0
                   if stats{i}{t}.x(1) < 128 % first arm
                      sesStats.x2D(i,t) = 0.1-0.05 + (0.1)*rand(1);
                      lin = linspace(0.9,0.1,128);
                      sesStats.y2D(i,t) = lin(stats{i}{t}.x(1));
                   elseif stats{i}{t}.x(1) >= 128 && stats{i}{t}.x(1) <= 192 % hztal arm
                      lin = linspace(0.1,0.5,64);
                      if stats{i}{t}.x(1) == 128;stats{i}{t}.x(1) = 129; end
                      sesStats.x2D(i,t) = lin(stats{i}{t}.x(1)-128);
                      sesStats.y2D(i,t) = 0.1-0.05 + (0.1)*rand(1);
                   elseif stats{i}{t}.x(1) > 192 % center arm
                      sesStats.x2D(i,t) = 0.5-0.05 + (0.1)*rand(1);
                      lin = linspace(0.1,0.9,128);
                      sesStats.y2D(i,t) = lin(stats{i}{t}.x(1)-192);
                   end
                else
                      sesStats.x2D(i,t) = NaN; sesStats.y2D(i,t) = NaN; 
                end
            end    
            for t = 3 % outb left
                if stats{i}{t}.peak(1) ~= 0
                   if stats{i}{t}.x(1) < 128 % first arm
                      sesStats.x2D(i,t) = 0.5-0.05 + (0.1)*rand(1);
                      lin = linspace(0.9,0.1,128);
                      sesStats.y2D(i,t) = lin(stats{i}{t}.x(1));
                   elseif stats{i}{t}.x(1) >= 128 && stats{i}{t}.x(1) <= 192 % hztal arm
                      lin = linspace(0.5,0.9,64);
                      if stats{i}{t}.x(1) == 128;stats{i}{t}.x(1) = 129; end
                      sesStats.x2D(i,t) = lin(stats{i}{t}.x(1)-128);
                      sesStats.y2D(i,t) = 0.1-0.05 + (0.1)*rand(1);
                   elseif stats{i}{t}.x(1) > 192 % center arm
                      sesStats.x2D(i,t) = 0.9-0.05 + (0.1)*rand(1);
                      lin = linspace(0.1,0.9,128);
                      sesStats.y2D(i,t) = lin(stats{i}{t}.x(1)-192);
                   end
                else
                      sesStats.x2D(i,t) = NaN; sesStats.y2D(i,t) = NaN; 
                end
            end            
            for t = 4 % outb right
                if stats{i}{t}.peak(1) ~= 0
                   if stats{i}{t}.x(1) < 128 % first arm
                      sesStats.x2D(i,t) = 0.5-0.05 + (0.1)*rand(1);
                      lin = linspace(0.9,0.1,128);
                      sesStats.y2D(i,t) = lin(stats{i}{t}.x(1));
                   elseif stats{i}{t}.x(1) >= 128 && stats{i}{t}.x(1) <= 192 % hztal arm
                      lin = linspace(0.1,0.5,64);
                      if stats{i}{t}.x(1) == 128;stats{i}{t}.x(1) = 129; end
                      sesStats.x2D(i,t) = lin(stats{i}{t}.x(1)-128);
                      sesStats.y2D(i,t) = 0.1-0.05 + (0.1)*rand(1);
                   elseif stats{i}{t}.x(1) > 192 % center arm
                      sesStats.x2D(i,t) = 0.1-0.05 + (0.1)*rand(1);
                      lin = linspace(0.1,0.9,128);
                      sesStats.y2D(i,t) = lin(stats{i}{t}.x(1)-192);
                   end
                else
                      sesStats.x2D(i,t) = NaN; sesStats.y2D(i,t) = NaN; 
                end
            end             
        end 
        save([basename '.PCsStatsNew.mat'],'sesStats');
       
        clear ripSpkT curve stats sesStats spikes 
        disp([animals{rat} sessions{rat}{ses} ' DONE']); 
        end
    end
end

%% save isPC for each session
for rat = 1:length(animals)
    for ses = 1:length(sessions{rat})
        cd([dirData animals{rat} sessions{rat}{ses}]);
        basename = bz_BasenameFromBasepath(pwd); 
        load([basename '.PCsStatsNew.mat']);
        load([basename '.cellMetrics.mat']);
        
        for i = 1:length(sesStats.isPC)
            if sum(sesStats.isPC(i,:)) > 0
               cellMetrics.isPC(i,1) = 1;
            else
               cellMetrics.isPC(i,1) = 0;
            end
        end
        save([basename '.cellMetrics.mat'],'cellMetrics');
        clear sesStats cellMetrics
    end
end

%% Pool all PCs
 sesPool = [2:5 10:12 16];
 sesDPool = [];
 count = 0;
for rat = 1:length(animals)
    for ses = 1:length(sessions{rat})
        cd([dirData animals{rat} sessions{rat}{ses}]);
        basename = bz_BasenameFromBasepath(pwd);     
        count = count +1;
        
        load([basename '.PCsStatsNew.mat'],'sesStats');
        sesPCstats{count} = orderfields(sesStats);
    end
end
   countD = 0;
   for i = 1:length(sesD)
       ii=sesD(i);
       countD = countD +1;
       sesDPCstats{countD} = orderfields(sesPCstats{ii});
   end
sesDPCstats = cell2mat(sesDPCstats);
allDPCstats = bz_CollapseStruct(sesDPCstats,'match');  
sesPCstats = cell2mat(sesPCstats(sesPool));
allPCstats = bz_CollapseStruct(sesPCstats,'match');  

for i = 1:length(allPCstats.isPC)
    if sum(allPCstats.isPC(i,:)) > 0
        isPC(i) = 1;
    else 
        isPC(i) = 0;
    end
end
for i = 1:length(allDPCstats.isPC)
    if sum(allDPCstats.isPC(i,:)) > 0
        isPCD(i) = 1;
    else 
        isPCD(i) = 0;
    end
end

%save('F:\ORproject\allPCstats.mat','allPCstats','allDPCstats','isPC','isPCD');
%%

