function peri_event_spikePSTH(basepath, varargin)


% INPUTS:
%       basepath          (required) -  path to data
%       events            (required) -  events timestamps (in sec)
%                                       containing events.peaks (in sec)
%

%% Defaults and Params
isevents = @(x) sum(strcmp(x,{'ripple','IED','DS2'}));

p = inputParser;
addParameter(p,'eventType','ripple',isevents)

parse(p,varargin{:});
eventType = p.Results.eventType;

basename = bz_BasenameFromBasepath(basepath);
% load .Tracking.Behavior.mat and .ripples.events.mat from session folder, and then run
if ~isempty(dir('*.cell_metrics.cellinfo.mat'))
    cellMeticsFile = dir('*.cell_metrics.cellinfo.mat');
    load(cellMeticsFile.name);
else
    error('Needs Cell_metrics file! Run ProcessCellMetrics.m')
end

switch eventType
    case 'ripple'
        if ~isempty(dir('*.ripples.events.mat'))
            ripplesFile = dir('*.ripples.events.mat');
            load(ripplesFile.name);
        else
            error('Needs ripples file!')
        end
        events = ripples.peaks;
    case 'IED'
        if ~isempty(dir('*.IEDclean.events.mat'))
            ripplesFile = dir('*.IEDclean.events.mat');
            load(ripplesFile.name);
        else
            error('Needs IED file!')
        end
        events = IED.peaks;
    case 'DS2'
        if ~isempty(dir('*.DS2.events.mat'))
            dentateSpikesFile = dir('*.DS2.events.mat');
            load(dentateSpikesFile.name);
        else
            error('Needs DS2 file!')
        end
        events = DS2.peaks;
end


if ~isempty(dir('*.spikes.cellinfo.mat'))
    spikesFile = dir('*.spikes.cellinfo.mat');
    load(spikesFile.name);
else
    error('Needs spikes.cellinfo file! Run ProcessCellMetrics.m')
end

%% Example peri-SPW-R firing histograms of CA1-CA2-CA3a, b, and c pyramidal cells

Pyr_ID = find(contains(cell_metrics.putativeCellType,'Pyramidal Cell') == 1);
CA3_Pyr = intersect(cell_metrics.tags.CA3,Pyr_ID);  %find pyramidal cells in CA3
%[~,id] = intersect(CA3_Pyr,cell_metrics.tags.Bad);      %find bad cells among CA3 pyramidal cells
%CA3_Pyr(id) = [];   %delete bad CA3 pyramidal cells
CA2_Pyr = intersect(cell_metrics.tags.CA2,Pyr_ID);  %find pyramidal cells in CA2
%[~,id] = intersect(CA2_Pyr,cell_metrics.tags.Bad);      %find bad cells among CA2 pyramidal cells
%CA2_Pyr(id) = [];   %delete bad CA2 pyramidal cells
CA1_Pyr = intersect(cell_metrics.tags.CA1,Pyr_ID);  %find pyramidal cells in CA1
%[~,id] = intersect(CA1_Pyr,cell_metrics.tags.Bad);      %find bad cells among CA1 pyramidal cells
%CA1_Pyr(id) = [];   %delete bad CA1 pyramidal cells
DG_Pyr = [intersect(cell_metrics.tags.DG,Pyr_ID) intersect(cell_metrics.tags.hilus,Pyr_ID)];  %find pyramidal cells in DG
%[~,id] = intersect(DG_Pyr,cell_metrics.tags.Bad);      %find bad cells among DG pyramidal cells
%DG_Pyr(id) = [];   %delete bad CA1 pyramidal cells
ctx_Pyr = intersect(cell_metrics.tags.cortex,Pyr_ID);  %find pyramidal cells in ctx
%[~,id] = intersect(ctx_Pyr,cell_metrics.tags.Bad);      %find bad cells among ctx pyramidal cells
%ctx_Pyr(id) = [];   %delete bad ctx pyramidal cells


PSTH_out_CA3 = [];
PSTH_out_CA2 = [];
PSTH_out_CA1 = [];
PSTH_out_DG = [];
PSTH_out_ctx = [];

event_times = events;
duration = 0.2;
padding = 0;
binCount = 100;
binSize = max(duration/binCount*1000,0.5)/1000; % minimum binsize is 0.5ms.
binsToKeep = 1:duration*2/binSize;

for ii = 1:length(CA3_Pyr)
    spikesTemp = cell_metrics.spikes.times{1,CA3_Pyr(ii)};
    % Calculating PSTH
    sr = 20000; % Sampling rate
    [spike_times,index] = sort([spikesTemp;event_times(:)]);
    spike_cluster_index = [ones(size(spikesTemp));2*ones(size(event_times(:)))];
    [ccg,~] = CCG(spike_times,spike_cluster_index(index),'binSize',binSize,'duration',(duration+padding)*2,'Fs',1/sr);
    PSTH_out_CA3(:,ii) = ccg(binsToKeep+1,2,1)./numel(event_times)/binSize;
end
clear ii

for ii = 1:length(CA2_Pyr)
    spikesTemp = cell_metrics.spikes.times{1,CA2_Pyr(ii)};
    % Calculating PSTH
    sr = 20000; % Sampling rate
    [spike_times,index] = sort([spikesTemp;event_times(:)]);
    spike_cluster_index = [ones(size(spikesTemp));2*ones(size(event_times(:)))];
    [ccg,~] = CCG(spike_times,spike_cluster_index(index),'binSize',binSize,'duration',(duration+padding)*2,'Fs',1/sr);
    PSTH_out_CA2(:,ii) = ccg(binsToKeep+1,2,1)./numel(event_times)/binSize;
end
clear ii

for ii = 1:length(CA1_Pyr)
    spikesTemp = cell_metrics.spikes.times{1,CA1_Pyr(ii)};
    % Calculating PSTH
    sr = 20000; % Sampling rate
    [spike_times,index] = sort([spikesTemp;event_times(:)]);
    spike_cluster_index = [ones(size(spikesTemp));2*ones(size(event_times(:)))];
    [ccg,~] = CCG(spike_times,spike_cluster_index(index),'binSize',binSize,'duration',(duration+padding)*2,'Fs',1/sr);
    PSTH_out_CA1(:,ii) = ccg(binsToKeep+1,2,1)./numel(event_times)/binSize;
end
clear ii

for ii = 1:length(DG_Pyr)
    spikesTemp = cell_metrics.spikes.times{1,DG_Pyr(ii)};
    % Calculating PSTH
    sr = 20000; % Sampling rate
    [spike_times,index] = sort([spikesTemp;event_times(:)]);
    spike_cluster_index = [ones(size(spikesTemp));2*ones(size(event_times(:)))];
    [ccg,~] = CCG(spike_times,spike_cluster_index(index),'binSize',binSize,'duration',(duration+padding)*2,'Fs',1/sr);
    PSTH_out_DG(:,ii) = ccg(binsToKeep+1,2,1)./numel(event_times)/binSize;
end
clear ii

for ii = 1:length(ctx_Pyr)
    spikesTemp = cell_metrics.spikes.times{1,ctx_Pyr(ii)};
    % Calculating PSTH
    sr = 20000; % Sampling rate
    [spike_times,index] = sort([spikesTemp;event_times(:)]);
    spike_cluster_index = [ones(size(spikesTemp));2*ones(size(event_times(:)))];
    [ccg,~] = CCG(spike_times,spike_cluster_index(index),'binSize',binSize,'duration',(duration+padding)*2,'Fs',1/sr);
    PSTH_out_ctx(:,ii) = ccg(binsToKeep+1,2,1)./numel(event_times)/binSize;
end
clear ii

figure
subplot(4,5,1)
bar(mean(PSTH_out_CA3,2))
ylimit = get(gca,'ylim');
hold on
plot([100 100],[ylimit(1) ylimit(2)],'k--')
hold off
box off
set(gca,'xtick',[0 100 200],'xticklabel',[-100 0 100])
%xlabel('Time lag (ms)')
title(['CA3 n=' num2str(size(PSTH_out_CA3,2))])

subplot(4,5,6)
plot(PSTH_out_CA3,'color',[0.7 0.7 0.7 0.4])
ylimit = get(gca,'ylim');
hold on
plot([100 100],[ylimit(1) ylimit(2)],'k--')
hold off
box off
set(gca,'xtick',[0 100 200],'xticklabel',[-100 0 100])
%xlabel('Time lag (ms)')
%title('CA3')

subplot(4,5,2)
bar(mean(PSTH_out_CA2,2))
ylimit = get(gca,'ylim');
hold on
plot([100 100],[ylimit(1) ylimit(2)],'k--')
hold off
box off
set(gca,'xtick',[0 100 200],'xticklabel',[-100 0 100])
%xlabel('Time lag (ms)')
title(['CA2 n=' num2str(size(PSTH_out_CA2,2))])

subplot(4,5,7)
plot(PSTH_out_CA2,'color',[0.7 0.7 0.7 0.4])
ylimit = get(gca,'ylim');
hold on
plot([100 100],[ylimit(1) ylimit(2)],'k--')
hold off
box off
set(gca,'xtick',[0 100 200],'xticklabel',[-100 0 100])
%xlabel('Time lag (ms)')
%title('CA2')

subplot(4,5,3)
bar(mean(PSTH_out_CA1,2))
ylimit = get(gca,'ylim');
hold on
plot([100 100],[ylimit(1) ylimit(2)],'k--')
hold off
box off
set(gca,'xtick',[0 100 200],'xticklabel',[-100 0 100])
%xlabel('Time lag (ms)')
title(['CA1 n=' num2str(size(PSTH_out_CA1,2))])

subplot(4,5,8)
plot(PSTH_out_CA1,'color',[0.7 0.7 0.7 0.4])
ylimit = get(gca,'ylim');
hold on
plot([100 100],[ylimit(1) ylimit(2)],'k--')
hold off
box off
set(gca,'xtick',[0 100 200],'xticklabel',[-100 0 100])
%xlabel('Time lag (ms)')
%title('CA1')

subplot(4,5,4)
bar(mean(PSTH_out_DG,2))
ylimit = get(gca,'ylim');
hold on
plot([100 100],[ylimit(1) ylimit(2)],'k--')
hold off
box off
set(gca,'xtick',[0 100 200],'xticklabel',[-100 0 100])
%xlabel('Time lag (ms)')
title(['DG n=' num2str(size(PSTH_out_DG,2))])

subplot(4,5,9)
plot(PSTH_out_DG,'color',[0.7 0.7 0.7 0.4])
ylimit = get(gca,'ylim');
hold on
plot([100 100],[ylimit(1) ylimit(2)],'k--')
hold off
box off
set(gca,'xtick',[0 100 200],'xticklabel',[-100 0 100])
%xlabel('Time lag (ms)')
%title('DG')

subplot(4,5,5)
bar(mean(PSTH_out_ctx,2))
ylimit = get(gca,'ylim');
hold on
plot([100 100],[ylimit(1) ylimit(2)],'k--')
hold off
box off
set(gca,'xtick',[0 100 200],'xticklabel',[-100 0 100])
%xlabel('Time lag (ms)')
title(['Ctx n=' num2str(size(PSTH_out_ctx,2))])

subplot(4,5,10)
plot(PSTH_out_ctx,'color',[0.7 0.7 0.7 0.4])
ylimit = get(gca,'ylim');
hold on
plot([100 100],[ylimit(1) ylimit(2)],'k--')
hold off
box off
set(gca,'xtick',[0 100 200],'xticklabel',[-100 0 100])
xlabel('Time lag (ms)')
%title('Ctx')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
INT_ID = find(contains(cell_metrics.putativeCellType,'Int') == 1);
CA3_INT = intersect(cell_metrics.tags.CA3,INT_ID);  %find pyramidal cells in CA3
%[~,id] = intersect(CA3_INT,cell_metrics.tags.Bad);      %find bad cells among CA3 pyramidal cells
%CA3_INT(id) = [];   %delete bad CA3 pyramidal cells
CA2_INT = intersect(cell_metrics.tags.CA2,INT_ID);  %find pyramidal cells in CA2
%[~,id] = intersect(CA2_INT,cell_metrics.tags.Bad);      %find bad cells among CA2 pyramidal cells
%CA2_INT(id) = [];   %delete bad CA2 pyramidal cells
CA1_INT = intersect(cell_metrics.tags.CA1,INT_ID);  %find pyramidal cells in CA1
%[~,id] = intersect(CA1_INT,cell_metrics.tags.Bad);      %find bad cells among CA1 pyramidal cells
%CA1_INT(id) = [];   %delete bad CA1 pyramidal cells
DG_INT = [intersect(cell_metrics.tags.DG,INT_ID) intersect(cell_metrics.tags.hilus,INT_ID)];  %find pyramidal cells in DG
%[~,id] = intersect(DG_INT,cell_metrics.tags.Bad);      %find bad cells among DG pyramidal cells
%DG_INT(id) = [];   %delete bad CA1 pyramidal cells
ctx_INT = intersect(cell_metrics.tags.cortex,INT_ID);  %find pyramidal cells in ctx
%[~,id] = intersect(ctx_INT,cell_metrics.tags.Bad);      %find bad cells among ctx pyramidal cells
%ctx_INT(id) = [];   %delete bad ctx pyramidal cells


PSTH_out_CA3_INT = [];
PSTH_out_CA2_INT = [];
PSTH_out_CA1_INT = [];
PSTH_out_DG_INT = [];
PSTH_out_ctx_INT = [];

for ii = 1:length(CA3_INT)
    spikesTemp = cell_metrics.spikes.times{1,CA3_INT(ii)};
    % Calculating PSTH
    sr = 20000; % Sampling rate
    [spike_times,index] = sort([spikesTemp;event_times(:)]);
    spike_cluster_index = [ones(size(spikesTemp));2*ones(size(event_times(:)))];
    [ccg,~] = CCG(spike_times,spike_cluster_index(index),'binSize',binSize,'duration',(duration+padding)*2,'Fs',1/sr);
    PSTH_out_CA3_INT(:,ii) = ccg(binsToKeep+1,2,1)./numel(event_times)/binSize;
end
clear ii

for ii = 1:length(CA2_INT)
    spikesTemp = cell_metrics.spikes.times{1,CA2_INT(ii)};
    % Calculating PSTH
    sr = 20000; % Sampling rate
    [spike_times,index] = sort([spikesTemp;event_times(:)]);
    spike_cluster_index = [ones(size(spikesTemp));2*ones(size(event_times(:)))];
    [ccg,~] = CCG(spike_times,spike_cluster_index(index),'binSize',binSize,'duration',(duration+padding)*2,'Fs',1/sr);
    PSTH_out_CA2_INT(:,ii) = ccg(binsToKeep+1,2,1)./numel(event_times)/binSize;
end
clear ii

for ii = 1:length(CA1_INT)
    spikesTemp = cell_metrics.spikes.times{1,CA1_INT(ii)};
    % Calculating PSTH
    sr = 20000; % Sampling rate
    [spike_times,index] = sort([spikesTemp;event_times(:)]);
    spike_cluster_index = [ones(size(spikesTemp));2*ones(size(event_times(:)))];
    [ccg,~] = CCG(spike_times,spike_cluster_index(index),'binSize',binSize,'duration',(duration+padding)*2,'Fs',1/sr);
    PSTH_out_CA1_INT(:,ii) = ccg(binsToKeep+1,2,1)./numel(event_times)/binSize;
end
clear ii

for ii = 1:length(DG_INT)
    spikesTemp = cell_metrics.spikes.times{1,DG_INT(ii)};
    % Calculating PSTH
    sr = 20000; % Sampling rate
    [spike_times,index] = sort([spikesTemp;event_times(:)]);
    spike_cluster_index = [ones(size(spikesTemp));2*ones(size(event_times(:)))];
    [ccg,~] = CCG(spike_times,spike_cluster_index(index),'binSize',binSize,'duration',(duration+padding)*2,'Fs',1/sr);
    PSTH_out_DG_INT(:,ii) = ccg(binsToKeep+1,2,1)./numel(event_times)/binSize;
end
clear ii

for ii = 1:length(ctx_INT)
    spikesTemp = cell_metrics.spikes.times{1,ctx_INT(ii)};
    % Calculating PSTH
    sr = 20000; % Sampling rate
    [spike_times,index] = sort([spikesTemp;event_times(:)]);
    spike_cluster_index = [ones(size(spikesTemp));2*ones(size(event_times(:)))];
    [ccg,~] = CCG(spike_times,spike_cluster_index(index),'binSize',binSize,'duration',(duration+padding)*2,'Fs',1/sr);
    PSTH_out_ctx_INT(:,ii) = ccg(binsToKeep+1,2,1)./numel(event_times)/binSize;
end
clear ii

subplot(4,5,11)
bar(mean(PSTH_out_CA3_INT,2))
ylimit = get(gca,'ylim');
hold on
plot([100 100],[ylimit(1) ylimit(2)],'k--')
hold off
box off
set(gca,'xtick',[0 100 200],'xticklabel',[-100 0 100])
xlabel('Time lag (ms)')
title(['CA3 n=' num2str(size(PSTH_out_CA3_INT,2))])

subplot(4,5,16)
plot(PSTH_out_CA3_INT,'color',[0.7 0.7 0.7 0.4])
ylimit = get(gca,'ylim');
hold on
plot([100 100],[ylimit(1) ylimit(2)],'k--')
hold off
box off
set(gca,'xtick',[0 100 200],'xticklabel',[-100 0 100])
xlabel('Time lag (ms)')
title('CA3')

subplot(4,5,12)
bar(mean(PSTH_out_CA2_INT,2))
ylimit = get(gca,'ylim');
hold on
plot([100 100],[ylimit(1) ylimit(2)],'k--')
hold off
box off
set(gca,'xtick',[0 100 200],'xticklabel',[-100 0 100])
xlabel('Time lag (ms)')
title(['CA2 n=' num2str(size(PSTH_out_CA2_INT,2))])

subplot(4,5,17)
plot(PSTH_out_CA2_INT,'color',[0.7 0.7 0.7 0.4])
ylimit = get(gca,'ylim');
hold on
plot([100 100],[ylimit(1) ylimit(2)],'k--')
hold off
box off
set(gca,'xtick',[0 100 200],'xticklabel',[-100 0 100])
xlabel('Time lag (ms)')
title('CA2')

subplot(4,5,13)
bar(mean(PSTH_out_CA1_INT,2))
ylimit = get(gca,'ylim');
hold on
plot([100 100],[ylimit(1) ylimit(2)],'k--')
hold off
box off
set(gca,'xtick',[0 100 200],'xticklabel',[-100 0 100])
xlabel('Time lag (ms)')
title(['CA1 n=' num2str(size(PSTH_out_CA1_INT,2))])

subplot(4,5,18)
plot(PSTH_out_CA1_INT,'color',[0.7 0.7 0.7 0.4])
ylimit = get(gca,'ylim');
hold on
plot([100 100],[ylimit(1) ylimit(2)],'k--')
hold off
box off
set(gca,'xtick',[0 100 200],'xticklabel',[-100 0 100])
xlabel('Time lag (ms)')
title('CA1')

subplot(4,5,14)
bar(mean(PSTH_out_DG_INT,2))
ylimit = get(gca,'ylim');
hold on
plot([100 100],[ylimit(1) ylimit(2)],'k--')
hold off
box off
set(gca,'xtick',[0 100 200],'xticklabel',[-100 0 100])
xlabel('Time lag (ms)')
title(['DG n=' num2str(size(PSTH_out_DG_INT,2))])

subplot(4,5,19)
plot(PSTH_out_DG_INT,'color',[0.7 0.7 0.7 0.4])
ylimit = get(gca,'ylim');
hold on
plot([100 100],[ylimit(1) ylimit(2)],'k--')
hold off
box off
set(gca,'xtick',[0 100 200],'xticklabel',[-100 0 100])
xlabel('Time lag (ms)')
title('DG')

subplot(4,5,15)
bar(mean(PSTH_out_ctx_INT,2))
ylimit = get(gca,'ylim');
hold on
plot([100 100],[ylimit(1) ylimit(2)],'k--')
hold off
box off
set(gca,'xtick',[0 100 200],'xticklabel',[-100 0 100])
xlabel('Time lag (ms)')
title(['Ctx n=' num2str(size(PSTH_out_ctx_INT,2))])

subplot(4,5,20)
plot(PSTH_out_ctx_INT,'color',[0.7 0.7 0.7 0.4])
ylimit = get(gca,'ylim');
hold on
plot([100 100],[ylimit(1) ylimit(2)],'k--')
hold off
box off
set(gca,'xtick',[0 100 200],'xticklabel',[-100 0 100])
xlabel('Time lag (ms)')
title('Ctx')
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if exist('PSTH')==7
    cd('PSTH')
else
    [status, msg, msgID] = mkdir('PSTH');
    cd('PSTH')
end

set(gcf,'renderer','painters');
saveas(gcf,[basename '-peri-' eventType '-spiking-histogram.pdf'])
saveas(gcf,[basename '-peri-' eventType '-spiking-histogram.png'])
close all
cd ..

end