% LFP ANALYSIS AND LAYER DEFINITION


% ndnf lfp responses
analogChannelsList = 1;
pulses = bz_getAnalogPulses('analogCh',analogChannelsList,'manualThr',false);
xml = LoadParameters;       
st = pulses.timestamps(pulses.analogChannel == analogChannelsList(1),1);

% CSD FOR THE LFP RESPONSE
mm = 1;
twin = [0.1]; % time windows
excludeShanks = [];
shanks = xml.AnatGrps;
shanks(excludeShanks) = [];
figure
set(gcf,'Position',[100 100 1400 600])
for jj = 1:length(shanks)
    lfp = bz_GetLFP(shanks(jj).Channels,'noPrompts', true);
    
    [csd,lfpAvg] = bz_eventCSD(lfp,st,'twin',[twin twin],'plotLFP',false,'plotCSD',false);
    taxis = linspace(-twin,twin,size(csd.data,1));
    cmax = max(max(csd.data)); 
    subplot(1,size(xml.AnatGrps,2),jj);
    contourf(taxis,1:size(csd.data,2),csd.data',40,'LineColor','none');hold on;
    set(gca,'YDir','reverse'); xlabel('time (s)'); ylabel('channel'); title(strcat('AnalogCh #',num2str(analogChannelsList(mm))),'FontWeight','normal'); 
    colormap jet; try caxis([-cmax cmax]); end
    hold on
    for kk = 1:size(lfpAvg.data,2)
        plot(taxis,(lfpAvg.data(:,kk)/100)+kk-1,'k')
    end
end
csd_pulses = csd;
csd_pulses.taxis = taxis;
xlim([-.01 0.1]);
caxis([-3 3])
saveas(gcf,['SummaryFigures\analogPulsesCSD_ch',num2str(analogChannelsList(mm)), '.png']);

% CSD FOR SLOWWAVES
disp('Slow-waves CSD and PSTH...');
UDStates = detectUD;
xml = LoadParameters;
shanks = xml.AnatGrps;
shanks(excludeShanks) = [];
% CSD
twin = 0.2;
evs = UDStates.timestamps.DOWN(:,1);
figure
set(gcf,'Position',[100 100 1400 600])
for jj = 1:size(shanks,2)
    lfp = bz_GetLFP(shanks(jj).Channels,'noPrompts', true);
    [csd,lfpAvg] = bz_eventCSD(lfp,evs,'twin',[twin twin],'plotLFP',false,'plotCSD',false);
    taxis = linspace(-twin,twin,size(csd.data,1));
    cmax = max(max(csd.data)); 
    subplot(1,size(shanks,2),jj);
    contourf(taxis,1:size(csd.data,2),csd.data',40,'LineColor','none');hold on;
    set(gca,'YDir','reverse'); xlabel('time (s)'); ylabel('channel'); title(strcat('DOWN-UP, Shank #',num2str(jj)),'FontWeight','normal'); 
    colormap jet; caxis([-cmax cmax]);
    hold on
    for kk = 1:size(lfpAvg.data,2)
        plot(taxis,(lfpAvg.data(:,kk)/100)+kk-1,'k')
    end
end
csd_sw = csd;
csd_sw.taxis = taxis;
saveas(gcf,'SummaryFigures\downUpStatesCSD.png');

% CELL SPATIAL DISTRIBUTION
targetFile = dir('*cell_metrics.cellinfo.mat'); load(targetFile.name);
figure
subplot(1,3,1)
contourf(csd_pulses.taxis,1:size(csd_pulses.data,2),csd_pulses.data',40,'LineColor','none');
set(gca,'YDir','reverse'); xlabel('time (s)'); 
colormap jet; caxis([-cmax cmax]);
xlim([-.01 0.1]);
ylabel('Channels');
title('Light pulses');

subplot(1,3,2)
contourf(csd_sw.taxis,1:size(csd_sw.data,2),csd_sw.data',40,'LineColor','none');
set(gca,'YDir','reverse'); xlabel('time (s)'); 
colormap jet; caxis([-cmax cmax]);
ylabel('Channels');
title('Slow waves');

subplot(1,3,3)
[N, edges] = histcounts(cell_metrics.maxWaveformChannelOrder,[1:5:64]);
centers = edges(1:end-1) + mean(diff(edges))/2;
plot(N, centers,'k');
set(gca,'YDir','reverse'); xlabel('time (s)'); 
ylim([1 64]);
title('Cell distribution');
cellDephDist.N = N;
cellDephDist.centers = centers;

disp('Specify layer 3/4 border... some description');
[pos] = ginput(1);
hold on
plot([0 max(N)], pos([2 2]),'-r');
border34 = pos(2);

disp('Specify layer 5/6 border... some description');
[pos] = ginput(1);
hold on
plot([0 max(N)], pos([2 2]),'-g');
border56 = pos(2);

disp('Specify layer 1/2 border... some description');
[pos] = ginput(1);
hold on
plot([0 max(N)], pos([2 2]),'-b');
border12 = pos(2);

disp('Specify end of layer 6... some description');
[pos] = ginput(1);
hold on
plot([0 max(N)], pos([2 2]),'--k');
borderEnd = pos(2);

layer1Cells = find(cell_metrics.maxWaveformChannelOrder<border12);
layer23Cells = find(cell_metrics.maxWaveformChannelOrder>border12 & cell_metrics.maxWaveformChannelOrder<border34);
layer5Cells = find(cell_metrics.maxWaveformChannelOrder>border34 & cell_metrics.maxWaveformChannelOrder<border56);
layer6Cells = find(cell_metrics.maxWaveformChannelOrder>border56 & cell_metrics.maxWaveformChannelOrder<borderEnd);

layerDefinition.layer1Cells = layer1Cells;
layerDefinition.layer23Cells = layer23Cells;
layerDefinition.layer5Cells = layer5Cells;
layerDefinition.layer6Cells = layer6Cells;
layerDefinition.border34 = border34;
layerDefinition.border56 = border56;
layerDefinition.border12 = border12;
layerDefinition.borderEnd = borderEnd;
layerDefinition.csd_slowWaves = csd_sw;
layerDefinition.csd_pulses = csd_pulses;
layerDefinition.cellDephDist = cellDephDist;

filename = split(pwd,filesep); filename = filename{end};
save([filename '.layerDefinition.cellinfo.mat'],'layerDefinition');








