diff --git a/image_analysis/geom_refinement/XtalOptimizer.cpp b/image_analysis/geom_refinement/XtalOptimizer.cpp
index 2d9b35ff..a6d030b3 100644
--- a/image_analysis/geom_refinement/XtalOptimizer.cpp
+++ b/image_analysis/geom_refinement/XtalOptimizer.cpp
@@ -9,7 +9,6 @@ struct XtalResidual {
     XtalResidual(double x, double y,
                  double lambda,
                  double pixel_size,
-                 double rot1, double rot2,
                  const Eigen::Matrix3d &gonio_back_rot,
                  double exp_h, double exp_k,
                  double exp_l,
@@ -20,8 +19,6 @@ struct XtalResidual {
           exp_h(exp_h),
           exp_k(exp_k),
           exp_l(exp_l),
-          rot1(rot1),
-          rot2(rot2),
           gonio_back_rot_(gonio_back_rot),
           symmetry(symmetry) {
     }
@@ -30,14 +27,15 @@ struct XtalResidual {
     bool operator()(const T *const beam_x,
                     const T *const beam_y,
                     const T *const distance_mm,
+                    const T *const detector_rot,
                     const T *const p0,
                     const T *const p1,
                     const T *const p2,
                     T *residual) const {
-        T c1 = ceres::cos(T(rot1));
-        T c2 = ceres::cos(T(rot2));
-        T s1 = ceres::sin(T(rot1));
-        T s2 = ceres::sin(T(rot2));
+        T c1 = ceres::cos(T(detector_rot[0]));
+        T c2 = ceres::cos(T(detector_rot[1]));
+        T s1 = ceres::sin(T(detector_rot[0]));
+        T s2 = ceres::sin(T(detector_rot[1]));
 
         // x_lab in mm
         T x_lab = (T(obs_x) - beam_x[0]) * T(pixel_size);
@@ -141,7 +139,6 @@ struct XtalResidual {
     const double exp_h;
     const double exp_k;
     const double exp_l;
-    const double rot1, rot2;
     const Eigen::Matrix3d gonio_back_rot_;
     gemmi::CrystalSystem symmetry;
 };
@@ -379,6 +376,8 @@ bool XtalOptimizerInternal(XtalOptimizerData &data,
         double beam_y = data.geom.GetBeamY_pxl();
         double distance_mm = data.geom.GetDetectorDistance_mm();
 
+        double detector_rot[2] = {data.geom.GetPoniRot1_rad(), data.geom.GetPoniRot2_rad()};
+
         ceres::Problem problem;
 
         double latt_vec0[3], latt_vec1[3], latt_vec2[3];
@@ -447,16 +446,11 @@ bool XtalOptimizerInternal(XtalOptimizerData &data,
                 continue;
 
 
-
-
-
             problem.AddResidualBlock(
-                new ceres::AutoDiffCostFunction<XtalResidual, 3, 1, 1, 1, 3, 3, 3>(
+                new ceres::AutoDiffCostFunction<XtalResidual, 3, 1, 1, 1, 2, 3, 3, 3>(
                     new XtalResidual(pt.x, pt.y,
                                      data.geom.GetWavelength_A(),
                                      data.geom.GetPixelSize_mm(),
-                                     data.geom.GetPoniRot1_rad(),
-                                     data.geom.GetPoniRot2_rad(),
                                      gonio_back_rot,
                                      h, k, l,
                                      data.crystal_system)),
@@ -464,6 +458,7 @@ bool XtalOptimizerInternal(XtalOptimizerData &data,
                 &beam_x,
                 &beam_y,
                 &distance_mm,
+                detector_rot,
                 latt_vec0,
                 latt_vec1,
                 latt_vec2
@@ -475,12 +470,27 @@ bool XtalOptimizerInternal(XtalOptimizerData &data,
 
         if (!data.refine_distance_mm)
             problem.SetParameterBlockConstant(&distance_mm);
+        else {
+            const double dist_range = 0.1;
+            problem.SetParameterLowerBound(&distance_mm, 0, distance_mm * (1.0 - dist_range));
+            problem.SetParameterLowerBound(&distance_mm, 0, distance_mm * (1.0 + dist_range));
+        }
 
         if (!data.refine_beam_center) {
             problem.SetParameterBlockConstant(&beam_x);
             problem.SetParameterBlockConstant(&beam_y);
         }
 
+        if (!data.refine_detector_angles) {
+            problem.SetParameterBlockConstant(detector_rot);
+        } else {
+            const double rot_range = 3.0 / 180.0 * M_PI;
+            for (int i = 0; i < 2; ++i) {
+                problem.SetParameterLowerBound(detector_rot, i, detector_rot[i] - rot_range);
+                problem.SetParameterUpperBound(detector_rot, i, detector_rot[i] + rot_range);
+            }
+        }
+
         // Parameter bounds
         // Lengths
         for (int i = 0; i < 3; ++i) {
@@ -522,6 +532,9 @@ bool XtalOptimizerInternal(XtalOptimizerData &data,
         if (data.refine_distance_mm)
             data.geom.DetectorDistance_mm(distance_mm);
 
+        if (data.refine_detector_angles)
+            data.geom.PoniRot1_rad(detector_rot[0]).PoniRot2_rad(detector_rot[1]);
+
         if (data.crystal_system == gemmi::CrystalSystem::Orthorhombic)
             data.latt = AngleAxisAndLengthsToLattice(latt_vec0, latt_vec1, false);
         else if (data.crystal_system == gemmi::CrystalSystem::Tetragonal) {
@@ -551,18 +564,22 @@ bool XtalOptimizerInternal(XtalOptimizerData &data,
             const double sg = std::sin(latt_vec2[2]);
 
             // a along x, b in x-y, c general
-            B(0,0) = 1.0;  B(1,0) = 0.0;  B(2,0) = 0.0;
-            B(0,1) = cg;   B(1,1) = sg;   B(2,1) = 0.0;
+            B(0, 0) = 1.0;
+            B(1, 0) = 0.0;
+            B(2, 0) = 0.0;
+            B(0, 1) = cg;
+            B(1, 1) = sg;
+            B(2, 1) = 0.0;
 
             const double cx = cb;
-            const double cy = (ca - cb*cg) / sg;
-            const double cz = std::sqrt(std::max(0.0, 1.0 - cx*cx - cy*cy));
+            const double cy = (ca - cb * cg) / sg;
+            const double cz = std::sqrt(std::max(0.0, 1.0 - cx * cx - cy * cy));
 
-            B(0,2) = cx;
-            B(1,2) = cy;
-            B(2,2) = cz;
+            B(0, 2) = cx;
+            B(1, 2) = cy;
+            B(2, 2) = cz;
 
-            Eigen::DiagonalMatrix<double,3> D(latt_vec1[0], latt_vec1[1], latt_vec1[2]);
+            Eigen::DiagonalMatrix<double, 3> D(latt_vec1[0], latt_vec1[1], latt_vec1[2]);
             Eigen::Matrix3d latt = R * D * B;
 
             data.latt = CrystalLattice(Coord(latt(0, 0), latt(1, 0), latt(2, 0)),
diff --git a/image_analysis/geom_refinement/XtalOptimizer.h b/image_analysis/geom_refinement/XtalOptimizer.h
index 5cc0d09e..155128f4 100644
--- a/image_analysis/geom_refinement/XtalOptimizer.h
+++ b/image_analysis/geom_refinement/XtalOptimizer.h
@@ -25,6 +25,7 @@ struct XtalOptimizerData {
 
     bool refine_beam_center = true;
     bool refine_distance_mm = false;
+    bool refine_detector_angles = false;
 
     std::optional<GoniometerAxis> axis;