function figure_neurotrace_prepare(fld,monkey,do_plots)

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

% settings 
snrthres = 2.5;
daythres = 3;
green = [23, 105, 13]./255;
red = [201, 0, 34]./255;
cols = [green; 0.3 0.3 0.3; red];

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

%% initialize a figure if necessary 
if do_plots
    figure; set(gcf,'Position',[43 476 1197 622]);
    rws = ceil(sqrt(length(chans)));
    cls = ceil(length(chans)/rws);
end

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

% first loop all channels 
numtrls = [];
traces = [];
tracesLUT = [];
grandtraces = [];
grandtracesLUT = [];
cnt=1;
for ch = chans'
    disp(['analysing channel ' num2str(ch) ...
        ', ' num2str(cnt) '/' num2str(size(chans,1))]);
    cnt=cnt+1;
    chidx = find(chans==ch);
    
    incl_sessions = MP.chan_id==ch & MP.snr>snrthres;
    cur_mp = MP(incl_sessions,:);
    numtrls(chidx) = sum(cellfun(@sum,cur_mp.goodtrials));
    
    % 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)
            continue
        end
        
        % load the data for this session for this channel
        load(chanfile); % 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
        
        env = [env; data'];
        clut = [clut; lut];
    end
    
    lut = clut;
    
    % include only correct trials, in which the respective stimulus was in
    % the RF 
    correct = strcmp(lut.Status,'Target');
    goodtrials = [cur_mp.sample_wise_goodtrials{:}]';
    targetInRF = lut.TarPos==MP.rfpos(ch);
    nontargInRF = lut.TarPos~=MP.rfpos(ch) & lut.DistPos~=MP.rfpos(ch);
    distInRF = lut.DistPos==MP.rfpos(ch);

    condincl = {targetInRF, nontargInRF, distInRF};
    
    % store in traces, this is the concatenated average
    for cond = 1:3
        incl = correct & goodtrials & condincl{cond};
        mn = nanmean(env(incl,:));
        
        traces(end+1,:) = mn;
        tracesLUT(end+1,:) = [ch,cond,sum(incl)];
    end
    
    % store in grandtraces, this is the average of the session average    
    do_sessions = unique(lut.session);
    for sid = 1:length(do_sessions)
        this_session = strcmp(lut.session,do_sessions{sid});
        for cond = 1:3
            incl = this_session & correct & goodtrials & condincl{cond};
            mn = nanmean(env(incl,:));

            grandtraces(end+1,:) = mn;
            grandtracesLUT(end+1,:) = [ch,cond,sum(incl),sid];
        end
    end
    
    % do analysis for individual channels
    if do_plots
        subplot(rws,cls,find(chans==ch));
        anova_data = [];
        anova_group = [];
        for cond = 1:3
            incl = correct & goodtrials & condincl{cond};

            win_idx = find(env_t>0.15 & env_t<0.20);
            anova_data = [anova_data; nanmean(env(incl,win_idx),2)]; % window average for anova
            anova_group = [anova_group; repmat(cond,sum(incl),1)];
            
            mn = nanmean(env(incl,:));
            sem = nanstd(env(incl,:))./sqrt(sum(incl));

            t2 = [env_t, fliplr(env_t)];
            inBetween = [mn+sem, fliplr(mn-sem)];
            fill(t2, inBetween, 'g','LineStyle','none','FaceColor',...
                cols(cond,:),'FaceAlpha',0.3); hold all
            plot(env_t,mn,'Color',cols(cond,:));
        end
        [p,tbl,stats] = anova1(anova_data,anova_group,'off');
        xlim([0 .25]);
        title(['ch' num2str(ch) ', p=' num2str(round(p*100)/100) ...
            ', ' num2str(length(do_sessions)) ' days'],'FontSize',8);
    end
end

% save output to make a combined plot of all channels
savedir = fullfile(fld.basedir,'results','figure_neurotrace');
save(fullfile(savedir, [monkey '_averages_snr' num2str(snrthres) ...
    '_mindays' num2str(daythres) '.mat']));

if do_plots
    savefig(fullfile(savedir, [monkey '_channelwise']));
end