diff --git a/5_reconstruction/gtForwardSimulate_v2.m b/5_reconstruction/gtForwardSimulate_v2.m
index 6f7be299587073acf91b95a3e5b2e63b3f0e42a2..ab609103060f0daace20a14a9ffd139914438738 100644
--- a/5_reconstruction/gtForwardSimulate_v2.m
+++ b/5_reconstruction/gtForwardSimulate_v2.m
@@ -137,20 +137,20 @@ for ii = 1: length(difspotIDs)
 end
 
     
-mymcmd    = sprintf('select ldirX, ldirY, ldirZ, omega, omegaB, difAID, difBID from %s ',acq.pair_tablename);
+mymcmd    = sprintf('select ldirX, ldirY, ldirZ, omega, omegaB, difAID, difBID from %s ', acq.pair_tablename);
 [allpairs.ldirX, allpairs.ldirY, allpairs.ldirZ, allpairs.omegaA, allpairs.omegaB, allpairs.difAID, allpairs.difBID] = mym(mymcmd);
 
 difspot_num = length(difspotIDs);   % number of difspots in database
 spot2grain  = cell(difspot_num, 1); % cell array linking difspots with grainids
 spot2phase  = cell(difspot_num, 1);
 
-filetable   = sprintf('%s_filetable',acq.name);
+filetable   = sprintf('%s_filetable', acq.name);
 filename    = fullfile(acq.dir,'4_grains','spot2grain.mat');
 spotTable   = GtMergeTable(filetable, filename, {'spot2grain', 'spot2phase'});
 
 samplefile  = fullfile(acq.dir,'4_grains','sample.mat');
-sample      = GtSample(phases_num);
-sample.fileTable = sprintf('%s_filetable',parameters.acq.name);
+sample      = GtSample.loadFromFile(samplefile);
+%sample.fileTable = sprintf('%s_filetable',parameters.acq.name);
 % initialise sample structure
 phase = GtPhase(parameters.cryst(phaseID).name, grain_num);
 sample.phases{phaseID} = phase;
diff --git a/5_reconstruction/gtSetupForwardSimulation.m b/5_reconstruction/gtSetupForwardSimulation.m
index 008f68328a8f0700e1642084969699144cd2b8a7..65e63d6ae368e0dba3c7f181b81fe39cda7c24f5 100644
--- a/5_reconstruction/gtSetupForwardSimulation.m
+++ b/5_reconstruction/gtSetupForwardSimulation.m
@@ -26,11 +26,13 @@ disp(['Setting up Forward Simulation for phase ' phaseID_str ' ' parameters.crys
 filetable = sprintf('%s_filetable', parameters.acq.name);
 % get the latest version of indexter output
 grain = [];
-filename = fullfile(parameters.acq.dir, '4_grains', ['phase_' phaseID_str], 'index.mat');
 
-load(filename, 'grain');
-% how many grains in this dataset?
-ngrains = length(grain);
+for ii = 1 : nof_phases
+    filename = fullfile(parameters.acq.dir, '4_grains', sprintf('phase_%02d', ii), 'index.mat');
+    load(filename, 'grain');
+    % how many grains in this dataset?
+    ngrains(ii) = length(grain);
+end    
 
 
 %create / reset sample.mat spot2grain.mat
@@ -40,14 +42,17 @@ if exist(filename, 'file') && (nof_phases > 1)
 else
     sample = GtSample(nof_phases);
     sample.fileTable = sprintf('%s_filetable',parameters.acq.name);
+    for ii = 1 : nof_phases
+        phase = GtPhase(parameters.cryst(ii).name, ngrains(ii));
+        sample.phases{ii} = phase;
+    end
+    save(filename, 'sample');
 end
 
-phase = GtPhase(parameters.cryst(phaseID).name, ngrains);
-sample.phases{phaseID} = phase;
-save(filename, 'sample');
+
 disp(['Reset phase number ' num2str(phaseID) ' in the sample.mat:'])
 disp(['Name : ' parameters.cryst(phaseID).name])
-disp(['nof_grains: ' num2str(ngrains)])
+disp(['nof_grains: ' num2str(ngrains(phaseID))])
 
 
 filename    = fullfile(parameters.acq.dir,'4_grains','spot2grain.mat');
@@ -74,8 +79,8 @@ check=inputwdefault('Launch Forward Simulation with OAR? [y/n]', 'y');
 if strcmpi(check, 'y')
     disp('Launching several separate OAR jobs')
     
-    oar.njobs    = min(ceil(ngrains/20), 30); % max 30 jobs
-    oar.walltime = ngrains / oar.njobs * 5;  % allow 5 minutes per grain
+    oar.njobs    = min(ceil(ngrains(phaseID)/20), 30); % max 30 jobs
+    oar.walltime = ngrains(phaseID) / oar.njobs * 5;  % allow 5 minutes per grain
     
     OARtmp   = ceil(oar.walltime); % total minutes (could be > 60)
     OARm     = mod(OARtmp, 60); % minutes
@@ -89,12 +94,12 @@ if strcmpi(check, 'y')
     OAR_parameters=gtModifyStructure(OAR_parameters, OAR_list, 'Parameters for OAR');
     
 
-    OAR_make('gtForwardSimulate_v2', 1, ngrains, OAR_parameters.njobs, ...
+    OAR_make('gtForwardSimulate_v2', 1, ngrains(phaseID), OAR_parameters.njobs, ...
         [parameters.acq.dir ' ' num2str(phaseID)], true, ...
         'walltime', OAR_parameters.walltime);
     
 else
-    gtForwardSimulate_v2(1, ngrains, pwd, phaseID);
+    gtForwardSimulate_v2(1, ngrains(phaseID), pwd, phaseID);
 end
 
 end % end of function