% % Jon Newman
% % Look at stats of Unit firing during OU stimulation
% 
% fid = '2012-02-21_ollie-363_OU-stim_tau0.05sec-u-2.5-s-0.625_num-real-2_00-a';
% s = loadspk([fid '.spk']);
% SS = SqueakySpk(fid,25000,1,s);
% % load([fid '.SS'],'-mat');
% d = loaddig([fid '.dig']);
% [a at] = loadraw([fid '.aux']);
% 
% % Clean and sort
% SS.RemoveUnit(0);
% % SS.RemoveChannel([1 3 8 11 33 37 41 43 46 50 54 57 61 63 64]);
% % SS.HardThreshold(1000,25);
% % SS.WaveClus;
% 
% % Decode digital information
% b = dec2bin(d.state);
% t_r = d.time(d.state ~= 0);
% t_r_1 = d.time(d.state == 1);
% t_r_2 = d.time(d.state == 2);
% 
% % ASDR
% SS.ASDR(0.1);
% SS.Save();
%% Big trial raster

% Data storage
data = struct();
trialnames = {};
for t = 1:length(t_r)
    trialnames{t} = ['trial' num2str(t)];
end

figure()
tot = 0;
maxu = 0;
col = hsv(length(t_r));
height_light = 1000;
hold on
plot([0 10],[0 0],'k-')
plotted_1 = 0;
plotted_2 = 0;
for t = 1:length(t_r)
    
    % Parse out data sections within a given OU repeat
    s_t = SS.time(SS.time > tot & SS.time <= t_r(t));
    s_u = SS.unit(SS.time > tot & SS.time <= t_r(t));
    s_w = SS.waveform(:,SS.time > tot & SS.time <= t_r(t));
    s_c = SS.channel(SS.time > tot & SS.time <= t_r(t));
    a_t = at(at > tot & at <= t_r(t));
    a_i = a(1, at > tot & at <= t_r(t));
    a_l = a(3, at > tot & at <= t_r(t));
    
    
    % Update data struct
    data.(trialnames{t}).spk.time = s_t;
    data.(trialnames{t}).spk.unit = s_u;
    data.(trialnames{t}).spk.channel = s_c;
    data.(trialnames{t}).spk.waveform = s_w';
    
    data.(trialnames{t}).stm.time = a_t;
    data.(trialnames{t}).stm.currentAmps = a_i;
    data.(trialnames{t}).stm.lightArb = a_l;
    
    % Increment where we are in the plot
    maxu = maxu + max(s_u);

    % Plot the light signal
    if (plotted_1 == 0 && ismember(t_r(t),t_r_1) && t_r(t) > 310)
        plot(a_t - tot ,height_light*((a_l-min(a_l))/max(a_l)) ,'r')
        axis tight
        xlim([0 10])
        ylabel('\propto light power','FontSize',16)
        plotted_1 = 1;
    elseif (plotted_2 == 0 && ismember(t_r(t),t_r_2)&& t_r(t) > 310)
        plot(a_t - tot ,height_light*((a_l-min(a_l))/max(a_l)) ,'b')
        axis tight
        xlim([0 10])
        ylabel('\propto light power','FontSize',16)
        plotted_2 = 1;
    end
    
    % Raster Plot
    if ismember(t_r(t),t_r_1)
        col = [1 0 0];
    elseif ismember(t_r(t),t_r_2)
        col = [0 0 1];
    end
    plot(s_t - tot ,s_u - maxu,'.','MarkerSize',2,'Color',col)
    plot([0 10],[-maxu  -maxu],'k-')
    axis tight
    xlim([0 10])
    
    % Increment temporal position
    tot = t_r(t);
    
    %Plot raster
%     drawnow;
    
end

set(gca,'YTick',[]);
ylabel('trial #/ unit # & \propto light power (black)','FontSize',16)
xlabel('time (sec)','FontSize',16)

%% Spont Data

% Pre spont
s_t = SS.time(SS.time > 0 & SS.time <= t_r(1));
s_u = SS.unit(SS.time > 0 & SS.time <= t_r(1));
s_w = SS.waveform(:,SS.time > 0 & SS.time <= t_r(1));
s_c = SS.channel(SS.time > 0 & SS.time <= t_r(1));
data.pre_spont.spk.time = s_t;
data.pre_spont.spk.unit = s_u;
data.pre_spont.spk.channel = s_c;
data.pre_spont.spk.waveform = s_w';

%Post spont
s_t = SS.time(SS.time > t_r(end));
s_u = SS.unit(SS.time > t_r(end));
s_w = SS.waveform(:,SS.time > t_r(end));
s_c = SS.channel(SS.time > t_r(end));
data.post_spont.spk.time = s_t;
data.post_spont.spk.unit = s_u;
data.post_spont.spk.channel = s_c;
data.post_spont.spk.waveform = s_w';
%% save data and plot
% export_fig 2012-02-21_ollie-363_OU-stim_tau0.05sec-u-2.5-s-0.625_num-real-2_00-a -pdf -transparent
% save([fid '.mat'],'data');
% saveas(gca,[fid '.fig']);
