
global DATA
import CMBHOME.*
if ~isstruct(DATA)
    [filename,path] = uigetfile('*.xml','Please select a parameter file for this session');
    SetCurrentSession([path filename])
    [xml, rxml] = LoadXml(filename);
    
end
path = DATA.session.path;

sl = regexp(path,'/');
pardir = path(1:sl(end));
%get analogin

fils = getAllDatFiles(pardir);

fils = fils(cellfun(@(a) any(regexp(a,'auxiliary')),fils));
cmd = [];
if length(fils)>1
    for i = 1:length(fils)
    cmd = [cmd fils{i} ' '];
    end
    cmd = ['cat ' cmd ' > ' filename(1:end-4) '_auxiliary.dat'];
    system(cmd)
    
    auxfil  = [filename(1:end-4) '_auxiliary.dat'];
elseif length(fils)>0
    auxfil = fils{1};
end


%concat all dat files, downsample and save

%%
%allocate memory
if length(fils)>0
d = LoadBinary(auxfil,'nChannels',3,'channels',[1:3]);
  
d = diff([0 0;double(d(1:16:end,[1 3]))]);
outfile = [DATA.session.basename '_CMB.mat'];
vel = sqrt(d(:,1).^2 + d(:,2).^2);
ts = 1/1250:1/1250:length(d)/1250;

else
    
    
end


b_ts = 1/39.125: 1/39.125:max(ts);
x = nan(length(b_ts),1);
y = nan(length(b_ts),1);
hd = nan(length(b_ts),1);
vel = avghist(ts,vel,b_ts);
event = [];
root = CMBHOME.Session('name', outfile, 'b_ts', b_ts,...
    'b_x',x, 'b_y', y, 'b_headdir', hd,'b_vel',vel,...
    'event', {}, 'raw_pos', 1, 'raw_headdir', 1, 'date_created', now, ...
    'epoch', ts([1 end]), 'fs_video', 39.0625, 'path_raw_data', path);

%%



units = GetUnits;
spikes = GetSpikeTimes(units,'output','numbered');


for ii = 1:size(units,1)
    
    spk_ts = spikes(spikes(:,2)==ii,1);
    spk_i = time2ind(root.b_ts,spk_ts);
    spk_ts = spk_ts(spk_i>0);
    spk_i = spk_i(spk_i>0);
    
    root.spike(units(ii,1),units(ii,2)) = CMBHOME.Spike('i',spk_i,'ts',spk_ts,'tet_name',num2str(units(ii,:)));
end

%%

root.path_lfp = repmat({[path '/' filename]},128,1);




%get lfp
%%

fils = getAllExtFiles(pardir,'eeg');


%concat all dat files, downsample and save



    lfp = LoadBinary(fils{1},'nChannels',32,'channels',23);
   
 
lfp = double(lfp);

%%



root.path_lfp = repmat({[path '/' filename]},5,1);

%%
    tet = unique(root.cells(:,1));
fs = 1250;
%fs = DATA.rates.lfp;
channelMap = xml.SpkGrps(tet);




for i = 1:4 
    
    if exist(([root.path_lfp{1}(1:end-3) 'eeg']),'file')
        suf = 'eeg';
    elseif exist(([root.path_lfp{1}(1:end-3) 'lfp']),'file')
        suf = 'lfp';
    end
    %select channel with highest gamma
    [alldat,ch] = PYR_LFP([root.path_lfp{1}(1:end-3) suf],channelMap(i).Channels);
    
    
    signal = alldat;
    ts = 0:(1/fs):((length(signal)-1)/fs);

    channel_name=['tet ' num2str(i) ' ch ' num2str(ch)];
    
    root.b_lfp(i) = CMBHOME.LFP(signal, ts, fs, channel_name);
    root.b_lfp(i) = root.b_lfp(i).AppendTheta;
    root.b_lfp(i) = root.b_lfp(i).AppendThetaAmplitude;
    root.b_lfp(i) = root.b_lfp(i).AppendThetaPhase;
    
    if ~isfield( root.lfp_manager,'theta_freq')
        root.lfp_manager.theta_freq=0;
    end
    root.b_lfp(i) = root.b_lfp(i).AppendThetaFreq;
    
    
    root.b_lfp(i).channel_name = channel_name;
end



for ii = 1:size(root.cells,1)
    
    spk_lfp = Utils.SpkInds(ts,root.spike(root.cells(ii,1),root.cells(ii,2)).ts);
    root.spike(root.cells(ii,1),root.cells(ii,2)).i_lfp = spk_lfp;
end

%%

%events
events = [];idx = 1;
for i = 1:4
temp = getAllExtFiles(pwd,['on' num2str(i)]);
evs1 = LoadEvents(temp{1});

temp = getAllExtFiles(pwd,['of' num2str(i)]);
evs2 = LoadEvents(temp{1});


for j = 1:length(evs1.time)
    events{idx,1} = evs1.description{j};
    events{idx,2} = evs1.time(j);
    idx = idx+1;
    
    events{idx,1} = evs2.description{j};
    events{idx,2} = evs2.time(j);
    idx = idx+1;
    
    
end
end

[~,b] = sort(cell2mat(events(:,2)));
events = events(b,:);
root.event = events;
%%
save([path '/' outfile],'root','-v7.3')

%%
ts1 = root.Label2Time('anin_on 1');
ts2 = root.Label2Time('anin_on 2');
ts3 = root.Label2Time('anin_on 3');
ts4 = root.Label2Time('anin_on 4');

binnedPop=populationMatrix(root,.25,.25,50,ts1);
kp = (max(mean(binnedPop(:,15:30,:),3),[],2)./mean(mean(binnedPop(:,1:15,:),3),2))>2.5;

stim = [ts4(1) ts4(end)];
figure
ax = tight_subplot(3,1);
k = gaussian2Dfilter([200 1],2);
for i =  1:3
    
    switch i
        case 1
            ts = ts1(ts1<stim(1));
            
        case 2
            ts = ts1(ts1>stim(1)& ts1<stim(2));
        case 3
            ts = ts1(ts1>stim(2));
    end

binnedPop=populationMatrix(root,.25,.25,500,ts,k',kp);
axes(ax(i));
imagesc(-.25:.01:.25,[],zscore(nanmean(binnedPop,3),[],2),[-5 5])
hold on

for j = 1:4
    dt = [-.05:.001:.05];
    switch j
        case 1
            off = 0;
           
        case 2
            off = .025;
        case 3
            off = .05;
        case 4
                off = .075;

    end
%    dt = dt +off;
     rang = sum(root.cells(kp,1)==j);
    offset = find(root.cells(kp,1)==j,1,'last')+.5;
    if i~=2
        if ismember(j,1:3)
            plot(dt,-rang*cos([0:.001:.1]* 2*pi*5 - pi/2) + offset,'w','linewidth',2)
        end
    else
         plot(dt,-rang*cos([0:.001:.1]* 2*pi*5 - pi/2) + offset,'w','linewidth',2)
    end
end
colormap('jet')
xlim([-.25 .25])
if i~=3
   set(gca,'xticklabel',[]) 
end
set(gca,'yticklabel',[])
end
%%


