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

Conflicts:
	templates/librarylink.cpp.in
	templates/standard_model_high_scale_spectrum_generator.hpp.in
	templates/standard_model_low_scale_spectrum_generator.hpp.in
	templates/two_scale_high_scale_spectrum_generator.cpp.in
	templates/two_scale_low_scale_spectrum_generator.cpp.in

ChangeLog

@@ -115,6 +115,9 @@ FlexibleSUSY-1.7.3 [not released yet]
  * Bugfix: Implement missing limits of threshold correction functions
    from arXiv:1407.4081.
 
+ * Bugfix (commit 581080f): Catch further NaNs from inside the MSSM 2L
+   Higgs mass routines of Pietro Slavich.
+
 FlexibleSUSY-1.7.2 [December, 15 2016]
 
  * Feature (commit b052e35): New flag FlexibleSUSY[23] to disable the

README

@@ -39,8 +39,8 @@ following setup:
 
    $ cd ~/.Mathematica/Applications/
    $ wget https://www.hepforge.org/archive/sarah/SARAH-4.7.0.tar.gz
-   $ tar -xf SARAH-4.9.1.tar.gz
-   $ ln -s $PWD/SARAH-4.9.1/ SARAH
+   $ tar -xf SARAH-4.9.3.tar.gz
+   $ ln -s $PWD/SARAH-4.9.3/ SARAH
 
    $ cd ~/.Mathematica/Kernel/
    $ echo "AppendTo[\$Path, \"$HOME/.Mathematica/Applications/SARAH/\"];" \

doc/install.dox

@@ -8,7 +8,7 @@ FlexibleSUSY requires SARAH to be installed and to be loadable with
 the ``Needs["SARAH`"]`` command from inside Mathematica.  We recommend
 the following setup:
 
-    VERSION=4.9.1
+    VERSION=4.9.3
     cd ~/.Mathematica/Applications/
     wget https://www.hepforge.org/archive/sarah/SARAH-${VERSION}.tar.gz
     tar -xf SARAH-${VERSION}.tar.gz

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.1"
+sarah_version="4.9.3"
 
 # determine Mathematica user directory
 case "$operating_system" in

meta/FlexibleSUSY.m

@@ -1563,7 +1563,7 @@ corresponding tadpole is real or imaginary (only in models with CP
                   _, Print["Error: invalid BVP solver requested: ", solver];
                      Quit[1];
                  ];
-           body = "data.reset(new " <> FlexibleSUSY`FSModelName <> "_model_data<" <> class <> ">());\nbreak;\n";
+           body = "data.reset(new Model_data_impl<" <> class <> ">());\nbreak;\n";
            result = "case " <> key <> ":\n" <> IndentText[body];
            EnableForBVPSolver[solver, IndentText[result]] <> "\n"
           ];

meta/Parameters.m

@@ -1137,12 +1137,11 @@
                                           allModelParameters, allOutputParameters]];
 
 IncreaseIndexLiterals[expr_, num_Integer, heads_List] :=
-    Module[{indexedSymbols, rules, decrExpr, allHeads},
+    Module[{indexedSymbols, rules, allHeads},
            allHeads = Join[heads /. FlexibleSUSY`M -> Identity, {SARAH`Delta, SARAH`ThetaStep}];
-           indexedSymbols = Cases[{expr}, s_[__] /; MemberQ[allHeads, s], Infinity];
+           indexedSymbols = Extract[{expr}, Position[{expr}, s_[__] /; MemberQ[allHeads, s], Infinity]];
            rules = Rule[#, IncreaseIndices[#,num]] & /@ indexedSymbols;
-           decrExpr = expr /. rules;
-           Return[decrExpr]
+           expr /. rules
           ];
 
 DecreaseIndexLiterals[expr_] :=

meta/SelfEnergies.m

@@ -253,6 +253,10 @@
            symbol[Sequence @@ indices]
           ];
 
+indexCount = 0;
+MakeUniqueIdx[] :=
+    Symbol["id" <> ToString[indexCount++]];
+
 (* creates a C++ function that calculates a coupling
  *
  * Return: {prototypes_String, definitions_String, rules_List}
@@ -333,10 +337,12 @@
 
 CreateVertexExpressions[vertexRules_List] :=
     Module[{k, prototypes = "", defs = "", rules, coupling, expr,
-            p, d, r},
+            p, d, r, MakeIndex},
+           MakeIndex[i_Integer] := MakeUniqueIdx[];
+           MakeIndex[i_] := i;
            rules = Table[0, {Length[vertexRules]}];
            For[k = 1, k <= Length[vertexRules], k++,
-               coupling = Vertices`ToCp[vertexRules[[k,1]]];
+               coupling = Vertices`ToCp[vertexRules[[k,1]]] /. p_[{idx__}] :> p[MakeIndex /@ {idx}];
                expr = vertexRules[[k,2]];
                WriteString["stdout", "."];
                If[Mod[k, 50] == 0, WriteString["stdout","\n"]];

meta/TreeMasses.m

@@ -90,6 +90,9 @@
 GetParticleIndices::usage = "returns list of particle indices with
  names";
 
+ParticleQ::usage = "returns True if argument is a particle, False
+ otherwise."
+
 FindMixingMatrixSymbolFor::usage="returns the mixing matrix symbol for
 a given field";
 
@@ -226,6 +229,9 @@
 GetSMParticles[states_:FlexibleSUSY`FSEigenstates] :=
     Select[GetParticles[states], (SARAH`SMQ[#])&];
 
+ParticleQ[p_, states_:FlexibleSUSY`FSEigenstates] :=
+    MemberQ[GetParticles[states], p];
+
 IsOfType[sym_Symbol, type_Symbol, states_:FlexibleSUSY`FSEigenstates] :=
     SARAH`getType[sym, False, states] === type;
 

meta/TwoLoopSM.m

@@ -0,0 +1,85 @@
+BeginPackage["TwoLoopSM`"];
+EndPackage[];
+
+GetSMHiggsMass::usage = "Returns the loop corrections to the Higgs
+ mass in the Standard Model.  The loop-corrections are taken from
+ arxiv:1205.6497 and arxiv:1504.05200.
+
+Note: The return value contains the contributions from tadpole
+ diagrams.
+
+Usage: GetSMHiggsMass[loopOrder -> {1,1,1}, corrections -> {1,1}, simplifications -> { p -> 0 }]
+
+Parameters:
+
+- loopOrder (optional): List of factors multiplied by each loop order.
+  #1: tree-level
+  #2: 1-loop level
+  #3: 2-loop level
+  (default: {1,1,1})
+
+- corrections (optional): List of factors multiplied by each 2-loop
+  correction. (default: {1, 1})
+  #1: alpha_t * alpha_s (arxiv:1205.6497, Eq. (20))
+  #2: alpha_t^2         (arxiv:1504.05200, Eq. (20))
+
+- simplifications (optional): List of replacement rules to simplify
+  the result (default: { p -> 0, B0[0,m12_,m22_,Q2_] :>
+     B0zero[Sqrt[m12],Sqrt[m22],Sqrt[Q2]] }).
+  Example: simplifications -> {} (* no simplifications *)
+";
+
+(* MS-bar parameters *)
+{ p, yt, mt, g3, Q };
+
+(* loop functions *)
+{ B0 };
+
+(* options *)
+{ loopOrder, corrections };
+
+Begin["TwoLoopSM`Private`"];
+
+(* B0 for p = 0 *)
+B0zero[m1_, m2_, mu_] :=
+    Which[PossibleZeroQ[m1 - m2],
+          Log[mu^2/m2^2],
+          PossibleZeroQ[m1],
+          1 + Log[mu^2/m2^2],
+          PossibleZeroQ[m2],
+          1 + Log[mu^2/m1^2],
+          True,
+          1 + (m1^2 Log[mu^2/m1^2] - m2^2 Log[mu^2/m2^2])/(m1^2 - m2^2)
+   ];
+
+Options[GetSMHiggsMass] = {
+    loopOrder -> {1,1,1},
+    corrections -> {1,1},
+    simplifications -> { p -> 0,
+                         B0[0,m12_,m22_,Q2_] :> B0zero[Sqrt[m12],Sqrt[m22],Sqrt[Q2]] }
+};
+
+GetSMHiggsMass1LAlphaT[] :=
+    Module[{h = 1/(4 Pi)^2},
+           h 3 yt^2 (4 mt^2 - p^2) B0[p^2,mt^2,mt^2,Q^2]
+          ];
+
+GetSMHiggsMass2LAlphaTAlphaS[] :=
+    Module[{LogT = Log[mt^2/Q^2], h = 1/(4 Pi)^2},
+           h^2 2 mt^2 16 g3^2 yt^2 (3 LogT^2 + LogT)
+          ];
+
+GetSMHiggsMass2LAlphaTAlphaT[] :=
+    Module[{LogT = Log[mt^2/Q^2], h = 1/(4 Pi)^2},
+           h^2 2 mt^2 (-3 yt^4 (3 LogT^2 - 7 LogT + 2 + Pi^2/3))
+          ];
+
+GetSMHiggsMass[OptionsPattern[]] :=
+    (
+        OptionValue[loopOrder][[1]] 0 +
+        OptionValue[loopOrder][[2]] GetSMHiggsMass1LAlphaT[] +
+        OptionValue[loopOrder][[3]] OptionValue[corrections][[1]] GetSMHiggsMass2LAlphaTAlphaS[] +
+        OptionValue[loopOrder][[3]] OptionValue[corrections][[2]] GetSMHiggsMass2LAlphaTAlphaT[]
+    ) //. OptionValue[simplifications];
+
+End[];