Merge branch 'development' into feature-loop-vertices

FlexibleSUSY · Feb 9, 2016 · 633cb25 · 633cb25
2 parents 10c9e08 + 2211abb
commit 633cb25
Show file tree

Hide file tree

Showing 7 changed files with 82 additions and 23 deletions.
diff --git a/ChangeLog b/ChangeLog
@@ -73,6 +73,9 @@ FlexibleSUSY-1.4.0 [not released yet]
    commit, the phase of Dirac fermion singlets is set to e^(i Pi) if
    their mass is less than zero.
 
+ * Bugfix (commits 060b492, a6f7741, 306385b): Implement massless
+   limits in C0, D0 and D27 functions.
+
 FlexibleSUSY-1.3.2 [January, 10 2016]
 
  * Bugfix (commit d76ca79): Fix compilation error with Eigen

diff --git a/...DiracGauginos/LesHouches.in.Diracgauginos → ...DiracGauginos/LesHouches.in.DiracGauginos b/...DiracGauginos/LesHouches.in.Diracgauginos → ...DiracGauginos/LesHouches.in.DiracGauginos
diff --git a/src/numerics.cpp b/src/numerics.cpp
@@ -503,9 +503,18 @@ double d0(double m1, double m2, double m3, double m4) {
   if (close(m1, m2, EPSTOL)) {
     double m2sq = sqr(m2), m3sq = sqr(m3), m4sq = sqr(m4);
 
-    if (close(m2, m3, EPSTOL) && close(m2, m4, EPSTOL)) 
+    if (close(m2,0.,EPSTOL)) {
+       // d0 is undefined for m1 == m2 == 0
+       return 0.;
+    } else if (close(m3,0.,EPSTOL)) {
+       return (-sqr(m2) + sqr(m4) - sqr(m2) * log(sqr(m4/m2)))/
+          sqr(m2 * sqr(m2) - m2 * sqr(m4));
+    } else if (close(m4,0.,EPSTOL)) {
+       return (-sqr(m2) + sqr(m3) - sqr(m2) * log(sqr(m3/m2)))/
+          sqr(m2 * sqr(m2) - m2 * sqr(m3));
+    } else if (close(m2, m3, EPSTOL) && close(m2, m4, EPSTOL)) {
       return 1.0 / (6.0 * sqr(m2sq));
-    else if (close(m2, m3, EPSTOL)) {
+    } else if (close(m2, m3, EPSTOL)) {
       return (sqr(m2sq) - sqr(m4sq) + 2.0 * m4sq * m2sq * log(m4sq / m2sq)) / 
 	(2.0 * m2sq * sqr(m2sq - m4sq) * (m2sq - m4sq));
     } else if (close(m2, m4, EPSTOL)) {
@@ -546,34 +555,58 @@ double c0(double m1, double m2, double m3) {
   double c0l = C0(psq, psq, psq, m1*m1, m2*m2, m3*m3).real();
 #endif
 
-  double ans;
+  double ans = 0.;
 
-  if (close(m2, m3, EPSTOL)) {
+  if (close(m1,0.,EPSTOL) && close(m2,0.,EPSTOL) && close(m3,0.,EPSTOL)) {
+     // c0 is undefined for m1 == m2 == m3 == 0
+     ans = 0.;
+  } else if (close(m2,0.,EPSTOL) && close(m3,0.,EPSTOL)) {
+     // c0 is undefined for m2 == m3 == 0
+     ans = 0.;
+  } else if (close(m1,0.,EPSTOL) && close(m3,0.,EPSTOL)) {
+     // c0 is undefined for m1 == m3 == 0
+     ans = 0.;
+  } else if (close(m1,0.,EPSTOL) && close(m2,0.,EPSTOL)) {
+     // c0 is undefined for m1 == m2 == 0
+     ans= 0.;
+  } else if (close(m1,0.,EPSTOL)) {
+     if (close(m2,m3,EPSTOL)) {
+        ans = -1./sqr(m2);
+     } else {
+        ans = (-log(sqr(m2)) + log(sqr(m3)))/(sqr(m2) - sqr(m3));
+     }
+  } else if (close(m2,0.,EPSTOL)) {
+     if (close(m1,m3,EPSTOL)) {
+        ans = -1./sqr(m1);
+     } else {
+        ans = log(sqr(m3/m1))/(sqr(m1) - sqr(m3));
+     }
+  } else if (close(m3,0.,EPSTOL)) {
+     if (close(m1,m2,EPSTOL)) {
+        ans = -1./sqr(m1);
+     } else {
+        ans = log(sqr(m2/m1))/(sqr(m1) - sqr(m2));
+     }
+  } else if (close(m2, m3, EPSTOL)) {
     if (close(m1, m2, EPSTOL)) {
       ans = ( - 0.5 / sqr(m2) ); // checked 14.10.02
-    }
-    else {
+    } else {
       ans = ( sqr(m1) / sqr(sqr(m1)-sqr(m2) ) * log(sqr(m2)/sqr(m1))
                + 1.0 / (sqr(m1) - sqr(m2)) ) ; // checked 14.10.02
     }
+  } else if (close(m1, m2, EPSTOL)) {
+     ans = ( - ( 1.0 + sqr(m3) / (sqr(m2)-sqr(m3)) * log(sqr(m3)/sqr(m2)) )
+             / (sqr(m2)-sqr(m3)) ) ; // checked 14.10.02
+  } else if (close(m1, m3, EPSTOL)) {
+     ans = ( - (1.0 + sqr(m2) / (sqr(m3)-sqr(m2)) * log(sqr(m2)/sqr(m3)))
+             / (sqr(m3)-sqr(m2)) ); // checked 14.10.02
+  } else {
+     ans = (1.0 / (sqr(m2) - sqr(m3)) *
+            (sqr(m2) / (sqr(m1) - sqr(m2)) *
+             log(sqr(m2) / sqr(m1)) -
+             sqr(m3) / (sqr(m1) - sqr(m3)) *
+             log(sqr(m3) / sqr(m1))) );
   }
-  else
-    if (close(m1, m2, EPSTOL)) {
-      ans = ( - ( 1.0 + sqr(m3) / (sqr(m2)-sqr(m3)) * log(sqr(m3)/sqr(m2)) )
-               / (sqr(m2)-sqr(m3)) ) ; // checked 14.10.02
-    }
-    else
-      if (close(m1, m3, EPSTOL)) {
-        ans = ( - (1.0 + sqr(m2) / (sqr(m3)-sqr(m2)) * log(sqr(m2)/sqr(m3))) 
-                 / (sqr(m3)-sqr(m2)) ); // checked 14.10.02
-      }
-      else {
-	ans = (1.0 / (sqr(m2) - sqr(m3)) * 
-		   (sqr(m2) / (sqr(m1) - sqr(m2)) *
-		    log(sqr(m2) / sqr(m1)) -
-		    sqr(m3) / (sqr(m1) - sqr(m3)) *
-		    log(sqr(m3) / sqr(m1))) );
-      }
 
 #ifdef USE_LOOPTOOLS
   if (!close(c0l, ans, 1.0e-3)) {

diff --git a/src/weinberg_angle.cpp b/src/weinberg_angle.cpp
@@ -25,6 +25,7 @@
 #include "numerics.h"
 
 #include <limits>
+#include <cmath>
 
 #define WARN_IF_ZERO(p,fun)                     \
    if (is_zero(p))                              \
@@ -172,7 +173,16 @@ int Weinberg_angle::calculate(double rho_start, double sin_start)
                              number_of_loops);
 
       if (deltaR > 1.) {
+#if defined(ENABLE_VERBOSE) || defined(ENABLE_DEBUG)
          WARNING("delta_r > 1");
+#endif
+         deltaR = 0.;
+      }
+
+      if (!std::isfinite(deltaR)) {
+#if defined(ENABLE_VERBOSE) || defined(ENABLE_DEBUG)
+         WARNING("delta_r non-finite");
+#endif
          deltaR = 0.;
       }
 

diff --git a/templates/two_scale_initial_guesser_low_scale_model.cpp.in b/templates/two_scale_initial_guesser_low_scale_model.cpp.in
@@ -183,10 +183,16 @@ void @ModelName@_initial_guesser<Two_scale>::guess_soft_parameters()
 
    model->run_to(susy_scale_guess, running_precision);
 
+   // set EWSB loop order to 0 temporarily
+   const unsigned lo = model->get_ewsb_loop_order();
+   model->set_ewsb_loop_order(0);
+
    // apply susy-scale constraint
    susy_constraint.set_model(model);
    susy_constraint.apply();
 
+   model->set_ewsb_loop_order(lo);
+
    model->run_to(low_scale_guess, running_precision);
 
    // apply EWSB constraint

diff --git a/templates/two_scale_low_scale_constraint.cpp.in b/templates/two_scale_low_scale_constraint.cpp.in
@@ -235,6 +235,12 @@ void @ModelName@_low_scale_constraint<Two_scale>::calculate_threshold_correction
    AlphaS = alpha_s_drbar;
    EDRbar = e_drbar;
    ThetaWDRbar = calculate_theta_w(alpha_em_drbar);
+
+   if (!IsFinite(ThetaWDRbar)) {
+      model->get_problems().flag_non_perturbative_parameter(
+         "sin(theta_W)", ThetaWDRbar, model->get_scale(), 0);
+      ThetaWDRbar = ArcSin(Electroweak_constants::sinThetaW);
+   }
 }
 
 double @ModelName@_low_scale_constraint<Two_scale>::calculate_theta_w(double alpha_em_drbar)

diff --git a/test/test_run_all_spectrum_generators.sh b/test/test_run_all_spectrum_generators.sh
@@ -92,6 +92,7 @@ HGTHDMIIMSSMBC,_DEFAULT_,0
 TMSSM,_DEFAULT_,0
 UMSSM,_DEFAULT_,0
 YukawaCMSSM,${DEFAULT_CMSSM_INPUT},0
+DiracGauginos,_DEFAULT_,0
 "
 
 TMP_FILE="${BASEDIR}/test_spectrum_generator.sh.tmp"