%
% This is a batch script to get the median (mean, trimmed mean) 
%  of each subject's RT data, plot the grand mean and standard error 
%  for selected conditions. 
% 
% by Jason Taylor (17/11/2008)
% + updated (jt 17/11/2008): added error bars
% + updated (jt 17/11/2008): added plotvar, options, conds loop
% + updated (jt 13/01/2011): added extract raw RTs 
% + updated (jt 20/11/2011): some small cosmetic improvements
% + updated (jt 19/11/2012): changed axis scaling; made windows compatible

%% (0) Define options:

% Plot format:
barcolor  = [.5 .5 .5];
ebarcolor = [0 0 0];
ebarsize  = 3;
plotfont  = 12;

% Processing options:
plotvar = 'median';  % 'median', 'mean', 'trim<N>' (N%-trimmed mean)
dosave  = 0;         % save grandmean data?
doplot  = 1;         % plot grandmean data?

% Data options:
conds     = [1 2];
condlabs  = {'word', 'nonword'};
Nevents   = [240 240];


%% (1) Define directory, filename, subject parameters:

% Project directory:
if ispc
    projdir = '\\cbsu\data\imaging\jt03\demo\rtdata\subjects';
else
    projdir = '/imaging/jt03/demo/rtdata/subjects';
end

if dosave
    % Working directory: (CHANGE TO A DIR YOU ARE PERMITTED TO WRITE IN!)
    %wkdir = '/imaging/jt03/demo/rtdata/subjects/ga15';
    %wkdir = '/imaging/<user>/jt_rt_demo/ga15';
    wkdir = uigetdir(projdir,'Select directory to save group results in');
end

% Filename (input):
fname = 'word_nonword.mat';

% Subjects:
subjects = [1:15];

 
%% (2) Get each subject's median/mean/trimmed mean RT:

% Initialise variable to collect md/m/tm RTs:
mrt = zeros(length(subjects),length(conds));
RT  = zeros(Nevents(1),length(subjects),length(conds));

% Loop over subjects:
for i = 1:length(subjects)

	% Loop over conditions
	for j = 1:length(conds)
	
		% Get current subject label:
		subj = sprintf('s%02d',subjects(i)); 
	
		% Go to subject's directory, load data:
		cd(fullfile(projdir,subj));
		load(fname);
		rt = D.rt(D.event==conds(j));
    
        % Store raw RT:
        RT(:,i,j) = rt;
    
		switch plotvar
    		case 'median'
        		mrt(i,j) = median(rt);
            case 'mean'
                mrt(i,j) = mean(rt);
    		otherwise  % trim<N>
        		trimpct = str2double(plotvar(5:end));
            	rt = rt(rt>prctile(rt,trimpct/2) & rt<prctile(rt,100-(trimpct/2)));
    			mrt(i,j) = mean(rt);
		end % switch

	end % j in conds
	
end % i in subjects

cd(projdir)


%% (3) Compute grand mean, standard error:

% Compute mean (collapsing over rows):
gm = mean(mrt,1);

% Get standard error:
se = std(mrt)/sqrt(size(mrt,1));

% Save it as a .mat file in working directory:
if dosave
	cd(wkdir)
	save rtdata.mat mrt gm se
end


%% (4) Plot:

if doplot
	% Open a figure, make background white:
	fig = figure('color',[1 1 1]);
	
	% Plot means:
	bar(gm, 'facecolor', barcolor);
	
	% Rescale:
	ymax = ceil(max(gm+se));
	ymin = floor(min(gm-se));
    yrng = ymax-ymin;
	set(gca, 'ylim', [ymin-(.25*yrng) ymax+(.25*yrng)]);
	
	% Plot and format error bars:
	for j = 1:length(conds)
		ebar1 = line([j j],[gm(j) gm(j)+se(j)]);
		ebar2 = line([j j],[gm(j) gm(j)-se(j)]);
		set([ebar1 ebar2], 'linewidth', ebarsize, 'color', ebarcolor);
	end
	
	% Apply title, labels, etc.:
	ttl = title(sprintf('Grand Mean of %s RTs',plotvar));
	xl  = xlabel('Stimulus Type');
	yl  = ylabel('RT + SEM (ms)');
	set(gca,'xticklab',condlabs);
    set(gca,'fontsize',plotfont);
	set([ttl xl yl],'fontsize',plotfont+.20*plotfont);
    box off
    
end

% End gracefully:
fprintf(1,'\n++ done! ++\n\n');
