diff --git a/common/GoniometerAxis.cpp b/common/GoniometerAxis.cpp
index 163316d9..f2e95ad7 100644
--- a/common/GoniometerAxis.cpp
+++ b/common/GoniometerAxis.cpp
@@ -37,10 +37,11 @@ GoniometerAxis &GoniometerAxis::ScreeningWedge(const std::optional<float> &input
 }
 
 GoniometerAxis &GoniometerAxis::Axis(const Coord &input) {
-    if (input.Length() == 0.0f)
+    float len = input.Length();
+    if (len == 0.0f)
         throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                               "Rotation axis cannot have 0 length");
-
+    // increment *= len;
     axis = input.Normalize();
     return *this;
 }
diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp
index 7ce56467..66d57c39 100644
--- a/image_analysis/IndexAndRefine.cpp
+++ b/image_analysis/IndexAndRefine.cpp
@@ -38,7 +38,8 @@ IndexAndRefine::IndexingOutcome IndexAndRefine::DetermineLatticeAndSymmetry(Data
                     .BeamY_pxl(result->geom.GetBeamY_pxl())
                     .DetectorDistance_mm(result->geom.GetDetectorDistance_mm())
                     .PoniRot1_rad(result->geom.GetPoniRot1_rad())
-                    .PoniRot2_rad(result->geom.GetPoniRot2_rad());
+                    .PoniRot2_rad(result->geom.GetPoniRot2_rad())
+            .Goniometer(result->axis);
             outcome.symmetry.centering = result->search_result.centering;
             outcome.symmetry.niggli_class = result->search_result.niggli_class;
             outcome.symmetry.crystal_system = result->search_result.system;
diff --git a/image_analysis/RotationIndexer.cpp b/image_analysis/RotationIndexer.cpp
index dc2c1254..7affc00d 100644
--- a/image_analysis/RotationIndexer.cpp
+++ b/image_analysis/RotationIndexer.cpp
@@ -87,6 +87,7 @@ void RotationIndexer::TryIndex() {
             .refine_beam_center = true,
             .refine_distance_mm =  false,
             .refine_detector_angles = true,
+            .refine_rotation_axis = true,
             .index_ice_rings = experiment.GetIndexingSettings().GetIndexIceRings(),
             .axis = axis_
         };
@@ -97,6 +98,7 @@ void RotationIndexer::TryIndex() {
         if (XtalOptimizer(data, v_sel)) {
             indexed_lattice = data.latt;
             updated_geom_ = data.geom;
+            axis_ = data.axis;
         }
         PredictionRotationSettings pred_settings{
             .high_res_A = 1.0,
@@ -139,7 +141,8 @@ std::optional<RotationIndexerResult> RotationIndexer::ProcessImage(int64_t image
     return RotationIndexerResult{
         .lattice = indexed_lattice.value(),
         .search_result = search_result_,
-        .geom = updated_geom_
+        .geom = updated_geom_,
+        .axis = axis_
     };
 }
 
@@ -149,7 +152,8 @@ std::optional<RotationIndexerResult> RotationIndexer::GetLattice() {
     return RotationIndexerResult{
         .lattice = indexed_lattice.value(),
         .search_result = search_result_,
-        .geom = updated_geom_
+        .geom = updated_geom_,
+        .axis = axis_
     };
 }
 
diff --git a/image_analysis/RotationIndexer.h b/image_analysis/RotationIndexer.h
index 14b2efd5..56d55768 100644
--- a/image_analysis/RotationIndexer.h
+++ b/image_analysis/RotationIndexer.h
@@ -23,6 +23,7 @@ struct RotationIndexerResult {
     CrystalLattice lattice;
     LatticeSearchResult search_result;
     DiffractionGeometry geom;
+    std::optional<GoniometerAxis> axis;
 };
 
 class RotationIndexer {
diff --git a/image_analysis/geom_refinement/XtalOptimizer.cpp b/image_analysis/geom_refinement/XtalOptimizer.cpp
index 5e000cbd..c2d9b1ba 100644
--- a/image_analysis/geom_refinement/XtalOptimizer.cpp
+++ b/image_analysis/geom_refinement/XtalOptimizer.cpp
@@ -564,18 +564,8 @@ bool XtalOptimizerInternal(XtalOptimizerData &data,
         if (data.refine_detector_angles)
             data.geom.PoniRot1_rad(detector_rot[0]).PoniRot2_rad(detector_rot[1]);
 
-        if (!data.refine_unit_cell) {
-            Coord rot_vector{
-                static_cast<float>(latt_vec0[0]), static_cast<float>(latt_vec0[1]), static_cast<float>(latt_vec0[2])
-            };
-            Coord rot_vector_norm = rot_vector.Normalize();
-            std::cout << "X " << rot_vector.Length() * 180.0 / M_PI << " " << rot_vector_norm << std::endl;
-        }
-
-        if (data.axis && data.refine_rotation_axis) {
+        if (data.axis && data.refine_rotation_axis)
             data.axis.value().Axis(Coord(rot_vec[0], rot_vec[1], rot_vec[2]));
-            std::cout << "Refined rotation axis " << rot_vec[0] << " " << rot_vec[1] << " " << rot_vec[2] << std::endl;
-        }
 
         if (data.crystal_system == gemmi::CrystalSystem::Orthorhombic)
             data.latt = AngleAxisAndLengthsToLattice(latt_vec0, latt_vec1, false);