function RTs = figure_behavior(fld)

fprintf('\n======================================================\n');
fprintf('-- Creating Figure for behavior --\n');
fprintf('======================================================\n');

%% Get sessions ---------------------------------------------------------
datainfo = session_info();
datadir = fld.procdatadir;
mkdir(fullfile(fld.basedir,'results'));
    
%% Settings -------------------------------------------------------------
green = [23, 105, 13]./255;
red = [201, 0, 34]./255;
cols = [green; 0.3 0.3 0.3; red];

%% Plots ----------------------------------------------------------------
% load traces
monkeys = {'M1','M2'};
f1 = figure;
f2 = figure;
f3 = figure;
for mi = 1:2
    monkey = monkeys{mi};
    
    info = datainfo(mi);
    sessions = info.dates;
    
    clut = [];
    for sid = 1:length(sessions)
        session = sessions{sid};
        % session info
        q = strsplit(session,'_');
        cdate = q{2};
        block = q{3}(4);
        
        % load the data 
        session_id = [monkey '_' cdate '_' num2str(block)];
        sessiondir = fullfile(datadir, monkey, session_id);
        
        % load the data for this session for this channel
        load(fullfile(sessiondir, [session_id '_LUT.mat'])); % loads 'lut'
        clut = [clut; lut];
    end
    
    lut = clut;
    
    disp(['for monkey ' monkey ' we recorded ' num2str(size(lut,1)) ...
        ' trials in ' num2str(length(sessions)) ' days']);
    
    % calculate concatenated means
    choices = {'Target','Nontarget','Distractor'};
    q = cellfun(@(x) sum(strcmp(lut.Status,x)),choices);
    response_ratio = q./sum(q);
    
    RT = cellfun(@(x) mean(lut.reactTime2(strcmp(lut.Status,x))),choices);
    for ic=1:3
        collectRT{mi,ic} = lut.reactTime2(strcmp(lut.Status,choices{ic}));
    end
    RTstd = cellfun(@(x) std(lut.reactTime2(strcmp(lut.Status,x))),choices);
    RTntrls = cellfun(@(x) sum(strcmp(lut.Status,x)),choices);
    RTsem = RTstd./sqrt(RTntrls);
    
    % calculate response ratio per session to provide SEM estimate
    incl_sessions = unique(lut.session);
    all_ratios = [];
    for sid = 1:length(incl_sessions)
        clut = lut(strcmp(lut.session,incl_sessions{sid}),:);
        q = cellfun(@(x) sum(strcmp(clut.Status,x)),choices);
        all_ratios(sid,:) = q./sum(q);
    end
    
    sem = std(all_ratios);%./sqrt(size(all_ratios,1)); 
    % NB! SEMs are so small, plot SD over sessions instead
    
    % collect choice proportions for stats
    choiceprop{mi} = all_ratios;
    
    % behavioral responses
    figure(f1);
    subplot(2,2,mi);
    mn = response_ratio;
    b = bar(mn); hold all
    b.FaceColor = 'flat';
    b.CData = cols;
    errorbar(mn,sem,'k.');
    ylim([0 1]);
    title(monkey);
    if mi==1
        ylabel('Ratio');
    end
    set(gca,'XTickLabel',{'T','ND','SD'});
    
    % same, corrected for different chance levels
    figure(f3);
    subplot(2,2,mi);
    mn = response_ratio;
    mn(2) = mn(2)/4;
    sem(2) = sem(2)/4;
    b = bar(mn); hold all
    b.FaceColor = 'flat';
    b.CData = cols;
    errorbar(mn,sem,'k.');
    ylim([0 1]);
    title([monkey ', corrected for chance level']);
    if mi==1
        ylabel('Ratio');
    end
    set(gca,'XTickLabel',{'T','ND','SD'});
    
    % reaction time
    figure(f1)
    subplot(2,2,mi+2);
    b2 = bar(RT); hold all
    b2.FaceColor = 'flat';
    b2.CData = cols;
    errorbar(RT,RTsem,'k.');
    ylim([190 290])
    if mi==1
        ylabel('Reaction time (ms)');
    end
    xlabel('Chosen target');
    set(gca,'XTickLabel',{'T','ND','SD'});
    fprintf(['Mean RTs for ' monkey ': T ' num2str(RT(1)) ...
        'ms, ND ' num2str(RT(2)) , 'ms , SD ' num2str(RT(3)) 'ms\n'])
    RTs{mi}=RT;
    
    %% RT distribution
    figure(f2);
    subplot(2,2,mi);
    histogram(lut.reactTime2);
    xlabel('Reaction Time (ms)');
    xlim([0 500]);
    title(monkey);
    
    subplot(2,2,mi+2); hold off
    for c = 3:-1:1
        rts = lut.reactTime2(strcmp(lut.Status,choices{c}));
        rts(rts>500) = [];
        histogram(rts,70,'FaceColor',cols(c,:),'Normalization',...
            'probability'); hold all
    end
    xlim([100 400]);

    %% RT distribution collect
    for c = 3:-1:1
        RTcoll{mi,c} = lut.reactTime2(strcmp(lut.Status,choices{c}));
    end
end
% save all RTs
save(fullfile(fld.basedir, 'results','figure_behavior','RTcoll.mat'),'RTcoll');

% Generate some plots
RT_SlidingWindow;
% RT_Distributions;

% add the pooled average RT
RTs{mi+1} = [mean([collectRT{1,1};collectRT{2,1}]) ...
    mean([collectRT{1,2};collectRT{2,2}]) ...
    mean([collectRT{1,3};collectRT{2,3}]) ];
fprintf(['Mean pooled RTs : T ' num2str(RTs{3}(1)) ...
        'ms, D ' num2str(RTs{3}(2)) , 'ms , SD ' num2str(RTs{3}(3)) 'ms\n'])

% add the pooled choice proportions
choiceprop{3} = [choiceprop{1}; choiceprop{2}];

%% Stats ----------------------------------------------------------------
fprintf ('===============================================\n');
fprintf (' Statistics CHOICE PROPORTIONS\n');
fprintf ('===============================================\n');
for mi=1:2
    % choices
    cp = choiceprop{mi};
    % against chance
    [h,p,~,stats]=ttest(cp(:,1),1/6,'tail','right');
    fprintf([monkeys{mi} ' pTARG > chance: p = ' num2str(p) ...
        ', t = ' num2str(stats.tstat) ', df = ' num2str(stats.df) '\n']);
    [h,p,~,stats]=ttest(cp(:,2),4/6,'tail','left');
    fprintf([monkeys{mi} ' pNDIST < chance: p = ' num2str(p) ...
        ', t = ' num2str(stats.tstat) ', df = ' num2str(stats.df) '\n']);
    [h,p,~,stats]=ttest(cp(:,3),1/6,'tail','left');
    fprintf([monkeys{mi} ' pSDIST < chance: p = ' num2str(p) ...
        ', t = ' num2str(stats.tstat) ', df = ' num2str(stats.df) '\n']);
    % pSD<pND
    [h,p,~,stats]=ttest(cp(:,3),cp(:,2),'tail','left');
    fprintf([monkeys{mi} ' pSD < pND: p = ' num2str(p) ...
        ', t = ' num2str(stats.tstat) ', df = ' num2str(stats.df) '\n']);
end
% choices
cp = choiceprop{3};
% against chance
[h,p,~,stats]=ttest(cp(:,1),1/6,'tail','right');
fprintf(['POOLED pTARG > chance: p = ' num2str(p) ...
    ', t = ' num2str(stats.tstat) ', df = ' num2str(stats.df) '\n']);
[h,p,~,stats]=ttest(cp(:,2),4/6,'tail','left');
fprintf(['POOLED pNDIST < chance: p = ' num2str(p) ...
    ', t = ' num2str(stats.tstat) ', df = ' num2str(stats.df) '\n']);
[h,p,~,stats]=ttest(cp(:,3),1/6,'tail','left');
fprintf(['POOLED pSDIST < chance: p = ' num2str(p) ...
    ', t = ' num2str(stats.tstat) ', df = ' num2str(stats.df) '\n']);
% pSD<pND
[h,p,~,stats]=ttest(cp(:,3),cp(:,2),'tail','left');
fprintf(['POOLED pSD < pND: p = ' num2str(p) ...
    ', t = ' num2str(stats.tstat) ', df = ' num2str(stats.df) '\n']);

fprintf ('===============================================\n');
fprintf (' Statistics REACTION TIME\n');
fprintf ('===============================================\n');
for mi=1:2
    % choices
    fprintf(['Stats for ' monkeys{mi} '\n']);
    rt = [collectRT{mi,1};collectRT{mi,2};collectRT{mi,3}];
    s = categorical([ones(size(collectRT{mi,1}));...
        2*ones(size(collectRT{mi,2})); ...
        3*ones(size(collectRT{mi,3}))]);
    [p,tbl,stats] = kruskalwallis(rt,s);
    fprintf(['Kruskal Wallis -- H = ' num2str(tbl{2,5}(1)) ...
        ', df = ' num2str(tbl{2,3}(1)) ', p = ' num2str(p) '\n']);
    [c,m,h,gnames] = multcompare(stats);
    fprintf('Tukey HSD -----\n')
    fprintf('T=1, ND=2, SD=3 \n')
    for i=1:size(c,1)
        fprintf([gnames{c(i,1)} ' vs ' gnames{c(i,2)} ...
            ', p = ' num2str(c(i,6))  '\n'])
    end
end

%% Save figures ---------------------------------------------------------    
figure(f1)
savefig(fullfile(fld.basedir,'results','figure_behavior',...
    'figure_behavior.fig'));
figure(f2)
savefig(fullfile(fld.basedir, 'results', 'figure_behavior',...
    'figure_rt_distribution.fig'));

%% Save RTs -------------------------------------------------------------
save(fullfile(fld.basedir, 'results','figure_behavior','RTs.mat'),'RTs');