diff --git a/4_grains/gtCalculateGrain.m b/4_grains/gtCalculateGrain.m
index 71a71c2405d42f9a542313f274ed4b859c610bc2..f092911597aa2acf0bdf0050af43db0188ba0268 100644
--- a/4_grains/gtCalculateGrain.m
+++ b/4_grains/gtCalculateGrain.m
@@ -283,17 +283,14 @@ end
 % % Computing the remaining physical quantities
 
 % Diffraction vectors % takes coloumn vectors
-dveclab = gtFedPredictDiffVecMultiple(pllab, labgeo.beamdir');
-
-% Normalize vectors
-dvec_norms = sqrt(sum(dveclab .* dveclab, 1));
-dveclab = (dveclab ./ dvec_norms([1 1 1], :))';
+dveclab = gtFedPredictDiffVecMultiple(pllab, labgeo.beamdir')';
+dveclab = gtMathsNormalizeVectorsList(dveclab);
 
 % Eta angles % takes row vector
 eta = gtGeoEtaFromDiffVec(dveclab, labgeo);
 
 % Diffraction vector in Sample reference
-dvecsam = gtGeoLab2Sam(dveclab, om, labgeo, samgeo, true);
+dvecsam = gtGeoLab2Sam(dveclab, rot_l2s, labgeo, samgeo, true);
 
 % Let's convert grain center Sam -> Lab
 rot_s2l = permute(rot_l2s, [2 1 3]);
diff --git a/4_grains/gtCalculateGrain_p.m b/4_grains/gtCalculateGrain_p.m
index 785f6704a67270d862677d142d3b34bc77dc3951..5cbe8b94addd426da486df8413fce87edfa1c755 100644
--- a/4_grains/gtCalculateGrain_p.m
+++ b/4_grains/gtCalculateGrain_p.m
@@ -338,17 +338,14 @@ hklsp = permute(hklsp, [2 1 3]);
 % % Computing the remaining physical quantities
 
 % Diffraction vectors % takes coloumn vectors
-dveclab = gtFedPredictDiffVecMultiple(pllab, labgeo.beamdir');
-
-% Normalize vectors
-dvec_norms = sqrt(sum(dveclab .* dveclab, 1));
-dveclab = (dveclab ./ dvec_norms([1 1 1], :))';
+dveclab = gtFedPredictDiffVecMultiple(pllab, labgeo.beamdir')';
+dveclab = gtMathsNormalizeVectorsList(dveclab);
 
 % Eta angles % takes row vector
 eta = gtGeoEtaFromDiffVec(dveclab, labgeo);
 
 % Diffraction vector in Sample reference
-dvecsam = gtGeoLab2Sam(dveclab, om, labgeo, samgeo, true);
+dvecsam = gtGeoLab2Sam(dveclab, rot_l2s, labgeo, samgeo, true);
 
 % Let's convert grain center Sam -> Lab
 rot_s2l = permute(rot_l2s, [2 1 3]);
diff --git a/zUtil_Deformation/GtOrientationSampling.m b/zUtil_Deformation/GtOrientationSampling.m
index e8344cbad2d5bd166bf024ebecf92eb629e2c651..edcd263055d46b2b45be0378b0f940f8675d320c 100644
--- a/zUtil_Deformation/GtOrientationSampling.m
+++ b/zUtil_Deformation/GtOrientationSampling.m
@@ -814,14 +814,13 @@ classdef GtOrientationSampling < handle
                 num_chars_gr = fprintf('computing projection geometries..');
                 bbpos_oi = repmat(bbpos, [num_orients, 1]);
 
-                allblobs = cellfun(@(x){x.allblobs}, self.lattice(o_ii).gr);
-                allblobs = [allblobs{:}];
+                allblobs = cellfun(@(x)(x.allblobs), self.lattice(o_ii).gr);
                 dvecs = cat(1, allblobs(:).dvecsam);
-                omegas = cat(1, allblobs(:).omega);
+                rot_w_l2s = cat(3, allblobs(:).srot);
 
                 for n = 1:numel(self.parameters.detgeo)
                     proj_geom = gtGeoProjForReconstruction(...
-                        dvecs, omegas, grain_center, bbpos_oi, [], ...
+                        dvecs, rot_w_l2s, grain_center, bbpos_oi, [], ...
                         self.parameters.detgeo(n), ...
                         self.parameters.labgeo, ...
                         self.parameters.samgeo, ...
@@ -856,6 +855,9 @@ classdef GtOrientationSampling < handle
                 bbpos(b_ii, :) = [sbl.bbuim(1), sbl.bbvim(1), sbl.bbsize(1:2)];
             end
 
+            num_ss_orients = prod(self.ss_factor);
+            bbpos = repmat(bbpos, [num_ss_orients, 1]);
+
             ref_omind = self.ref_gr.allblobs.omind;
             ref_included = self.ondet(self.included);
 
@@ -872,18 +874,13 @@ classdef GtOrientationSampling < handle
                         self.lattice_ss(o_ii).gr{g_ii}, self.parameters, ...
                         'ref_omind', ref_omind, 'included', ref_included );
 
-                    num_ss_orients = numel(self.lattice_ss(o_ii).gr{g_ii});
-                    bbpos_oi = repmat(bbpos, [num_ss_orients, 1]);
-
-                    allblobs = cellfun(@(x){x.allblobs}, self.lattice_ss(o_ii).gr{g_ii});
-                    allblobs = [allblobs{:}];
+                    allblobs = cellfun(@(x)(x.allblobs), self.lattice_ss(o_ii).gr{g_ii});
                     dvecs = cat(1, allblobs(:).dvecsam);
-                    omegas = cat(1, allblobs(:).omega);
+                    rot_w_l2s = cat(3, allblobs(:).srot);
 
-                    num_ss_orient = numel(self.lattice_ss(o_ii).gr{g_ii});
                     for n = 1:numel(self.parameters.detgeo)
                         proj_geom = gtGeoProjForReconstruction(...
-                            dvecs, omegas, grain_center, bbpos_oi, [], ...
+                            dvecs, rot_w_l2s, grain_center, bbpos, [], ...
                             self.parameters.detgeo(n), ...
                             self.parameters.labgeo, ...
                             self.parameters.samgeo, ...
@@ -891,7 +888,7 @@ classdef GtOrientationSampling < handle
                             'ASTRA_grain');
 
                         count = 1;
-                        for ss_ii = 1:num_ss_orient
+                        for ss_ii = 1:num_ss_orients
                             num_vecs = size(self.lattice_ss(o_ii).gr{g_ii}{ss_ii}.allblobs.dvecsam, 1);
                             self.lattice_ss(o_ii).gr{g_ii}{ss_ii}.allblobs.detector(n).proj_geom ...
                                 = proj_geom(count:(count+num_vecs-1), :);
diff --git a/zUtil_Geo/gtGeoProjForReconstruction.m b/zUtil_Geo/gtGeoProjForReconstruction.m
index 1be1691721e77bdaae048537482d9c3b7119e071..9be68697d547fabfcb0af0843be3fdab0f148f39 100644
--- a/zUtil_Geo/gtGeoProjForReconstruction.m
+++ b/zUtil_Geo/gtGeoProjForReconstruction.m
@@ -53,12 +53,20 @@ function proj_geom = gtGeoProjForReconstruction(projvec, omega, graincenter,...
     end
 
     % Bounding box center U,V coordinates
-    bbcent_det  = [ bbpos(:, 1) + (bbpos(:, 3) - 1) / 2 + bboff(:, 1) , ...
-                    bbpos(:, 2) + (bbpos(:, 4) - 1) / 2 + bboff(:, 2) ];
+    bbcent_det  = [ ...
+        bbpos(:, 1) + (bbpos(:, 3) - 1) / 2 + bboff(:, 1) , ...
+        bbpos(:, 2) + (bbpos(:, 4) - 1) / 2 + bboff(:, 2) ];
 
-    % All functions that use the same omegas, need the same rotation tensors:
-    rot_comp_w = gtMathsRotationMatrixComp(labgeo.rotdir', 'col');
-    rot_omega = gtMathsRotationTensor(omega, rot_comp_w);
+    if (all([size(omega, 1), size(omega, 2)] == [3 3]))
+        % Rotation tensors were given already
+        rot_omega = omega;
+        ones_omega = ones(size(omega, 3), 1);
+    else
+        % All functions that use the same omegas, need the same rotation tensors:
+        rot_comp_w = gtMathsRotationMatrixComp(labgeo.rotdir', 'col');
+        rot_omega = gtMathsRotationTensor(omega, rot_comp_w);
+        ones_omega = ones(length(omega), 1);
+    end
 
     % Bounding box center X,Y,Z coordinates in LAB and REC reference
     bbcent_lab  = gtGeoDet2Lab(bbcent_det, detgeo, 0);
@@ -68,32 +76,29 @@ function proj_geom = gtGeoProjForReconstruction(projvec, omega, graincenter,...
     detdiru = detgeo.detdiru ./ sqrt(sum(detgeo.detdiru .^ 2)) * detgeo.pixelsizeu;
     detdirv = detgeo.detdirv ./ sqrt(sum(detgeo.detdirv .^ 2)) * detgeo.pixelsizev;
 
-    ones_omega = ones(length(omega), 1);
-    detdiru_rec = gtGeoLab2Sam(detdiru(ones_omega, :), rot_omega, labgeo, recgeo, 1, 0);
-    detdirv_rec = gtGeoLab2Sam(detdirv(ones_omega, :), rot_omega, labgeo, recgeo, 1, 0);
+    detdiru_rec = gtGeoLab2Sam(detdiru(ones_omega, :), rot_omega, labgeo, recgeo, true, false);
+    detdirv_rec = gtGeoLab2Sam(detdirv(ones_omega, :), rot_omega, labgeo, recgeo, true, false);
 
     if (strcmpi(rectype, 'ASTRA_absorption'))
         projvec = labgeo.beamdir(ones_omega, :);
-        projvec = gtGeoLab2Sam(projvec, rot_omega, labgeo, recgeo, 1, 0);
+        projvec = gtGeoLab2Sam(projvec, rot_omega, labgeo, recgeo, true, false);
     end
 
     % Normalised projection vector in REC reference
-    projvec_rec = gtGeoSam2Sam(projvec, samgeo, recgeo, 1, 0);
+    projvec_rec = gtGeoSam2Sam(projvec, samgeo, recgeo, true, false);
     projvec_rec = gtMathsNormalizeVectorsList(projvec_rec);
 
-    switch (rectype)
-        case {'ASTRA_grain'}
-            % Grain center in REC reference
-            graincenter_rec = gtGeoSam2Sam(graincenter, samgeo, recgeo, 0, 0);
-            graincenter_rec = graincenter_rec(ones(size(bbcent_rec, 1), 1), :);
-
-            bbcent_rec = bbcent_rec - graincenter_rec;
+    if (strcmpi(rectype, 'ASTRA_grain'))
+        % Grain center in REC reference
+        graincenter_rec = gtGeoSam2Sam(graincenter, samgeo, recgeo, false, false);
+        graincenter_rec = graincenter_rec(ones_omega, :);
 
-            proj_geom = [projvec_rec, bbcent_rec, detdiru_rec, detdirv_rec];
-
-        case {'ASTRA_full', 'ASTRA_absorption'}
-            proj_geom = [projvec_rec, bbcent_rec, detdiru_rec, detdirv_rec];
+        bbcent_rec = bbcent_rec - graincenter_rec;
+    else
+        % case {'ASTRA_full', 'ASTRA_absorption'}
     end
+
+    proj_geom = [projvec_rec, bbcent_rec, detdiru_rec, detdirv_rec];
 end