Merge branch 'feature-2.0' into feature-2.0-semianalytic-solver

Conflicts: meta/FlexibleSUSY.m src/module.mk
FlexibleSUSY · Apr 2, 2017 · c0a56c2 · c0a56c2
2 parents 9f3aaf4 + e68ba03
commit c0a56c2
Show file tree

Hide file tree

Showing 88 changed files with 5,923 additions and 8,439 deletions.
diff --git a/.gitattributes b/.gitattributes
@@ -14,6 +14,8 @@ examples/lattice_numerical_fmssm.cpp        export-ignore
 examples/lattice_numerical_fmssm_fmssmn.cpp export-ignore
 examples/switch_MSSM.cpp                    export-ignore
 
+meta/SM/HSSUSY_corrections.m    export-ignore
+
 model_files/BetaSM              export-ignore
 model_files/cCMSSM              export-ignore
 model_files/complexMSSM         export-ignore
@@ -31,6 +33,7 @@ model_files/CMSSMMediumPrecision export-ignore
 model_files/CMSSMHighPrecision  export-ignore
 model_files/CMSSMMassWInput     export-ignore
 model_files/CMSSMConvergenceTester export-ignore
+model_files/CMSSMYt2L           export-ignore
 model_files/cSM                 export-ignore
 model_files/BLSMlightZp         export-ignore
 model_files/InertMSSM           export-ignore

diff --git a/ChangeLog b/ChangeLog
@@ -31,6 +31,15 @@ FlexibleSUSY 2.0.0 [not released yet]
 
  * Feature: Support for SLHA-2 input block IMEXTPAR.
 
+ * Feature: The full 2-loop O(alpha_s^2) corrections to the DR-bar top
+   Yukawa coupling in the MSSM [arxiv:hep-ph/0210258,
+   arxiv:hep-ph/0507139] can be added by setting
+
+     UseMSSMYukawa3LoopSQCD = True
+
+   in the model file.
+   Thanks to Alexander Bednyakov.
+
  * Change: The interface for adding constraints and matching
    conditions to the Two_scale_solver class has been simplified.
    Constrains and matching conditions are now added using the add()
@@ -85,6 +94,16 @@ FlexibleSUSY 2.0.0 [not released yet]
 
    Thanks to Peter Athron and Pat Scott.
 
+FlexibleSUSY-1.7.4 [not released yet]
+
+ * Bugfix (commit f434e30): Rename internal IndexSum symbol which
+   conflicts with SARAH version 4.11.0 and higher.
+
+ * Bugfix (commit b8d5dcf): Correcting gauge-dependent term in 2-loop
+   beta function of SM vacuum expectation value after a corresponding
+   bugfix in SARAH 4.11.0.  This bugfix affects the Higgs mass
+   prediction with FlexibleEFTHiggs by around 10 MeV.
+
 FlexibleSUSY-1.7.3 [February, 27 2017]
 
  * Change (commit 43bb03a): FlexibleSUSY now aborts the code

diff --git a/README b/README
@@ -19,7 +19,7 @@ Requirements
  * C++ compiler (g++ >= 4.7.2 or clang++ >= 3.1 or icpc >= 12.1)
  * Fortran compiler (gfortran, ifort)
  * Mathematica (version 7.0 or higher)
- * SARAH (version 4.10.2 or higher)   http://sarah.hepforge.org
+ * SARAH (version 4.11.0 or higher)   http://sarah.hepforge.org
  * Boost (version 1.37.0 or higher)   http://www.boost.org
  * Eigen 3 (version 3.1 or higher)    http://eigen.tuxfamily.org
  * GNU scientific library             http://www.gnu.org/software/gsl/

diff --git a/addons/GM2Calc/MSSMNoFV_onshell.cpp b/addons/GM2Calc/MSSMNoFV_onshell.cpp
@@ -522,12 +522,22 @@ void MSSMNoFV_onshell::convert_ml2()
    const double g12 = sqr(get_g1());
    const double g22 = sqr(get_g2());
 
+   if (verbose_output) {
+      std::cout << "Converting msl(2,2) to on-shell scheme ...\n"
+         "   Using MSvm_pole = " << MSvmL_pole << '\n';
+   }
+
    // calculate ml2(1,1) from muon sneutrino pole mass
    const double ml211
       = sqr(MSvmL_pole) + 0.125*(0.6*g12*(vu2 - vd2) + g22*(vu2 - vd2));
 
    set_ml2(1,1,ml211);
    calculate_MSvmL();
+
+   if (verbose_output) {
+      std::cout << "   New msl(2,2) = " << ml211
+                << ", new MSvmL = " << get_MSvmL() << '\n';
+   }
 }
 
 /**
@@ -656,9 +666,6 @@ double MSSMNoFV_onshell::convert_me2_fpi(
       std::cout << "Converting mse(2,2) to on-shell scheme with FPI ...\n"
                    "   Goal: MSm(" << right_index << ") = "
                 << MSm_goal(right_index) << '\n';
-      std::cout << "      "
-                << "MSm_pole = " << get_physical().MSm.transpose()
-                << ", MSm_pole_sorted = " << MSm_pole_sorted.transpose() << '\n';
    }
 
    bool accuracy_goal_reached =

diff --git a/configure b/configure
@@ -239,8 +239,8 @@ required_mathematica_version="7"
 
 # required SARAH version
 required_sarah_major="4"
-required_sarah_minor="10"
-required_sarah_patch="2"
+required_sarah_minor="11"
+required_sarah_patch="0"
 required_sarah_version="${required_sarah_major}.${required_sarah_minor}.${required_sarah_patch}"
 SARAH_VERSION="unknown"
 SARAH_MAJOR="0"

diff --git a/doc/building.dox b/doc/building.dox
@@ -9,7 +9,7 @@
  - C++ compiler (`g++` >= 4.7.2 or `clang++` >= 3.1 or `icpc` >= 12.1)
  - Fortran compiler (`gfortran`, `ifort`)
  - Mathematica (version 7.0 or higher)
- - SARAH (version 4.0.4 or higher)    http://sarah.hepforge.org
+ - SARAH (version 4.11.0 or higher)   http://sarah.hepforge.org
  - Boost (version 1.37.0 or higher)   http://www.boost.org
  - Eigen 3 (version 3.1 or higher)    http://eigen.tuxfamily.org
  - GNU scientific library             http://www.gnu.org/software/gsl/

diff --git a/doc/model_file.dox b/doc/model_file.dox
@@ -1333,4 +1333,16 @@ in the model file
 UseMSSM3LoopRGEs = True; (* use three-loop MSSM RGEs *)
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
+### Two-loop threshold corrections ###
+
+#### MSSM ####
+
+In the MSSM the known two-loop relation between the top pole mass and
+the DR-bar top mass from [arxiv:hep-ph/0210258,arxiv:hep-ph/0507139]
+can be used by setting in the model file
+
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.m}
+UseMSSMYukawa3LoopSQCD = True; (* use two-loop threshold for yt *)
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
 */
diff --git a/install-sarah b/install-sarah
@@ -13,7 +13,7 @@ ABSBASEDIR=$(cd $BASEDIR; pwd)
 operating_system="`(uname -s) 2>/dev/null || echo unknown`"
 
 # default SARAH version
-sarah_version="4.9.3"
+sarah_version="4.11.0"
 
 # determine Mathematica user directory
 case "$operating_system" in

diff --git a/meta/EffectiveCouplings.m b/meta/EffectiveCouplings.m
@@ -476,13 +476,13 @@ previous results (e.g. define along the lines of f[p] := f[p] = ...)    *)
 
            result = "const double scale = model.get_scale();\nconst Eigen::ArrayXd saved_parameters(model.get());\n\n"
                     <> "const double saved_mt = PHYSICAL("
-                    <> CConversion`RValueToCFormString[TreeMasses`GetThirdGenerationMass[SARAH`TopQuark]]
+                    <> CConversion`RValueToCFormString[TreeMasses`GetThirdGenerationMass[TreeMasses`GetSMTopQuarkMultiplet[]]]
                     <> ");\nPHYSICAL("
-                    <> CConversion`RValueToCFormString[TreeMasses`GetThirdGenerationMass[SARAH`TopQuark]]
+                    <> CConversion`RValueToCFormString[TreeMasses`GetThirdGenerationMass[TreeMasses`GetSMTopQuarkMultiplet[]]]
                     <> ") = qedqcd.displayPoleMt();\n\n"
                     <> result;
            result = result <> "PHYSICAL("
-                           <> CConversion`RValueToCFormString[TreeMasses`GetThirdGenerationMass[SARAH`TopQuark]]
+                           <> CConversion`RValueToCFormString[TreeMasses`GetThirdGenerationMass[TreeMasses`GetSMTopQuarkMultiplet[]]]
                            <> ") = saved_mt;\n";
            result = result <> "model.set_scale(scale);\nmodel.set(saved_parameters);\n";
 

diff --git a/meta/FlexibleSUSY.m b/meta/FlexibleSUSY.m
@@ -135,6 +135,7 @@ FlexibleSUSY model file (FlexibleSUSY.m).
 EffectiveMASqr;
 UseSM3LoopRGEs = False;
 UseMSSM3LoopRGEs = False;
+UseMSSMYukawa3LoopSQCD = False;
 UseHiggs2LoopSM = False;
 UseHiggs3LoopSplit = False;
 UseYukawa3LoopQCD = Automatic;
@@ -1379,6 +1380,7 @@ corresponding tadpole is real or imaginary (only in models with CP
             reorderDRbarMasses = "", reorderPoleMasses = "",
             checkPoleMassesForTachyons = "",
             twoLoopHiggsHeaders = "", threeLoopHiggsHeaders = "",
+            twoLoopThresholdHeaders = "",
             lspGetters = "", lspFunctions = "",
             convertMixingsToSLHAConvention = "",
             convertMixingsToHKConvention = "",
@@ -1417,7 +1419,6 @@ corresponding tadpole is real or imaginary (only in models with CP
            divideTadpoleByVEV           = SelfEnergies`DivideTadpoleByVEV[Parameters`DecreaseIndexLiterals @ CreateVEVToTadpoleAssociation[], "tadpole"];
            If[SARAH`UseHiggs2LoopMSSM === True || FlexibleSUSY`UseHiggs2LoopNMSSM === True,
               calculateTwoLoopTadpoles  = SelfEnergies`FillArrayWithTwoLoopTadpoles[SARAH`HiggsBoson, "tadpole", "-"];
-              {thirdGenerationHelperPrototypes, thirdGenerationHelperFunctions} = TreeMasses`CreateThirdGenerationHelpers[];
              ];
            If[FlexibleSUSY`UseHiggs2LoopSM === True,
               {twoLoopSelfEnergyPrototypes, twoLoopSelfEnergyFunctions} = SelfEnergies`CreateTwoLoopSelfEnergiesSM[{SARAH`HiggsBoson}];
@@ -1437,6 +1438,14 @@ corresponding tadpole is real or imaginary (only in models with CP
               {twoLoopSelfEnergyPrototypes, twoLoopSelfEnergyFunctions} = SelfEnergies`CreateTwoLoopSelfEnergiesNMSSM[{SARAH`HiggsBoson, SARAH`PseudoScalar}];
               twoLoopHiggsHeaders = "#include \"sfermions.hpp\"\n#include \"mssm_twoloophiggs.hpp\"\n#include \"nmssm_twoloophiggs.hpp\"\n";
              ];
+           If[FlexibleSUSY`UseMSSMYukawa3LoopSQCD === True,
+              twoLoopThresholdHeaders = "#include \"mssm_twoloop_mt.hpp\"";
+             ];
+           If[SARAH`UseHiggs2LoopMSSM === True ||
+              FlexibleSUSY`UseHiggs2LoopNMSSM === True ||
+              FlexibleSUSY`UseMSSMYukawa3LoopSQCD === True,
+              {thirdGenerationHelperPrototypes, thirdGenerationHelperFunctions} = TreeMasses`CreateThirdGenerationHelpers[];
+             ];
            {selfEnergyPrototypes, selfEnergyFunctions} = SelfEnergies`CreateNPointFunctions[nPointFunctions, vertexRules];
            phasesDefinition             = Phases`CreatePhasesDefinition[phases];
            phasesGetterSetters          = Phases`CreatePhasesGetterSetters[phases];
@@ -1541,6 +1550,7 @@ corresponding tadpole is real or imaginary (only in models with CP
                             "@twoLoopSelfEnergyPrototypes@" -> IndentText[twoLoopSelfEnergyPrototypes],
                             "@twoLoopSelfEnergyFunctions@"  -> twoLoopSelfEnergyFunctions,
                             "@twoLoopHiggsHeaders@"       -> twoLoopHiggsHeaders,
+                            "@twoLoopThresholdHeaders@"   -> twoLoopThresholdHeaders,
                             "@threeLoopSelfEnergyPrototypes@" -> IndentText[threeLoopSelfEnergyPrototypes],
                             "@threeLoopSelfEnergyFunctions@"  -> threeLoopSelfEnergyFunctions,
                             "@threeLoopHiggsHeaders@"         -> threeLoopHiggsHeaders,

diff --git a/meta/GMuonMinus2.m b/meta/GMuonMinus2.m
@@ -32,7 +32,7 @@ If you add new kinds of vertices (e.g for new diagram types):
 (************* Begin public interface *******************)
 
 GMuonMinus2IsSMParticle[particle_] :=
-    TreeMasses`IsSMParticle[particle] || TreeMasses`IsSMGoldstone[particle];
+    TreeMasses`IsSMParticle[particle];
 
 CreateSMParticleFlags[particle_] :=
     Module[{result = "", i,
@@ -273,14 +273,14 @@ If you add new kinds of vertices (e.g for new diagram types):
 (************************ Begin helper routines *******************************)
 
 (* Return the name of the SARAH particle family containing the muon *)
-GetMuonFamily[] := If[TreeMasses`GetDimension[SARAH`Electron] =!= 1,
-                        SARAH`Electron,
+GetMuonFamily[] := If[TreeMasses`GetDimension[TreeMasses`GetSMMuonLeptonMultiplet[]] =!= 1,
+                        TreeMasses`GetSMMuonLeptonMultiplet[],
                         Cases[SARAH`ParticleDefinitions[FlexibleSUSY`FSEigenstates],
                               {p_, {Description -> "Muon", ___}} -> p, 1][[1]]
                         ];
 (* If the muon has a generation index, return it, otherwise return Null.
  e.g. CMSSMNoFV has no muon generation index *)
-GetMuonIndex[] := If[TreeMasses`GetDimension[SARAH`Electron] =!= 1, 2, Null];
+GetMuonIndex[] := If[TreeMasses`GetDimension[TreeMasses`GetSMMuonLeptonMultiplet[]] =!= 1, 2, Null];
 
 GetPhoton[] := SARAH`Photon;
 

diff --git a/meta/LoopMasses.m b/meta/LoopMasses.m
@@ -56,23 +56,64 @@
            Return[result];
           ];
 
-AddMtPoleQCDCorrections[1, expr_] := "\
+FillMt2LStruct[] := "\
+double mst_1, mst_2, theta_t;
+" <> TreeMasses`CallThirdGenerationHelperFunctionName[SARAH`TopSquark, "mst_1", "mst_2", "theta_t"] <> ";
+
+mssm_twoloop_mt::Parameters pars;
+pars.g3 = " <> CConversion`RValueToCFormString[SARAH`strongCoupling /. Parameters`ApplyGUTNormalization[]] <> ";
+pars.mt = " <> CConversion`RValueToCFormString[TreeMasses`GetThirdGenerationMass[TreeMasses`GetSMTopQuarkMultiplet[]]] <> ";
+pars.mg = " <> CConversion`RValueToCFormString[FlexibleSUSY`M[SARAH`Gluino]] <> ";
+pars.mst1 = mst_1;
+pars.mst2 = mst_2;
+pars.msusy = " <> CConversion`RValueToCFormString[Sqrt[Sqrt[SARAH`SoftSquark[2,2] SARAH`SoftDown[2,2]]]] <> ";
+pars.xt = Sin(2*theta_t) * (Sqr(mst_1) - Sqr(mst_2)) / (2. * pars.mt);
+pars.Q = get_scale();";
+
+
+AddMtPoleQCDCorrections[1, expr_] :=
+    If[FlexibleSUSY`UseMSSMYukawa3LoopSQCD === True,
+"\
+double qcd_1l = 0.;
+
+{
+" <> IndentText[FillMt2LStruct[]] <> "
+
+   qcd_1l = - mssm_twoloop_mt::dMt_over_mt_1loop(pars);
+}
+"
+       ,
+"\
 double qcd_1l = 0.;
 
 {
    const double currentScale = get_scale();
    qcd_1l = " <> CConversion`RValueToCFormString[expr] <> ";
 }
-";
+"
+    ];
+
+AddMtPoleQCDCorrections[2, expr_] :=
+    If[FlexibleSUSY`UseMSSMYukawa3LoopSQCD === True,
+"\
+double qcd_2l = 0.;
+
+if (pole_mass_loop_order > 1 && TOP_POLE_QCD_CORRECTION > 0) {
+" <> IndentText[FillMt2LStruct[]] <> "
 
-AddMtPoleQCDCorrections[2, expr_] := "\
+   qcd_2l = - mssm_twoloop_mt::dMt_over_mt_2loop(pars);
+}
+"
+       ,
+ "\
 double qcd_2l = 0.;
 
 if (pole_mass_loop_order > 1 && TOP_POLE_QCD_CORRECTION > 0) {
    const double currentScale = get_scale();
    qcd_2l = " <> CConversion`RValueToCFormString[expr] <> ";
 }
-";
+"
+    ];
 
 AddMtPoleQCDCorrections[3, expr_] := "\
 double qcd_3l = 0.;
@@ -108,7 +149,7 @@
     "PHYSICAL(" <> ToValidCSymbolString[FlexibleSUSY`M[particle]] <> ") = " <>
     "calculate_singlet_mass(M_loop);\n";
 
-Do1DimFermion[particle_ /; particle === SARAH`TopQuark, massMatrixName_String,
+Do1DimFermion[particle_ /; particle === TreeMasses`GetSMTopQuarkMultiplet[], massMatrixName_String,
               _String, _String, _String, momentum_String, type_] :=
     Module[{massName,
             topSelfEnergyFunctionS, topSelfEnergyFunctionPL, topSelfEnergyFunctionPR,
@@ -467,7 +508,7 @@ be set to the (positive) tree-level mass.  M_tree
            selfEnergyMatrixType = TreeMasses`GetMassMatrixType[particle];
            selfEnergyMatrixCType = CreateCType[selfEnergyMatrixType];
            eigenArrayType = CreateCType[CConversion`ArrayType[CConversion`realScalarCType, dim]];
-           topTwoLoop = particle === SARAH`TopQuark;
+           topTwoLoop = particle === TreeMasses`GetSMTopQuarkMultiplet[];
            If[topTwoLoop,
               topSelfEnergyFunctionS  = SelfEnergies`CreateHeavySelfEnergyFunctionName[particle[1], 1];
               topSelfEnergyFunctionPL = SelfEnergies`CreateHeavySelfEnergyFunctionName[particle[PL], 1];
@@ -872,7 +913,7 @@ be set to the (positive) tree-level mass.  M_tree
            Return[result];
           ];
 
-CreateRunningDRbarMassFunction[particle_ /; particle === SARAH`BottomQuark, renormalizationScheme_] :=
+CreateRunningDRbarMassFunction[particle_ /; particle === TreeMasses`GetSMBottomQuarkMultiplet[], renormalizationScheme_] :=
     Module[{result, body, selfEnergyFunctionS, selfEnergyFunctionPL,
             selfEnergyFunctionPR, name, alphaS, drbarConversion, gPrime,
             dimParticle, treeLevelMass},
@@ -959,7 +1000,21 @@ be set to the (positive) tree-level mass.  M_tree
            Return[result <> IndentText[body] <> "}\n\n"];
           ];
 
-CreateRunningDRbarMassFunction[particle_ /; particle === SARAH`TopQuark, _] :=
+CreateMSSM1LoopSQCDContributions[result_:"qcd_1l"] := "\
+{
+" <> IndentText[FillMt2LStruct[]] <> "
+
+   " <> result <> " = - mssm_twoloop_mt::dMt_over_mt_1loop(pars);
+}";
+
+CreateMSSM2LoopSQCDContributions[result_:"qcd_2l"] :=
+    FillMt2LStruct[] <> "
+
+const double q_2l = mssm_twoloop_mt::dMt_over_mt_2loop(pars);
+
+" <> result <> " = -q_2l + qcd_1l * qcd_1l;";
+
+CreateRunningDRbarMassFunction[particle_ /; particle === TreeMasses`GetSMTopQuarkMultiplet[], _] :=
     Module[{result, body, selfEnergyFunctionS, selfEnergyFunctionPL,
             selfEnergyFunctionPR, name, qcdOneLoop, qcdTwoLoop, qcdThreeLoop = 0,
             dimParticle, treeLevelMass},
@@ -1004,10 +1059,19 @@ be set to the (positive) tree-level mass.  M_tree
                 ];
               body = body <>
               "const double currentScale = get_scale();\n" <>
-              "const double qcd_1l = " <> CConversion`RValueToCFormString[qcdOneLoop /. FlexibleSUSY`M[particle] -> treeLevelMass] <> ";\n" <>
-              "double qcd_2l = 0., qcd_3l = 0.;\n\n" <>
+              "double qcd_1l = 0., qcd_2l = 0., qcd_3l = 0.;\n\n" <>
+                  If[FlexibleSUSY`UseMSSMYukawa3LoopSQCD === True,
+                     CreateMSSM1LoopSQCDContributions[],
+                     "qcd_1l = " <> CConversion`RValueToCFormString[qcdOneLoop /. FlexibleSUSY`M[particle] -> treeLevelMass] <> ";"
+                    ] <>
+              "\n\n" <>
               "if (get_thresholds() > 1) {\n" <>
-              IndentText["qcd_2l = " <> CConversion`RValueToCFormString[qcdTwoLoop /. FlexibleSUSY`M[particle] -> treeLevelMass] <> ";"] <> "\n" <>
+              IndentText[
+                  If[FlexibleSUSY`UseMSSMYukawa3LoopSQCD === True,
+                     CreateMSSM2LoopSQCDContributions[],
+                     "qcd_2l = " <> CConversion`RValueToCFormString[qcdTwoLoop /. FlexibleSUSY`M[particle] -> treeLevelMass] <> ";"
+                    ]
+              ] <> "\n" <>
               "}\n\n" <>
               If[qcdThreeLoop =!= 0,
               "if (get_thresholds() > 2) {\n" <>

diff --git a/meta/MSSM/dmtas2.expr b/meta/MSSM/dmtas2.expr