Merge remote-tracking branch 'organization/development' into feature-2.0

FlexibleSUSY · Nov 8, 2016 · fb93c43 · fb93c43
2 parents aef27cd + 627502c
commit fb93c43
Show file tree

Hide file tree

Showing 5 changed files with 43 additions and 23 deletions.
diff --git a/meta/ThreeLoopQCD.m b/meta/ThreeLoopQCD.m
@@ -50,7 +50,7 @@
            colorPosition = Position[SARAH`Gauge, SARAH`color][[1,1]];
            CF = SA`Casimir[quark, colorPosition];
            as = SARAH`strongCoupling^2 / (4 Pi);
-           CF as / Pi (-1 + 3/4 Log[FlexibleSUSY`M[quark]^2/Q^2])
+           CF as / Pi (-1 + 3/4 Log[FlexibleSUSY`Pole[FlexibleSUSY`M[quark]]^2/Q^2])
           ];
 
 (* 2-loop coefficients *)
@@ -109,7 +109,7 @@
            as = SARAH`strongCoupling^2 / (4 Pi);
            CF (as/Pi)^2 (
                CF d12 + CA d22 + TR NL d32 + TR NH d42 +
-               Get2LLogs[FlexibleSUSY`M[quark], Q, CF, CA]
+               Get2LLogs[FlexibleSUSY`Pole[FlexibleSUSY`M[quark]], Q, CF, CA]
            )
           ];
 
@@ -230,13 +230,26 @@
           ];
 
 GetMTopMSbarOverMTopPole[loopOrder_List:{1,1,1,1}, quark_:SARAH`TopQuark, Q_:Q, NH_:1, NL_:5] :=
-    Module[{},
+    Module[{result, Mpole, h},
            Assert[Length[loopOrder] == 4];
+           result = (
+               h^0 GetMTopMSbarOverMTopPole0L[] +
+               h^1 GetMTopMSbarOverMTopPole1L[Q, quark] +
+               h^2 GetMTopMSbarOverMTopPole2L[Q, quark, NH, NL] +
+               h^3 GetMTopMSbarOverMTopPole3L[Q, quark, NH, NL]
+           );
+           (* rewrite in terms of the running mass *)
+           Mpole = FlexibleSUSY`M[quark] / result;
+           result = Normal @ Series[
+               result /. FlexibleSUSY`Pole[FlexibleSUSY`M[quark]] -> Mpole /.
+                         FlexibleSUSY`Pole[FlexibleSUSY`M[quark]] -> Mpole /.
+                         FlexibleSUSY`Pole[FlexibleSUSY`M[quark]] -> FlexibleSUSY`M[quark],
+               {h,0,3}];
            (
-               loopOrder[[1]] GetMTopMSbarOverMTopPole0L[] +
-               loopOrder[[2]] GetMTopMSbarOverMTopPole1L[Q, quark] +
-               loopOrder[[3]] GetMTopMSbarOverMTopPole2L[Q, quark, NH, NL] +
-               loopOrder[[4]] GetMTopMSbarOverMTopPole3L[Q, quark, NH, NL]
+               loopOrder[[1]] Coefficient[result, h, 0] +
+               loopOrder[[2]] Coefficient[result, h, 1] +
+               loopOrder[[3]] Coefficient[result, h, 2] +
+               loopOrder[[4]] Coefficient[result, h, 3]
            )
           ];
 

diff --git a/meta/TwoLoopQCD.m b/meta/TwoLoopQCD.m
@@ -22,7 +22,11 @@
 
 GetDeltaMPoleOverMRunningQCDTwoLoopMSbar::usage="Returns two-loop QCD
 contributions to Delta M_f/m_f in the MS-bar scheme.  Taken from
-hep-ph/9803493 Eq. (17).  (M_f = pole mass, m_f = MS-bar mass)";
+hep-ph/9803493 Eq. (17).  (M_f = pole mass, m_f = MS-bar mass)
+
+Note: Eq. (17) of [hep-ph/9803493] is expressed in terms of
+Log[Q^2/M_f^2].  This function, however, returns Delta M_f/m_f as a
+function of Log[Q^2/m_f^2].  We have accounted for the difference.";
 
 Begin["`Private`"];
 
@@ -91,7 +95,7 @@
            result = (alphaStrong / (4 Pi))^2 (
                CF CA (1111/24 - 8 Zeta[2] - 4 I31 + 185/6 log + 11/2 log^2)
                - CF t Nf (71/6 + 8 Zeta[2] + 26/3 log + 2 log^2)
-               + CF^2 (121/8 + 30 Zeta[2] + 8 I31 + 27/2 log + 9/2 log^2)
+               + CF^2 (-71/8 + 30 Zeta[2] + 8 I31 - 9/2 log + 9/2 log^2)
                - 12 CF t (1 - 2 Zeta[2]));
            result
           ];

diff --git a/src/standard_model.cpp b/src/standard_model.cpp
@@ -4268,8 +4268,8 @@ void Standard_model::calculate_MFu_pole()
 
    if (pole_mass_loop_order > 1 && TOP_POLE_QCD_CORRECTION > 0) {
       const double currentScale = get_scale();
-      qcd_2l = -0.006995771808874528*Power(g3,4) -
-         0.004518101565212638*Power(g3,4)*Log(Sqr(currentScale)/Sqr(MFu(2))) -
+      qcd_2l = -0.005284774766427138*Power(g3,4) -
+         0.0032348537833770956*Power(g3,4)*Log(Sqr(currentScale)/Sqr(MFu(2))) -
          0.0008822328500119351*Power(g3,4)*Sqr(Log(Power(currentScale,2)/Sqr(
          MFu(2))));
    }
@@ -4278,10 +4278,10 @@ void Standard_model::calculate_MFu_pole()
 
    if (pole_mass_loop_order > 2 && TOP_POLE_QCD_CORRECTION > 1) {
       const double currentScale = get_scale();
-      qcd_3l = Power(g3,6)*(-0.0017408026847411467 -
-         0.000984413176263005*Log(Sqr(currentScale)/Sqr(MFu(2))) -
-         0.00003352082872926087*Power(Log(Sqr(currentScale)/Sqr(MFu(2))),3) -
-         0.00029813221915266867*Sqr(Log(Power(currentScale,2)/Sqr(MFu(2)))));
+      qcd_3l = -0.00003352082872926087*Power(g3,6)*(35.702577217116016
+         + 15.387410814884797*Log(Sqr(currentScale)/Sqr(MFu(2))) + 1.*Power(
+         Log(Sqr(currentScale)/Sqr(MFu(2))),3) + 5.378787878787879*Sqr(Log(
+         Power(currentScale,2)/Sqr(MFu(2)))));
    }
 
    Eigen::Matrix<double,3,3> self_energy_1;
@@ -4458,18 +4458,20 @@ double Standard_model::calculate_MFu_DRbar(double m_pole, int idx) const
    double qcd_2l = 0., qcd_3l = 0.;
 
    if (get_thresholds() > 1) {
-      qcd_2l = -0.005855107113909601*Power(g3,4) -
-         0.0028071045227652486*Power(g3,4)*Log(Sqr(currentScale)/Sqr(MFu(idx)))
+      qcd_2l = -0.0041441100714622115*Power(g3,4) -
+         0.0015238567409297061*Power(g3,4)*Log(Sqr(currentScale)/Sqr(MFu(idx)))
          - 0.00024060895909416413*Power(g3,4)*Sqr(Log(Power(currentScale,2)
          /Sqr(MFu(idx))));
    }
 
    if (get_thresholds() > 2) {
-      qcd_3l = -0.0013067805969741943*Power(g3,6) -
+      qcd_3l = -0.0008783313853540776*Power(g3,6) -
          0.0004114970933517977*Power(g3,6)*Log(Sqr(currentScale)/Sqr(MFu(idx)))
          - 5.078913443827405e-6*Power(g3,6)*Power(Log(Sqr(currentScale)/Sqr(
-         MFu(idx))),3) - 0.00007466002762426286*Power(g3,6)*Sqr(Log(Power(
-         currentScale,2)/Sqr(MFu(idx))));
+         MFu(idx))),3) - 0.0002952541682011665*Power(g3,6)*Log(Sqr(MFu(idx))
+         /Sqr(currentScale)) + 0.00005282069981580501*Power(g3,6)*Sqr(Log(Power
+         (MFu(idx),2)/Sqr(currentScale))) - 0.00007466002762426286*Power(g3,6)*
+         Sqr(Log(Power(currentScale,2)/Sqr(MFu(idx))));
    }
 
    const double m_susy_drbar = m_pole + self_energy_1 + m_pole * (

diff --git a/test/test_SM_three_loop_spectrum.cpp b/test/test_SM_three_loop_spectrum.cpp
@@ -47,7 +47,7 @@ BOOST_AUTO_TEST_CASE( test_SM_two_loop_top_pole_mass )
 
    unsigned iterations = 100;
 
-   // calculate top DR-bar mass from top pole mass using two-loop
+   // calculate top DR-bar mass from top pole mass using three-loop
    // corrections
    do {
       Eigen::Matrix<double,3,3> mt_drbar_3loop(Eigen::Matrix<double,3,3>::Zero());
@@ -71,5 +71,5 @@ BOOST_AUTO_TEST_CASE( test_SM_two_loop_top_pole_mass )
 
    BOOST_TEST_MESSAGE("mt_pole(3-loop) = " << mt_pole_3loop);
 
-   BOOST_CHECK_CLOSE_FRACTION(mt_pole_input, mt_pole_3loop, 3.0e-4);
+   BOOST_CHECK_CLOSE_FRACTION(mt_pole_input, mt_pole_3loop, 2.0e-4);
 }
diff --git a/test/test_ThreeLoopQCD.m b/test/test_ThreeLoopQCD.m
@@ -52,7 +52,8 @@
     Yu[__] -> 0,
     Yd[__] -> 0,
     g1 -> 0,
-    g2 -> 0
+    g2 -> 0,
+    FlexibleSUSY`M[Fu] -> m[Q]
 };
 
 ass = {Q > 0, FlexibleSUSY`M[Fu] > 0, m[Q] > 0};