function figure_neurobar(fld, limit)

fprintf('\n======================================================\n');
fprintf('-- Creating Figure neurobar --\n');
fprintf('======================================================\n');

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

% NOTE: there is an embedded function in the bottom that has some more
% settings.
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

%% plots
% load traces
monkeys = {'M1','M2'};
alltraces = [];
alltracesLUT = [];
figure; set(gcf,'Position',[ 100 100 1000 1200]);

for mi = 1:2
    monkey = monkeys{mi};
    savedir = fullfile(fld.basedir,'results','figure_neurobar');
    [~,~]=mkdir(savedir);
    
    load(fullfile(savedir, [monkey '_averages_snr' num2str(snrthres) ...
        '_mindays' num2str(mindays) limitRT C2T C2SD '.mat']),...
        'traces','tracesLUT','env_t');

    subj = monkey;
    bar_data = make_bar_plot(subj,traces,tracesLUT,env_t,3,mi);

    %% do some statistics
    % plot 1
    D = bar_data{1}; % from the first plot, correct trials
    [~,p,~,st] = ttest(D(:,1),D(:,2));
    fprintf('\n');
    disp('--- correct trials ---');
    disp([monkey ', target vs non-target, p=' num2str(p) ...
        ', t=' num2str(st.tstat) ', df=' num2str(st.df)]);
    [~,p,~,st] = ttest(D(:,3),D(:,2));
    disp([monkey ', distractor vs non-target, p=' num2str(p) ...
        ', t=' num2str(st.tstat) ', df=' num2str(st.df)]);
    disp(' ');
    
    disp('--- monkey chose non-target ---');
    D = bar_data{2};
    [~,p,~,st] = ttest(D(:,1),D(:,3));
    disp([monkey ', chosen non-target vs target, p=' num2str(p) ...
        ', t=' num2str(st.tstat) ', df=' num2str(st.df)]);
    [~,p,~,st] = ttest(D(:,4),D(:,3));
    disp([monkey ', chosen non-target vs distractor, p=' num2str(p) ...
        ', t=' num2str(st.tstat) ', df=' num2str(st.df)]);
    [~,p,~,st] = ttest(D(:,2),D(:,3));
    disp([monkey ', chosen non-target vs non-chosen non-target, p=' ...
        num2str(p) ', t=' num2str(st.tstat) ', df=' num2str(st.df)]);
    disp(' ');
    [~,p,~,st] = ttest(D(:,1),D(:,4));
    disp([monkey ', target vs distractor, p=' num2str(p) ', t=' ...
        num2str(st.tstat) ', df=' num2str(st.df)]);
    [~,p,~,st] = ttest(D(:,1),D(:,2));
    disp([monkey ', target vs non-chosen non-target, p=' num2str(p) ...
        ', t=' num2str(st.tstat) ', df=' num2str(st.df)]);
    [~,p,~,st] = ttest(D(:,4),D(:,2));
    disp([monkey ', distractor vs non-chosen non-target, p=' num2str(p) ...
        ', t=' num2str(st.tstat) ', df=' num2str(st.df)]);
    disp(' ');
    
    disp('--- monkey chose distractor ---');
    D = bar_data{3};
    [~,p,~,st] = ttest(D(:,1),D(:,2));
    disp([monkey ', target vs non-target, p=' num2str(p) ', t=' ...
        num2str(st.tstat) ', df=' num2str(st.df)]);
    [~,p,~,st] = ttest(D(:,3),D(:,2));
    disp([monkey ', distractor vs non-target, p=' num2str(p) ', t=' ...
        num2str(st.tstat) ', df=' num2str(st.df)]);
    [~,p,~,st] = ttest(D(:,1),D(:,3));
    disp([monkey ', target vs distractor, p=' num2str(p) ', t=' ...
        num2str(st.tstat) ', df=' num2str(st.df)]);
    
    %% concatenate traces for combined plot
    alltraces = [alltraces; traces]; %#ok<*AGROW> 
    alltracesLUT = [alltracesLUT; tracesLUT];
end
subj = 'BOTH';
make_bar_plot(subj,alltraces,alltracesLUT,env_t,3,3);
st=suptitle(['SELECT: ' limitRT C2T C2SD ]);
set(st,'interpreter','none');
% save figure
savefig(gcf,fullfile(savedir, ['figure_neurobar' limitRT C2T C2SD]));
print(gcf,fullfile(savedir, ['figure_neurobar' limitRT C2T C2SD]),'-dpng');

end

function bar_data = make_bar_plot(subj,traces,tracesLUT,env_t,nrows,rownr)

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

smoothfact = 10;
plot_type = 2; % 1 for default (bar), 2 for scatter (showing all channels)

%% make include matrix
tl = tracesLUT;
choices = unique(tracesLUT.choice_id,'stable');
rfstims = unique(tracesLUT.stimulus_in_rf,'stable');

include = cell(0);
% for first figure, vary the stimulus in the RF
for i = 1:3
    include{1}(:,i) = strcmp(tl.choice_id,'Target') & ...
        strcmp(tl.stimulus_in_rf,rfstims{i}) & strcmp(tl.choice_in_rf,'Both');
end

% for second figure, we split the case when a nontarget is chosen
include{2}(:,1) = strcmp(tl.choice_id,'NDist') & ...
    strcmp(tl.stimulus_in_rf,'Target') & strcmp(tl.choice_in_rf,'Both');
include{2}(:,2) = strcmp(tl.choice_id,'NDist') & ...
    strcmp(tl.stimulus_in_rf,'NDist') & strcmp(tl.choice_in_rf,'Out');
include{2}(:,3) = strcmp(tl.choice_id,'NDist') & ...
    strcmp(tl.stimulus_in_rf,'NDist') & strcmp(tl.choice_in_rf,'In');
include{2}(:,4) = strcmp(tl.choice_id,'NDist') & ...
    strcmp(tl.stimulus_in_rf,'SDist') & strcmp(tl.choice_in_rf,'Both');

for i = 1:3
    include{3}(:,i) = strcmp(tl.choice_id,'SDist') & ...
        strcmp(tl.stimulus_in_rf,rfstims{i}) & strcmp(tl.choice_in_rf,'Both');
end

%% plot figures in a loop
use_colors{1} = cols;
use_colors{2} = [green; 0.5 0.5 0.5; 0.1 0.1 0.1; red];
use_colors{3} = cols;
tick_labels{1} = {'T','NT','D'};
tick_labels{2} = {'T','NT-','NT+','D'};
tick_labels{3} = {'T','NT','D'};

%figure; set(gcf,'Position',[ 92 621 1096 371]);
bar_data = cell(0);
for plt = 1:3
    % subplot(1,3,plt);
    subplot(nrows,3,(3*(rownr-1))+plt);
    
    % make window_means
    window_means = [];
    ntrls = [];
    for l = 1:size(include{plt},2) % loop l for level
        incl = include{plt}(:,l);
        ntrls(l) = sum(tl.numtrls(incl));

        % get average within window
        win_idx = find(env_t>0.15 & env_t<0.20);
        get_traces = traces(incl,:);
        s_traces = [];
        for t = 1:size(get_traces,1)
            s_traces(t,:) = smooth(get_traces(t,:),smoothfact);
        end
        window_means(:,l) = mean(s_traces(:,win_idx),2);
    end

    % store window_means in bar_data
    bar_data{plt} = window_means;

    % plot
    mn = mean(window_means);
    sem = std(window_means)./sqrt(size(window_means,1));

    % make a bar plot (DEFAULT)
    b = bar(mean(window_means)); hold all;
    b.FaceColor = 'flat';
    b.CData = use_colors{plt};
    errorbar(mn,sem,'k+');
    ylim([0.4 0.7]);

    if plot_type==2
        % to get info per channel, do a scatter or line plot
        b.FaceAlpha = 0.5;
        x = ones(size(window_means,1),1);
        x = x + (randi(10,length(x),1)-5)/30;
        cx = [];
        for l=1:size(window_means,2)
            cx(:,l) = x + 1*(l-1);
        end

        % find lines that are increasing, only in third plot
        if plt==3
            q = window_means;
            % find increasing lines, i.e. lines that are in line with
            % what we expect
            makebold = q(:,3)>q(:,1);
            plot(cx(makebold,:)',q(makebold,:)','Color',[1 .5 .5]); hold all;
            plot(cx(~makebold,:)',q(~makebold,:)','Color',[.8 .8 .8]);
        else
            % plot lines
            plot(cx',window_means','Color',[.8 .8 .8]); hold all;
        end

        for l=1:size(window_means,2)
            scatter(cx(:,l),window_means(:,l),20,use_colors{plt}(l,:),'f');
        end
        ylim([0 1]);
    end
    set(gca,'XTickLabel',tick_labels{plt});
    title([subj ': ' choices{plt}]);
    if plot_type==2
        txtheight = ones(size(include{plt},2),1)*0.1;
    else
        txtheight = mns+0.01;
    end
    t = text([1:size(include{plt},2)],txtheight,strsplit(num2str(ntrls),' '),...
        'FontSize',8,'HorizontalAlignment','center');
    if plt==2
        xlabel('Stimulus in RF');
    end
    if plt==1
        ylabel('Response magnitude');
    end
end
end