load('buster359_FullRecNoWaveformNoClean.mat');



%%
tc = sortrows([time,channel],1);
dt = 0.001;

pre.time = tc(tc(:,1)<=e1(1),1);
pre.channel = tc(tc(:,1)<=e1(1),2);
pre.ISI = [NaN; diff(pre.time)];
pre.bins = 0.5*dt:dt:max(pre.ISI(pre.ISI<60));
pre.hist = hist(pre.ISI,pre.bins);
pre.smooth = smooth(pre.hist.*pre.bins,round(5*dt*1000));

[xmax imax xmin imin] = extrema(pre.smooth);
pre.ISIcut = imin(dt*imin>=0.01);
pre.ISIcut = min(pre.ISIcut);

figure;
bar(pre.bins,pre.hist.*pre.bins);
hold on;
plot(pre.bins,pre.smooth,'--r');
xlim([-0.001 0.1]);
ylim([0 200])
plot(pre.bins(pre.ISIcut),pre.smooth(pre.ISIcut(1)),'*r');

pre.isburst = zeros(size(pre.ISI));
pre.isburst(pre.ISI<pre.bins(pre.ISIcut)) = 1;

figure;
plot(pre.time(pre.isburst==1),pre.channel(pre.isburst==1)+0.1*randn(size(pre.channel(pre.isburst==1))),'.r','MarkerSize',1);
hold on; plot(pre.time(pre.isburst==0),pre.channel(pre.isburst==0)+0.1*randn(size(pre.channel(pre.isburst==0))),'.k','MarkerSize',1);

%%
tc = sortrows([time,channel],1);
dt = 0.001;

CNQX.time = tc(tc(:,1)>e1(1),1);
CNQX.channel = tc(tc(:,1)>e1(1),2);
CNQX.ISI = [NaN; diff(CNQX.time)];
CNQX.bins = 0.5*dt:dt:max(CNQX.ISI(CNQX.ISI<60));
CNQX.hist = hist(CNQX.ISI,CNQX.bins);
CNQX.smooth = smooth(CNQX.hist.*CNQX.bins,round(5*dt*1000));

[xmax imax xmin imin] = extrema(CNQX.smooth);
CNQX.ISIcut = imin(dt*imin>=0.01);
CNQX.ISIcut = min(CNQX.ISIcut);

figure;
bar(CNQX.bins,CNQX.hist.*CNQX.bins);
hold on;
plot(CNQX.bins,CNQX.smooth,'--r');
xlim([-0.001 0.1]);
ylim([0 200])
plot(CNQX.bins(CNQX.ISIcut),CNQX.smooth(CNQX.ISIcut(1)),'*r');

CNQX.isburst = zeros(size(CNQX.ISI));
CNQX.isburst(CNQX.ISI<CNQX.bins(CNQX.ISIcut)) = 1;

figure;
plot(CNQX.time(CNQX.isburst==1),CNQX.channel(CNQX.isburst==1)+0.1*randn(size(CNQX.channel(CNQX.isburst==1))),'.r','MarkerSize',1);
hold on; plot(CNQX.time(CNQX.isburst==0),CNQX.channel(CNQX.isburst==0)+0.1*randn(size(CNQX.channel(CNQX.isburst==0))),'.k','MarkerSize',1);


%%
disp(['FR in bursts before drug: ' num2str(sum(pre.isburst)/e1(1)) 'Hz']);
disp(['FR outside of bursts before drug: ' num2str(sum(abs(pre.isburst-1))/e1(1)) 'Hz']);
disp(['FR in bursts after CNQX: ' num2str(sum(CNQX.isburst)/(max(time)-e1(1))) 'Hz']);
disp(['FR outside of after CNQX: ' num2str(sum(abs(CNQX.isburst-1))/(max(time)-e1(1))) 'Hz']);


