diff --git a/image_analysis/bragg_integration/BraggPrediction.cpp b/image_analysis/bragg_integration/BraggPrediction.cpp
index 7073be7a..2bb78a74 100644
--- a/image_analysis/bragg_integration/BraggPrediction.cpp
+++ b/image_analysis/bragg_integration/BraggPrediction.cpp
@@ -8,7 +8,6 @@ std::multimap<float, Reflection> CalcBraggPredictions(const DiffractionExperimen
                                                       float high_res_A,
                                                       float ewald_dist_cutoff,
                                                       int max_hkl) {
-    // Assumption: Scattering vector is (0,0,1)
     std::multimap<float, Reflection> reflections;
 
     auto geom = experiment.GetDiffractionGeometry();
@@ -18,11 +17,12 @@ std::multimap<float, Reflection> CalcBraggPredictions(const DiffractionExperimen
     float one_over_dmax = 1.0f / high_res_A;
     float one_over_dmax_sq = one_over_dmax * one_over_dmax;
 
-    float one_over_wavelength = 1.0f/ geom.GetWavelength_A();
+    float one_over_wavelength = 1.0f / geom.GetWavelength_A();
 
     Coord Astar = lattice.Astar();
     Coord Bstar = lattice.Bstar();
     Coord Cstar = lattice.Cstar();
+    Coord S0 = geom.GetScatteringVector();
 
     for (int h = -max_hkl; h < max_hkl; h++) {
         // Precompute A* h contribution
@@ -45,38 +45,35 @@ std::multimap<float, Reflection> CalcBraggPredictions(const DiffractionExperimen
                 float recip_y = AhBk_y + Cstar.y * l;
                 float recip_z = AhBk_z + Cstar.z * l;
 
-                // Check resolution (actually through 1/d^2)
-                float dot_xy = recip_y * recip_y + recip_z * recip_z;
-                float dot = dot_xy + recip_x * recip_x;
+                float dot = recip_x * recip_x + recip_y * recip_y + recip_z * recip_z;
                 if (dot > one_over_dmax_sq)
                     continue;
 
-                // Get Distance from Ewald sphere, to check if reflection is in diffracting condition
-                float S_z = recip_z + one_over_wavelength;
-                float S_len = sqrtf(dot_xy + S_z * S_z);
-
+                float S_x = recip_x + S0.x;
+                float S_y = recip_y + S0.y;
+                float S_z = recip_z + S0.z;
+                float S_len = sqrtf(S_x * S_x + S_y * S_y + S_z * S_z);
                 float dist_ewald_sphere = std::fabs(S_len - one_over_wavelength);
 
-                if (!std::isfinite(dist_ewald_sphere) || dist_ewald_sphere > ewald_dist_cutoff )
-                    continue;
+                if (dist_ewald_sphere <= ewald_dist_cutoff ) {
+                    Coord recip(recip_x, recip_y, recip_z);
 
-                Coord recip(recip_x, recip_y, recip_z);
+                    auto [x,y] = geom.RecipToDector(recip);
 
-                auto [x,y] = geom.RecipToDector(recip);
+                    if ((x < 0) || (x >= det_width_pxl) || (y < 0) || (y >= det_height_pxl))
+                        continue;
 
-                if ((x < 0) || (x >= det_width_pxl) || (y < 0) || (y >= det_height_pxl))
-                    continue;
+                    float d = 1.0f / sqrtf(dot);
 
-                float d = 1.0f / sqrtf(dot);
-
-                reflections.insert(std::make_pair(dist_ewald_sphere, Reflection{
-                        .h = h,
-                        .k = k,
-                        .l = l,
-                        .predicted_x = x,
-                        .predicted_y = y,
-                        .d = d
-                }));
+                    reflections.insert(std::make_pair(dist_ewald_sphere, Reflection{
+                            .h = h,
+                            .k = k,
+                            .l = l,
+                            .predicted_x = x,
+                            .predicted_y = y,
+                            .d = d
+                    }));
+                }
             }
         }
     }