function check_snr_over_sessions(fld,monkey,plot_responses)
%
% this function explores some information about the channels and their SNR.
% it first concatenates MP, which is created in
% analyse_individual_channels.m. it chucks out channels that never see a
% stimulus, and it plots the distribution of SNR over sessions. it can plot
% the average responses for individual sessions and individual channels,
% but this is pretty time consuming, because data needs to be loaded. 

%% find which data to include 

fprintf('\n======================================================\n');
fprintf(['-- Running check_snr_over_sessions for ' monkey ' --\n']);
fprintf('======================================================\n');

datainfo = session_info();

if strcmp(monkey,'M1')
    info = datainfo(1);
else 
    info = datainfo(2);
end
sessions = info.dates;

datadir = fld.procdatadir;
MP = get_mp(fld.basedir,monkey,sessions);

%% inspect SNR

chan = findgroups(MP.chan_id);
snr = splitapply(@nanmean,MP.snr,chan);
snr_std = splitapply(@nanstd,MP.snr,chan);
numses = splitapply(@length,MP.snr,chan);

disp(['out of ' num2str(size(MP,1)) ' rows, ' num2str(sum(MP.snr>2.5)) ...
    ' have SNR bigger than 2.5, ' num2str(sum(MP.snr>1)) ' bigger than 1']);

figure; set(gcf,'Position',[129         547        1097         500]);
boxplot(MP.snr,MP.chan_id);hold all
p2 = plot([1 length(chan)],[2.5 2.5],'r');
p1 = plot([1 length(chan)],[1 1],'g');
set(gca,'YScale','log');
title(['SNR per channel for ' monkey ', ' ...
    num2str(length(unique(MP.session_id))) ' sessions']);

chans = unique(MP.chan_id);

%% let's try to get the mean response per channel per day
if plot_responses
    cnt = 1;
    mns = []; mns_lut = [];
    for ch = chans'
        disp(['loading channel ' num2str(ch)]);
        for sid = 1:length(sessions)
            % session info
            session = sessions{sid};
            q = strsplit(session,'_');
            cdate = q{2};
            block = str2num(q{3}(end-4));

            % load the data 
            session_id = [monkey '_' cdate '_' num2str(block)];
            sessiondir = fullfile(datadir, monkey, session_id);
            channeldir = fullfile(sessiondir, 'individual_channels');
            chanfile = fullfile(channeldir, ...
                ['Env_preprocessed_ch' num2str(ch) '_' session_id '.mat']);
            if ~exist(chanfile)
                continue
            end
            load(chanfile);
            load(fullfile(sessiondir, [session_id '_LUT.mat']));

            % calculate the average response to when the target is in the RF
            inclTrls = lut.TarPos==1 & strcmp(lut.Status,'Target');

            % store in a matrix
            mns(cnt,:) = nanmean(data(:,inclTrls),2)';
            mns_lut(cnt,:) = [ch, sid];
            cnt = cnt + 1;
        end
    end

    save([monkey '_mns_per_channel_per_day.mat'],'mns','mns_lut','MP');

    %% plot average response per session
    cols = ceil(sqrt(length(chans)));
    rows = ceil(length(chans)/cols);
    for sid = 1:length(sessions)
        session = sessions{sid};
        q = strsplit(session,'_');
        cdate = q{2};
        block = str2num(q{3}(end-4));
        figure; set(gcf,'Position',[92 475 1055 620])
        cnt = 1;
        for ch = chans'
            subplot(rows,cols,cnt);

            q = mns_lut(:,1)==ch & mns_lut(:,2)==sid;
            mn = mns(q,:);

            % find snr in this session
            q = MP.session_id==sid & MP.chan_id==ch;
            snr = MP.snr(q);

            % plot and set title
            if snr>2.5
                plot(mn,'r');
            elseif snr>1
                plot(mn,'g');
            else
                plot(mn,'k');
            end
            title(['snr:' num2str(snr)]);

            cnt = cnt + 1;
        end
        mtit([monkey ' ' cdate '\_B' num2str(block)],'xoff',0,'yoff',.025);
    end

    %% plot average across sessions 
    figure; set(gcf,'Position',[92 475 1055 620])
    cnt = 1;
    for ch = chans'

        % get means from sessions higher than snr 2.5
        inclSessions = MP.session_id(MP.chan_id==ch & MP.snr>2.5);
        snrs = MP.snr(MP.chan_id==ch & MP.snr>1);
        q = mns_lut(:,1)==ch & ismember(mns_lut(:,2),inclSessions);
        mn = nanmean(mns(q,:));

        subplot(rows,cols,cnt);
        plot(mn);
        title(['snr:' num2str(round(mean(snrs)*100)/100) ' #' ...
            num2str(length(inclSessions))]);

        ylim([0 1]);

        cnt = cnt + 1;
    end
end

%% make a bar plot to indicate how many sessions would be included per channel
% for different levels of SNR. 
threshes = [1 2 2.5];
figure; set(gcf,'Position',[83 500 1180 452]);
for ti = 1:length(threshes)
    thres = threshes(ti);
    incl = MP.snr>thres;
    q = MP(incl,:);
    [chid, inclch] = findgroups(q.chan_id);
    session_per_chan = splitapply(@length,q.session_id,chid);
    subplot(1,length(threshes),ti);
    b = bar(session_per_chan);
    set(gca,'XTick',[1:length(inclch)],'XTickLabel',inclch);
    xtickangle(90);
    title(['SNR>' num2str(thres) ', #' ...
        num2str(sum(session_per_chan>5)) ' chans>5 days'],'fontsize',9);
    ylabel('#sessions that meet SNR criterion');
    xlabel('channel id');
end
mtit([monkey],'xoff',0,'yoff',0.035);