Fix pyramid form factor

Core/FormFactors/src/FormFactorPyramid.cpp

@@ -101,14 +101,18 @@ complex_t FormFactorPyramid::fullPyramidPrimitive(complex_t a, complex_t b, comp
 {
     const complex_t im(0, 1);
     if (std::norm(a*z) > Numeric::double_epsilon && std::norm(b*z) > Numeric::double_epsilon) {
-        complex_t phase = std::exp(im * c * z);
-        complex_t numerator = std::sin(a * z) * (b * (a * a - b * b + c * c) * std::cos(b * z)
-                                                 + im * c * (a * a + b * b - c * c) * std::sin(b * z))
-                              + a * std::cos(a * z) * ((-a * a + b * b + c * c) * std::sin(b * z)
-                                                       + 2.0 * im * b * c * std::cos(b * z));
-        complex_t denominator = a * b * (a - b - c) * (a + b - c) * (a - b + c) * (a + b + c);
-//        return -4.0 * phase * numerator / denominator;
-        return 2.0 * (g(a-b, c, z) - g(a+b, c, z)) / (a*b);
+        if (std::abs((a-b)*(a-b)-c*c)*z*z > Numeric::double_epsilon &&
+            std::abs((a+b)*(a+b)-c*c)*z*z > Numeric::double_epsilon) {
+            complex_t phase = std::exp(im * c * z);
+            complex_t numerator = std::sin(a * z) * (b * (a * a - b * b + c * c) * std::cos(b * z)
+                                                     + im * c * (a * a + b * b - c * c) * std::sin(b * z))
+                    + a * std::cos(a * z) * ((-a * a + b * b + c * c) * std::sin(b * z)
+                                             + 2.0 * im * b * c * std::cos(b * z));
+            complex_t denominator = a * b * (a - b - c) * (a + b - c) * (a - b + c) * (a + b + c);
+            return -4.0 * phase * numerator / denominator;
+        } else {
+            return 2.0 * (g(a-b, c, z) - g(a+b, c, z)) / (a*b);
+        }
     } else if (std::norm(a*z) <= Numeric::double_epsilon
                && std::norm(b*z) <= Numeric::double_epsilon) {
         if (std::norm(c*z) <= Numeric::double_epsilon) {