function [all,channelID,count,ht,day,deltas,session] = BatchCheckGammaDuringNonsilentDeltas(session,varargin)


%% 
if strcmp(varargin{8},'Sampling'),day=1; elseif strcmp(varargin{8},'Test'),day=2; end

[~,sessionName] = fileparts(session);

global BASEPATH
deltas = dlmread(fullfile(BASEPATH, 'data', [sessionName '-deltas-lfp-deltagroups']));;

groups = round(deltas(1,4:end));
deltas(1,:) = [];
occ = deltas(:,4:end);
deltas = deltas(:,1:3);

%% Load session
SetCurrentSession(session,'verbose','off');
units = GetUnits;

[PFCgroups,~,~,~,~,~,HPCgroups] = varargin{:};

pfcUnits = units(ismember(units(:,1),PFCgroups),:);
hpcUnits = units(ismember(units(:,1),HPCgroups),:);

hpc = GetSpikeTimes(hpcUnits,'output','numbered');
pfc = GetSpikeTimes(pfcUnits,'output','numbered');


%% count spikes within local spikes

% for each electrode group, count delta spikes only within locally detected deltas
spikes = pfc;
group = round(pfcUnits(spikes(:,2)));
groups = round(groups);
intervals = bsxfun(@plus,deltas,[-1 1]*0.015);

[~,group] = ismember(group,groups); % which of the columns of 'occ' does the electrode group correspond to
noDeltas = group==0;
spikes(noDeltas,:) = []; group(noDeltas) = [];

count0 = zeros(size(intervals,1),1);
for g=1:size(occ,2),
    ok = occ(:,g)>0.5;
    count(ok) = count0(ok)+CountInIntervals(spikes(group==g),intervals(ok,:));
end

channels = dlmread(fullfile(BASEPATH, 'data', [sessionName '-good-PFC-channels']));;nChannels = length(channels);

nDeltas = size(deltas,1);
all = nan(nChannels*nDeltas,100);
channelID = repelem([1:length(channels)]',nDeltas);
for i=1:length(channels),
    lfp = GetLFP(channels(i));
    gamma = FilterLFP(lfp,'passband',[40 100]);
    gamma(:,2) = zscore(gamma(:,2));
    [~,gamma] = Phase(gamma);
    [h,ht] = PlotSync(gamma,deltas,'durations',[-1 1]*2);
    all(channelID==i,:) = h;
end
