%% Script to analyse individual data from the study examining the stopping of
% planned versus ongoing actions in the same stop signal task.
% Task involves reaching actions performed with a computer mouse.

%% Code is suitable for use with task code "StopWhatYoureDoing15"

% Some variables have been removed from csv file, so requires new import
% options

%% 
clc; clear; close all;
set(0,'DefaultFigureWindowStyle','docked');

% All subs in first pass run before filtering
subjs = {'81ut_Pilot15', 'k7kb_Pilot15', 'u75q_Pilot15', 'dj19_Pilot15'};

% Good subs (see filter in _allSubs.m script)
% subjs = { };


% get the directory of your input files:
pathnamedat = strcat(pwd, '/processed');
pathnamefigs = strcat(pwd, '/indie_figs');

%% Setup the Import Options
opts = delimitedTextImportOptions("NumVariables", 50);

% Specify range and delimiter
opts.DataLines = [2, Inf];
opts.Delimiter = ",";

% Specify column names and types
opts.VariableNames = ["view_history", "rt", "CueType", "trial_type", "trial_index", "time_elapsed", "internal_node_id", "subject", "browser_name", "browser_version", "os_name", "os_version", "screen_resolution", "window_resolution", "Targ1Ycent", "Targ1XLcent", "Targ1XRcent", "TargetWidth", "TargetHeight", "HomeXcent", "HomeYcent", "HomeHeight", "HomeWidth", "responses", "question_order", "response_type", "key_press", "avg_frame_time", "center_x", "center_y", "click_x", "click_y", "CueDirection", "TrialType", "Staircase", "BlockName", "BlockNo", "Xpos", "Ypos", "IsHome", "TooSoon", "Jitter", "ReactTime", "MoveTime", "TargetHit", "Outcome", "StopSignalDelay", "MouseX", "MouseY", "MouseTime"];
opts.VariableTypes = ["string", "categorical", "categorical", "categorical", "double", "double", "string", "categorical", "categorical", "double", "categorical", "categorical", "categorical", "categorical", "double", "double", "double", "double", "double", "double", "double", "double", "double", "string", "double", "categorical", "categorical", "double", "double", "double", "double", "double", "categorical", "categorical", "string", "categorical", "double", "double", "double", "categorical", "categorical", "double", "double", "double", "categorical", "categorical", "double", "double", "double", "double"];
opts = setvaropts(opts, [1, 7, 24, 35], "WhitespaceRule", "preserve");
opts = setvaropts(opts, [25, 48, 49, 50], "TrimNonNumeric", true);
opts = setvaropts(opts, [25, 48, 49, 50], "ThousandsSeparator", ",");
opts = setvaropts(opts, [1, 2, 3, 4, 7, 8, 9, 11, 12, 13, 14, 24, 26, 27, 33, 34, 35, 36, 40, 41, 45, 46], "EmptyFieldRule", "auto");
opts.ExtraColumnsRule = "ignore";
opts.EmptyLineRule = "read";


%% Loop through importing and analysing individual subject data
for i = 4%1:length(subjs) % -----------------------------------------------------------------------
    
    subjID  = string(subjs(i));
    disp(subjID);
    fileName = strcat("IndieData_", subjID, ".mat");
    folderName = strcat('/Users/riccihannah/Google Drive/Online_Stop_Signal/Online_Stop_Task/Task_replication/data/indie_figs/',num2str(i),'_',subjID);
    mkdir(folderName);
    
    fileDat = strcat("/Users/riccihannah/Google Drive/Online_Stop_Signal/Online_Stop_Task/Task_replication/data/Blockwise_PvO_Jstop_", subjID ,".csv");
    dataTable = readtable(fileDat, opts);
    
    % Get demographics
    prolificInfo = dataTable.responses(4);
    prolificID = extractBetween(prolificInfo,17,40);
    age = extractBetween(prolificInfo,50,51);
    
    prolificDemog = dataTable.responses(5);
    sex = extractBetween(prolificDemog,10,10);
    handedness = extractBetween(prolificDemog,"handedness","hispanic");
    handedness = extractBetween(handedness,5,5);
    hispanicYN = extractBetween(prolificDemog, "hispanic" ,"race");
    hispanicYN = extractBetween(hispanicYN,5,5);
    race = extractBetween(prolificDemog, "race" ,"]");
    if length(race) < 3
        race = 'NDisc';
    else
        race = extractBetween(race,5,9);
    end
    
    demographics = string([subjID, prolificID, age, sex, handedness, hispanicYN, race]);
    data.Demographics = array2table(demographics, 'VariableNames', {'SubID', ...
        'ProlificID', 'Age', 'Sex', 'Handedness', 'Hispanic', 'Race'});
    
    % Now get rid of rows that don't pertain to the Go/Stop cues
    trl_GoStop = dataTable.BlockName == 'MainStopExp' & dataTable.CueType == 'Go/Stop';
    dataTable = dataTable(trl_GoStop,:);
    
    % Get monitor, browser and OS data
    frameRate = nanmedian(dataTable.avg_frame_time);
    browser = dataTable.browser_name(1);
    OSname = CatStr(dataTable.os_name(1));
    OSversion = dataTable.os_version(1);
    screenRes = CatStr(dataTable.screen_resolution(1));
    screenResWidth = (extractBefore(screenRes,' x'));
    screenResHeight = (extractAfter(screenRes,'x '));
    windowRes = CatStr(dataTable.window_resolution(1));
    windowResWidth = (extractBefore(windowRes,' x'));
    windowResHeight = (extractAfter(windowRes,'x '));
    centerSreen_x = dataTable.center_x(1);
    centerSreen_y = dataTable.center_y(1);
    targetYcenter = dataTable.Targ1Ycent(1);
    targetXLcenter = dataTable.Targ1XLcent(1);
    targetXRcenter = dataTable.Targ1XRcent(1);
    
    % Work out distance from home to target
    homeXpos =  str2double(windowResWidth) / 2;
    homeYpos = str2double(windowResHeight) * 0.85;
    homeWidth = str2double(windowResWidth) * 0.006;
    homeHeight = homeWidth;
    leaveHomeHeightY = homeYpos - homeHeight;
    
    targWidth = str2double(windowResWidth) * 0.025;
    targHeight = targWidth;
    enterTargetY = targetYcenter + (targHeight / 2);
    
    homeToTarget_Ydist = leaveHomeHeightY - enterTargetY;
    
    computer_setup = ([string(frameRate), string(browser),OSname ,string(OSversion) ,...
        screenResWidth,screenResHeight ,windowResWidth , ...
        windowResHeight , string(centerSreen_x), string(centerSreen_y),string(targetYcenter),...
        string(targetXLcenter), string(targetXRcenter), ...
        leaveHomeHeightY, enterTargetY, homeToTarget_Ydist]);
    data.Computer_setup = array2table(computer_setup, 'VariableNames', {'frameRate',...
        'Browser','OSname' ,'OSversion' ,  ...
        'screenResWidth','screenResHeight' ,'windowResWidth' , ...
        'windowResHeight' , 'centerSreen_x', 'centerSreen_y','targetYcenter',...
        'targetXLcenter', 'targetXRcenter',...
        'leaveHomeYlim','enterTargetYlim','home2targetYdist'});
    
    %% Remove unused variables from the table
    dataTable = removevars(dataTable, {'rt','response_type','key_press','trial_type','trial_index',...
        'time_elapsed', 'internal_node_id','view_history','responses','question_order','Xpos','Ypos',...
        'click_x','click_y','IsHome','MouseX','MouseY','MouseTime',...
        'avg_frame_time','browser_name','browser_version','os_name','os_version','screen_resolution','window_resolution',...
        'Targ1XLcent', 'Targ1XRcent', 'center_x','center_y'});
    
    %% Bit to get a matrix of mouseX and mouseY coordinates and times
    csv = readtable(fileDat);
    names = {'MouseX', 'MouseY', 'MouseTime'};
    array_outcome = csv{:,'Outcome'};
    array_blockName = csv{:, 'BlockName'};
    name_mainExp = 'MainStopExp';
    %.................................................... problem is that now
    % have practice go, practice stop and main exp all in same file, so need to
    % only focus on main exp here
    for n = 1: length(names)
        currentCol = csv.Properties.VariableNames{names(n)};
        array = csv{:,{currentCol}};
        %     array_cell = array(cellfun(@(x) ~isempty(x), array_outcome)); old
        
        array_cell = array(cellfun(@(x) strcmp(x,name_mainExp), array_blockName));
        array_outcome_new = array_outcome(cellfun(@(x) strcmp(x,name_mainExp), array_blockName));
        array_cell = array_cell(cellfun(@(x) ~isempty(x), array_outcome_new));
        
        array_cell = cellfun(@(x) str2double(strsplit(x,',')),array_cell,'uni',0);
        mouseXYT(:,:,n) = padcat(array_cell{:,:}); % This function pads the array with NaNs as it concatenates
    end
    
    clear csv array_outcome;
    
    %% Index trials
    % NaN RT and MT trials == 0
    dataTable.ReactTime(dataTable.ReactTime == 0) = NaN;
    dataTable.MoveTime(dataTable.MoveTime == 0) = NaN;
    
    % Go trials
    trials.GoAll = dataTable.TrialType == 'Go' ;
    trials.GoCorr =  dataTable.Outcome == 'Correct_Go';
    trials.GoErr =  dataTable.TrialType == 'Go' & dataTable.TargetHit == 'Miss' | ...
        dataTable.TrialType == 'Go' & dataTable.Outcome == 'False_Alarm';
    
    trials.pcntErr_go = round(sum(trials.GoErr) / sum(trials.GoAll) * 100);
    
    % False alarms
    trials.FalseAlarm_allTrl = dataTable.Outcome == 'False_Alarm' & dataTable.TargetHit ~= 'Miss';
    sum(trials.FalseAlarm_allTrl);
    trials.pcnt_FalseAlarm_allTrl = round(sum(trials.FalseAlarm_allTrl) / size(dataTable,1) * 100);
    
    % Omitted responses or missed target in go trials
    trials.Omit_go = dataTable.TrialType == 'Go' & isnan(dataTable.ReactTime);
    trials.Miss_go = dataTable.TrialType == 'Go' & dataTable.Outcome == 'Miss';
    trials.pcnt_omit_go = round(sum(trials.Omit_go) / sum (trials.GoAll) * 100);
    trials.pcnt_miss_go = round(sum(trials.Miss_go) / sum (trials.GoAll) * 100);
    
    % Overall omission or miss - classified simply as a miss
    trials.Miss = dataTable.Outcome == 'Miss';
    sum(trials.Miss)
    
    % Stop trials with the Homepad staircase
    trials.StopAll_HomeStr = dataTable.Staircase == 'HomeStair' ;
    trials.StopCorr_HomeStr = dataTable.Staircase == 'HomeStair' & dataTable.Outcome == 'Correct_Stop_Home' ;
    trials.StopFail_HomeStr = dataTable.Staircase == 'HomeStair' & dataTable.Outcome == 'Failed_Stop_Home' | ...
        dataTable.Staircase == 'HomeStair' & dataTable.Outcome == 'Failed_Stop_Target' ;
    
    % Stop trials with the Target staircase
    trials.StopAll_TargetStr = dataTable.Staircase == 'TargetStair' & ~isnan(dataTable.StopSignalDelay);
    trials.StopCorr_TargetStr = dataTable.Staircase == 'TargetStair' & dataTable.Outcome == 'Correct_Stop_Home' | ...
        dataTable.Staircase == 'TargetStair' & dataTable.Outcome == 'Failed_Stop_Home' ;
    trials.StopFail_TargetStr = dataTable.Staircase == 'TargetStair' & dataTable.Outcome == 'Failed_Stop_Target' ;
    
    % Omitted and missed targets in stop trials
    trials.Omit_home = dataTable.Staircase == 'HomeStair' & isnan(dataTable.ReactTime);
    trials.Miss_home = dataTable.Staircase == 'HomeStair' & dataTable.Outcome == 'Miss';
    trials.pcnt_omit_home = sum(trials.Omit_home) / sum (trials.StopAll_HomeStr) * 100;
    trials.pcnt_miss_home = sum(trials.Miss_home) / sum (trials.StopAll_HomeStr) * 100;
    
    trials.Omit_target = dataTable.Staircase == 'TargetStair' & isnan(dataTable.ReactTime);
    trials.Miss_target = dataTable.Staircase == 'TargetStair' & dataTable.Outcome == 'Miss';
    trials.pcnt_omit_target = sum(trials.Omit_target) / sum (trials.StopAll_TargetStr) * 100;
    trials.pcnt_miss_target = sum(trials.pcnt_miss_go) / sum (trials.StopAll_TargetStr) * 100;
    
    % How many times do they hit the target on home staircase? i.e. a true trigger failure
    trials.Trigger_failure  = dataTable.Staircase == 'HomeStair' & dataTable.MoveTime > 0;
    n_trigger_failure = sum(trials.Trigger_failure);
    
    % Sanity check - display number of Go and Stop trials, and proportion
    disp(strcat('Total trials =  ', num2str(size(dataTable,1))));
    disp(strcat('Go trials =  ', num2str(sum(trials.GoAll))));
    nStoptrls = sum(trials.StopAll_HomeStr) + sum(trials.StopAll_HomeStr);
    GoStopRatio = nStoptrls / size(dataTable,1);
    disp(strcat('Stop trials =  ', num2str(nStoptrls)));
    disp(strcat('Proportion stop trials =  ', num2str(GoStopRatio)));
    disp(strcat('Stop Home trials =  ', num2str(sum(trials.StopAll_HomeStr))));
    disp(strcat('Stop Target trials =  ', num2str(sum(trials.StopAll_HomeStr)))); 


    %% Filter and Get basic data
    
    % Remove "miss" trials
    dataTable.ReactTime(dataTable.Outcome == 'Miss') = NaN;
    dataTable.MoveTime(dataTable.Outcome == 'Miss') = NaN;
    
    % All RTs and MT1s
    RT_all = dataTable.ReactTime;
    MT_all = dataTable.MoveTime;
    
    % Remove outliers from RTs -
    prctile_RT          = prctile(dataTable.ReactTime(trials.GoAll),[25 75]);
    trials.RT_outlier = (dataTable.ReactTime < (prctile_RT(1) - 1.5 * iqr(dataTable.ReactTime(trials.GoAll))) | ...
        dataTable.ReactTime >(prctile_RT(2)+ 1.5 * iqr(dataTable.ReactTime(trials.GoAll))) | ...
        dataTable.ReactTime < 150);
    dataTable.ReactTime(dataTable.ReactTime < (prctile_RT(1) - 1.5 * iqr(dataTable.ReactTime(trials.GoAll))) | ...
        dataTable.ReactTime >(prctile_RT(2)+ 1.5 * iqr(dataTable.ReactTime(trials.GoAll))) | ...
        dataTable.ReactTime < 150) = NaN;
    % Remove outliers from MTs - gets complicated because MTs vary if Go and
    % Jump
    prctile_MT          = prctile(dataTable.MoveTime(trials.GoAll),[25 75]);
    trials.MT_outlier = (dataTable.MoveTime < (prctile_MT(1) - 1.5 * iqr(dataTable.MoveTime(trials.GoAll))) | ...
        dataTable.MoveTime>(prctile_MT(2) + 1.5 * iqr(dataTable.MoveTime(trials.GoAll))) | ...
        dataTable.MoveTime < 150);
    dataTable.MoveTime(dataTable.MoveTime < (prctile_MT(1) - 1.5 * iqr(dataTable.MoveTime(trials.GoAll))) | ...
        dataTable.MoveTime>(prctile_MT(2) + 1.5 * iqr(dataTable.MoveTime(trials.GoAll))) | ...
        dataTable.MoveTime < 150) = NaN;
    
    % All RTs and MT1s
    RT_GoAll = dataTable.ReactTime(trials.GoAll);
    MT_GoAll = dataTable.MoveTime(trials.GoAll);
    
    % Go RT and MT
    RT_GoCorr = dataTable.ReactTime(trials.GoCorr);
    RT_GoErr = dataTable.ReactTime(trials.GoErr);
    MT_GoCorr = dataTable.MoveTime(trials.GoCorr);
    MT_GoErr = dataTable.MoveTime(trials.GoErr);
    
    % Correct Stop RT (Target staircase)
    RT_StopCorr_TargetStair = dataTable.ReactTime(trials.StopCorr_TargetStr);
    MT_StopCorr_TargetStair = dataTable.MoveTime(trials.StopCorr_TargetStr);
    
    % Failed Stop RT and MT
    RT_StopFail_HomeStair = dataTable.ReactTime(trials.StopFail_HomeStr);
    MT_StopFail_HomeStair = dataTable.MoveTime(trials.StopFail_HomeStr);
    RT_StopFail_TargetStair = dataTable.ReactTime(trials.StopFail_TargetStr);
    MT_StopFail_TargetStair = dataTable.MoveTime(trials.StopFail_TargetStr);
    
    % Total response times
    TotRespT_GoCorr = RT_GoAll + MT_GoAll;
    TotRespT_StopFail_HomeStair = RT_StopFail_HomeStair + MT_StopFail_HomeStair;
    TotRespT_StopFail_TargetStair = RT_StopFail_TargetStair + MT_StopFail_TargetStair;
    
    % Mean Go RT and MT
    RT_GoCorr_mn = round(nanmedian(RT_GoCorr));
    MT_GoCorr_mn = round(nanmedian(MT_GoCorr));
    
    % Mean Correct Stop RT (Target staircase)
    RT_StopCorr_TargetStair_mn = round(nanmedian(RT_StopCorr_TargetStair));
    MT_StopCorr_TargetStair_mn = round(nanmedian(MT_StopCorr_TargetStair));
    
    % Mean Failed Stop Change RT and MT
    RT_StopFail_HomeStair_mn = round(nanmedian(RT_StopFail_HomeStair));
    MT_StopFail_HomeStair_mn = round(nanmedian(MT_StopFail_HomeStair));
    RT_StopFail_TargetStair_mn = round(nanmedian(RT_StopFail_TargetStair));
    MT_StopFail_TargetStair_mn = round(nanmedian(MT_StopFail_TargetStair));
    
    % Mean Total response times
    TotRespT_GoCorr_mn = round(nanmedian(TotRespT_GoCorr));
    TotRespT_StopFail_HomeStair_mn = round(nanmedian(TotRespT_StopFail_HomeStair));
    TotRespT_StopFail_TargetStair_mn = round(nanmedian(TotRespT_StopFail_TargetStair));
    
    % Standard deviation of RT, MT and TT on Go trials
    RT_std_GoCorr = round(nanstd(RT_GoCorr));
    MT_std_GoCorr = round(nanstd(MT_GoCorr));
    TT_std_GoCorr = round(nanstd(TotRespT_GoCorr));
    
    % All SSDs
    SSD_Stop_HomeStair = round(dataTable.StopSignalDelay(trials.StopAll_HomeStr));
    SSD_Stop_TargetStair = round(dataTable.StopSignalDelay(trials.StopAll_TargetStr));
    
    % Mean SSD
    SSD_Stop_HomeStair_mn = round(nanmedian(SSD_Stop_HomeStair));
    SSD_Stop_TargetStair_mn = round(nanmedian(SSD_Stop_TargetStair));
    
    % Stop success
    pStop_HomeStair = round((sum(trials.StopCorr_HomeStr) / sum(sum([trials.StopCorr_HomeStr,trials.StopFail_HomeStr]))) *100);
    pStop_TargetStair = round((sum(trials.StopCorr_TargetStr) / sum(sum([trials.StopCorr_TargetStr,trials.StopFail_TargetStr]))) *100);
    
    
    %% Plot RTs and MTs
    figure()
    
    subplot(2,2,1);
    histogram(RT_GoCorr,'BinWidth',25, 'Normalization','probability', 'FaceColor','g');
    hold on
    %     histogram((0:25:200),'BinWidth',25, 'Normalization','probability');
    histogram(RT_StopCorr_TargetStair,'BinWidth',25, 'Normalization','probability','FaceColor','b');
    title('RTs: Go & Corr. Stop');
    xlabel('Reaction Times (ms)');
    ylabel('Probability');
    xlim([0,round(max(RT_GoCorr) + 100)]);
    ylim([0,0.6]);
    tx1 = (['Go: ' num2str(round(RT_GoCorr_mn)) 'ms']);
    tx2 = ('SS-H: N/A ms');
    tx3 = (['SS-T: ' num2str(round(RT_StopCorr_TargetStair_mn)) 'ms']);
    text(40, 0.38,tx1);
    text(40, 0.32,tx2);
    text(40, 0.26,tx3);
    legend('Go','SS-T');
    
    subplot(2,2,2)
    histogram(RT_GoCorr,'BinWidth',25, 'Normalization','probability','FaceColor','g');
    hold on
    histogram(RT_StopFail_HomeStair,'BinWidth',25, 'Normalization','probability','FaceColor','r','LineStyle','-.');
    histogram(RT_StopFail_TargetStair,'BinWidth',25, 'Normalization','probability','FaceColor','m');
    title('RTs: Go & Fail. Stop');
    xlabel('Reaction Times (ms)');
    ylabel('Probability');
    xlim([0,round(max(RT_GoCorr) + 100)]);
    ylim([0,0.6]);
    tx1 = (['Go: ' num2str(round(RT_GoCorr_mn)) 'ms']);
    tx2 = (['FS-H: ' num2str(round(RT_StopFail_HomeStair_mn)) 'ms']);
    tx3 = (['FS-T: ' num2str(round(RT_StopFail_TargetStair_mn)) 'ms']);
    text(40, 0.38,tx1);
    text(40, 0.32,tx2);
    text(40, 0.26,tx3);
    legend('Go','FS-H','FS-T');
    
    subplot(2,2,3)
    histogram(MT_GoCorr,'BinWidth',25, 'Normalization','probability','FaceColor','g')
    hold on
    histogram(MT_StopFail_HomeStair,'BinWidth',25, 'Normalization','probability','FaceColor','r','LineStyle','-.');
    histogram(MT_StopFail_TargetStair,'BinWidth',25, 'Normalization','probability','FaceColor','m');
    title('MTs: Go & Fail. Stopp');
    xlabel('Movement Times (ms)');
    ylabel('Probability');
    xlim([0,round(max(MT_GoCorr) + 100)]);
    ylim([0,0.6]);
    tx1 = (['Go MT: ' num2str(round(MT_GoCorr_mn)) 'ms']);
    tx2 = ('SS-H: N/A ms');
    tx3 = (['FS-F MT: ' num2str(round(MT_StopFail_TargetStair_mn)) 'ms']);
    text(20, 0.38,tx1);
    text(20, 0.32,tx2);
    text(20, 0.26,tx3);
    legend('Go','FS-H','FS-T');
    
    subplot(2,2,4);
    histogram(TotRespT_GoCorr,'BinWidth',20, 'Normalization','probability','FaceColor','g')
    hold on
    %     histogram((0:25:200),'BinWidth',25, 'Normalization','probability');
    histogram(TotRespT_StopFail_TargetStair,'BinWidth',20, 'Normalization','probability','FaceColor','m')
    title('Tot Time: Go & Fail. Stopp');
    xlabel('Total Times (ms)');
    ylabel('Probability');
    xlim([0,round(max(TotRespT_GoCorr) + 100)]);
    ylim([0,0.3]);
    tx1 = (['Go: ' num2str(round(nanmedian(TotRespT_GoCorr))) 'ms']);
    tx2 = (['FS-H: ' num2str(round(nanmedian(TotRespT_StopFail_HomeStair))) 'ms']);
    tx3 = (['FS-T: ' num2str(round(nanmedian(TotRespT_StopFail_TargetStair))) 'ms']);
    text(40, 0.28,tx1);
    text(40, 0.22,tx2);
    text(40, 0.16,tx3);
    legend('Go','FS-T');
    
    print(gcf,strcat(folderName,'/',subjID,'_PvO_Stop_RT_MT.pdf'),'-dpdf','-r400');
    
    % print(gcf,strcat('indie_figs/',subjID,'_PvO_Stop_RT_MT.pdf'),'-dpdf','-r400');
    
    %% Calculate SSRT for RTs and MTs
    % Home
    if sum(trials.Miss) > 0
        new_RT = dataTable.ReactTime;
        new_RT(trials.Miss) = nanmax(RT_GoCorr);
        new_RT = new_RT(trials.GoAll);
        
        new_MT = dataTable.MoveTime;
        new_MT(trials.Miss) = nanmax(MT_GoCorr);
        new_MT = new_MT(trials.GoAll);
        
        new_TotRespT_GoCorr = new_RT + new_MT;
    else
        new_RT = dataTable.ReactTime(trials.GoAll);
        new_MT = dataTable.MoveTime(trials.GoAll);
        new_TotRespT_GoCorr = new_RT + new_MT;
    end
    
    % Home stair case - using RTs
    sort_RT = sort(new_RT(~isnan(new_RT)));
    nth_RT_index_HomeStair  = round((pStop_HomeStair/100).*length(new_RT));
    if nth_RT_index_HomeStair > 0
        nth_RT_HomeStair        = sort_RT(nth_RT_index_HomeStair);
        SSRT_HomeStair          = nth_RT_HomeStair - SSD_Stop_HomeStair_mn;
    else
        SSRT_HomeStair = NaN;
    end
    
    % Target stair case using movement time *** You get the very similar
    % answer as with TT because you still have to include RT in the
    % calc
    sort_MT         = sort(new_MT(~isnan(new_MT)));
    nth_MT_index_TargetStair = round((pStop_TargetStair/100).*length(new_MT));
    if nth_MT_index_TargetStair > 0
        nth_MT_TargetStair        = sort_MT(nth_MT_index_TargetStair);
        SSRT_TargetStair_MT         = nth_MT_TargetStair - SSD_Stop_TargetStair_mn + nanmedian(new_RT);
    else
        SSRT_TargetStair_MT = NaN;
    end
    
    
    % Target stair case using Total response time
    sort_TotRespT            = sort(new_TotRespT_GoCorr(~isnan(new_TotRespT_GoCorr)));
    nth_TRT_index_TargetStair = round((pStop_TargetStair/100).*length(sort_TotRespT));
    if nth_TRT_index_TargetStair > 0
        nth_TRT_TargetStair        = sort_TotRespT(nth_TRT_index_TargetStair);
        SSRT_TargetStair_TT         = nth_TRT_TargetStair - SSD_Stop_TargetStair_mn;
    else
        SSRT_TargetStair_TT = NaN;
    end
    
    
    
    %% Plot stopping staircases
    % Need to first calculate pStop at each SSD
    % how many SSDs?
    SSD_uniq_HomeStair = unique(SSD_Stop_HomeStair);
    SSD_uniq_TargetStair = unique(SSD_Stop_TargetStair);
    max_ntrls = max(length(SSD_uniq_HomeStair), length(SSD_uniq_TargetStair));
    
    pStop_SSD = NaN(max_ntrls,2);
    FailedStop_RT_ToT_SSD = NaN(max_ntrls,2);
    ZRFT_SSD = NaN(max_ntrls,2);
    trlwise_FSRT_GoRT_pcnt = NaN(max_ntrls,2);
    % -----------------------------------------------------
    % ZRFT - normalise the functions according to the variability in response times
    % ZRFT=(mean[RTgo]-SSD-SSRT)/SD[RTgo]
    % -----------------------------------------------------
    
    % -----------------------------------------------------
    % Should maybe consider removing SSDs with < 1 or 2 trial when comparing function?
    % Or just see what Logan has done in the past
    % -----------------------------------------------------
    trials.RTisNaN = find(isnan(dataTable.ReactTime));
    trials.MTisNaN = find(isnan(dataTable.MoveTime));
    
    % Home stair case
    for s = 1:length(SSD_uniq_HomeStair)
        % for each SSD
        currentSSD = SSD_uniq_HomeStair(s);
        % index and count how many trials for that SSD in that stair case
        idx_all = dataTable.StopSignalDelay(trials.StopAll_HomeStr) == currentSSD;
        nTrls_all = sum(idx_all);
        % index and count the the successful stops
        idx_success = dataTable.Staircase == 'HomeStair' & ...
            dataTable.StopSignalDelay == currentSSD & dataTable.Outcome == 'Correct_Stop_Home';
        nTrls_correct = sum(idx_success);
        % Calculate pStop at that SSD
        pStop_SSD(s,1 ) = (nTrls_correct / nTrls_all) * 100;
        
        % Index and get the failed stop RTs
        idx_fail = dataTable.Staircase == 'HomeStair' & ...
            dataTable.StopSignalDelay == currentSSD & dataTable.Outcome == 'Failed_Stop_Home' | ...
            dataTable.Staircase == 'HomeStair' & ...
            dataTable.StopSignalDelay == currentSSD & dataTable.Outcome == 'Failed_Stop_Target';
        find_fail = find(dataTable.Staircase == 'HomeStair' & ...
            dataTable.StopSignalDelay == currentSSD & dataTable.Outcome == 'Failed_Stop_Home' | ...
            dataTable.Staircase == 'HomeStair' & ...
            dataTable.StopSignalDelay == currentSSD & dataTable.Outcome == 'Failed_Stop_Target');
        FailedStop_RT_ToT_SSD(s,1) = round(nanmedian(dataTable.ReactTime(idx_fail)));
        
        % Index the n-1 or n-2 go trials with respect to the failed stops
        idx_n_min1_isGo = find(idx_fail(2:end) + trials.GoAll(1:end-1) == 2);
        if ~isempty(idx_n_min1_isGo)
            idx_n_min2_isGo = find(idx_fail(2:end) - trials.GoAll(1:end-1) == 1);
            if isempty(idx_n_min2_isGo)
                idx_n_minx_isGo = idx_n_min1_isGo;
            else
                idx_n_minx_isGo = sort([idx_n_min1_isGo; idx_n_min2_isGo]);
            end
            % Get arrays of FS RTs and the N-1 or N-2 Go RTs
            FSRT_array = dataTable.ReactTime(idx_fail);
            GoRT_N_minx_isGo_S = dataTable.ReactTime(idx_n_minx_isGo);
            % Remove values where Go RT has been NaN'd
            idx_nans = isnan(FSRT_array) | isnan(GoRT_N_minx_isGo_S);
            FSRT_array(idx_nans) = [];
            GoRT_N_minx_isGo_S(idx_nans) = [];
            % Check if FSRT < GoRT
            FSRT_lessthen_GoRT = FSRT_array < GoRT_N_minx_isGo_S;
            % Number of FS trials at this SSD =
            n_FStrial = length(FSRT_lessthen_GoRT);
            n_FSRT_less_GoRT = sum(FSRT_lessthen_GoRT);
            % Calculate % trials where FSRT < GoRT
            pcnt_FSRT_lessthan_GoRT = round(n_FSRT_less_GoRT / n_FStrial * 100);
            % Put value into an array for each SSD
            trlwise_FSRT_GoRT_pcnt(s,1) = n_FStrial;
            trlwise_FSRT_GoRT_pcnt(s,2) = n_FSRT_less_GoRT;
            trlwise_FSRT_GoRT_pcnt(s,3) = pcnt_FSRT_lessthan_GoRT;
        else
            trlwise_FSRT_GoRT_pcnt(s,1) = NaN;
            trlwise_FSRT_GoRT_pcnt(s,2) = NaN;
            trlwise_FSRT_GoRT_pcnt(s,3) = NaN;
        end
        
        % Calculate the relative finishing times normalised with respect to
        % the distribution of go response times (i.e. the std)
        ZRFT_SSD(s,1) = (nanmedian(RT_GoCorr) - currentSSD - SSRT_HomeStair) / ...
            nanstd(RT_GoCorr);
    end
    
    % Target stair case
    for s = 1:length(SSD_uniq_TargetStair)
        currentSSD = SSD_uniq_TargetStair(s);
        %     count how many trials for that SSD in that stair case
        idx_all = dataTable.StopSignalDelay(trials.StopAll_TargetStr) == currentSSD;
        nTrls_all = sum(idx_all);
        % index and count the the successful stops
        idx_success = dataTable.Staircase == 'TargetStair' & ...
            dataTable.StopSignalDelay == currentSSD & dataTable.Outcome == 'Correct_Stop_Home' | ...
            dataTable.Staircase == 'TargetStair' & ...
            dataTable.StopSignalDelay == currentSSD & dataTable.Outcome == 'Failed_Stop_Home';
        nTrls_correct = sum(idx_success);
        % Calculate pStop at that SSD
        pStop_SSD(s,2) = (nTrls_correct / nTrls_all) * 100;
        
        % Index and get the failed stopTTs
        idx_fail = dataTable.Staircase == 'TargetStair' & ...
            dataTable.StopSignalDelay == currentSSD & dataTable.Outcome == 'Failed_Stop_Target';
        FailedStop_RT_ToT_SSD(s,2 ) = round(nanmedian(dataTable.ReactTime(idx_fail) + dataTable.MoveTime(idx_fail)));
        
        % Index the n-1 or n-2 go trials with respect to the failed stops
        idx_n_min1_isGo = find(idx_fail(2:end) + trials.GoAll(1:end-1) == 2);
        if ~isempty(idx_n_min1_isGo)
            idx_n_min2_isGo = find(idx_fail(2:end) - trials.GoAll(1:end-1) == 1);
            if isempty(idx_n_min2_isGo);
                idx_n_minx_isGo = idx_n_min1_isGo;
            else
                idx_n_minx_isGo = sort([idx_n_min1_isGo; idx_n_min2_isGo]);
            end
            % Get arrays of FS RTs and the N-1 or N-2 Go RTs
            FSRT_array = dataTable.ReactTime(idx_fail);
            GoRT_N_minx_isGo_S = dataTable.ReactTime(idx_n_minx_isGo);
            % Remove values where Go RT has been NaN'd
            idx_nans = isnan(FSRT_array) | isnan(GoRT_N_minx_isGo_S);
            FSRT_array(idx_nans) = [];
            GoRT_N_minx_isGo_S(idx_nans) = [];
            % Check if FSRT < GoRT
            FSRT_lessthen_GoRT = FSRT_array < GoRT_N_minx_isGo_S;
            % Number of FS trials at this SSD =
            n_FStrial = length(FSRT_lessthen_GoRT);
            n_FSRT_less_GoRT = sum(FSRT_lessthen_GoRT);
            % Calculate % trials where FSRT < GoRT
            pcnt_FSRT_lessthan_GoRT = round(n_FSRT_less_GoRT / n_FStrial * 100);
            % Put value into an array for each SSD
            trlwise_FSRT_GoRT_pcnt(s,4) = n_FStrial;
            trlwise_FSRT_GoRT_pcnt(s,5) = n_FSRT_less_GoRT;
            trlwise_FSRT_GoRT_pcnt(s,6) = pcnt_FSRT_lessthan_GoRT;
        else
            trlwise_FSRT_GoRT_pcnt(s,4) = NaN;
            trlwise_FSRT_GoRT_pcnt(s,5) = NaN;
            trlwise_FSRT_GoRT_pcnt(s,6) = NaN;
        end
        
        ZRFT_SSD(s,2) = (nanmedian(TotRespT_GoCorr) - currentSSD - SSRT_TargetStair_TT) / ...
            nanstd(TotRespT_GoCorr);
    end
    
    idx_HomeStair = dataTable.Staircase == 'HomeStair' & ~isnan(dataTable.StopSignalDelay);
    idx_TargetStair = dataTable.Staircase == 'TargetStair' & ~isnan(dataTable.StopSignalDelay);
    
    figure()
    % Stopping stair case
    subplot(2,2,1)
    plot(dataTable.StopSignalDelay(idx_HomeStair),  'b');
    hold on
    plot(dataTable.StopSignalDelay(idx_TargetStair),  '--b');
    title('Staircasing','FontSize',14);
    xlabel('Trial','FontSize',14);
    ylabel('Stop signal delay (ms)','FontSize',14);
    % xlim([0,max(SSD_uniq_TargetStair)]);
    ylim([0,max([SSD_uniq_HomeStair;SSD_uniq_TargetStair])]);
    legend('Home','Target','FontSize',14, 'Location', 'East')
    legend boxoff
    % print(gcf,strcat(folderName,'/',subjID,'_PvO_Stop_Staircase.pdf'),'-dpdf','-r400');
    
    %FS RT as function of SSD
    subplot(2,2,2)
    scatter(SSD_uniq_HomeStair, FailedStop_RT_ToT_SSD(1:length(SSD_uniq_HomeStair),1), 50, 'filled', 'b');
    hold on
    scatter(SSD_uniq_TargetStair, FailedStop_RT_ToT_SSD(1:length(SSD_uniq_TargetStair),2), 50, 'b');
    plot([0, max(SSD_uniq_TargetStair)],[nanmedian(RT_GoCorr_mn), nanmedian(RT_GoCorr_mn)],'--k');
    plot([0, max(SSD_uniq_TargetStair)],[nanmedian(TotRespT_GoCorr), nanmedian(TotRespT_GoCorr)],'--k');
    title('Failed Stop Move Time','FontSize',14);
    xlabel('Stop signal delay (ms)','FontSize',14);
    ylabel('Failed Stop response time (ms)','FontSize',14);
    xlim([0,max(SSD_uniq_TargetStair)]);
    ylim([0,nanmedian(TotRespT_GoCorr)+100]);
    legend('Home','Target','RT','TT','FontSize',14, 'location', 'SouthEast')
    legend boxoff
    % print(gcf,strcat(folderName,'/',subjID,'_PvO_Stop_SSD_FSMT.pdf'),'-dpdf','-r400');
    
    SSD_FSRT = NaN(22, 4);
    SSD_FSRT(1:length(SSD_uniq_HomeStair),1) = SSD_uniq_HomeStair;
    SSD_FSRT(1:length(SSD_uniq_HomeStair),2) = FailedStop_RT_ToT_SSD(1:length(SSD_uniq_HomeStair), 1);
    SSD_FSRT(1:length(SSD_uniq_TargetStair),3) = SSD_uniq_TargetStair;
    SSD_FSRT(1:length(SSD_uniq_TargetStair),4) = FailedStop_RT_ToT_SSD(1:length(SSD_uniq_TargetStair),2);
    data.SSD_FSRT = array2table(SSD_FSRT, ...
        'VariableNames',{'SSDh','FSRTh','SSDt','FSTTt'});
    
    % Stop function plot
    SSD_pStop = NaN(22, 4);
    SSD_pStop(1:length(SSD_uniq_HomeStair),1) = SSD_uniq_HomeStair;
    SSD_pStop(1:length(SSD_uniq_HomeStair),2) = pStop_SSD(1:length(SSD_uniq_HomeStair), 1);
    SSD_pStop(1:length(SSD_uniq_TargetStair),3) = SSD_uniq_TargetStair;
    SSD_pStop(1:length(SSD_uniq_TargetStair),4) = pStop_SSD(1:length(SSD_uniq_TargetStair),2);
    data.SSD_pStop = array2table(SSD_pStop, ...
        'VariableNames',{'SSDh','pStoph','SSDt','pStopt'});
    
    % want to plot raw data and the robust fit regression lines
%     idx_nan_home = isnan(data.SSD_pStop.SSDh) |isnan(data.SSD_pStop.pStoph);
%     
%     corr_result_home = fitlm(data.SSD_pStop.SSDh(~idx_nan_home),data.SSD_pStop.pStoph(~idx_nan_home), 'RobustOpts','on');
%     rh = sqrt(corr_result_home.Rsquared.Adjusted);
%     ph = corr_result_home.Coefficients.pValue(2);
%     interch = corr_result_home.Coefficients.Estimate(1);
%     slopeh = corr_result_home.Coefficients.Estimate(2);
%     
%     idx_nan_target = isnan(data.SSD_pStop.SSDt) |isnan(data.SSD_pStop.pStopt);
%     
%     corr_result_target = fitlm(data.SSD_pStop.SSDt(~idx_nan_target),data.SSD_pStop.pStopt(~idx_nan_target), 'RobustOpts','on');
%     rt = sqrt(corr_result_target.Rsquared.Adjusted);
%     pt = corr_result_target.Coefficients.pValue(2);
%     interct = corr_result_target.Coefficients.Estimate(1);
%     slopet = corr_result_target.Coefficients.Estimate(2);
%     
%     %     y = a + bX (where a is intercept, i.e. value of y when x = 0, and b is the slope)
%     minSSDh =  round((100 - interch )/ slopeh); % i.e. value of x when y is 100
%     maxSSDh =  round((0 - interch )/ slopeh);% i.e. value of x when y is 0
%     meanSSDh = round(mean([maxSSDh, minSSDh])); % Should be value of x when y is 50, i.e. mean SSD
%     rangeSSDh = maxSSDh - minSSDh;
%     
%     %     y = a + bX (where a is intercept, i.e. value of y when x = 0, and b is the slope)
%     minSSDt =  round((100 - interct )/ slopet); % i.e. value of x wten y is 100
%     maxSSDt =  round((0 - interct )/ slopet);% i.e. value of x when y is 0
%     meanSSDt = round(mean([maxSSDt, minSSDt])); % Should be value of x when y is 50, i.e. mean SSD
%     rangeSSDt = maxSSDt - minSSDt;
    
%     data.StopFunction_linReg = array2table([rh, ph, interch, slopeh, minSSDh, maxSSDh,meanSSDh, rangeSSDh,...
%         rt, pt, interct,slopet,minSSDt,maxSSDt,meanSSDt,rangeSSDt],...
%         'VariableNames',{'rh','ph','intercepth','slopeh','minSSDh','maxSSDh','meanSSDh','rangeSSDh',...
%         'rt','pt','interceptt','slopet','minSSDt','maxSSDt','meanSSDt','rangeSSDt'});
%     
    subplot(2,2,3)
    scatter(SSD_uniq_HomeStair, pStop_SSD(1:length(SSD_uniq_HomeStair),1), 50, 'filled', 'b');
    hold on
    scatter(SSD_uniq_TargetStair, pStop_SSD(1:length(SSD_uniq_TargetStair),2), 50, 'b');
    plot([0, max(SSD_uniq_TargetStair)],[50,50],'--k');
%     plot([minSSDh,maxSSDh],[100,0],'k--x')
%     plot([minSSDt,maxSSDt],[100,0],'k--x')
    title('Stop functions','FontSize',14);
    xlabel('Stop signal delay (ms)','FontSize',14);
    ylabel('Prob. stopping (%)','FontSize',14);
    xlim([0,max(SSD_uniq_TargetStair)]);
    ylim([0,100]);
    legend('Home','Target','robFit_h','robFit_t','FontSize',14, 'Location', 'South')
    legend boxoff
    % print(gcf,strcat(folderName,'/',subjID,'_PvO_Stop_StopFunction.pdf'),'-dpdf','-r400');
    
    
    % ZRFT
    subplot(2,2,4)
    scatter(ZRFT_SSD(1:length(SSD_uniq_HomeStair),1), pStop_SSD(1:length(SSD_uniq_HomeStair),1), 50, 'filled', 'b');
    hold on
    scatter(ZRFT_SSD(1:length(SSD_uniq_TargetStair),2), pStop_SSD(1:length(SSD_uniq_TargetStair),2), 50, 'b');
    plot([0, max(SSD_uniq_TargetStair)],[50,50],'--k');
    lsline;
    title('Stop functions','FontSize',14);
    xlabel('ZRFT','FontSize',14); % NB positive/negative ZRFT = short/long SSD respectively
    ylabel('Norm. Prob. stopping (%)','FontSize',14);
    xlim([- 2, 2]);
    ylim([0,100]);
    legend('Home','Target','FontSize',14, 'Location', 'NorthWest')
    legend boxoff
    % print(gcf,strcat(folderName,'/',subjID,'_PvO_Stop_NormStopFunction.pdf'),'-dpdf','-r400');
    
    ZRFT_pStop = NaN(22, 4);
    ZRFT_pStop(1:length(SSD_uniq_HomeStair),1) = ZRFT_SSD(1:length(SSD_uniq_HomeStair),1);
    ZRFT_pStop(1:length(SSD_uniq_HomeStair),2) = pStop_SSD(1:length(SSD_uniq_HomeStair), 1);
    ZRFT_pStop(1:length(SSD_uniq_TargetStair),3) = ZRFT_SSD(1:length(SSD_uniq_TargetStair),1);
    ZRFT_pStop(1:length(SSD_uniq_TargetStair),4) = pStop_SSD(1:length(SSD_uniq_TargetStair),2);
    data.ZRFT_pStop = array2table(ZRFT_pStop, ...
        'VariableNames',{'SSDh','pStoph','SSDt','pStopt'});
    
    print(gcf,strcat(folderName,'/',subjID,'_PvO_Stop_StopPerf.pdf'),'-dpdf','-r400');
    
    figure()
    scatter(SSD_uniq_HomeStair(1:length(SSD_uniq_HomeStair),1), trlwise_FSRT_GoRT_pcnt(1:length(SSD_uniq_HomeStair),3), 50, 'filled', 'b');
    hold on
    scatter(SSD_uniq_TargetStair(1:length(SSD_uniq_TargetStair),1), trlwise_FSRT_GoRT_pcnt(1:length(SSD_uniq_TargetStair),6), 50, 'b');
    title('Trlwise violation: context independence','FontSize',14);
    xlabel('SSD (ms)','FontSize',14); % NB positive/negative ZRFT = short/long SSD respectively
    ylabel('FSRT<GoRT (%)','FontSize',14);
    ylim([0,100]);
    xlim([0, max(SSD_uniq_TargetStair)+50]);
    legend('Home','Target','FontSize',14, 'Location', 'NorthWest')
    legend boxoff
    
    Trlwise_CntxInd_pcnt = NaN(22, 8);
    Trlwise_CntxInd_pcnt(1:length(SSD_uniq_HomeStair),1) = SSD_uniq_HomeStair;
    Trlwise_CntxInd_pcnt(1:length(SSD_uniq_HomeStair),2) = trlwise_FSRT_GoRT_pcnt(1:length(SSD_uniq_HomeStair), 1);
    Trlwise_CntxInd_pcnt(1:length(SSD_uniq_HomeStair),3) = trlwise_FSRT_GoRT_pcnt(1:length(SSD_uniq_HomeStair), 2);
    Trlwise_CntxInd_pcnt(1:length(SSD_uniq_HomeStair),4) = trlwise_FSRT_GoRT_pcnt(1:length(SSD_uniq_HomeStair), 3);
    Trlwise_CntxInd_pcnt(1:length(SSD_uniq_TargetStair),5) = SSD_uniq_TargetStair;
    Trlwise_CntxInd_pcnt(1:length(SSD_uniq_TargetStair),6) = trlwise_FSRT_GoRT_pcnt(1:length(SSD_uniq_TargetStair),4);
    Trlwise_CntxInd_pcnt(1:length(SSD_uniq_TargetStair),7) = trlwise_FSRT_GoRT_pcnt(1:length(SSD_uniq_TargetStair),5);
    Trlwise_CntxInd_pcnt(1:length(SSD_uniq_TargetStair),8) = trlwise_FSRT_GoRT_pcnt(1:length(SSD_uniq_TargetStair),6);
    data.Trlwise_CntxInd_pcnt = array2table(Trlwise_CntxInd_pcnt, ...
        'VariableNames',{'SSDh','nFStrl_h','nFSlessGo_h','pcntFSlessGo_h',...
        'SSDt','nFStrl_t','nFSlessGo_t','pcntFSlessGo_t'});
    
    %% Plot RT MT TotT SSRT and mean SSD
    figure()
    subplot(2,3,1)
    bar([RT_GoCorr_mn,0, 0,0,0], 'k');
    hold on
    bar([0, 0, 0,0,0],'b');
    bar([0,0, RT_StopCorr_TargetStair_mn,0,0], 'b','LineStyle','--');
    bar([0,0,0,RT_StopFail_HomeStair_mn, 0],'r');
    bar([0,0,0,0 RT_StopFail_TargetStair_mn], 'r','LineStyle','--');
    er = errorbar([1,2,3,4,5],[RT_GoCorr_mn,0, RT_StopCorr_TargetStair_mn, RT_StopFail_HomeStair_mn,RT_StopFail_TargetStair_mn],...
        [nanstd(RT_GoCorr),0, nanstd(RT_StopCorr_TargetStair),nanstd(RT_StopFail_HomeStair), nanstd(RT_StopFail_TargetStair)],'k');
    er.LineStyle = 'none';
    title('React. Time','FontSize',14);
    xlabel('Staircase','FontSize',14);
    xticks([1 2 3 4 5])
    xticklabels({'Go','SSH','SST','FSH','FST'})
    ylabel('RT (ms)','FontSize',14);
    ylim([200,round(RT_GoCorr_mn + 200)]);
    
    subplot(2,3,2)
    bar([MT_GoCorr_mn,0, 0,0,0], 'k');
    hold on
    bar([0,0, 0,0,0],'b');
    bar([0,0, 0,0,0],'b' ,'LineStyle','--');
    bar([0,0,0,MT_StopFail_HomeStair_mn,0],'r');
    bar([0,0,0,0, MT_StopFail_TargetStair_mn], 'r','LineStyle','--');
    er = errorbar([1,2,3,4,5],[MT_GoCorr_mn, 0, 0,MT_StopFail_HomeStair_mn,MT_StopFail_TargetStair_mn],...
        [nanstd(MT_GoCorr), 0, 0, nanstd(MT_StopFail_HomeStair), nanstd(MT_StopFail_HomeStair)],'k');
    er.LineStyle = 'none';
    title('Move. Time','FontSize',14);
    xlabel('Staircase','FontSize',14);
    xticks([1 2 3 4 5])
    xticklabels({'Go','SSH','SST','FSH','FST'})
    ylabel('MT (ms)','FontSize',14);
    ylim([200,round(MT_GoCorr_mn + 200)]);
    
    
    subplot(2,3,3)
    bar([nanmedian(TotRespT_GoCorr),0, 0], 'k');
    hold on
    bar([0,0, 0,0,0],'b');
    bar([0,0, 0,0,0], 'b','LineStyle','--');
    bar([0,0,0, 0, 0],'r');
    bar([0,0,0,0,  nanmedian(TotRespT_StopFail_TargetStair)], 'r','LineStyle','--');
    er = errorbar([1,2,3, 4, 5],[nanmedian(TotRespT_GoCorr),0,0,0,nanmedian(TotRespT_StopFail_TargetStair)],...
        [nanstd(TotRespT_GoCorr),0, 0,0, nanstd(TotRespT_StopFail_TargetStair)],'k');
    er.LineStyle = 'none';
    title('Total Time','FontSize',14);
    xlabel('Staircase','FontSize',14);
    xticks([1 2 3 4 5])
    xticklabels({'Go','SSH','SST','FSH','FST'})
    ylabel('Total Time (ms)','FontSize',14);
    ylim([500,round(nanmedian(TotRespT_GoCorr + 200))]);
    
    subplot(2,3,4)
    bar([pStop_HomeStair,0], 'b');
    hold on
    bar([0,pStop_TargetStair], 'b','LineStyle','--');
    title('pStop','FontSize',14);
    xlabel('Staircase','FontSize',14);
    xticks([1 2])
    xticklabels({'Home','Target'})
    ylabel('pStop (%)','FontSize',14);
    ylim([0,100]);
    
    subplot(2,3,5)
    bar([SSD_Stop_HomeStair_mn,0], 'b');
    hold on
    bar([0,SSD_Stop_TargetStair_mn], 'b','LineStyle','--');
    title('SSD','FontSize',14);
    xlabel('Staircase','FontSize',14);
    xticks([1 2])
    xticklabels({'Home','Target'})
    ylabel('SSD (ms)','FontSize',14);
    ylim([0,round(max([SSD_Stop_HomeStair_mn;SSD_Stop_TargetStair_mn]) + 50)]);
    
    subplot(2,3,6)
    bar([SSRT_HomeStair,0, 0], 'b');
    hold on
    bar([0,SSRT_TargetStair_MT,0], 'b','LineStyle','--');
    bar([0,0,SSRT_TargetStair_TT], 'b','LineStyle','--');
    title('SSRT','FontSize',14);
    xlabel('Staircase','FontSize',14);
    xticks([1 2 3])
    xticklabels({'hRT','tMT','tTT'})
    ylabel('SSRT (ms)','FontSize',14);
    ylim([0,round(max([SSRT_HomeStair,SSRT_TargetStair_MT,SSRT_TargetStair_TT]) + 75)]);
    lsline
    
    print(gcf,strcat(folderName,'/',subjID,'_PvO_Stop_MeanBehav.pdf'),'-dpdf','-r400');
    
    
    
    %% Analyse mouse movements
    % Qu. Successful stop - velocity change useful measure of stopping latency?
    % Qu. Also want to simply visualise movements
    
    % Want to index the time point/column number when the mouse entered the box
    % so that we're not looking at all the junk after.
    % Need a way that captures all trials and the different maximums and times
    % that they move for.
    
    % For go trials, can use the min Y value, but this won't work for
    % Target/reversal trials since it will cut it off
    
    mouseXYT(trials.Miss,:,1) = NaN; % NaN miss trials (applies to Go only)
    mouseXYT(trials.Miss,:,2) = NaN;
    mouseXYT(trials.Miss,:,3) = NaN;
    
    mouseXYT(trials.RT_outlier,:,1) = NaN;
    mouseXYT(trials.RT_outlier,:,2) = NaN;
    mouseXYT(trials.RT_outlier,:,3) = NaN;
    
    mouseXYT(trials.MT_outlier,:,1) = NaN;
    mouseXYT(trials.MT_outlier,:,2) = NaN;
    mouseXYT(trials.MT_outlier,:,3) = NaN;
    
    data.mouseX = mouseXYT(:,:,1);
    data.mouseY = mouseXYT(:,:,2);
    data.mouseT = mouseXYT(:,:,3);
    
    % Get reaction times, and NaN values less than RT. This is because
    % sometimes there are small mouse movements that do not leave the box and
    % thus generate an RT initially. This messes up correspondence between
    % data.mouseT and RTs.
    col_underRT = data.mouseT < round(dataTable.ReactTime);
    data.mouseX(col_underRT)= NaN;
    data.mouseY(col_underRT)= NaN;
    data.mouseT(col_underRT)= NaN;
    
    % Get total movement time and NaN anything that is above that
    % This has no effect on trials without a movement time
    % since they are NaNs in TotRespT_T1_all so nothing to compare to.
    %This therefore leaves successful stops unaffected, which is good.
    TotRespT_T1_all = round(dataTable.ReactTime + dataTable.MoveTime);
    col_overTRT = data.mouseT > TotRespT_T1_all + round(frameRate);
    data.mouseX(col_overTRT)= NaN;
    data.mouseY(col_overTRT)= NaN;
    data.mouseT(col_overTRT)= NaN;
    
    % New time variable so that RT is now 0 ms, so can just look at movement
    % times
    data.mouseT_alignRT = data.mouseT;
    min_mouseT = min(data.mouseT_alignRT');
    max_mouseT = max(data.mouseT_alignRT') - min_mouseT;
    data.mouseT_alignRT = bsxfun(@minus, data.mouseT_alignRT, min_mouseT');
    
    % Shift X, Y and T so first value is in first column
    idx_mouseX_NaN = isnan(data.mouseX);
    
    for j = 1:size(data.mouseX,1)
        data.mouseX(j,:)  = [data.mouseX(j,~idx_mouseX_NaN(j,:)) data.mouseX(j,idx_mouseX_NaN(j,:))];
        data.mouseY(j,:)  = [data.mouseY(j,~idx_mouseX_NaN(j,:)) data.mouseY(j,idx_mouseX_NaN(j,:))];
        data.mouseT_alignRT(j,:)  = [data.mouseT_alignRT(j,~idx_mouseX_NaN(j,:)) data.mouseT_alignRT(j,idx_mouseX_NaN(j,:))];
    end
    
    %  Remember actually that movements were from high to low in terms of pixels, so
    % have corrected Y by making negative and then adding the minimum Y to all
    % values
    data.mouseY = -data.mouseY; % make positive
    data.mouseY = data.mouseY - data.mouseY(:,1);
    
    % Find the location in the time array at which time ends
    [TmaxVal,TmaxLoc] = max(data.mouseT_alignRT');
    
    % Find the location in the Y array at which Y is at it's maximum (now
    % Y-change data has been made positive (see above)
    [data.YmaxVal,data.YmaxLoc] = max(data.mouseY');
    
    figure()
    for m = 1:size(data.mouseX,1)
        if dataTable.Outcome(m) == 'Correct_Go'
            final_datpoint = TmaxLoc(m);
            color = 'g';
        elseif dataTable.Outcome(m) == 'Correct_Stop_Home' & dataTable.Staircase(m) == 'HomeStair'
            %         final_datpoint = TmaxLoc(m);
            final_datpoint = data.YmaxLoc(m);
            color = 'b';
        elseif dataTable.Outcome(m) == 'Failed_Stop_Home' && dataTable.Staircase(m) == 'HomeStair'
            %         final_datpoint = size(data.mouseX,2);
            final_datpoint = data.YmaxLoc(m);
            color = 'r';
        elseif dataTable.Outcome(m) == 'Failed_Stop_Home' && dataTable.Staircase(m) == 'TargetStair'
            %         final_datpoint = TmaxLoc(m);
            final_datpoint = TmaxLoc(m);
            color = '-.b';
        elseif dataTable.Outcome(m) == 'Failed_Stop_Target' && dataTable.Staircase(m) == 'TargetStair'
            %         final_datpoint = size(data.mouseX,2);
            final_datpoint = TmaxLoc(m);
            color = '-.r';
        elseif dataTable.Outcome(m) == 'Miss'
            final_datpoint = TmaxLoc(m);
            color = '.g';
        end
        
        plot(data.mouseX(m,1:final_datpoint), data.mouseY(m,1:final_datpoint),color)
        
        hold on
    end
    
    xlabel('Horizontal (px)','FontSize',14);
    ylabel('Vertical (px)','FontSize',14);
    title('Mouse movement trajectories');
    
    print(gcf,strcat(folderName,'/',subjID,'_PvO_Stop_Trajectory2dXY.pdf'),'-dpdf','-r400');
    
    
    %% Test whether some "failed" home stops really include trigger failures
    % i.e. a TF followed by a error-triggered stop, which should be slow
    % and perhaps also associated with late/slow rising RT, since without
    % that it might become a full-blown TF and hit the target
    
    % Plot Y positions relative to the time of movement onset for Go, FSh
    % and SSt
    figure()
    data.Ypos_GoCorr = data.mouseY(trials.GoCorr,:);
    T_GoCorr = data.mouseT_alignRT(trials.GoCorr,:);
    for trl = 1:size(data.Ypos_GoCorr,1)
        plot(T_GoCorr(trl,:), data.Ypos_GoCorr(trl,:), '-g');
        hold on
    end
    
    data.Ypos_StopFail_Home = data.mouseY(trials.StopFail_HomeStr,:);
    T_StopFail_Home = data.mouseT_alignRT(trials.StopFail_HomeStr,:);
    for trl = 1:size(data.Ypos_StopFail_Home,1)
        plot(T_StopFail_Home(trl,:),data.Ypos_StopFail_Home(trl,:),'-r');
    end
    
    data.Ypos_StopCorr_Targ = data.mouseY(trials.StopCorr_TargetStr,:);
    T_StopCorr_Targ = data.mouseT_alignRT(trials.StopCorr_TargetStr,:);
    for trl = 1:size(data.Ypos_StopCorr_Targ,1)
        plot(T_StopCorr_Targ(trl,:),data.Ypos_StopCorr_Targ(trl,:),'--b');
    end
    xlabel('T.move onset (ms)','FontSize',14);
    ylabel('Y pos (px)','FontSize',14);
    legend('Go (green)','FSh (red)','SSt (blue)','location','NorthEast');
    title('Y pos vs. time from movement onset','FontSize',14);
    print(gcf,strcat(folderName,'/',subjID,'_PvO_TriggerFailure_FailStopH_Ypos_wrtMoveOn.pdf'),'-dpdf','-r400');
    
    %---------------------------------------------------------------------
    % Plot figure of Y position relative to time from Stop Signal
    data.mouseT_SSDnorm_Home = data.mouseT_alignRT;
    
    % Align Go to mean SSD of target trial, all home to home SSD and all target to target SSD
    data.mouseT_SSDnorm_Home(trials.GoCorr,:) = data.mouseT_SSDnorm_Home(trials.GoCorr,:) - ...
        nanmedian(dataTable.StopSignalDelay(trials.StopFail_HomeStr)) + RT_GoCorr_mn;
    
    data.mouseT_SSDnorm_Home(trials.StopFail_HomeStr, :) = data.mouseT_SSDnorm_Home(trials.StopFail_HomeStr, :) - ...
        dataTable.StopSignalDelay(trials.StopFail_HomeStr) +  dataTable.ReactTime(trials.StopFail_HomeStr) ;
    
    data.mouseT_SSDnorm_Home(trials.StopCorr_TargetStr, :) = data.mouseT_SSDnorm_Home(trials.StopCorr_TargetStr, :) - ...
        dataTable.StopSignalDelay(trials.StopCorr_TargetStr) +  dataTable.ReactTime(trials.StopCorr_TargetStr) ;
    
    % Failed stops, home staircase
    % Separate TFs from regular Failed stops
    dummy_table = dataTable(trials.StopFail_HomeStr,:);
    find_TF =  find(dummy_table.MoveTime > 0);
    
    data.Ypos_GoCorr = data.mouseY(trials.GoCorr,:);
    data.T_GoCorr_SSDnorm_Home = data.mouseT_SSDnorm_Home(trials.GoCorr,:);
    figure()
    subplot(2,1,1)
    % plot go relative to home SSD
    for trl = 1:size(data.Ypos_GoCorr,1)
        plot(data.T_GoCorr_SSDnorm_Home(trl,:), data.Ypos_GoCorr(trl,:), '-g');
        hold on
    end
    % plot FSh relative to home SSD
    data.Ypos_StopFail_Home = data.mouseY(trials.StopFail_HomeStr,:);
    data.T_StopFail_Home_SSDnorm = data.mouseT_SSDnorm_Home(trials.StopFail_HomeStr,:);
    for trl = 1:size(data.Ypos_StopFail_Home,1)
        plot(data.T_StopFail_Home_SSDnorm(trl,:),data.Ypos_StopFail_Home(trl,:),'-r');
    end
    if ~isempty(find_TF)
        % plot go relative to home SSD
        for trl = 1:length(find_TF)
            
            plot(data.T_StopFail_Home_SSDnorm(find_TF(trl),:),data.Ypos_StopFail_Home(find_TF(trl),:),'-k');
        end
    end
    
    %    [val,loc] = max(data.Ypos_StopFail_Home);
    xlabel('T.stop (ms)','FontSize',14);
    ylabel('Y pos (px)','FontSize',14);
    xlim([-300, 1200]);
    legend('...');
    text(-200,550,'Go (green)');
    text(-200,500,'FSh (red)');
    text(-200,450,'TF (black)');
    title('Peak Y pos as function of RT wrt stop','FontSize',14);
    
    %---------------------------------------------------------------------
    % Re-align data again so now relative to Target SSD
    data.mouseT_SSDnorm_Target = data.mouseT_alignRT;
    
    data.mouseT_SSDnorm_Target(trials.GoCorr,:) = data.mouseT_SSDnorm_Target(trials.GoCorr,:) - ...
        nanmedian(dataTable.StopSignalDelay(trials.StopCorr_TargetStr)) + RT_GoCorr_mn;
    T_GoCorr_SSDnorm_Target = data.mouseT_SSDnorm_Target(trials.GoCorr,:);
    
    data.mouseT_SSDnorm_Target(trials.StopCorr_TargetStr, :) = data.mouseT_SSDnorm_Target(trials.StopCorr_TargetStr, :) - ...
        dataTable.StopSignalDelay(trials.StopCorr_TargetStr) +  dataTable.ReactTime(trials.StopCorr_TargetStr) ;
    
    subplot(2,1,2)
    for trl = 1:size(data.Ypos_GoCorr,1)
        plot(T_GoCorr_SSDnorm_Target(trl,:), data.Ypos_GoCorr(trl,:), '-g');
        hold on
    end
    
    data.Ypos_StopCorr_Targ = data.mouseY(trials.StopCorr_TargetStr,:);
    data.T_StopCorr_Targ_SSDnorm = data.mouseT_SSDnorm_Target(trials.StopCorr_TargetStr,:);
    for trl = 1:size(data.Ypos_StopCorr_Targ,1)
        plot(data.T_StopCorr_Targ_SSDnorm(trl,:),data.Ypos_StopCorr_Targ(trl,:),'--b');
    end
    xlabel('T.stop (ms)','FontSize',14);
    ylabel('Y pos (px)','FontSize',14);
    xlim([-300, 1200]);
    legend('...');
    text(-200,550,'Go (green)');
    text(-200,500,'SSt (blue)');
    title('Peak Y pos as function of RT wrt stop','FontSize',14);
    
    print(gcf,strcat(folderName,'/',subjID,'_PvO_TriggerFailure_FailStopH_Ypos_wrtStop.pdf'),'-dpdf','-r400');
    
    %---------------------------------------------------------------------
    % Max Y failed stop home as % of average max Y on Go Correct trials
    maxYpos_StopFail_Home = round((data.YmaxVal(trials.StopFail_HomeStr) ./ str2double(data.Computer_setup.home2targetYdist)) * 100);
    RangeYpos_StopFail_Home(1,1) = round(nanmin(maxYpos_StopFail_Home));
    RangeYpos_StopFail_Home(1,2) = round(nanmax(maxYpos_StopFail_Home));
    RangeYpos_StopFail_Home(1,3) = round(nanmedian(maxYpos_StopFail_Home));
    RangeYpos_StopFail_Home(1,4) = round(nanstd(maxYpos_StopFail_Home));
    
    % Max Y correct stop target as % of average max Y on Go Correct trials
    maxYpos_StopCorr_Target = round((data.YmaxVal(trials.StopCorr_TargetStr) ./ str2double(data.Computer_setup.home2targetYdist)) * 100);
    RangeYpos_StopCorr_Target(1,1) = round(nanmin(maxYpos_StopCorr_Target));
    RangeYpos_StopCorr_Target(1,2) = round(nanmax(maxYpos_StopCorr_Target));
    RangeYpos_StopCorr_Target(1,3) = round(nanmedian(maxYpos_StopCorr_Target));
    RangeYpos_StopCorr_Target(1,4) = round(nanstd(maxYpos_StopCorr_Target));
    
    data.maxYpos_pcntTargDist = array2table([ RangeYpos_StopFail_Home(1,3), RangeYpos_StopCorr_Target(1,3)],...
        'VariableNames',{'Ymaxpos_h','Ymaxpos_t'});
    
    %---------------------------------------------------------------------
    % Get time of peak relative to the time of the stop signal for Home
    % Staircase
    % Go Home
    [valGho,locGoh] = max(data.Ypos_GoCorr,[],2); % find max Y on go trials
    row = (1:1:size(data.Ypos_GoCorr,1))';
    col = locGoh;
    sz = [size(data.Ypos_GoCorr,1) size(data.Ypos_GoCorr,2)];
    indGo = sub2ind(sz,row,col); % function used for linear indexing
    time_maxY_GoCorr_Home = data.T_GoCorr_SSDnorm_Home(indGo);
    
    prctile_time_maxY_GoCorr_Home        = prctile(time_maxY_GoCorr_Home,[25 75]);
    idx_outliers = time_maxY_GoCorr_Home < (prctile_time_maxY_GoCorr_Home(1) - 1.5 *...
        iqr(time_maxY_GoCorr_Home)) | ...
        time_maxY_GoCorr_Home > (prctile_time_maxY_GoCorr_Home(2) + 1.5 * iqr(time_maxY_GoCorr_Home));
%     time_maxY_GoCorr_Home(idx_outliers) = NaN;
    
    time_maxY_GoCorr_Home_mn = round(nanmedian(time_maxY_GoCorr_Home(~idx_outliers)));
    time_maxY_GoCorr_Home_std = round(nanstd(time_maxY_GoCorr_Home(~idx_outliers)));
    
    % Failed stop home
    [valFSh,locFSh] = max(data.Ypos_StopFail_Home,[],2); % find max Y on go trials
    row = (1:1:size(data.Ypos_StopFail_Home,1))';
    col = locFSh;
    sz = [size(data.Ypos_StopFail_Home,1) size(data.Ypos_StopFail_Home,2)];
    indFS = sub2ind(sz,row,col); % function used for linear indexing
    time_maxY_StopFail_Home = data.T_StopFail_Home_SSDnorm(indFS);
    
    prctile_time_maxY_StopFail_Home        = prctile(time_maxY_StopFail_Home,[25 75]);
    idx_outliers = time_maxY_StopFail_Home < (prctile_time_maxY_StopFail_Home(1) - 1.5 *...
        iqr(time_maxY_StopFail_Home)) | ...
        time_maxY_StopFail_Home > (prctile_time_maxY_StopFail_Home(2) + 1.5 * iqr(time_maxY_StopFail_Home));
%     time_maxY_StopFail_Home(idx_outliers) = NaN;
    
    time_maxY_StopFail_Home_mn = round(nanmedian(time_maxY_StopFail_Home(~idx_outliers)));
    time_maxY_StopFail_Home_std = round(nanstd(time_maxY_StopFail_Home(~idx_outliers)));
    time_maxY_StopFail_Home_min = nanmin(time_maxY_StopFail_Home(~idx_outliers));
    time_maxY_StopFail_Home_max = nanmax(time_maxY_StopFail_Home(~idx_outliers));
    if ~isempty(find_TF)
        time_maxY_StopFail_TF_Home_mn = round(nanmedian(time_maxY_StopFail_Home(find_TF)));
    else
        time_maxY_StopFail_TF_Home_mn = NaN;
    end
    
    % How many outliers and what was the time of Y max for them?
    n_outliers_time_maxY_StopFail_Home = sum(idx_outliers);
    outliers_time_maxY_StopFail_Home_mn = nanmedian(time_maxY_StopFail_Home(idx_outliers,1));
    
    % Onset times for FSh and FSh outliers in terms of Y max time
    OnsetTime_StopFail_Home_mn = round(nanmedian(data.T_StopFail_Home_SSDnorm(:,1)));
    outliers_OnsetTime_StopFail_Home_mn = round(nanmedian(data.T_StopFail_Home_SSDnorm(idx_outliers,1)));
    
    %---------------------------------------------------------------------
    % Get time of peak relative to the time of the stop signal for Target
    % Staircase
    % Go target
    [valGot,locGot] = max(data.Ypos_GoCorr,[],2); % find max Y on go trials
    row = (1:1:size(data.Ypos_GoCorr,1))';
    col = locGot;
    sz = [size(data.Ypos_GoCorr,1) size(data.Ypos_GoCorr,2)];
    indGo = sub2ind(sz,row,col); % function used for linear indexing
    time_maxY_GoCorr_Targ = T_GoCorr_SSDnorm_Target(indGo);
    
    prctile_time_maxY_GoCorr_Targ        = prctile(time_maxY_GoCorr_Targ,[25 75]);
    idx_outliers = time_maxY_GoCorr_Targ < (prctile_time_maxY_GoCorr_Targ(1) - 1.5 *...
        iqr(time_maxY_GoCorr_Targ)) | ...
        time_maxY_GoCorr_Targ > (prctile_time_maxY_GoCorr_Targ(2) + 1.5 * iqr(time_maxY_GoCorr_Targ));
%     time_maxY_GoCorr_Targ(idx_outliers) = NaN;
    
    time_maxY_GoCorr_Targ_mn = round(nanmedian(time_maxY_GoCorr_Targ(~idx_outliers)));
    time_maxY_GoCorr_Targ_std = round(nanstd(time_maxY_GoCorr_Targ(~idx_outliers)));
    
    
    % Successful stop target.......///
    [valSSt,locSSt] = max(data.Ypos_StopCorr_Targ,[],2); % find max Y on go trials
    row = (1:1:size(data.Ypos_StopCorr_Targ,1))';
    col = locSSt;
    sz = [size(data.Ypos_StopCorr_Targ,1) size(data.Ypos_StopCorr_Targ,2)];
    indSS = sub2ind(sz,row,col); % function used for linear indexing
    time_maxY_StopCorr_Targ = data.T_StopCorr_Targ_SSDnorm(indSS);
    
    prctile_time_maxY_StopCorr_Targ        = prctile(time_maxY_StopCorr_Targ,[25 75]);
    idx_outliers = time_maxY_StopCorr_Targ < (prctile_time_maxY_StopCorr_Targ(1) - 1.5 *...
        iqr(time_maxY_StopCorr_Targ)) | ...
        time_maxY_StopCorr_Targ > (prctile_time_maxY_StopCorr_Targ(2) + 1.5 * iqr(time_maxY_StopCorr_Targ));
%     time_maxY_StopCorr_Targ(idx_outliers) = NaN;
    
    time_maxY_StopCorr_Targ_mn = round(nanmedian(time_maxY_StopCorr_Targ(~idx_outliers)));
    time_maxY_StopCorr_Targ_std = round(nanstd(time_maxY_StopCorr_Targ(~idx_outliers)));
    time_maxY_StopCorr_Targ_min = nanmin(time_maxY_StopCorr_Targ(~idx_outliers));
    time_maxY_StopCorr_Targ_max = nanmax(time_maxY_StopCorr_Targ(~idx_outliers));
    
    %---------------------------------------------------------------------
    % Plot distribtions of Y max times with respect to the Stop signal
    figure()
    subplot(3,1,1)
    histogram(time_maxY_GoCorr_Home,'Normalization','probability','BinWidth',25, 'FaceColor','g','BinWidth',50); % NB this is using mean SSD
    hold on;
    histogram(time_maxY_StopFail_Home,'Normalization','probability','BinWidth',25, 'FaceColor','r','BinWidth',50);
    histogram(time_maxY_StopFail_Home(find_TF),'Normalization','probability','BinWidth',25, 'FaceColor','k','BinWidth',50);
    xlabel('Ymax times wrt stop (ms)');
    ylabel('Probability');
    xlim([0,1200]);
    legend('Go','FSh','TF-FSh','location','NorthWest');
    title('Home Stair: Go-FSh');
    
    subplot(3,1,2)
    histogram(time_maxY_GoCorr_Targ, 'Normalization','probability','BinWidth',25, 'FaceColor','g','BinWidth',50); % NB this is using mean SSD
    hold on;
    histogram(time_maxY_StopCorr_Targ, 'Normalization','probability','BinWidth',25, 'FaceColor','b','BinWidth',50);
    xlabel('Ymax times wrt stop (ms)');
    ylabel('Probability');
    xlim([0,1200]);
    legend('Go','SSt','location','NorthWest');
    title('Target Stair: Go-SSt');
    
    subplot(3,1,3)
    histogram(time_maxY_StopFail_Home,'Normalization','probability','BinWidth',25,  'FaceColor','r','BinWidth',50); % NB this is using mean SSD
    hold on;
    histogram(time_maxY_StopCorr_Targ, 'Normalization','probability','BinWidth',25, 'FaceColor','b','BinWidth',50);
    histogram(time_maxY_StopFail_Home(find_TF),'Normalization','probability','BinWidth',25, 'FaceColor','k','BinWidth',50);
    xlabel('Ymax times wrt stop (ms)');
    ylabel('Probability');
    xlim([0,1200]);
    legend('Go','FSh','TF-FSh','location','NorthWest');
    title('Home vs Target Stair: FSh-SSt');
    print(gcf,strcat(folderName,'/',subjID,'_PvO_TriggerFailure_Ymaxtime_Hist.pdf'),'-dpdf','-r400');
    
    % Put data in table
    data.time_maxY = array2table([time_maxY_GoCorr_Home_mn, time_maxY_GoCorr_Home_std,...
        time_maxY_StopFail_Home_mn, time_maxY_StopFail_Home_std, time_maxY_StopFail_Home_min, time_maxY_StopFail_Home_max,...
        time_maxY_StopFail_TF_Home_mn,...
        time_maxY_GoCorr_Targ_mn, time_maxY_GoCorr_Targ_std, ...
        time_maxY_StopCorr_Targ_mn,time_maxY_StopCorr_Targ_std , time_maxY_StopCorr_Targ_min, time_maxY_StopCorr_Targ_max,...
        n_outliers_time_maxY_StopFail_Home, outliers_time_maxY_StopFail_Home_mn,...
        OnsetTime_StopFail_Home_mn, outliers_OnsetTime_StopFail_Home_mn],...
        'VariableNames', {'TimeYmaxGoh_mn','TimeYmaxGoh_std',...
        'TimeYmaxFSh_mn','TimeYmaxFSh_std','TimeYmaxFSh_Min','TimeYmaxFSh_Max',...
        'TimeYmaxTFFSh_mn',...
        'TimeYmaxGot_mn','TimeYmaxGot_std',...
        'TimeYmaxSSt_mn','TimeYmaxSSt_std','TimeYmaxSSt_Min','TimeYmaxSSt_Max',...
        'Nh_outliers','TimeYmaxFSh_Out','TimeonFSh_mn','TimeonFSh_Outlier'});
    
    
    %% What happens in trigger failures? Are they fast go, and therefore
    % stop failure could still be the result of relatively slower stop,
    % rather than a failure to trigger?
    if sum(trials.Trigger_failure) > 0
        RT_TF = dataTable.ReactTime(trials.Trigger_failure);
        MT_TF = dataTable.MoveTime(trials.Trigger_failure);
        TotRespT_TF = RT_TF + MT_TF;
        SSD_TF = dataTable.StopSignalDelay(trials.Trigger_failure);
    else
        RT_TF = NaN;
        MT_TF = NaN;
        TotRespT_TF = NaN;
        SSD_TF = NaN;
    end
    
    figure();
    subplot(3,1,1);
    histogram(RT_GoCorr,'BinWidth',20, 'Normalization','probability','FaceColor','g');
    hold on
    histogram(RT_TF,'BinWidth',20, 'Normalization','probability','FaceColor','r');
    ylabel('Probability', 'FontSize',14);
    xlabel('RT (ms)', 'FontSize',14);
    xlim([0, round(max(RT_GoCorr)+200)]);
    legend('Go','TF', 'FontSize',14);
    title('Are TFs fast or slow?', 'FontSize',14);
    
    subplot(3,1,2);
    histogram(MT_GoCorr,'BinWidth',20, 'Normalization','probability','FaceColor','g');
    hold on
    histogram(MT_TF,'BinWidth',20, 'Normalization','probability','FaceColor','r');
    ylabel('Probability', 'FontSize',14);
    xlabel('MT (ms)', 'FontSize',14);
    xlim([0, round(max(MT_GoCorr)+200)]);
    legend('Go','TF', 'FontSize',14);
    %     title('Are TFs fast or slow?', 'FontSize',14);
    
    subplot(3,1,3);
    histogram(TotRespT_GoCorr,'BinWidth',20, 'Normalization','probability','FaceColor','g');
    hold on
    histogram(TotRespT_TF,'BinWidth',20, 'Normalization','probability','FaceColor','r');
    ylabel('Probability', 'FontSize',14);
    xlabel('TT (ms)', 'FontSize',14);
    xlim([0, round(max(TotRespT_GoCorr)+200)]);
    legend('Go','TF', 'FontSize',14);
    %     title('Are TFs fast or slow?', 'FontSize',14);
    print(gcf,strcat(folderName,'/',subjID,'_PvO_TriggerFailure_ResponseTimes_vs_Go.pdf'),'-dpdf','-r400');
    
    
    %% Averaged 3d plot
    % ---------------------------------------------------------
    % Separate out positions and correct for different X directions
    startPos_mouseX_mn = nanmedian(data.mouseX(:,1));
    
    % Go
    trials.GoCorr_left =  dataTable.Outcome == 'Correct_Go' ...
        & dataTable.TooSoon == 'false' & dataTable.CueDirection == 'Left';
    trials.GoCorr_right =  dataTable.Outcome == 'Correct_Go' ...
        & dataTable.TooSoon == 'false' & dataTable.CueDirection == 'Right';
    
    % Home
    trials.StopCorr_HomeStr_left = dataTable.Staircase == 'HomeStair'...
        & dataTable.Outcome == 'Correct_Stop_Home' & dataTable.CueDirection == 'Left';
    trials.StopCorr_HomeStr_right = dataTable.Staircase == 'HomeStair'...
        & dataTable.Outcome == 'Correct_Stop_Home' & dataTable.CueDirection == 'Right';
    trials.StopFail_HomeStr_left = dataTable.Staircase == 'HomeStair'...
        & dataTable.Outcome == 'Failed_Stop_Home' & dataTable.CueDirection == 'Left' | ...
        dataTable.Staircase == 'HomeStair'...
        & dataTable.Outcome == 'Failed_Stop_Target' & dataTable.CueDirection == 'Left';
    trials.StopFail_HomeStr_right = dataTable.Staircase == 'HomeStair'...
        & dataTable.Outcome == 'Failed_Stop' & dataTable.CueDirection == 'Right' | ...
        dataTable.Staircase == 'HomeStair'...
        & dataTable.Outcome == 'Failed_Stop_Home' & dataTable.CueDirection == 'Right';
    % Target
    trials.StopCorr_TargetStr_left = dataTable.Staircase == 'TargetStair' & ...
        dataTable.Outcome == 'Correct_Stop_Home' & dataTable.CueDirection == 'Left' | ...
        dataTable.Staircase == 'TargetStair' & ...
        dataTable.Outcome == 'Failed_Stop_Home' & dataTable.CueDirection == 'Left';
    trials.StopCorr_TargetStr_right = dataTable.Staircase == 'TargetStair' &...
        dataTable.Outcome == 'Correct_Stop_Home' & dataTable.CueDirection == 'Right' |...
        dataTable.Staircase == 'TargetStair' &...
        dataTable.Outcome == 'Failed_Stop_Home' & dataTable.CueDirection == 'Right';
    trials.StopFail_TargetStr_left = dataTable.Staircase == 'TargetStair'...
        & dataTable.Outcome == 'Failed_Stop_Target' & dataTable.CueDirection == 'Left';
    trials.StopFail_TargetStr_right = dataTable.Staircase == 'TargetStair'...
        & dataTable.Outcome == 'Failed_Stop_Target' & dataTable.CueDirection == 'Right';
    
    % Get positions and times for each direction of movement
    % x pos
    mouseX_GoCorr_left = nanmedian(data.mouseX(trials.GoCorr_left,:));
    mouseX_GoCorr_right = nanmedian(data.mouseX(trials.GoCorr_right,:));
    
    mouseX_StopCorr_HomeStair_left = nanmedian(data.mouseX(trials.StopCorr_HomeStr_left,:));
    mouseX_StopCorr_HomeStair_right = nanmedian(data.mouseX(trials.StopCorr_HomeStr_right,:));
    
    mouseX_StopCorr_TargetStair_left = nanmedian(data.mouseX(trials.StopCorr_TargetStr_left,:));
    mouseX_StopCorr_TargetStair_right = nanmedian(data.mouseX(trials.StopCorr_TargetStr_right,:));
    
    mouseX_StopFail_HomeStair_left = nanmedian(data.mouseX(trials.StopFail_HomeStr_left,:));
    mouseX_StopFail_HomeStair_right = nanmedian(data.mouseX(trials.StopFail_HomeStr_right,:));
    
    mouseX_StopFail_TargetStair_left = nanmedian(data.mouseX(trials.StopFail_TargetStr_left,:));
    mouseX_StopFail_TargetStair_right = nanmedian(data.mouseX(trials.StopFail_TargetStr_right,:));
    
    % y pos
    mouseY_GoCorr_left = nanmedian(data.mouseY(trials.GoCorr_left,:));
    mouseY_GoCorr_right = nanmedian(data.mouseY(trials.GoCorr_right,:));
    
    mouseY_StopCorr_HomeStair_left = nanmedian(data.mouseY(trials.StopCorr_HomeStr_left,:));
    mouseY_StopCorr_HomeStair_right = nanmedian(data.mouseY(trials.StopCorr_HomeStr_right,:));
    
    mouseY_StopCorr_TargetStair_left = nanmedian(data.mouseY(trials.StopCorr_TargetStr_left,:));
    mouseY_StopCorr_TargetStair_right = nanmedian(data.mouseY(trials.StopCorr_TargetStr_right,:));
    
    mouseY_StopFail_HomeStair_left = nanmedian(data.mouseY(trials.StopFail_HomeStr_left,:));
    mouseY_StopFail_HomeStair_right = nanmedian(data.mouseY(trials.StopFail_HomeStr_right,:));
    
    mouseY_StopFail_TargetStair_left = nanmedian(data.mouseY(trials.StopFail_TargetStr_left,:));
    mouseY_StopFail_TargetStair_right = nanmedian(data.mouseY(trials.StopFail_TargetStr_right,:));
    
    % time
    mouseT_GoCorr_left = nanmedian(data.mouseT_alignRT(trials.GoCorr_left,:));
    mouseT_GoCorr_right = nanmedian(data.mouseT_alignRT(trials.GoCorr_right,:));
    
    mouseT_StopCorr_HomeStair_left = nanmedian(data.mouseT_alignRT(trials.StopCorr_HomeStr_left,:));
    mouseT_StopCorr_HomeStair_right = nanmedian(data.mouseT_alignRT(trials.StopCorr_HomeStr_right,:));
    
    mouseT_StopCorr_TargetStair_left = nanmedian(data.mouseT_alignRT(trials.StopCorr_TargetStr_left,:));
    mouseT_StopCorr_TargetStair_right = nanmedian(data.mouseT_alignRT(trials.StopCorr_TargetStr_right,:));
    
    mouseT_StopFail_HomeStair_left = nanmedian(data.mouseT_alignRT(trials.StopFail_HomeStr_left,:));
    mouseT_StopFail_HomeStair_right = nanmedian(data.mouseT_alignRT(trials.StopFail_HomeStr_right,:));
    
    mouseT_StopFail_TargetStair_left = nanmedian(data.mouseT_alignRT(trials.StopFail_TargetStr_left,:));
    mouseT_StopFail_TargetStair_right = nanmedian(data.mouseT_alignRT(trials.StopFail_TargetStr_right,:));
    
    % Plot figure
    figure()
    plot3(mouseX_GoCorr_left, ...
        mouseY_GoCorr_left,...
        mouseT_GoCorr_left,'g','LineWidth',2);
    hold on
    plot3(mouseX_GoCorr_right,...
        mouseY_GoCorr_right,...
        mouseT_GoCorr_right,'g','LineWidth',2);
    plot3(mouseX_StopCorr_TargetStair_left, ... % + 10 - need to sort our precisely when
        mouseY_StopCorr_TargetStair_left,...
        mouseT_StopCorr_TargetStair_left,'--b','LineWidth',2);
    plot3(mouseX_StopCorr_TargetStair_right,...
        mouseY_StopCorr_TargetStair_right,...
        mouseT_StopCorr_TargetStair_right,'--b','LineWidth',2);
    plot3(mouseX_StopFail_TargetStair_left, ...
        mouseY_StopFail_TargetStair_left,...
        mouseT_StopFail_TargetStair_left,'--r','LineWidth',2);
    plot3(mouseX_StopFail_TargetStair_right,...
        mouseY_StopFail_TargetStair_right,...
        mouseT_StopFail_TargetStair_right,'--r','LineWidth',2);
    
    xlabel('Horizontal (px)','FontSize',14);
    ylabel('Vertical (px)','FontSize',14);
    zlabel('Time from go (ms)','FontSize',14);
    title('Average Movement Trajectories', 'FontSize',14);
    legend('GoL','GoR','SSt','FSt','TL','TR','location','SouthEast','FontSize',14)
    view(0,90)
    %     set(gca, 'YDir','reverse')
    print(gcf,strcat(folderName,'/',subjID,'_PvO_Stop_Trajectory2dXZ.pdf'),'-dpdf','-r400');
    
    % view(90,0)
    % print(gcf,strcat(folderName,'/',subjID,'_PvO_Stop_Trajectory2dYZ.pdf'),'-dpdf','-r400');
    
    
    %% Get velocities for all conditions
    % Calculate lateral (X) velocity in px.s
    mouseX_GoCorr = data.mouseX(trials.GoCorr,:);
    mouseX_StopCorr_HomeStair = data.mouseX(trials.StopCorr_HomeStr,:);
    mouseX_StopFail_HomeStair = data.mouseX(trials.StopFail_HomeStr,:);
    mouseX_StopCorr_TargetStair = data.mouseX(trials.StopCorr_TargetStr,:);
    mouseX_StopFail_TargetStair = data.mouseX(trials.StopFail_TargetStr,:);
    
    mouseXdiff_GoCorr = diff(mouseX_GoCorr,1,2);
    mouseXvel_GoCorr = abs(mouseXdiff_GoCorr./(frameRate/1000));
    [xVelvalCorrGo, xVellocCorrGo] = max(mouseXvel_GoCorr');
    
    mouseXdiff_StopCorr_HomeStair = diff(mouseX_StopCorr_HomeStair,1,2);
    mouseXvel_StopCorr_HomeStair = abs(mouseXdiff_StopCorr_HomeStair./(frameRate/1000));
    % mouseXvel_StopCorr_HomeStair = (mouseXdiff_StopCorr_HomeStair./(frameRate/1000));
    [xVelvalCorrHome, xVellocCorrHome] = max(mouseXvel_StopCorr_HomeStair');
    
    mouseXdiff_StopFail_HomeStair = diff(mouseX_StopFail_HomeStair,1,2);
    mouseXvel_StopFail_HomeStair = abs(mouseXdiff_StopFail_HomeStair./(frameRate/1000));
    [xVelvalFailHome, xVellocFailHome] = max(mouseXvel_StopFail_HomeStair');
    
    mouseXdiff_StopCorr_TargetStair = diff(mouseX_StopCorr_TargetStair,1,2);
    mouseXvel_StopCorr_TargetStair = abs(mouseXdiff_StopCorr_TargetStair./(frameRate/1000));
    [xVelvalCorrTarget, xVellocCorrTarget] = max(mouseXvel_StopCorr_TargetStair');
    
    mouseXdiff_StopFail_TargetStair = diff(mouseX_StopFail_TargetStair,1,2);
    mouseXvel_StopFail_TargetStair = abs(mouseXdiff_StopFail_TargetStair./(frameRate/1000));
    [xVelvalFailTarget, xVellocFailTarget] = max(mouseXvel_StopFail_TargetStair');
    
    % average X velocity and time of the peak velocity? expect later time of peak
    % velocity for failed stops reflecting later change of direction
    meanStopCue_target =  MT_GoCorr_mn - mean([SSRT_HomeStair,SSRT_TargetStair_TT]);
    
    figure();
    plot((frameRate/2:frameRate:frameRate*size(mouseXvel_GoCorr,2) - (frameRate/2)),nanmedian(mouseXvel_GoCorr),'-g')
    hold on
    plot([meanStopCue_target,meanStopCue_target],[0, 500],'--k')
    plot([0,0],[0, 500],'-k')
    plot([MT_GoCorr_mn,MT_GoCorr_mn],[0,200],'.-k')
    plot((frameRate/2:frameRate:frameRate*size(mouseXvel_StopCorr_HomeStair,2) - (frameRate/2)),nanmedian(mouseXvel_StopCorr_HomeStair),'-b')
    plot((frameRate/2:frameRate:frameRate*size(mouseXvel_StopCorr_TargetStair,2) - (frameRate/2)),nanmedian(mouseXvel_StopCorr_TargetStair),'--b')
    plot((frameRate/2:frameRate:frameRate*size(mouseXvel_StopFail_HomeStair,2) - (frameRate/2)),nanmedian(mouseXvel_StopFail_HomeStair),'-r')
    plot((frameRate/2:frameRate:frameRate*size(mouseXvel_StopFail_TargetStair,2) - (frameRate/2)),nanmedian(mouseXvel_StopFail_TargetStair),'--r')
    plot([meanStopCue_target,meanStopCue_target],[0,500],'-k')
    nanmedian(TotRespT_GoCorr)
    legend('Go','Stop','MoveOn','MoveEnd','SSH','SST','FSH','FST','location','NorthEast','FontSize',14);
    xlabel('Time from movement onset (ms)','FontSize',14);
    ylabel('X velocity (px.s)','FontSize',14);
    xlim([-100, MT_GoCorr_mn + 100]);
    ylim([0, nanmax(nanmedian(mouseXvel_GoCorr))+250]);
    text(meanStopCue_target +10, 150, 'Stop (Targ.)','FontSize',14,'Rotation', -90);
    text(MT_GoCorr_mn +10, 150, 'Target Hit','FontSize',14,'Rotation', -90);
    print(gcf,strcat(folderName,'/',subjID,'_PvO_Stop_VelocityAbsX.pdf'),'-dpdf','-r400');
    
    % Calculate vertical (Y) velocity in px.s
    mouseY_GoCorr = data.mouseY(trials.GoCorr,:);
    mouseY_StopCorr_HomeStair = data.mouseY(trials.StopCorr_HomeStr,:);
    mouseY_StopFail_HomeStair = data.mouseY(trials.StopFail_HomeStr,:);
    mouseY_StopCorr_TargetStair = data.mouseY(trials.StopCorr_TargetStr,:);
    mouseY_StopFail_TargetStair = data.mouseY(trials.StopFail_TargetStr,:);
    
    mouseYdiff_GoCorr = diff(mouseY_GoCorr,1,2);
    mouseYvel_GoCorr = (mouseYdiff_GoCorr./(frameRate/1000));
    [yVelvalCorrGo, yVellocCorrGo] = max(mouseYvel_GoCorr');
    
    mouseYdiff_StopCorr_HomeStair = diff(mouseY_StopCorr_HomeStair,1,2);
    mouseYvel_StopCorr_HomeStair = (mouseYdiff_StopCorr_HomeStair./(frameRate/1000));
    [yVelvalCorrHome, yVellocCorrHome] = max(mouseYvel_StopCorr_HomeStair');
    
    mouseYdiff_StopFail_HomeStair = diff(mouseY_StopFail_HomeStair,1,2);
    mouseYvel_StopFail_HomeStair = (mouseYdiff_StopFail_HomeStair./(frameRate/1000));
    [yVelvalFailHome, yVellocFailHome] = max(mouseYvel_StopFail_HomeStair');
    
    mouseYdiff_StopCorr_TargetStair = diff(mouseY_StopCorr_TargetStair,1,2);
    mouseYvel_StopCorr_TargetStair = (mouseYdiff_StopCorr_TargetStair./(frameRate/1000));
    [yVelvalCorrTarget, yVellocCorrTarget] = max(mouseYvel_StopCorr_TargetStair');
    
    mouseYdiff_StopFail_TargetStair = diff(mouseY_StopFail_TargetStair,1,2);
    mouseYvel_StopFail_TargetStair = (mouseYdiff_StopFail_TargetStair./(frameRate/1000));
    [yVelvalFailTarget, yVellocFailTarget] = max(mouseYvel_StopFail_TargetStair');
    
    figure();
    plot(((frameRate/2:frameRate:frameRate*size(mouseYvel_GoCorr,2))- (frameRate/2)),nanmedian(mouseYvel_GoCorr),'-g');
    hold on
    plot([meanStopCue_target,meanStopCue_target],[0, 1000],'--k')
    plot([0,0],[-1000, 1000],'-k')
    plot([MT_GoCorr_mn,MT_GoCorr_mn],[0,1000],'.-k')
    plot((frameRate/2:frameRate:frameRate*size(mouseYvel_StopCorr_HomeStair,2)- (frameRate/2)),nanmedian(mouseYvel_StopCorr_HomeStair),'-b');
    plot((frameRate/2:frameRate:frameRate*size(mouseYvel_StopCorr_TargetStair,2)- (frameRate/2)),nanmedian(mouseYvel_StopCorr_TargetStair),'--b');
    plot((frameRate/2:frameRate:frameRate*size(mouseYvel_StopFail_HomeStair,2)- (frameRate/2)),nanmedian(mouseYvel_StopFail_HomeStair),'-r');
    plot((frameRate/2:frameRate:frameRate*size(mouseYvel_StopFail_TargetStair,2)- (frameRate/2)),nanmedian(mouseYvel_StopFail_TargetStair),'--r');
    legend('Go','Stop','MoveOn','MoveEnd','SSH','SST','FSH','FST','location','NorthEast','FontSize',14);xlabel('Time from movement onset (ms)','FontSize',14);
    ylabel('Y velocity (px.s)','FontSize',14);
    xlim([-100, MT_GoCorr_mn + 100]);
    ylim([min(nanmedian(mouseYvel_GoCorr)),round(nanmax(nanmedian(mouseYvel_GoCorr))+300)]);
    text(meanStopCue_target +10, 500, 'Stop (Targ.)','FontSize',14,'Rotation', -90);
    text(MT_GoCorr_mn +10, 500, 'Target Hit','FontSize',14,'Rotation', -90);
    print(gcf,strcat(folderName,'/',subjID,'_PvO_Stop_VelocityAbsY.pdf'),'-dpdf','-r400');
    
    
    %% Put key variables in an array
    
    summary_behaviour = [trials.pcnt_omit_go, trials.pcnt_miss_go,  trials.pcntErr_go, trials.pcnt_FalseAlarm_allTrl,...
        RT_GoCorr_mn, MT_GoCorr_mn, TotRespT_GoCorr_mn, ...
        RT_std_GoCorr, MT_std_GoCorr, TT_std_GoCorr, ...
        RT_StopFail_HomeStair_mn, MT_StopFail_HomeStair_mn,TotRespT_StopFail_HomeStair_mn, ...
        RT_StopCorr_TargetStair_mn, ...
        RT_StopFail_TargetStair_mn, MT_StopFail_TargetStair_mn, TotRespT_StopFail_TargetStair_mn, ...
        SSD_Stop_HomeStair_mn, SSD_Stop_TargetStair_mn, ...
        pStop_HomeStair, pStop_TargetStair, ...
        SSRT_HomeStair, SSRT_TargetStair_MT, SSRT_TargetStair_TT, ...
        n_trigger_failure, nanmedian(RT_TF), nanmedian(MT_TF), ...
        nanmedian(TotRespT_TF), nanmedian(SSD_TF)];
    data.Summary_Behavior = array2table(summary_behaviour, ...
        'VariableNames',...
        {'GoOmit', 'GoMiss', 'GoErr', 'FalseAlarm', ...
        'Go_RT','Go_MT', 'Go_TT', ...
        'Go_RTstd', 'Go_MTstd', 'Go_TTstd', ...
        'FS_RTh', 'FS_MTh', 'FS_TTh',...
        'SS_RTt', ...
        'FS_RTt', 'FS_MTt', 'FS_TTt', ...
        'SSDh','SSDt', ...
        'pStoph','pStopt', ...
        'SSRTh_RT','SSRTt_MT' 'SSRTt_TT', ...
        'N_TFh','TF_RT', 'TF_MT', ...
        'TF_TT', 'TF_SSD'});
    
    
    %% Save data
    
    matfile = fullfile(pathnamedat, strcat('IndieData_', string(subjID), '.mat'));
    save(matfile, 'data', 'dataTable', 'trials');
    
    
    %     loop but comment out if running single person to check
    clearvars -except opts subjs pathnamedat;
    close all
    
end

clearvars
close all

