diff --git a/zUtil_Fit/gtFITParameters.m b/zUtil_Fit/gtFITParameters.m
index bc54275895a6d4907324f091264c67831c6e2d58..f554e2426d7600fbfc6d092bff870d8f96d0436f 100644
--- a/zUtil_Fit/gtFITParameters.m
+++ b/zUtil_Fit/gtFITParameters.m
@@ -154,7 +154,12 @@ if isempty(par.energy)
 end
 
 if isempty(par.labgeoini)
-    par.labgeoini = parameters.labgeo;
+    disp('Using initial labgeo parameters from the parameters file.')
+    if (isfield(parameters, 'detgeo'))
+        par.labgeoini = gtGeoConvertDetgeo2LegacyLabgeo(parameters.detgeo, parameters.labgeo);
+    else
+        par.labgeoini = parameters.labgeo;
+    end
 end
 
 par.labgeoini.detorig = gtGeoDetOrig(par.labgeoini);
@@ -498,19 +503,19 @@ for ii = 1:length(gr)
 	end
 end
 
-
 % Global
+[detgeo, labgeo] = gtGeoConvertLegacyLabgeo2Detgeo(par.labgeoini);
 if any(activeglobal)
-    [labgeo, energy] = gtFitUpdateLabgeo(varsopt(1:sum(activeglobal)), ...
-                       activeglobal, par.iniglobal, par.labgeoini);
+    [labgeo, detgeo, energy] = gtFitUpdateLabgeo(varsopt(1:sum(activeglobal)), ...
+                       activeglobal, par.iniglobal, labgeo, detgeo);
+    labgeo = rmfield(labgeo, 'detnorm');
+    labgeo = rmfield(labgeo, 'detorig');
 else
-    labgeo = par.labgeoini;
     energy = par.iniglobal(10);
 end
 labgeo = rmfield(labgeo,'detnorm');
 labgeo = rmfield(labgeo,'detorig');
 
-
 % Drifts
 drift = par.drift;
 if ~isempty(par.ind.fit_drift_posX); drift.posx = [0, varsopt(par.ind.fit_drift_posX)']; end
@@ -521,8 +526,6 @@ if ~isempty(par.ind.fit_drift_rotY); drift.roty = [0, varsopt(par.ind.fit_drift_
 if ~isempty(par.ind.fit_drift_rotZ); drift.rotz = [0, varsopt(par.ind.fit_drift_rotZ)']; end
 if ~isempty(par.ind.fit_drift_energy); drift.energy = [0, varsopt(par.ind.fit_drift_energy)']; end
 
-
-
 % Plot strain
 meanstrain = mean(grainfit.strain, 2);
 
@@ -550,6 +553,7 @@ end
 % end
 
 % Output structure
+opt.detgeo        = detgeo;
 opt.labgeo        = labgeo;
 opt.energy        = energy;
 opt.grain         = grain;
@@ -587,7 +591,6 @@ sfPlotRes(opt)
 fprintf('\nTotal elapsed time (sec): %g\n', toc(tStart))
 beep, pause(0.5), beep, pause(0.5), beep
 
-
 end % of function
 
 
diff --git a/zUtil_Fit/gtFitFwdSimGrains.m b/zUtil_Fit/gtFitFwdSimGrains.m
index 688c872088150e1c0b936bbdc16d4bc3e5320e42..f1eb90f91819b145ef28b7bb631824b4be1a8138 100644
--- a/zUtil_Fit/gtFitFwdSimGrains.m
+++ b/zUtil_Fit/gtFitFwdSimGrains.m
@@ -43,28 +43,23 @@ end
 
 
 labgeo = parameters.labgeo;
+if (~isfield(labgeo, 'omstep') || isempty(labgeo.omstep))
+    labgeo.omstep = gtGetOmegaStepDeg(parameters);
+end
+if (~isfield(parameters, 'detgeo') || isempty(parameters.detgeo))
+    parameters.detgeo = gtGeoConvertLegacyLabgeo2Detgeo(labgeo);
+end
+detgeo = parameters.detgeo;
 samgeo = parameters.samgeo;
 cryst  = parameters.cryst(grain{1}.phaseid);
 energy = parameters.acq.energy;
-omstep = 180/parameters.acq.nproj;
 
 beamdir = labgeo.beamdir';
 rotdir  = labgeo.rotdir';
 
-% Scaling for Lab -> Detector transformation (mm -> pixels)
-detscaleu = 1/labgeo.pixelsizeu ;
-detscalev = 1/labgeo.pixelsizev ;
-
 % Rotation tensor components
 rotcomp = gtMathsRotationMatrixComp(labgeo.rotdir','col');
 
-% Detector projection matrix
-Qdet = gtFedDetectorProjectionTensor(labgeo.detdiru', labgeo.detdirv', ...
-                                     detscaleu, detscalev);
-
-% Detector normal
-detnorm = cross(labgeo.detdiru, labgeo.detdirv)';
-
 % B matrix for transformation from (hkl) to Cartesian Crystal reference 
 Bmat = gtCrystHKL2CartesianMatrix(cryst.latticepar);
 
@@ -139,8 +134,8 @@ for ii = 1:length(grain)
         gcsam = gcsam(:,ones(1,4*size(plsam,2)));
         
         % Absolute detector coordinates U,V in pixels
-        uv = gtFedPredictUVMultiple(rot4, dveclab, gcsam, labgeo.detrefpos', ...
-            detnorm, Qdet, [labgeo.detrefu; labgeo.detrefv]);
+        uv = gtFedPredictUVMultiple(rot4, dveclab, gcsam, detgeo.detrefpos', ...
+            detgeo.detnorm, detgeo.Qdet, [detgeo.detrefu; detgeo.detrefv]);
         
         
         % Absolute coordinates U,V,W (pixel,pixel,image); to account for
@@ -148,7 +143,7 @@ for ii = 1:length(grain)
         % to be mod(om,360)/omstep
         grain{ii}.fsim.cu        = uv(1,:);
         grain{ii}.fsim.cv        = uv(2,:);
-        grain{ii}.fsim.cw        = reshape(om4',[],1)'/omstep;
+        grain{ii}.fsim.cw        = reshape(om4',[],1)'/labgeo.omstep;
         grain{ii}.fsim.omind     = reshape(omind4',[],1)';
         grain{ii}.fsim.eta       = gtGeoEtaFromDiffVec(dveclab', labgeo)';
         grain{ii}.fsim.sinth     = [sinth, sinth, sinth, sinth];
@@ -252,8 +247,8 @@ for ii = 1:length(grain)
         %dveclab4 = gtFedPredictDiffVecMultiple(pllabd4, beamdir);
         
         % Absolute detector coordinates U,V in pixels
-        uv = gtFedPredictUVMultiple(rot, dveclab, gcsamd, labgeo.detrefpos', ...
-             detnorm, Qdet, [labgeo.detrefu; labgeo.detrefv]);
+        uv = gtFedPredictUVMultiple(rot, dveclab, gcsamd, detgeo.detrefpos', ...
+             detgeo.detnorm, detgeo.Qdet, [detgeo.detrefu; detgeo.detrefv]);
           
          %         uv2 = gtFedPredictUVMultiple(rot2, dveclab, gcsamd2, labgeo.detrefpos', ...
          %               detnorm, Qdet, [labgeo.detrefu; labgeo.detrefv]);
@@ -268,7 +263,7 @@ for ii = 1:length(grain)
         % to be mod(om,360)/omstep
         grain{ii}.fsim.cu        = uv(1,:);
         grain{ii}.fsim.cv        = uv(2,:);
-        grain{ii}.fsim.cw        = [om1, om2, om3, om4]/omstep;
+        grain{ii}.fsim.cw        = [om1, om2, om3, om4]/labgeo.omstep;
         grain{ii}.fsim.omind     = reshape(omind',[],1)';
         grain{ii}.fsim.eta       = gtGeoEtaFromDiffVec(dveclab', labgeo)';
         grain{ii}.fsim.sinth     = [sinthd1, sinthd2, sinthd3, sinthd4];
@@ -283,8 +278,8 @@ for ii = 1:length(grain)
 
     % Set NaN for spots that fall outside the detector area 
     if (~offdet)
-        offd = (uv(1,:) <= 0.5) | (labgeo.detsizeu+0.5 <= uv(1,:)) | ...
-               (uv(2,:) <= 0.5) | (labgeo.detsizev+0.5 <= uv(2,:));
+        offd = (uv(1,:) <= 0.5) | (detgeo.detsizeu+0.5 <= uv(1,:)) | ...
+               (uv(2,:) <= 0.5) | (detgeo.detsizev+0.5 <= uv(2,:));
 
         fn = fieldnames(grain{ii}.fsim);
         for jj = 1:length(fn)
diff --git a/zUtil_Fit/gtFitSpotPosDeviations.m b/zUtil_Fit/gtFitSpotPosDeviations.m
index 1438d1a526bbe065bc64cb0e90e69e54e3eb15f6..2f851e0968d92f4ab6136dbf8cf0e5597e9b209c 100644
--- a/zUtil_Fit/gtFitSpotPosDeviations.m
+++ b/zUtil_Fit/gtFitSpotPosDeviations.m
@@ -120,26 +120,18 @@ end
 %%% Update labgeo
 %%%%%%%%%%%%%%%%%%%%%%%%%
 
+[labgeo, detgeo] = gtGeoConvertLegacyLabgeo2Detgeo(par.labgeoini);
 % Recalculate labgeo
 if any(activeglobal)    
-    [labgeo, energy] = gtFitUpdateLabgeo(newvars(1:sum(activeglobal)), ...
-                       activeglobal, par.iniglobal, par.labgeoini);
+    [labgeo, detgeo, energy] = gtFitUpdateLabgeo( ...
+        newvars(1:sum(activeglobal)), activeglobal, par.iniglobal, labgeo, detgeo);
 else
-    labgeo = par.labgeoini;
     energy = par.iniglobal(10);
 end
 
 % Rotation tensor components
 rotcomp = gtMathsRotationMatrixComp(labgeo.rotdir','col');
 
-% Scaling for Lab -> Detector transformation (mm -> pixels)
-detscaleu = 1/labgeo.pixelsizeu ;
-detscalev = 1/labgeo.pixelsizev ;
-
-% Detector projection matrix
-Qdet = gtFedDetectorProjectionTensor(labgeo.detdiru', labgeo.detdirv', ...
-                                     detscaleu, detscalev);
-
 % In case the Sample axes are not aligned with the Lab axes at omega=0,
 % the plane normals need to be brought in the Lab reference where the
 % beam direction and rotation axis are defined.
@@ -233,8 +225,8 @@ sp.om = sp.om - (ind == 2)*par.numimages + (ind == 3)*par.numimages;
 sp.dveclab = gtFedPredictDiffVecMultiple(sp.pl, labgeo.beamdir');
 
 % Absolute detector coordinates U,V in pixels
-uv = gtFedPredictUVMultiple(sp.rot, sp.dveclab, sp.gcenter, labgeo.detrefpos', ...
-    labgeo.detnorm', Qdet, [labgeo.detrefu; labgeo.detrefv]);
+uv = gtFedPredictUVMultiple(sp.rot, sp.dveclab, sp.gcenter, detgeo.detrefpos', ...
+    detgeo.detnorm', detgeo.Qdet, [detgeo.detrefu; detgeo.detrefv]);
 
 % Absolute coordinates U,V,W (pixel,pixel,image); to account for
 % measurement errors, do not correct for the image number
diff --git a/zUtil_Indexter/gtINDEXDrawGrainUnitCells.m b/zUtil_Indexter/gtINDEXDrawGrainUnitCells.m
index 260fb651aca71d0fb366d290deb96ba0c06b9e4e..eb8829a93cf2a7b1ac737d193be9c82ccdc363be 100644
--- a/zUtil_Indexter/gtINDEXDrawGrainUnitCells.m
+++ b/zUtil_Indexter/gtINDEXDrawGrainUnitCells.m
@@ -59,15 +59,24 @@ sample.rad = parameters.labgeo.samenvrad;
 sample.bot = parameters.labgeo.samenvbot;
 sample.top = parameters.labgeo.samenvtop;
 
+labgeo = parameters.labgeo;
+if (~isfield(labgeo, 'omstep') || isempty(labgeo.omstep))
+    labgeo.omstep = gtGetOmegaStepDeg(parameters);
+end
+if (~isfield(parameters, 'detgeo') || isempty(parameters.detgeo))
+    parameters.detgeo = gtGeoConvertLegacyLabgeo2Detgeo(labgeo);
+end
+detgeo = parameters.detgeo;
+
 nof_grains = length(grains);
 
 fs = 20;
 samradmargin = 0.08;%0.05
 samtopmargin = 0.4;%0.2
-pixelsize = (parameters.labgeo.pixelsizeu + parameters.labgeo.pixelsizev)/2;
+pixelsize = (detgeo.pixelsizeu + detgeo.pixelsizev)/2;
 
-if isfield(parameters.labgeo, 'sourcepoint')
-    pixelsize = pixelsize * parameters.labgeo.sourcepoint(1)/(parameters.labgeo.detrefpos(1));
+if (isfield(parameters.labgeo, 'sourcepoint'))
+    pixelsize = pixelsize * parameters.labgeo.sourcepoint(1)/(detgeo.detrefpos(1));
     disp('applying projection magnification')
 end
 
diff --git a/zUtil_Indexter/gtINDEXFsimAllGrains.m b/zUtil_Indexter/gtINDEXFsimAllGrains.m
index b3ec142fbb019e45d810a3d74d9d320368b0073b..d7c4bfd84be9fff26eda5ddc4f534ef24d73855a 100644
--- a/zUtil_Indexter/gtINDEXFsimAllGrains.m
+++ b/zUtil_Indexter/gtINDEXFsimAllGrains.m
@@ -15,14 +15,22 @@ function gr = gtINDEXFsimAllGrains(phaseid, fsimID, saveFlag)
     load(fullfile('4_grains', sprintf('phase_%02d', phaseid), 'index.mat'), 'grain');
 
     parameters = gtLoadParameters();
+    labgeo = parameters.labgeo;
+    if (~isfield(labgeo, 'omstep') || isempty(labgeo.omstep))
+        labgeo.omstep = gtGetOmegaStepDeg(parameters);
+    end
+    if (~isfield(parameters, 'detgeo') || isempty(parameters.detgeo))
+        parameters.detgeo = gtGeoConvertLegacyLabgeo2Detgeo(labgeo);
+    end
+    detgeo = parameters.detgeo;
 
     grain = gtIndexUpdateGrains(grain, parameters.cryst(phaseid).latticepar);
     % add pairs info
     grain = gtIndexAddInfo(grain, parameters.acq, parameters.cryst(phaseid), phaseid);
     % fsim dct indexed spots
-    grain = gtIndexFwdSimPairs(grain, parameters.labgeo, ...
+    grain = gtIndexFwdSimPairs(grain, detgeo, labgeo, ...
         parameters.samgeo, [], parameters.cryst(phaseid).latticepar, ...
-        gtGetOmegaStepDeg(parameters), parameters.acq.energy);
+        parameters.acq.energy);
 
     % load DB and update dct grains
     difspotIDs = gtIndexAllGrainValues(grain, 'difspotID', [], [], [], false);
diff --git a/zUtil_Indexter/gtINDEXMatchGrains_2.m b/zUtil_Indexter/gtINDEXMatchGrains_2.m
index f14678987fa77afd78b683c2dc8ef92af1f8d092..0d39bb682786bef554980d05c3ed24d27c10bad0 100644
--- a/zUtil_Indexter/gtINDEXMatchGrains_2.m
+++ b/zUtil_Indexter/gtINDEXMatchGrains_2.m
@@ -130,10 +130,18 @@ nof_grains2 = length(grain2);
 
 out.odisp('Loading crystallographic info from parameters.mat ...')
 parameters = gtLoadParameters();
+labgeo = parameters.labgeo;
+if (~isfield(labgeo, 'omstep') || isempty(labgeo.omstep))
+    labgeo.omstep = gtGetOmegaStepDeg(parameters);
+end
+if (~isfield(parameters, 'detgeo') || isempty(parameters.detgeo))
+    parameters.detgeo = gtGeoConvertLegacyLabgeo2Detgeo(labgeo);
+end
+detgeo = parameters.detgeo;
 cryst = parameters.cryst(ph1);
 
-if ~exist('px1','var') || isempty(px1)
-    px1 = (parameters.labgeo.pixelsizeu+parameters.labgeo.pixelsizev)/2;
+if (~exist('px1', 'var') || isempty(px1))
+    px1 = (detgeo.pixelsizeu + detgeo.pixelsizev) / 2;
     disp(['pixelsize loaded from parameters: ' num2str(px1)])
 end
 if (~exist('px2', 'var') || isempty(px2))
diff --git a/zUtil_Indexter/gtIndexAddInfo.m b/zUtil_Indexter/gtIndexAddInfo.m
index 61bf573b7c70f7e138a424c38772d6a17053b684..06c73c62254aaa93e02e30d624fa9b6f32489686 100644
--- a/zUtil_Indexter/gtIndexAddInfo.m
+++ b/zUtil_Indexter/gtIndexAddInfo.m
@@ -152,19 +152,22 @@ grain = gtIndexUniqueReflections(grain, parameters.cryst(phaseid).dspacing,...
 
 % Forward simulation
 disp('Forward simulating reflections...')
-omstep = 180/parameters.acq.nproj;
 labgeo = parameters.labgeo;
-labgeo.detorig = gtGeoDetOrig(labgeo);
-labgeo.detnorm = cross(labgeo.detdiru, labgeo.detdirv);
+if (~isfield(labgeo, 'omstep') || isempty(labgeo.omstep))
+    labgeo.omstep = gtGetOmegaStepDeg(parameters);
+end
+if (~isfield(parameters, 'detgeo') || isempty(parameters.detgeo))
+    parameters.detgeo = gtGeoConvertLegacyLabgeo2Detgeo(labgeo);
+end
+detgeo = parameters.detgeo;
 
 % Fwd simulate pairs
-grain = gtIndexFwdSimPairs(grain, labgeo, parameters.samgeo,...
-        [], parameters.cryst(phaseid).latticepar, omstep, parameters.acq.energy);
+grain = gtIndexFwdSimPairs(grain, detgeo, labgeo, parameters.samgeo,...
+        [], parameters.cryst(phaseid).latticepar, parameters.acq.energy);
       
 % Fwd simulate individual spots    
 grain = gtFitFwdSimGrains(grain, parameters, 0, 0, []);
 
 toc
 
-
-end % of function
\ No newline at end of file
+end
\ No newline at end of file
diff --git a/zUtil_Indexter/gtIndexApplyDrifts.m b/zUtil_Indexter/gtIndexApplyDrifts.m
index 6adee79a1dcdd8d8e55346ed182f0a9eb2aca68b..9424cc850a8a7c8a7a4da590bd8dd94bd345b4f4 100644
--- a/zUtil_Indexter/gtIndexApplyDrifts.m
+++ b/zUtil_Indexter/gtIndexApplyDrifts.m
@@ -1,5 +1,5 @@
-function fsim = gtIndexApplyDrifts(pl, dsp, gc, om, drift, labgeo, ...
-                samgeo, rotcomp, Qdet, omstep, energy)
+function fsim = gtIndexApplyDrifts(pl, dsp, gc, om, drift, detgeo, labgeo, ...
+                samgeo, rotcomp, energy)
 % pl    - plane normals before drifts
 % sinth - sin(theta) in deformed state
 % om    - omega-s to interpolate the rotational drift
@@ -37,7 +37,8 @@ function fsim = gtIndexApplyDrifts(pl, dsp, gc, om, drift, labgeo, ...
     
     % Get fwd simulated data for the observed reflections
     fsim = gtIndexFwdSim(om, gcd, ...
-           sinthd, om4, pllab4, rot4, omind4, labgeo, samgeo, Qdet, omstep);
+           sinthd, om4, pllab4, rot4, omind4, ...
+           detgeo, labgeo, samgeo);
     
 end
 
diff --git a/zUtil_Indexter/gtIndexDifSpots_2.m b/zUtil_Indexter/gtIndexDifSpots_2.m
index 467ce8940a921f42a8efb4f4dd4d591b53532bcc..05c942214b8af85560b16ebf484b59158052a06c 100644
--- a/zUtil_Indexter/gtIndexDifSpots_2.m
+++ b/zUtil_Indexter/gtIndexDifSpots_2.m
@@ -30,145 +30,153 @@ function [gid, conf, sp] = gtIndexDifSpots_2(grain, sp, parameters, fsimID, lim)
 %       conf       = <double>  list of conflicts
 %       sp         = <struct>  updated sp info (like difspot table)
 
+    %%%%%%%%%%%%%%%%%%
+    % Prepare input
+    %%%%%%%%%%%%%%%%%%
 
-%%%%%%%%%%%%%%%%%%
-% Prepare input
-%%%%%%%%%%%%%%%%%%
-
-if ~exist('grain','var') || isempty(grain)
-    gtError('gtIndexDifspots:missingArgument','Argument ''grain'' must be given...Quitting')
-end
-if ~exist('sp','var') || isempty(sp)
-    gtError('gtIndexDifspots:missingArgument','Argument ''sp'' must be given...Quitting')
-end
-if ~exist('parameters','var') || isempty(parameters)
-    gtError('gtIndexDifspots:missingArgument','Argument ''parameters'' must be given...Quitting')
-end
-if ~exist('lim','var') || isempty(lim)
-    gtError('gtIndexDifspots:missingArgument','Argument ''lim'' must be given...Quitting')
-end
-
-if ~exist('fsimID','var') || isempty(fsimID)
-    fsimID = 1;
-end
+    if (~exist('grain', 'var') || isempty(grain))
+        gtError('gtIndexDifspots:missingArgument', ...
+            'Argument ''grain'' must be given...Quitting')
+    end
+    if (~exist('sp', 'var') || isempty(sp))
+        gtError('gtIndexDifspots:missingArgument', ...
+            'Argument ''sp'' must be given...Quitting')
+    end
+    if (~exist('parameters', 'var') || isempty(parameters))
+        gtError('gtIndexDifspots:missingArgument', ...
+            'Argument ''parameters'' must be given...Quitting')
+    end
+    if (~exist('lim', 'var') || isempty(lim))
+        gtError('gtIndexDifspots:missingArgument', ...
+            'Argument ''lim'' must be given...Quitting')
+    end
 
-omstep = 180/parameters.acq.nproj;
+    if (~exist('fsimID', 'var') || isempty(fsimID))
+        fsimID = 1;
+    end
 
-% Get spot and pair data
-if all(isfield(sp, {'difspotID','CentroidX','CentroidY','CentroidImage'}))
+    labgeo = parameters.labgeo;
+    if (~isfield(labgeo, 'omstep') || isempty(labgeo.omstep))
+        labgeo.omstep = gtGetOmegaStepDeg(parameters);
+    end
+    if (~isfield(parameters, 'detgeo') || isempty(parameters.detgeo))
+        parameters.detgeo = gtGeoConvertLegacyLabgeo2Detgeo(labgeo);
+    end
+    detgeo = parameters.detgeo;
 
-    sp.uvw  = [sp.CentroidX sp.CentroidY sp.CentroidImage];
-    sp.om   = sp.CentroidImage * omstep;
-    sp.id   = sp.difspotID;
+    % Get spot and pair data
+    if all(isfield(sp, {'difspotID', 'CentroidX', 'CentroidY', 'CentroidImage'}))
 
-    pairsinfo = true;
-elseif all(isfield(sp, {'difspotidA','difspotidB','omega','omegaB',...
-        'eta','etaB','theta','grainID',...
-        'samcentXA','samcentYA','samcentZA',...
-        'samcentXB','samcentYB','samcentZB'}))
+        sp.uvw  = [sp.CentroidX sp.CentroidY sp.CentroidImage];
+        sp.om   = sp.CentroidImage * labgeo.omstep;
+        sp.id   = sp.difspotID;
 
-    sp.eta   = [sp.eta, sp.etaB];
-    sp.id    = [sp.difspotidA, sp.difspotidB];
-    sp.om    = [sp.omega, sp.omegaB];
-    sp.uvw(:,3) = sp.om / omstep;
+        pairsinfo = true;
+    elseif all(isfield(sp, {'difspotidA', 'difspotidB', 'omega', 'omegaB', ...
+            'eta', 'etaB', 'theta', 'grainID', ...
+            'samcentXA', 'samcentYA', 'samcentZA', ...
+            'samcentXB', 'samcentYB', 'samcentZB'}))
 
-    xyzA  = [sp.samcentXA sp.samcentYA sp.samcentZA];
-    xyzB  = [sp.samcentXB sp.samcentYB sp.samcentZB];
-    xyz   = [xyzA; xyzB];
+        sp.eta   = [sp.eta, sp.etaB];
+        sp.id    = [sp.difspotidA, sp.difspotidB];
+        sp.om    = [sp.omega, sp.omegaB];
+        sp.uvw(:, 3) = sp.om / labgeo.omstep;
 
-    uvwLab = gtGeoSam2Lab(xyz, sp.om, parameters.labgeo, parameters.samgeo, 0, false);
-    sp.uvw(:,1:2) = gtGeoLab2Det(uvwLab, parameters.labgeo, 0);
+        xyzA  = [sp.samcentXA sp.samcentYA sp.samcentZA];
+        xyzB  = [sp.samcentXB sp.samcentYB sp.samcentZB];
+        xyz   = [xyzA; xyzB];
 
-    pairsinfo = true;
-% check if it is by chance a taper scan, so then it contains those fields
-elseif all(isfield(sp, {'spot3d_id','fc','sc','omega'}) & ~isfield(sp, 'cu'))
+        uvwLab = gtGeoSam2Lab(xyz, sp.om, labgeo, parameters.samgeo, 0, false);
+        sp.uvw(:, 1:2) = gtGeoLab2Det(uvwLab, detgeo, 0);
 
-    sp.uvw  = [sp.fc sp.sc sp.omega / omstep];
-    sp.om   = sp.omega;
-    sp.id   = sp.spot3d_id;
+        pairsinfo = true;
+    % check if it is by chance a taper scan, so then it contains those fields
+    elseif all(isfield(sp, {'spot3d_id', 'fc', 'sc', 'omega'}) & ~isfield(sp, 'cu'))
 
-    pairsinfo = false;
-else
-    gtError('gtIndexDifspots:wrong_input','Wrong input for ''sp''... Quitting!')
-end
+        sp.uvw  = [sp.fc sp.sc sp.omega / labgeo.omstep];
+        sp.om   = sp.omega;
+        sp.id   = sp.spot3d_id;
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% Vectorise fwd simulated grain spots
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-uvw = gtIndexAllGrainValues(grain, 'fsim', 'uvw', fsimID, [], false);
-uvw = [uvw{:}];
-gr.u = uvw(1,:);
-gr.v = uvw(2,:);
-gr.w = uvw(3,:);
-
-gr.om = gtIndexAllGrainValues(grain, 'fsim', 'om', fsimID, [], false);
-gr.om = [gr.om{:}];
-gr.id =[];
-for ii=1:length(grain)
-    gr.id = [gr.id; grain{ii}.id(ones(1,length(grain{ii}.fsim(fsimID).om)))'];
-end
-gr.id = gr.id';
+        pairsinfo = false;
+    else
+        gtError('gtIndexDifspots:wrong_input', 'Wrong input for ''sp''... Quitting!')
+    end
 
-disp(['current pixelsizeu: ' num2str(parameters.labgeo.pixelsizeu)])
-disp(['current pixelsizev: ' num2str(parameters.labgeo.pixelsizev)])
-disp(['omega step:         ' num2str(omstep)])
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    % Vectorise fwd simulated grain spots
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+    uvw = gtIndexAllGrainValues(grain, 'fsim', 'uvw', fsimID, [], false);
+    uvw = [uvw{:}];
+    gr.u = uvw(1, :);
+    gr.v = uvw(2, :);
+    gr.w = uvw(3, :);
+
+    gr.om = gtIndexAllGrainValues(grain, 'fsim', 'om', fsimID, [], false);
+    gr.om = [gr.om{:}];
+    gr.id =[];
+    for ii=1:length(grain)
+        gr.id = [gr.id; grain{ii}.id(ones(1, length(grain{ii}.fsim(fsimID).om)))'];
+    end
+    gr.id = gr.id';
 
-% pixelsize of fsim detector
-if (pairsinfo)
-    lim.u = lim.u / parameters.labgeo.pixelsizeu; % pixels
-    lim.v = lim.v / parameters.labgeo.pixelsizev; % pixels
-    lim.w = lim.w / omstep; % image number
-end
+    disp(['current pixelsizeu: ' num2str(detgeo.pixelsizeu)])
+    disp(['current pixelsizev: ' num2str(detgeo.pixelsizev)])
+    disp(['omega step:         ' num2str(labgeo.omstep)])
 
-print_structure(lim , 'lim')
+    % pixelsize of fsim detector
+    if (pairsinfo)
+        lim.u = lim.u / detgeo.pixelsizeu; % pixels
+        lim.v = lim.v / detgeo.pixelsizev; % pixels
+        lim.w = lim.w / labgeo.omstep; % image number
+    end
 
+    print_structure(lim, 'lim')
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%
-% Loop to find grain ID-s
-%%%%%%%%%%%%%%%%%%%%%%%%%%
+    %%%%%%%%%%%%%%%%%%%%%%%%%%
+    % Loop to find grain ID-s
+    %%%%%%%%%%%%%%%%%%%%%%%%%%
 
-gid  = zeros(1,length(sp.uvw));
-conf = zeros(1,length(sp.uvw));
-ok   = false(length(gr.id),3);
+    gid  = zeros(1, length(sp.uvw));
+    conf = zeros(1, length(sp.uvw));
+    ok   = false(length(gr.id), 3);
 
-gauge = GtGauge(length(sp.uvw), 'Processing spots: ');
+    gauge = GtGauge(length(sp.uvw), 'Processing spots: ');
 
-%tic
+    %tic
 
-% Loop through spots - could be a parallel parfor loop
-for ii = 1:length(sp.uvw)
-    gauge.incrementAndDisplay(100);
+    % Loop through spots - could be a parallel parfor loop
+    for ii = 1:length(sp.uvw)
+        gauge.incrementAndDisplay(100);
 
-    % Check if position is within limits
-    ok(:,1) = abs(gr.u - sp.uvw(ii,1)) < lim.u;
-    ok(:,2) = abs(gr.v - sp.uvw(ii,2)) < lim.v;
-    ok(:,3) = abs(gr.w - sp.uvw(ii,3)) < lim.w;
+        % Check if position is within limits
+        ok(:, 1) = abs(gr.u - sp.uvw(ii, 1)) < lim.u;
+        ok(:, 2) = abs(gr.v - sp.uvw(ii, 2)) < lim.v;
+        ok(:, 3) = abs(gr.w - sp.uvw(ii, 3)) < lim.w;
 
-    oks = all(ok, 2);
-    sumoks = sum(oks);
+        oks = all(ok, 2);
+        sumoks = sum(oks);
 
-    % Check metadata, if exactly one candidate grain found
-    if sumoks == 1
+        % Check metadata, if exactly one candidate grain found
+        if (sumoks == 1)
 
-        idg = gr.id(find(oks,1,'first'));
-        ids = gtIndexAllGrainValues(grain, 'id', [], [], []);
-        ind = find(ids == idg);
+            idg = gr.id(find(oks, 1, 'first'));
+            ids = gtIndexAllGrainValues(grain, 'id', [], [], []);
+            ind = find(ids == idg);
 
-        gid(ii) = idg;
+            gid(ii) = idg;
 
-    elseif sumoks > 1
-        % spot in conflict
-        conf(ii) = sum(oks);
+        elseif (sumoks > 1)
+            % spot in conflict
+            conf(ii) = sum(oks);
+        end
     end
-end
 
-gauge.delete();
+    gauge.delete();
 
-%toc;
-disp(' ')
-fprintf('No. of spots analysed   : %d\n', length(gid));
-fprintf('No. of spots indexed    : %d\n', sum(gid>0));
-fprintf('No. of spots in conflict: %d\n', sum(conf));
-
-end % end of function
+    %toc;
+    disp(' ')
+    fprintf('No. of spots analysed   : %d\n', length(gid));
+    fprintf('No. of spots indexed    : %d\n', sum(gid>0));
+    fprintf('No. of spots in conflict: %d\n', sum(conf));
+end
diff --git a/zUtil_Indexter/gtIndexFwdSim.m b/zUtil_Indexter/gtIndexFwdSim.m
index 8f67fbf5705ec5ef5b0e215d96a9b8e9f55786bd..4a35a068b586706cdf7f633a22791dff7dbb9375 100644
--- a/zUtil_Indexter/gtIndexFwdSim.m
+++ b/zUtil_Indexter/gtIndexFwdSim.m
@@ -1,18 +1,18 @@
-function fwdSim = gtIndexFwdSim(om_mes, gcsam, sinth, om4, pllab4, rot4, ...
-                omind4, labgeo, samgeo, Qdet, omstep)
+function fwdSim = gtIndexFwdSim(om_mes, gcsam, ~, om4, pllab4, rot4, ...
+                omind4, detgeo, labgeo, samgeo)
 % fwdSim = gtIndexFwdSim(om_mes, gcsam, sinth, om4, pllab4, rot4, ...
-%          omind4, labgeo, samgeo, Qdet, omstep)
+%          omind4, detgeo, labgeo, samgeo, omstep)
 % -------------------------------------------------------------------
 % Provides the fwd simulation data for a given omega index set.
 % /mntdirect/_data_id19_graintracking/DCT_Analysis/RD4_15N_taper_/workspace/grid/
-% peaks_t100_s_cleaned_smalls.flt.new 
+% peaks_t100_s_cleaned_smalls.flt.new
 
 % Absolute differences
-dom = om4 - om_mes([1 1 1 1],:);
+dom = om4 - om_mes([1 1 1 1], :);
 
-% Find the indices with the smallest absolute difference; consider 
+% Find the indices with the smallest absolute difference; consider
 % errors and periodicity of 360deg in omega
-[~, ind] = min( abs( [dom; dom - 360; dom + 360] ) , [], 1);
+[~, ind] = min( abs( [dom; dom - 360; dom + 360] ), [], 1);
 
 % Predicted omega is 360deg more
 omhigh = ( (5 <= ind) & (ind <= 8) );
@@ -26,33 +26,33 @@ ind(omlow)  = ind(omlow)  - 8;
 
 % Get final omega angles (can be <0) using the corresponding linear
 % indices of omind4
-linind            = ind + (0:length(ind)-1)*4;
+linind            = ind + (0:length(ind)-1) * 4;
 fwdSim.om         = om4(linind);
-fwdSim.om(omhigh) = fwdSim.om(omhigh) - 360; 
-fwdSim.om(omlow)  = fwdSim.om(omlow)  + 360; 
+fwdSim.om(omhigh) = fwdSim.om(omhigh) - 360;
+fwdSim.om(omlow)  = fwdSim.om(omlow)  + 360;
 
-% Omega index 
+% Omega index
 fwdSim.omind = omind4(linind);
 
 % Omegas in a column vector
-ind1 = (ind==1);
-ind2 = (ind==2);
-ind3 = (ind==3);
-ind4 = (ind==4);
+ind1 = (ind == 1);
+ind2 = (ind == 2);
+ind3 = (ind == 3);
+ind4 = (ind == 4);
 
 % Extract the correct pllab vectors according to omind
 fwdSim.pllab         = zeros(3, length(ind));
-fwdSim.pllab(:,ind1) = pllab4(:,ind1,1);
-fwdSim.pllab(:,ind2) = pllab4(:,ind2,2);
-fwdSim.pllab(:,ind3) = pllab4(:,ind3,3);
-fwdSim.pllab(:,ind4) = pllab4(:,ind4,4);
+fwdSim.pllab(:, ind1) = pllab4(:, ind1, 1);
+fwdSim.pllab(:, ind2) = pllab4(:, ind2, 2);
+fwdSim.pllab(:, ind3) = pllab4(:, ind3, 3);
+fwdSim.pllab(:, ind4) = pllab4(:, ind4, 4);
 
 % Extract the correct rotation matrices according to omind
 fwdSim.rot           = zeros(3, 3, length(ind));
-fwdSim.rot(:,:,ind1) = rot4(:,:,ind1,1);
-fwdSim.rot(:,:,ind2) = rot4(:,:,ind2,2);
-fwdSim.rot(:,:,ind3) = rot4(:,:,ind3,3);
-fwdSim.rot(:,:,ind4) = rot4(:,:,ind4,4);
+fwdSim.rot(:, :, ind1) = rot4(:, :, ind1, 1);
+fwdSim.rot(:, :, ind2) = rot4(:, :, ind2, 2);
+fwdSim.rot(:, :, ind3) = rot4(:, :, ind3, 3);
+fwdSim.rot(:, :, ind4) = rot4(:, :, ind4, 4);
 
 
 % Diffraction vector
@@ -60,25 +60,25 @@ fwdSim.dveclab = gtFedPredictDiffVecMultiple(fwdSim.pllab, labgeo.beamdir');
 fwdSim.dvecsam = gtGeoLab2Sam(fwdSim.dveclab', fwdSim.om, labgeo, samgeo, 1, 1)';
 
 % Make sure grain center is extended
-gcsam = gcsam(:,ones(size(fwdSim.dveclab,2),1));
+gcsam = gcsam(:, ones(size(fwdSim.dveclab, 2), 1));
 
-% Absolute detector coordinates U,V in pixels 
-uv = gtFedPredictUVMultiple(fwdSim.rot, fwdSim.dveclab, gcsam, labgeo.detrefpos', ...
-     labgeo.detnorm', Qdet, [labgeo.detrefu; labgeo.detrefv]);
+% Absolute detector coordinates U, V in pixels
+uv = gtFedPredictUVMultiple(fwdSim.rot, fwdSim.dveclab, gcsam, detgeo.detrefpos', ...
+     detgeo.detnorm', detgeo.Qdet, [detgeo.detrefu; detgeo.detrefv]);
 
-is_nan = isnan(uv(1,:));
+is_nan = isnan(uv(1, :));
 
-keep = (uv(1,:) >= 1) & (labgeo.detsizeu >= uv(1,:)) & ...
-       (uv(2,:) >= 1) & (labgeo.detsizev >= uv(2,:));
-% Absolute coordinates U,V,W (pixel,pixel,image); to account for 
+keep = (uv(1, :) >= 1) & (detgeo.detsizeu >= uv(1, :)) & ...
+       (uv(2, :) >= 1) & (detgeo.detsizev >= uv(2, :));
+% Absolute coordinates U, V, W (pixel, pixel, image); to account for
 % measurement errors, do not correct for the image number
-% to be mod(om,360)/omstep
+% to be mod(om, 360)/omstep
 index_original = ismember(om_mes, fwdSim.om);
 
-fwdSim.uvw        = [uv; fwdSim.om/omstep];
-fwdSim.uvw_nan    = [uv(:,is_nan); fwdSim.om(is_nan)/omstep];
+fwdSim.uvw        = [uv; fwdSim.om/labgeo.omstep];
+fwdSim.uvw_nan    = [uv(:, is_nan); fwdSim.om(is_nan)/labgeo.omstep];
 fwdSim.keep       = keep;
 fwdSim.is_nan     = is_nan;
-fwdSim.image      = mod(fwdSim.om,360)/omstep;
+fwdSim.image      = mod(fwdSim.om, 360)/labgeo.omstep;
 
 end % end of function
diff --git a/zUtil_Indexter/gtIndexFwdSimGrains.m b/zUtil_Indexter/gtIndexFwdSimGrains.m
index f79c9255df6460d2c7e89e7ba352c31bafd0c287..5b29a3b346f229f9092d99eed0237c5ec77fb240 100644
--- a/zUtil_Indexter/gtIndexFwdSimGrains.m
+++ b/zUtil_Indexter/gtIndexFwdSimGrains.m
@@ -74,12 +74,13 @@ else
 end
 
 labgeo = parameters.labgeo;
-if ~isfield(labgeo,'detorig')
-    labgeo.detorig = gtGeoDetOrig(labgeo);
+if (~isfield(labgeo, 'omstep') || isempty(labgeo.omstep))
+    labgeo.omstep = gtGetOmegaStepDeg(parameters);
 end
-if ~isfield(labgeo,'detnorm')
-    labgeo.detnorm = cross(labgeo.detdiru,labgeo.detidirv);
+if (~isfield(parameters, 'detgeo') || isempty(parameters.detgeo))
+    parameters.detgeo = gtGeoConvertLegacyLabgeo2Detgeo(labgeo);
 end
+detgeo = parameters.detgeo;
 
 energy = parameters.acq.energy;
 samgeo = parameters.samgeo;
@@ -104,14 +105,8 @@ end
 
 cryst.usedfamsp = ismember(cryst.thetatypesp, cryst.used_thetatype);
 
-% Scaling for Lab -> Detector transformation (mm -> pixels)
-detscaleu = 1/labgeo.pixelsizeu;
-detscalev = 1/labgeo.pixelsizev;
 % Rotation tensor components
 rotcomp = gtMathsRotationMatrixComp(labgeo.rotdir', 'col');
-% Detector projection matrix
-Qdet = gtFedDetectorProjectionTensor(labgeo.detdiru',labgeo.detdirv',...
-                                     detscaleu,detscalev);
 % B matrix for transformation from (hkl) to Cartesian Crystal reference
 Bmat = gtCrystHKL2CartesianMatrix(cryst.latticepar);
 
@@ -129,7 +124,7 @@ Rvec    = gtIndexAllGrainValues(grain, 'R_vector', [], 1, 1:3);
 cry2sam = gtFedRod2RotTensor(Rvec');
 
 % calc radius of u, v coordinate
-getradius = @(u,v) sqrt((u - labgeo.detrefu).^2 + (v - labgeo.detrefv).^2);
+getradius = @(u, v) sqrt((u - detgeo.detrefu).^2 + (v - detgeo.detrefv).^2);
 
 % loop through grains
 for ii = 1:length(grain)
diff --git a/zUtil_Indexter/gtIndexFwdSimLimits.m b/zUtil_Indexter/gtIndexFwdSimLimits.m
index 6e0364e3a00d916c74f09745b7b6b96aa937febd..efdc13ecb0d11bea5e3fddb7d8a16f05af827d61 100644
--- a/zUtil_Indexter/gtIndexFwdSimLimits.m
+++ b/zUtil_Indexter/gtIndexFwdSimLimits.m
@@ -15,56 +15,62 @@ function lim = gtIndexFwdSimLimits(grain, parameters, pairinfo)
 %          .v = <double>     Limits in mm for V coordinate
 %          .w = <double>     Limits in degrees for omega coordinate
 
-omstep = 180/parameters.acq.nproj;
+    labgeo = parameters.labgeo;
+    if (~isfield(labgeo, 'omstep') || isempty(labgeo.omstep))
+        labgeo.omstep = gtGetOmegaStepDeg(parameters);
+    end
+    if (~isfield(parameters, 'detgeo') || isempty(parameters.detgeo))
+        parameters.detgeo = gtGeoConvertLegacyLabgeo2Detgeo(labgeo);
+    end
+    detgeo = parameters.detgeo;
 
-if pairinfo
-    uvwA_sim = gtIndexAllGrainValues(grain,'fsimA','uvw',[],[],false); uvwA_sim = [uvwA_sim{:}];
-    uvwB_sim = gtIndexAllGrainValues(grain,'fsimB','uvw',[],[],false); uvwB_sim = [uvwB_sim{:}];
-    uvwA_exp = gtIndexAllGrainValues(grain,'uvwA',[],[],[],false); uvwA_exp = [uvwA_exp{:}];
-    uvwB_exp = gtIndexAllGrainValues(grain,'uvwB',[],[],[],false); uvwB_exp = [uvwB_exp{:}];
-    
-    uvw_sim = [uvwA_sim, uvwB_sim];
-    uvw_exp = [uvwA_exp, uvwB_exp];
-    
-     % Deviations between fwd simulation and measured spot positions
-    du = uvw_sim(1,:) - uvw_exp(1,:);
-    dv = uvw_sim(2,:) - uvw_exp(2,:);
-    dw = uvw_sim(3,:) - uvw_exp(3,:);
+    if pairinfo
+        uvwA_sim = gtIndexAllGrainValues(grain, 'fsimA', 'uvw', [], [], false); uvwA_sim = [uvwA_sim{:}];
+        uvwB_sim = gtIndexAllGrainValues(grain, 'fsimB', 'uvw', [], [], false); uvwB_sim = [uvwB_sim{:}];
+        uvwA_exp = gtIndexAllGrainValues(grain, 'uvwA', [], [], [], false); uvwA_exp = [uvwA_exp{:}];
+        uvwB_exp = gtIndexAllGrainValues(grain, 'uvwB', [], [], [], false); uvwB_exp = [uvwB_exp{:}];
 
-    lim.u = mean(du) + 3*std(du); % pixel
-    lim.v = mean(dv) + 3*std(dv); % pixel
-    lim.w = mean(dw) + 3*std(dw); % images
-    
-    % convert to mm and degrees
-    lim.u = lim.u*parameters.labgeo.pixelsizeu; % mm
-    lim.v = lim.v*parameters.labgeo.pixelsizev; % mm
-    lim.w = lim.w*omstep; % omega [degrees]
-else
-    % fsim only indexed peaks
+        uvw_sim = [uvwA_sim, uvwB_sim];
+        uvw_exp = [uvwA_exp, uvwB_exp];
 
-%     uvw = gtIndexAllGrainValues(grain,'fsimU','uvw',[],[],false);  uvw = [uvw{:}];
-%     ind = find(isnan(uvw(3,:)));
-%     % experimental positions
-%     fc = gtIndexAllGrainValues(grain,'fc',[],[],[],false);        fc = [fc{:}];
-%     sc = gtIndexAllGrainValues(grain,'sc',[],[],[],false);        sc = [sc{:}];
-%     omega = gtIndexAllGrainValues(grain,'omega',[],[],[],false);  omega = [omega{:}];
-% 
-%     uvw(:,ind) = [];
-% 
-%     fc(ind) = [];
-%     sc(ind) = [];
-%     omega(ind) = [];
-%     
-%     cent = [fc; sc];
-%     centim = omega/omstep;
-    
-    lim.u = 2; % pixel
-    lim.v = 2; % pixel
-    lim.w = 5; % image number
-end
+         % Deviations between fwd simulation and measured spot positions
+        du = uvw_sim(1, :) - uvw_exp(1, :);
+        dv = uvw_sim(2, :) - uvw_exp(2, :);
+        dw = uvw_sim(3, :) - uvw_exp(3, :);
+
+        lim.u = mean(du) + 3*std(du); % pixel
+        lim.v = mean(dv) + 3*std(dv); % pixel
+        lim.w = mean(dw) + 3*std(dw); % images
 
-disp(['current pixelsizeu: ' num2str(parameters.labgeo.pixelsizeu)])
-disp(['current pixelsizev: ' num2str(parameters.labgeo.pixelsizev)])
-disp(['omega step:         ' num2str(omstep)])
+        % convert to mm and degrees
+        lim.u = lim.u * detgeo.pixelsizeu; % mm
+        lim.v = lim.v * detgeo.pixelsizev; % mm
+        lim.w = lim.w * labgeo.omstep; % omega [degrees]
+    else
+        % fsim only indexed peaks
 
-end % end of function
+%         uvw = gtIndexAllGrainValues(grain, 'fsimU', 'uvw', [], [], false);  uvw = [uvw{:}];
+%         ind = find(isnan(uvw(3, :)));
+%         % experimental positions
+%         fc = gtIndexAllGrainValues(grain, 'fc', [], [], [], false);        fc = [fc{:}];
+%         sc = gtIndexAllGrainValues(grain, 'sc', [], [], [], false);        sc = [sc{:}];
+%         omega = gtIndexAllGrainValues(grain, 'omega', [], [], [], false);  omega = [omega{:}];
+%
+%         uvw(:, ind) = [];
+%
+%         fc(ind) = [];
+%         sc(ind) = [];
+%         omega(ind) = [];
+%         
+%         cent = [fc; sc];
+%         centim = omega/labgeo.omstep;
+
+        lim.u = 2; % pixel
+        lim.v = 2; % pixel
+        lim.w = 5; % image number
+    end
+
+    disp(['current pixelsizeu: ' num2str(detgeo.pixelsizeu)])
+    disp(['current pixelsizev: ' num2str(detgeo.pixelsizev)])
+    disp(['omega step:         ' num2str(labgeo.omstep)])
+end
diff --git a/zUtil_Indexter/gtIndexFwdSimPairs.m b/zUtil_Indexter/gtIndexFwdSimPairs.m
index 371d99f80cbdb6f9eb28bae6a958173be70c7002..288fc7c32ef426677139e100f2673036b4dab128 100644
--- a/zUtil_Indexter/gtIndexFwdSimPairs.m
+++ b/zUtil_Indexter/gtIndexFwdSimPairs.m
@@ -1,18 +1,18 @@
-function grain = gtIndexFwdSimPairs(grain, labgeo, samgeo, drift, ...
-                 latticepar, omstep, energy)
-% GTINDEXFWDSIMPAIRS  Forward simulates spot positions on the detector 
+function grain = gtIndexFwdSimPairs(grain, detgeo, labgeo, samgeo, drift, ...
+                 latticepar, energy)
+% GTINDEXFWDSIMPAIRS  Forward simulates spot positions on the detector
 %                     for all spots indexed as Friedel pairs.
-% 
-%     grain = gtIndexFwdSimPairs(grain, labgeo, samgeo, drift, latticepar, omstep, energy)
+%
+%     grain = gtIndexFwdSimPairs(grain, detgeo, labgeo, samgeo, drift, latticepar, energy)
 %     ------------------------------------------------------------------------------------
 %
 %     INPUT:
 %       grain      = <cell>       all the grains as in index.mat
+%       detgeo     = <struct>     as in parameters.detgeo
 %       labgeo     = <struct>     as in parameters.labgeo
 %       samgeo     = <struct>     as in parameters.samgeo
 %       drift      = <double>     sample drifts; if empty, no drifts are considered
 %       latticepar = <double>     lattice parameters (1x6)
-%       omstep     = <double>     omega rotational stepping (degrees/image)
 %       energy     = <double>     beam energy in keV
 %
 %     OUTPUT:
@@ -28,23 +28,11 @@ function grain = gtIndexFwdSimPairs(grain, labgeo, samgeo, drift, ...
 % Prepare input
 %%%%%%%%%%%%%%%%%%%
 
-% B matrix for transformation from (hkl) to Cartesian Crystal reference 
+% B matrix for transformation from (hkl) to Cartesian Crystal reference
 Bmat = gtCrystHKL2CartesianMatrix(latticepar);
 
-% Scaling for Lab -> Detector transformation (mm -> pixels)
-detscaleu = 1/labgeo.pixelsizeu ;
-detscalev = 1/labgeo.pixelsizev ;
-
-
-labgeo.detorig = gtGeoDetOrig(labgeo);
-labgeo.detnorm = cross(labgeo.detdiru, labgeo.detdirv);
-
 % Rotation tensor components
-rotcomp = gtMathsRotationMatrixComp(labgeo.rotdir','col');
-
-% Detector projection matrix
-Qdet = gtFedDetectorProjectionTensor(labgeo.detdiru', labgeo.detdirv', ...
-                                     detscaleu, detscalev);
+rotcomp = gtMathsRotationMatrixComp(labgeo.rotdir', 'col');
 
 % Rotation tensor describing crystal orientation; calculate all grains
 % at once for speed-up.
@@ -101,15 +89,15 @@ for ii = 1:length(grain)
 
         if isfield(grain{ii}, 'omegaA')
             grain{ii}.fsimA = gtIndexApplyDrifts(pl_samd, dsp_d, grain{ii}.center',...
-                grain{ii}.omegaA, drift, labgeo, samgeo, rotcomp, Qdet, omstep, energy);
+                grain{ii}.omegaA, drift, detgeo, labgeo, samgeo, rotcomp,  energy);
         end
         if isfield(grain{ii}, 'omegaB')
             grain{ii}.fsimB = gtIndexApplyDrifts(pl_samd, dsp_d, grain{ii}.center',...
-                grain{ii}.omegaB, drift, labgeo, samgeo, rotcomp, Qdet, omstep, energy);
+                grain{ii}.omegaB, drift, detgeo, labgeo, samgeo, rotcomp, energy);
         end
         if isfield(grain{ii}, 'omega') && ~isfield(grain{ii}, 'omegaA') && ~isfield(grain{ii}, 'omegaB')
             grain{ii}.fsimU = gtIndexApplyDrifts(pl_samd, dsp_d, grain{ii}.center',...
-                grain{ii}.omega, drift, labgeo, samgeo, rotcomp, Qdet, omstep, energy);
+                grain{ii}.omega, drift, detgeo, labgeo, samgeo, rotcomp, energy);
         end
 
     else
@@ -146,15 +134,15 @@ for ii = 1:length(grain)
         % or the one corrected for the pair).
         if isfield(grain{ii}, 'omegaA')
             grain{ii}.fsimA = gtIndexFwdSim(grain{ii}.omegaA, grain{ii}.center', ...
-                sinth_d, om4, pl_lab4, rot4, omind4, labgeo, samgeo, Qdet, omstep);
+                sinth_d, om4, pl_lab4, rot4, omind4, detgeo, labgeo, samgeo);
         end
         if isfield(grain{ii}, 'omegaB')
             grain{ii}.fsimB = gtIndexFwdSim(grain{ii}.omegaB, grain{ii}.center', ...
-                sinth_d, om4, pl_lab4, rot4, omind4, labgeo, samgeo, Qdet, omstep);
+                sinth_d, om4, pl_lab4, rot4, omind4, detgeo, labgeo, samgeo);
         end
         if isfield(grain{ii}, 'omega') && ~isfield(grain{ii}, 'omegaA') && ~isfield(grain{ii}, 'omegaB')
             grain{ii}.fsimU = gtIndexFwdSim(grain{ii}.omega, grain{ii}.center', ...
-                sinth_d, om4, pl_lab4, rot4, omind4, labgeo, samgeo, Qdet, omstep);
+                sinth_d, om4, pl_lab4, rot4, omind4, detgeo, labgeo, samgeo);
         end
     end
 end