function figure_neurobar_prepare(fld,monkey,limit)

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

% settings 
snrthres = 2.5;
daythres = 3;
do_cleantraces = true; % if true we throw out artefact traces

%% find which data to include 
datainfo = session_info();

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

MP = get_mp(fld,monkey,sessions);

%% set SNR threshold and find channels and sessions
% limit MP to snr threshold
incl = MP.snr>snrthres;
MP_ = MP; %#ok<*NASGU> % copy for backup

% only include channels that have at least 5 sessions
MP = MP(incl,:);
[chid, ch_list] = findgroups(MP.chan_id);
session_per_chan = splitapply(@length,MP.session_id,chid);

% make a list of channels that we will include
chans = ch_list(session_per_chan>daythres);

%% load and concatenate the data we need
datadir = fld.procdatadir;

% first loop all channels 
numtrls = [];
traces = [];
tracesLUT = [];
c = 1;
grandtraces = [];
grandtracesLUT = [];

chnum=1;
for ch = chans'
    disp(['analysing channel ' num2str(ch) ', ' num2str(chnum) ...
        '/' num2str(size(chans,1))]);
    
    chidx = find(chans==ch);
    chnum=chnum+1;
    
    incl_sessions = MP.chan_id==ch & MP.snr>snrthres;
    cur_mp = MP(incl_sessions,:);
    numtrls(chidx) = sum(cellfun(@sum,cur_mp.goodtrials)); %#ok<*FNDSB,*AGROW> 
    
    % look at each session for each channel individually
    env = []; clut = [];
    for sid = 1:size(cur_mp,1)

        % session info
        cdate = cur_mp.date{sid};
        block = cur_mp.block(sid);

        % 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,'file')
            continue
        end

        % load the data for this session for this channel
        load(chanfile); %#ok<*LOAD> % loads 'data'
        load(fullfile(sessiondir, [session_id '_LUT.mat'])); % loads 'lut'
        load(fullfile(sessiondir, [session_id '.mat']),'env_t');
        % this is overwritten every time but it's always the same so who cares

        if do_cleantraces % throw out artefact traces
            prewin = env_t > -0.1 & env_t < 0;
            postwin = env_t > 0 & env_t < 0.1;
            trace = mean(data');

            sd_pre = std(trace(prewin));
            m_post = mean(trace(postwin));

            if m_post < 3*sd_pre
                data = nan(size(data));
            end
        end

        env = [env; data'];
        clut = [clut; lut];
    end
    
    lut = clut;
    
    % we now have three independent variables, with 3 or 2 levels each: 
    % - choice_identity: 
    %   - target (correct)
    %   - nontarget
    %   - distractor
    % - stimulus_in_RF: 
    %   - target 
    %   - nontarget
    %   - distractor
    % - choice_in_rf: 
    %   - yes 
    %   - no
    
    % choice_identity
    ls = lut.Status;
    choice_identity = [strcmp(ls,'Target'), strcmp(ls,'Nontarget'), ...
        strcmp(ls,'Distractor')];
    choice_label = {'Target','NDist','SDist'};
    
    % get distances between chosen and target or salient distractor
    chdistT = [lut.TarPos lut.TarChoice abs(lut.TarChoice-lut.TarPos)];
    chdistSD = [lut.DistPos lut.TarChoice abs(lut.TarChoice-lut.DistPos)];

    chdistT(chdistT(:,3)<0,3)=chdistT(chdistT(:,3)<0,3)+6;
    chdistT(chdistT(:,3)>3,3) = 6-chdistT(chdistT(:,3)>3,3);
    chdistSD(chdistSD(:,3)<0,3)=chdistSD(chdistSD(:,3)<0,3)+6;
    chdistSD(chdistSD(:,3)>3,3) = 6-chdistSD(chdistSD(:,3)>3,3);

    for tr = 1:size(chdistT,1)
        if chdistT(tr,3)>3
            chdistT(tr,3) = mod(chdistT(tr,3),3);
        end
        if chdistSD(tr,3)>3
            chdistSD(tr,3) = mod(chdistSD(tr,3),3);
        end
    end

    % get RT
    RT = lut.reactTime2;
    % 30% fast/slow
    fastRT = RT<prctile(RT,limit.ReactimeTimePerc);
    slowRT = RT>prctile(RT,100-limit.ReactimeTimePerc);

    %% plot distance from target/dist on error trials
    if limit.PlotBehavior && chidx==1
        figure; set(gcf,'Position',[100 100 1000 1000])
        subplot(3,3,1);
        histogram(chdistT(choice_identity(:,2),3),0.5:1:3.5,...
            'Normalization','probability');
        title('Choice ND: dist. T ALL-RT');
        xlabel('#Stim from T');ylabel('Proportion');
        set(gca, 'ylim',[0 1]);

        subplot(3,3,2);
        histogram(chdistT(choice_identity(:,3),3),0.5:1:3.5,...
            'Normalization','probability');
        title('Choice SD: dist. T ALL-RT');
        xlabel('#Stim from T');ylabel('Proportion');
        set(gca, 'ylim',[0 1]);

        subplot(3,3,3);
        histogram(chdistSD(choice_identity(:,2),3),0.5:1:3.5,...
            'Normalization','probability');
        title('Choice ND: dist. SD ALL-RT');
        xlabel('#Stim from T');ylabel('Proportion');
        set(gca, 'ylim',[0 1]);

        subplot(3,3,4);
        histogram(chdistT(choice_identity(:,2) & fastRT,3),0.5:1:3.5,...
            'Normalization','probability');
        title('Choice ND: dist. T FAST-RT');
        xlabel('#Stim from T');ylabel('Proportion');
        set(gca, 'ylim',[0 1]);

        subplot(3,3,5);
        histogram(chdistT(choice_identity(:,3) & fastRT,3),0.5:1:3.5,...
            'Normalization','probability');
        title('Choice SD: dist. T FAST-RT');
        xlabel('#Stim from T');ylabel('Proportion');
        set(gca, 'ylim',[0 1]);

        subplot(3,3,6);
        histogram(chdistSD(choice_identity(:,2) & fastRT,3),0.5:1:3.5,...
            'Normalization','probability');
        title('Choice ND: dist. SD FAST-RT');
        xlabel('#Stim from T');ylabel('Proportion');
        set(gca, 'ylim',[0 1]);

        subplot(3,3,7);
        histogram(chdistT(choice_identity(:,2) & slowRT,3),0.5:1:3.5,...
            'Normalization','probability');
        title('Choice ND: dist. T SLOW-RT');
        xlabel('#Stim from T');ylabel('Proportion');
        set(gca, 'ylim',[0 1]);

        subplot(3,3,8);
        histogram(chdistT(choice_identity(:,3) & slowRT,3),0.5:1:3.5,...
            'Normalization','probability');
        title('Choice SD: dist. T SLOW-RT');
        xlabel('#Stim from T');ylabel('Proportion');
        set(gca, 'ylim',[0 1]);

        subplot(3,3,9);
        histogram(chdistSD(choice_identity(:,2) & slowRT,3),0.5:1:3.5,...
            'Normalization','probability');
        title('Choice ND: dist. SD SLOW-RT');
        xlabel('#Stim from T');ylabel('Proportion');
        set(gca, 'ylim',[0 1]);

        suptitle(['Subj: ' monkey ', ERRORS']);
    end

    % stimulus_in_RF
    targetInRF = lut.TarPos==MP.rfpos(ch);
    nontargInRF = lut.TarPos~=MP.rfpos(ch) & lut.DistPos~=MP.rfpos(ch);
    distInRF = lut.DistPos==MP.rfpos(ch);
    stimulus_in_RF = [targetInRF, nontargInRF, distInRF];
    stim_rf_label = {'Target','NDist','SDist'};
    
    % choice_in_RF; to make it easier, let's define three 
    rfchosen = lut.TarChoice==MP.rfpos(ch);
    choice_in_rf = [rfchosen ~rfchosen ones(length(rfchosen),1)];
    choice_rf_label = {'In','Out','Both'};
    
    % additional criteria that are not independent variables
    goodtrials = [cur_mp.sample_wise_goodtrials{:}]';
    
    inclRT = ones(size(goodtrials));
    if limit.ReactionTime~=0
        if limit.ReactionTime==-1
            pp = prctile(lut.reactTime2,limit.ReactimeTimePerc);
            inclRT = lut.reactTime2<pp;
        else
            pp = prctile(lut.reactTime2,100-limit.ReactimeTimePerc);
            inclRT = lut.reactTime2>pp;
        end
    end
    
    inclD_C2T = ones(size(goodtrials));
    if limit.DistChoiceToTarget
        inclD_C2T = (chdistT(:,3) >= limit.DistChoiceToTarget_range(1) & ...
            chdistT(:,3) <= limit.DistChoiceToTarget_range(2));
    end
    inclD_C2SD = ones(size(goodtrials));
    if limit.DistChoiceToSD
        inclD_C2SD = chdistSD(:,3) >= limit.DistChoiceToSD_range(1) & ...
            chdistSD(:,3) <= limit.DistChoiceToSD_range(2);
    end

    for choice_level = 1:3
        for stimulus_level = 1:3
            for choice_rf_level = 1:3
                incl = choice_identity(:,choice_level) & ...
                    stimulus_in_RF(:,stimulus_level) & ...
                    choice_in_rf(:,choice_rf_level) & ...
                    goodtrials & inclRT & inclD_C2T & inclD_C2SD;
                
                mn = mean(env(incl,:),'omitnan');
                traces(c,:) = mn;
                tracesLUT.ch(c,1) = ch;
                tracesLUT.choice_id{c,1} = choice_label{choice_level};
                tracesLUT.stimulus_in_rf{c,1} = stim_rf_label{stimulus_level};
                tracesLUT.choice_in_rf{c,1} = choice_rf_label{choice_rf_level};
                tracesLUT.numtrls(c,1) = sum(incl);
                c = c + 1;
            end
        end
    end
end
tracesLUT = struct2table(tracesLUT);

%% save output to make a combined plot of all channels
savedir = fullfile(fld.basedir,'results','figure_neurobar');
[~,~]=mkdir(savedir);
if limit.ReactionTime==0
    limitRT = '_allRT';
else
    limitRT = ['_limitRT_' num2str(limit.ReactionTime*limit.ReactimeTimePerc)];
end

if limit.DistChoiceToTarget
    C2T = ['_C2T-range_' num2str(limit.DistChoiceToTarget_range(1)) ...
        '-' num2str(limit.DistChoiceToTarget_range(2))];
else
    C2T = '_allC2T';
end

if limit.DistChoiceToSD
    C2SD = ['_C2SD-range_' num2str(limit.DistChoiceToSD_range(1)) ...
        '-' num2str(limit.DistChoiceToSD_range(2))];
else
    C2SD = '_allC2SD';
end

save(fullfile(savedir, [monkey '_averages_snr' num2str(snrthres) ...
    '_mindays' num2str(daythres) limitRT C2T C2SD '.mat']));