function figure_neurotrace(fld) 

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

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

snrthres = 2.5;
smoothfact = 15; mindays = 3;
falpha = 1; % set to 1 for converting to illustrator

do_cleantraces = true; % if true we throw out artefact traces
% citerion: avg emplitude in first 100ms after stimulus onset should be 
% at least 3 times the standard deviation in the preceding 100 ms
do_latency = true;

%% load RTs
if exist(fullfile(fld.basedir, 'results','figure_behavior','RTs.mat'),'file')
    load(fullfile(fld.basedir, 'results','figure_behavior','RTs.mat'),'RTs');  
else
    error('RTs missing, run figure_behavior.m first');
end

%% plots
% load traces
monkeys = {'M1','M2'};
wm = [];

%% plot
f1 = figure;
set(f1,'Position',[500 500 1600 800])
% individual animals
for mi = 1:2
    monkey = monkeys{mi};
    savedir = fullfile(fld.basedir,'results','figure_neurotrace');
    load(fullfile(savedir, [monkey '_averages_snr' num2str(snrthres) '_mindays' num2str(mindays) '.mat']),...
        'traces','tracesLUT','grandtraces','grandtracesLUT','env_t');

    traces_ALL{mi} = traces; %#ok<*AGROW,*NASGU> 
    tracesLUT_ALL{mi} = tracesLUT;
    grandtraces_ALL{mi} = grandtraces;
    grandtracesLUT_ALL{mi} = grandtracesLUT;
    env_t_ALL{mi} = env_t;
    
    figure(f1); subplot(2,3,mi+1); hold on;
    mns = []; sems = []; trcs = cell(0);
    for cond = 1:3
                    
        if do_cleantraces
            incl = grandtracesLUT(:,2)==cond;
            get_traces = grandtraces(incl,:);
            sel_LUT = grandtracesLUT(incl,:);
            incl_ALL{mi,cond} = incl;

            prewin = env_t > -0.1 & env_t < 0;
            postwin = env_t > 0 & env_t < 0.1;
            sd_pre = std(get_traces(:,prewin),0,2);
            m_post = mean(get_traces(:,postwin),2);
            snr_sel = m_post > 3*sd_pre;

            incl = incl(snr_sel');
            get_traces = get_traces(snr_sel',:);
            newLUT = sel_LUT(snr_sel',:);

            % now average sites over sessions
            get_traces_ = get_traces; % backup
            get_traces = [];
            for rs = unique(newLUT(:,1))'
                if sum(newLUT(:,1)==rs)>1
                    get_traces = [get_traces;...
                        mean(get_traces_(newLUT(:,1)==rs,:))];
                else
                    get_traces = [get_traces; get_traces_(newLUT(:,1)==rs,:)];
                end
            end
            recsites{mi,cond}=unique(newLUT(:,1));
        else 
            incl = tracesLUT(:,2)==cond;
            get_traces = traces(incl,:);
            sel_LUT = tracesLUT(incl,:);
            incl_ALL{mi,cond} = incl;

            newLUT = sel_LUT;
            recsites{mi,cond}=unique(newLUT(:,1));
        end

        % calculate the mean across channels        
        s_traces = [];
        for t = 1:size(get_traces,1)
            s_traces(t,:) = smooth(get_traces(t,:),smoothfact);
        end
        s_traces_ALL{mi,cond} = s_traces;

        mn = mean(s_traces);
        mns(cond,:) = mn; % for calculating the difference between conditions
        trcs{cond} = s_traces; % for calculating the modulation onset
        sem = std(s_traces)./sqrt(sum(incl));
        sems(cond,:) = sem;

        t2 = [env_t*1e3, fliplr(env_t*1e3)];
        inBetween = [mn+sem, fliplr(mn-sem)];
        fill(t2, inBetween, 'g','LineStyle','none','FaceColor',cols(cond,:),...
            'FaceAlpha',falpha); hold all
        plot(env_t*1e3,mn,'Color',cols(cond,:));
        
        RT=RTs{mi}(cond);
        plot([RT RT],[0 1],'--','Color',cols(cond,:),'LineWidth',2)
        
        xlim([0 0.25*1e3]); ylim([0 1]);
    end
    title(monkeys{mi}); xlabel('time(ms)'); ylabel('MUA'); 


    % only include sites with all conditions 
    if size(recsites{mi,1},1) <= size(recsites{mi,2},1) && ...
            size(recsites{mi,1},1) <= size(recsites{mi,3},1)
        rsidx = [1 2 3];
    elseif size(recsites{mi,2},1) < size(recsites{mi,1},1) && ...
            size(recsites{mi,2},1) < size(recsites{mi,3},1)
        rsidx = [2 1 3];
    elseif size(recsites{mi,3},1) < size(recsites{mi,1},1) && ...
            size(recsites{mi,3},1) < size(recsites{mi,2},1)
        rsidx = [3 1 2];
    end
    rsinc{mi,rsidx(1)} = ismember(recsites{mi,rsidx(1)},recsites{mi,rsidx(1)});
    rsinc{mi,rsidx(2)} = ismember(recsites{mi,rsidx(2)},recsites{mi,rsidx(1)});
    rsinc{mi,rsidx(3)} = ismember(recsites{mi,rsidx(3)},recsites{mi,rsidx(1)});
        
    inRF = {'target','','distractor'};
    for cond = [1 3 4]
        figure(f1)
        subplot(2,3,mi+4);
        if cond < 4
            df = trcs{cond}(rsinc{mi,cond},:) - ...
                trcs{2}(rsinc{mi,2},:); % difference for each channel separately
            mn = mean(df);
            csem = std(trcs{cond}(rsinc{mi,cond},:) - ...
                trcs{2}(rsinc{mi,2},:))./...
                sqrt(size(trcs{cond}(rsinc{mi,cond},:),1));
        else
            df = trcs{1}(rsinc{mi,1},:) - ...
                trcs{3}(rsinc{mi,3},:);
            mn = mean(df);
            csem = std(trcs{1}(rsinc{mi,1},:) - ...
                trcs{3}(rsinc{mi,3},:))./...
                sqrt(size(trcs{1}(rsinc{mi,1},:),1));
        end
            
        inBetween = [mn+csem, fliplr(mn-csem)];
        fill(t2, inBetween, 'g','LineStyle','none','FaceColor',cols(cond,:),...
            'FaceAlpha',falpha); hold all
        plot(env_t*1e3,mn,'Color',cols(cond,:));
        xlim([0 0.25*1e3]); ylim([-0.15 0.25]);
        
        % let's calculate the latency of modulation per channel
        if do_latency
            do_latfig=0;
            if mi==1
                until = find(env_t<0.25,1,'last');
            else
                until = find(env_t<0.21,1,'last');
            end
            lat = []; coeff = []; rs = []; converge = [];
            for i = 1:size(df,1) % loop channels
                cur_resp = df(i,1:until)';
                if cond==3
                    % the latency function cannot fit downward modulation
                    cur_resp = cur_resp*-1;
                end
                [lat(i),~] = latencyfit_jp(cur_resp,env_t(1:until)'*1000,do_latfig);
            end

            % latency on average
            cur_resp = mn(1:until)';
            if cond==3
                % the latency function cannot fit downward modulation
                cur_resp = cur_resp*-1;
            end
            lt = latencyfit_jp(smooth(cur_resp,20),env_t(1:until)'*1000,do_latfig);

            if cond < 4
                rng('default') % for reproducibility
                nsamp = 100;
                btstrp_LT = bootstrp(nsamp,@latencyfit_jp,smooth(cur_resp,20),env_t(1:until)'*1000,0);
                fprintf([monkey ' ' inRF{cond} ' - Bootstrapped latency [' num2str(nsamp) ' samples] =====\n'])
                fprintf(['mean: ' num2str(mean(btstrp_LT,'omitnan')) ', std:' num2str(std(btstrp_LT,'omitnan')) '\n'])
                BLT{mi,cond}=btstrp_LT;
            end

            if ~isnan(lt) % if lt = nan, latencyfit_jp gives no figure, so we shouldn't make a title either
                if cond < 4
                    title([monkey ', fit on average response, ' inRF{cond} ' in RF']);
                else
                    title([monkey ', fit on average response modulation']);
                end
            end

            if cond < 4
                disp([monkey ', latency of modulation onset when ' inRF{cond} ' is in RF: ' ...
                    num2str(round(mean(lat,'omitnan'))) 'ms (' num2str(round(lt)) ...
                    'ms as estimated on the average response)']);
            else
                disp([monkey ', latency of targ-distr modulation: ' ...
                    num2str(round(mean(lat,'omitnan'))) 'ms (' num2str(round(lt)) ...
                    'ms as estimated on the average response)']);
            end
            %disp(['individual channel estimates: ' num2str(lat)]);

            ModLatency{mi,cond} = lt;
            figure(f1);
            plot([lt lt],[-0.2 0.25],'--','Color',cols(cond,:),'LineWidth',2)
            if isnan(lt)
                plot([round(mean(lat,'omitnan')) round(mean(lat,'omitnan'))],[-0.2 0.25],...
                    ':','Color',cols(cond,:),'LineWidth',2)
            end
        end
    end
    figure(f1)
    plot([0 .25*1e3],[0 0],'k');
    title(monkeys{mi}); xlabel('time(ms)'); ylabel('Modulation(MUA)'); 

    % save the extra cleaned traces
    save(fullfile(savedir,'Traces_SNR_trialbased.mat'),...
        's_traces_ALL','recsites','rsinc');

    %% statistics
    window_means = [];
    win_idx = find(env_t>0.15 & env_t<0.20);

    for cond = 1:3
        if do_cleantraces
            incl = grandtracesLUT(:,2)==cond;
            get_traces = grandtraces(incl,win_idx); %#ok<*FNDSB> 
        else
            incl = tracesLUT(:,2)==cond;
            get_traces = traces(incl,win_idx);
        end
        incl = tracesLUT(:,2)==cond;

        % calculate the average in a window
        s_traces = [];
        for t = 1:size(get_traces,1)
            s_traces(t,:) = smooth(get_traces(t,:),smoothfact);
        end
        window_means(:,cond) = mean(s_traces,2);
    end
    anova1(window_means);
    
    wm = [wm; window_means];
    wm2{mi} = window_means;
end

figure(f1); subplot(2,3,1); hold on;
% pooled animals
for cond = 1:3
    traces = [s_traces_ALL{1,cond};s_traces_ALL{2,cond}];

    % calculate the mean across channels
    % get_traces = traces(incl,:);
    get_traces = traces;
    s_traces = [];
    for t = 1:size(get_traces,1)
        s_traces(t,:) = smooth(get_traces(t,:),smoothfact);
    end

    mn = mean(s_traces);
    mns(cond,:) = mn; % for calculating the difference between conditions
    trcs{cond} = s_traces; % for calculating the modulation onset
    sem = std(s_traces)./sqrt(sum(size(s_traces,1)));
    sems(cond,:) = sem;
    
    t2 = [env_t*1e3, fliplr(env_t*1e3)];
    inBetween = [mn+sem, fliplr(mn-sem)];
    fill(t2, inBetween, 'g','LineStyle','none','FaceColor',cols(cond,:),...
        'FaceAlpha',falpha); hold all
    
    RT=RTs{3}(cond);
    plot([RT RT],[0 1],'--','Color',cols(cond,:),'LineWidth',2)
    plot(env_t*1e3,mn,'Color',cols(cond,:));
    xlim([0 0.25*1e3]); ylim([0 1]);
end
title('BOTH animals pooled'); xlabel('time(ms)'); ylabel('MUA'); 

inRF = {'target','','distractor'};

for cond=1:3
    rsinc{3,cond} = [rsinc{1,cond};rsinc{2,cond}];
end

for cond = [1 3 4]
    figure(f1)
    subplot(2,3,4);

    if cond < 4
        df = trcs{cond}(rsinc{3,cond},:) - ...
            trcs{2}(rsinc{3,2},:); % difference for each channel separately
        mn = mean(df);
        csem = std(trcs{cond}(rsinc{3,cond},:) - ...
            trcs{2}(rsinc{3,2},:))./...
            sqrt(size(trcs{cond}(rsinc{3,cond},:),1));
    else
        df = trcs{1}(rsinc{3,1},:) - ...
            trcs{3}(rsinc{3,3},:);
        mn = mean(df);
        csem = std(trcs{1}(rsinc{3,1},:) - ...
            trcs{3}(rsinc{3,3},:))./...
            sqrt(size(trcs{1}(rsinc{3,1},:),1));
    end
    
    inBetween = [mn+csem, fliplr(mn-csem)];
    fill(t2, inBetween, 'g','LineStyle','none','FaceColor',cols(cond,:),...
        'FaceAlpha',falpha); hold all
    plot(env_t*1e3,mn,'Color',cols(cond,:));
    xlim([0 0.25*1e3]); ylim([-0.15 0.25]);
    
    % let's calculate the latency of modulation per channel
    if do_latency
        do_latfig=0;
        if mi==1
            until = find(env_t<0.25,1,'last');
        else
            until = find(env_t<0.21,1,'last');
        end
        lat = []; coeff = []; rs = []; converge = [];
        for i = 1:size(df,1) % loop channels
            cur_resp = df(i,1:until)';
            if cond==3
                % the latency function cannot fit downward modulation
                cur_resp = cur_resp*-1;
            end
            [lat(i),~] = latencyfit_jp(cur_resp,env_t(1:until)'*1000,do_latfig);
        end
        
        %collect latencty
        clat{cond} = lat;

        % latency on average
        do_latfig=1;
        cur_resp = mn(1:until)';
        if cond==3
            % the latency function cannot fit downward modulation
            cur_resp = cur_resp*-1;
        end
        lt = latencyfit_jp(smooth(cur_resp,20),env_t(1:until)'*1000,do_latfig);

        if cond < 4
            rng('default') % for reproducibility
            nsamp = 100;
            btstrp_LT = bootstrp(nsamp,@latencyfit_jp,smooth(cur_resp,20),env_t(1:until)'*1000,0);
            fprintf([inRF{cond} ' - Bootstrapped latency [' num2str(nsamp) ' samples] =====\n'])
            fprintf(['mean: ' num2str(mean(btstrp_LT,'omitnan')) ', std:' num2str(std(btstrp_LT,'omitnan')) '\n'])
            BLT{3,cond}=btstrp_LT;
        end

        if ~isnan(lt) % if lt = nan, latencyfit_jp gives no figure, so we shouldn't make a title either
            if cond < 4
                title(['BOTH monkeys, fit on average response, ' inRF{cond} ' in RF']);
            else
                title('BOTH monkeys, fit on average response modulation');
            end
        end
        if cond < 4
            disp(['BOTH monkeys, latency of modulation onset when ' inRF{cond} ' is in RF: ' ...
                num2str(round(mean(lat,'omitnan'))) 'ms (' num2str(round(lt)) ...
                'ms as estimated on the average response)']);
        else
            disp(['BOTH monkeys, latency of targ-distr modulation: ' ...
                num2str(round(mean(lat,'omitnan'))) 'ms (' num2str(round(lt)) ...
                'ms as estimated on the average response)']);
        end
        %disp(['individual channel estimates: ' num2str(lat)]);

        figure(f1)
        plot([lt lt],[-0.2 0.25],'--','Color',cols(cond,:),'LineWidth',2)
        if isnan(lt)
            plot([round(mean(lat,'omitnan')) round(mean(lat,'omitnan'))],[-0.2 0.25],...
                ':','Color',cols(cond,:),'LineWidth',2)
        end

    end
end
figure(f1)
plot([0 .25*1e3],[0 0],'k');
title('ALL animals pooled'); xlabel('time(ms)'); ylabel('Modulation (MUA)'); 

%% do a bootstrapped latency analysis to get some stats
rng('default') % for reproducibility
nsamp = 1000;
btstrp_LAT = bootstrp(nsamp,@nanmean,[clat{1}' clat{3}']);
[h,p,ci,stats]=ttest(btstrp_LAT(:,1)-btstrp_LAT(:,2));
fprintf('== BTSTRP comparison of latencies ===\n');
fprintf(['Mean difference in latency (T-SD): ' ...
    num2str(mean(btstrp_LAT(:,1)-btstrp_LAT(:,2))) ', SEM: ' ...
    num2str(std(btstrp_LAT(:,1)-btstrp_LAT(:,2))./sqrt(size(btstrp_LAT,1))) ...
    'ms\n']);
fprintf([num2str(sum((btstrp_LAT(:,1)-btstrp_LAT(:,2))<0)) '/' num2str(nsamp) ...
    ' T before SD, p = ' num2str(sum((btstrp_LAT(:,1)-btstrp_LAT(:,2))>0)/nsamp) ...
    ' that H0[T not before SD] is rejected \n']);
fprintf(['ttest t(' num2str(stats.df) ') = ' num2str(stats.tstat) ...
    ', p = ' num2str(p) '\n']);

%% use the bootstrapped latency estimates for stats
fprintf('== Comparison of bootstrapped T-SD latencies ===\n');
for mi=1:3
    [h,p,ci,stats]=ttest(BLT{mi,1}-BLT{mi,3});
    nn = sum(~isnan(BLT{mi,1}-BLT{mi,3}));
    if mi<3
        fprintf(['-- Monkey ' num2str(mi) ' --\n']);
    else
        fprintf('-- BOTH Monkeys --\n');
    end
    fprintf(['Mean difference in latency (T-SD): ' ...
        num2str(mean(BLT{mi,1}-BLT{mi,3},'omitnan')) ', SEM: ' ...
        num2str(std(BLT{mi,1}-BLT{mi,3},'omitnan')./sqrt(nn)) ...
        'ms\n']);
    fprintf([num2str(sum((BLT{mi,1}-BLT{mi,3})<0)) '/' num2str(nn) ...
        ' T before SD, p = ' num2str(sum((BLT{mi,1}-BLT{mi,3})>0)/nn) ...
        ' that H0[T not before SD] is rejected \n']);
    fprintf(['ttest t(' num2str(stats.df) ') = ' num2str(stats.tstat) ...
        ', p = ' num2str(p) '\n']);
end

%% statistics
% do a simple paired t-test for each monkey
fprintf('=============================\n');
fprintf(' Paired T-test \n');
fprintf('=============================\n');
clear h p stats
for mi = 1:2
    mns = wm2{mi};
    p = []; vp=[]; ci=1;
    for cond = [1 3]
        [vh,vp(end+1)]=vartest2(mns(:,cond),mns(:,2));
        [h,p(end+1),~,stats{mi,ci}] = ttest(mns(:,cond), mns(:,2));
        ci=ci+1;
    end
    [vh,vp(end+1)]=vartest2(mns(:,1),mns(:,3));
    [h,p(end+1),~,stats{mi,3}] = ttest(mns(:,1), mns(:,3));
        
    disp([monkeys{mi} ', target response increased: p=' num2str(p(1)) ...
        ' t = ' num2str(stats{mi,1}.tstat) ', df = ' num2str(stats{mi,1}.df) ... 
        ' (unequal var p = ' num2str(vp(1)) ')']);
    disp([monkeys{mi} ', distractor response decreased: p=' num2str(p(2)) ...
        ' t = ' num2str(stats{mi,2}.tstat) ', df = ' num2str(stats{mi,2}.df) ... 
        ' (unequal var p = ' num2str(vp(2)) ')']);
    disp([monkeys{mi} ', targ vs distractor: p=' num2str(p(3)) ...
        ' t = ' num2str(stats{mi,3}.tstat) ', df = ' num2str(stats{mi,3}.df) ... 
        ' (unequal var p = ' num2str(vp(3)) ')']);
end

%Both together
mns = [wm2{1}; wm2{2}];
p = [];vp=[]; ci=1;
for cond = [1 3]
    [vh,vp(end+1)]=vartest2(mns(:,cond),mns(:,2));
    [h,p(end+1),~,stats{3,ci}] = ttest(mns(:,cond), mns(:,2));
    ci=ci+1;
end
[vh,vp(end+1)]=vartest2(mns(:,1),mns(:,3));
[h,p(end+1),~,stats{3,3}] = ttest(mns(:,1), mns(:,3));

disp(['BOTH, target response increased: p=' num2str(p(1)) ...
    ' t = ' num2str(stats{3,1}.tstat) ', df = ' num2str(stats{3,1}.df) ... 
    ' (unequal var p = ' num2str(vp(1)) ')']);
disp(['BOTH, distractor response decreased: p=' num2str(p(2)) ...
    ' t = ' num2str(stats{3,2}.tstat) ', df = ' num2str(stats{3,2}.df) ... 
    ' (unequal var p = ' num2str(vp(2)) ')']);
disp(['BOTH, targ vs distractor: p=' num2str(p(3)) ...
    ' t = ' num2str(stats{3,3}.tstat) ', df = ' num2str(stats{3,3}.df) ... 
    ' (unequal var p = ' num2str(vp(3)) ')']);


%% save figure
figure(f1)
saveas(gcf,fullfile(savedir,'figure_neurotrace.svg'));