diff --git a/applyCARLARealCCEPs.m b/applyCARLARealCCEPs.m
index e179bfb..cda2ec2 100644
--- a/applyCARLARealCCEPs.m
+++ b/applyCARLARealCCEPs.m
@@ -30,10 +30,11 @@
 elecsPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_electrodes.tsv', sub));
 elecs = ieeg_readtableRmHyphens(elecsPath);
 
-mefPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_ieeg.mefd', sub, run));
+% known issue: Matlab crashing when readMef3 loads mef data from relative path. Temporary fix = fullfile with pwd. If you specified an absolute "dataPath", remove "pwd".
+mefPath = fullfile(pwd, dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_ieeg.mefd', sub, run));
 channelsPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_channels.tsv', sub, run));
 eventsPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_events.tsv', sub, run));
-annotPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_eventsAnnot.tsv', sub, run));
+annotPath = fullfile(dataPath, 'derivatives', 'event_annotations', sub, sprintf('%s_ses-ieeg01_task-ccep_run-%02d_eventsAnnot.tsv', sub, run));
 
 % Load CCEP data
 mef = ccep_PreprocessMef(mefPath, channelsPath, eventsPath);
@@ -41,7 +42,7 @@
 mef.highpass;
 mef.loadMefTrials([-1, 1]);
 
-mef.plotOutputs(1, [], 200);
+%mef.plotOutputs(1, [], 200);
 
 tt = mef.tt;
 srate = mef.srate;
@@ -54,7 +55,8 @@
 outdir = fullfile('output', 'realCCEPs');
 mkdir(outdir);
 
-try
+clear annots
+try % subject 4 has no annots because no interictal activity, hence the try
     annots = readtable(annotPath, 'FileType', 'text', 'Delimiter', '\t');
     annots = annots.status_description;
     assert(length(annots) == height(events), 'Error: length of event annotations does not match events table height');
@@ -67,7 +69,7 @@
 disp('SOZs:');
 disp(SOZ');
 
-clear *Path*
+clear elecsPath mefPath channelsPath eventsPath annotPath % cleaning
 
 %% Extract data from site of interest, NaN out individual bad trials
 
@@ -82,8 +84,10 @@
 chsSite = channels.name(goodChs);
 
 % remove bad trials and set channels to nan in individual trials, according to annotations
-annotsSite = annots(idxesSite); % get annotations corresponding to this stim site
-dataSite = rmBadTrialsAnnots(dataSite, chsSite, annotsSite);
+if exist('annots', 'var')
+    annotsSite = annots(idxesSite); % get annotations corresponding to this stim site
+    dataSite = rmBadTrialsAnnots(dataSite, chsSite, annotsSite);
+end
 
 figure('Position', [200, 200, 1200, 1200]);
 subplot(1, 2, 1); % first half of channels
diff --git a/applyCARLARealCCEPsLoop.m b/applyCARLARealCCEPsLoop.m
index 0bc346a..4da0acd 100644
--- a/applyCARLARealCCEPsLoop.m
+++ b/applyCARLARealCCEPsLoop.m
@@ -44,10 +44,11 @@
         
         run = runs{ss}(rr);
     
-        mefPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_ieeg.mefd', sub, run));
+        % known issue: Matlab crashing when readMef3 loads mef data from relative path. Temporary fix = fullfile with pwd. If you specified an absolute "dataPath", remove "pwd".
+        mefPath = fullfile(pwd, dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_ieeg.mefd', sub, run));
         channelsPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_channels.tsv', sub, run));
         eventsPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_events.tsv', sub, run));
-        annotPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_task-ccep_run-%02d_eventsAnnot.tsv', sub, run));
+        annotPath = fullfile(dataPath, 'derivatives', 'event_annotations', sub, sprintf('%s_ses-ieeg01_task-ccep_run-%02d_eventsAnnot.tsv', sub, run));
         elecsPath = fullfile(dataPath, sub, 'ses-ieeg01', 'ieeg', sprintf('%s_ses-ieeg01_electrodes.tsv', sub));
     
         % Load CCEP data from each subject
@@ -65,7 +66,8 @@
         sites = groupby(events.electrical_stimulation_site);
     
         % read manual annotations on events
-        try
+        clear annots
+        try % subject 4 has no annots because no interictal activity, hence the try
             annots = readtable(annotPath, 'FileType', 'text', 'Delimiter', '\t');
             annots = annots.status_description;
         catch
diff --git a/functions/groupby.m b/functions/groupby.m
new file mode 100644
index 0000000..89ea915
--- /dev/null
+++ b/functions/groupby.m
@@ -0,0 +1,32 @@
+%% Creates a 2-column cell array where the first column indicates unique occurreences in the input, and second col tracks indices of each occurrence
+% Modeled after Python groupby function
+%
+%   2024/02/14
+%
+%    If this code is used in a publication, please cite the manuscript:
+%    "CARLA: Adjusted common average referencing for cortico-cortical evoked potential data"
+%    by H Huang, G Ojeda Valencia, NM Gregg, GM Osman, MN Montoya,
+%    GA Worrell, KJ Miller, and D Hermes.
+%
+%    CARLA manuscript package.
+%    Copyright (C) 2023 Harvey Huang
+%
+%    This program is free software: you can redistribute it and/or modify
+%    it under the terms of the GNU General Public License as published by
+%    the Free Software Foundation, either version 3 of the License, or
+%    (at your option) any later version.
+%
+%    This program is distributed in the hope that it will be useful,
+%    but WITHOUT ANY WARRANTY; without even the implied warranty of
+%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+%    GNU General Public License for more details.
+%
+%    You should have received a copy of the GNU General Public License
+%    along with this program.  If not, see <https://www.gnu.org/licenses/>.
+%
+function groups = groupby(strs)
+    
+    gs = cellfun(@(x) find(strcmp(x, strs)), unique(strs), 'UniformOutput', false);
+    groups = [unique(strs), gs];
+    
+end
\ No newline at end of file
diff --git a/main.m b/main.m
index 345b886..c453e6a 100644
--- a/main.m
+++ b/main.m
@@ -37,7 +37,7 @@
 % Make sure you are currently working in the root manuscript directory (CARLA)
 addpath('functions');
 
-% Data should be downloaded from OpenNeuro and copied to a folder here (./data)
+% Data should be downloaded from OpenNeuro and copied to a folder here in the root directory (./data)
 dataPath = 'data';
 
 % Add paths to other dependencies. Change the dummy paths listed here to where they are located on your system.
@@ -140,7 +140,7 @@
 %% Figure S3: Average channel responses for sub 1, stim site 1 after applying CAR vs. CARLA
 
 % This supplemental figure is a direct addendum to Figure 5, showing the mean signal for the first 50 channels after rereferencing by CARLA or CAR
-% Run after the previous section, keeping all workspace variables
+% Run after the previous section, keeping all workspace variables. Output also saved to .\output\realCCEPs\.
 figS3_compareCARvCARLA
 
 %% Figure S4: Apply CARLA to real data at another 4 stim sites in subjects 1-4