XtalOptimizer: Work in progress to find out what is wrong with rotation axis

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;
 }

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;

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_
     };
 }
 

image_analysis/RotationIndexer.h

@@ -23,6 +23,7 @@ struct RotationIndexerResult {
     CrystalLattice lattice;
     LatticeSearchResult search_result;
     DiffractionGeometry geom;
+    std::optional<GoniometerAxis> axis;
 };
 
 class RotationIndexer {

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);