

eps = [0 4833.500;4833.500 7740.4];


%v = load('/data/SUBCA1/SM11/xml_chmap');
k = gaussian2Dfilter([100 1],5);


%v = load('/data/rat-40/xml_chmap');
%v = load('/data/rat-29/xml_chmap');
  import CMBHOME.*
            import CMBHOME.Utils.*


global DATA

if ~isstruct(DATA)
    [filename,path] = uigetfile('*.xml','Please select a parameter file for this session');
    

    SetCurrentSession([path filename])
    [xml, rxml] = LoadXml(filename);
    
end

coords = DATA.positions;
t2 = coords(:,1);
coords(coords <=0) = nan;



[~,x2,y2] = meanCoord(coords);

t2 = t2+eps(2,1);
x2 = x2*.51;
y2= y2*.51;
dt = mode(diff(t2));
t1 = [dt:dt:eps(2,1)]';
x1 = nan(length(t1),1);
y1 = nan(length(t1),1);

x = [x1;x2];
t = [t1;t2];
y = [y1;y2];


x(x==0) = nan;
y(y==0) = nan;

x = nanconvn(x,k,'nanout',0,'nonanout',1);
y = nanconvn(y,k,'nanout',0,'nonanout',1);

% Find the head direction based on the LEDs
if(exist('hd','var')~=1 && size(coords,2)>=5)
    hd = headdirectionLED(coords);
    hd(:,2) = hd(:,2).*(180/pi);
    %     [hd, ~, LEDloc] = headdirectionLED(coords);
else
    hd = nan(length(t),2);
end



%vel = cat(1, 0, sqrt(diff(x).^2+diff(y).^2)) / mode(diff(t));
% Unless the LEDs were mounted on the cup, the HD is inaccurate.
% if(strcmp(LEDloc,'headstage')); hd(:,2) = nan(size(hd,1),1); end

% Create an output file name by strapping '.CMB.mat' to the end of the
% existing file name.
path = DATA.session.path;
outfile = [DATA.session.basename '_CMB.mat'];





event = cell(0,2);
% Create CMBHOME object.
%%
root = CMBHOME.Session('name', outfile, 'b_ts', t,...
    'b_x',x, 'b_y', y, 'b_headdir', hd(:,2),...
    'event', event, 'raw_pos', 1, 'raw_headdir', 1, 'date_created', now, ...
    'epoch', t([1 end]), 'fs_video', 39.0625, 'path_raw_data', path);
%%
root = root.FixPos;




[len,pos,xp,yp,del] =linearizePositionTrack(root.b_x,root.b_y);
root.b_linpos = len;
root.spatial_scale = 1;



%%

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);

%%
tet = unique(root.cells(:,1));



channelMap = xml.SpkGrps(tet);


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

for i = 1:length(tet)
    
    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
    
    [alldat,ch] = PYR_LFP([root.path_lfp{1}(1:end-3) suf],channelMap(i).Channels);
    
    
    signal = alldat;
    ts = 0:(1/DATA.rates.lfp):((length(signal)-1)/DATA.rates.lfp);
    fs = DATA.rates.lfp;
    channel_name=['tet ' num2str(tet(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




%%
%root = defineEventsSes(root);

save([path '/' outfile],'root','-v7.3')

% 
% %%
% ts_13_22 = [root.spk_ts([12 3]);root.spk_ts([14 3])];
% ts_11_13 = [root.spk_ts([11 3]);root.spk_ts([13 3])];
% ch_13_22 = [90 83 92 87 93 91 94 95 71 78 70 79 72 76 75 74];
% ch = [80 89 81 88 82 86 84 85 77 68 73 66 69 67 65 64];
% FileName = '/data/SUBCA1/SM11/SM11_mar18/SM11_mar18.fil';
% for i = 1:16
% [data2 OrigIndex]= LoadBinary_bw(FileName, ch(i), [round(ts_11_13/(1/20000))-(.01/(1/20000)) round(ts_11_13/(1/20000))+(.01/(1/20000))]);
% data2 = reshape(data2',401,size(data2,2)/401)';
% dat{i} = data2;
% i
% end
% %%
% ts1 = round(1000*root.spk_ts([11 3]))/1000;
% ts2 = round(1000*root.spk_ts([13 3]))/1000;
% comp = round(ts_11_13*1000)/1000;
% BnotA = ~ismember(comp,ts1) & ismember(comp,ts2);
% AnotB = ismember(comp,ts1) & ~ismember(comp,ts2);
% AandB = ismember(comp,ts1) & ismember(comp,ts2);
% AnotB_dat = cell2mat(cellfun(@(a) nanmedian(a(AnotB,:)),dat,'uni',0)');
% BnotA_dat = cell2mat(cellfun(@(a) nanmedian(a(BnotA,:)),dat,'uni',0)');
% AandB_dat = cell2mat(cellfun(@(a) nanmedian(a(AandB,:)),dat,'uni',0)');
% 
% %%
% close all
% figure
% ax = tight_subplot(2,3);
% axes(ax(1))
% for i=1:8
% plot(-.01:(1/20000):.01,BnotA_dat(i,:)+80*i)
% hold on
% ylim([-50 650])
% axis off
% end
% 
% axes(ax(4))
% for i=9:16
% plot(-.01:(1/20000):.01,BnotA_dat(i,:)+80*(i-8))
% hold on
% ylim([-50 650])
% axis off
% end
% 
% 
% 
% 
% axes(ax(2))
% for i=1:8
% plot(-.01:(1/20000):.01,AandB_dat(i,:)+80*i)
% hold on
% ylim([-50 650])
% axis off
% end
% 
% axes(ax(5))
% for i=9:16
% plot(-.01:(1/20000):.01,AandB_dat(i,:)+80*(i-8))
% hold on
% ylim([-50 650])
% axis off
% end
% 
% 
% axes(ax(3))
% for i=1:8
% plot(-.01:(1/20000):.01,AnotB_dat(i,:)+80*i)
% hold on
% ylim([-50 650])
% axis off
% end
% 
% axes(ax(6))
% for i=9:16
% plot(-.01:(1/20000):.01,AnotB_dat(i,:)+80*(i-8))
% hold on
% ylim([-50 650])
% axis off
% end
% 
% %%
% 
% 
% figure
% spks = [root.spk_ts([11 3]);root.spk_ts([13 3])];
% gr = [ones(length(root.spk_ts([11 3])),1);2*ones(length(root.spk_ts([13 3])),1)];
% [spks,b] = sort(spks);
% gr = gr(b);
% 
% c = CCG(spks,gr,'binsize',1/20000);
% 
% bar(-1:1/20000:1,c(:,1,2)')
% xlim([-1 1])
