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

Conflicts: meta/FlexibleSUSY.m meta/Parameters.m
FlexibleSUSY · May 16, 2017 · 92599e6 · 92599e6
2 parents b7ebccf + 5192029
commit 92599e6
Show file tree

Hide file tree

Showing 15 changed files with 219 additions and 74 deletions.
diff --git a/AUTHORS b/AUTHORS
@@ -1,14 +1,12 @@
-Core developers
-===============
+FlexibleSUSY authors
+====================
 
 Peter Athron            peter.athron@coepp.org.au
+Markus Bach             Markus.Bach1@mailbox.tu-dresden.de
+Dylan Harries           dylan.harries@adelaide.edu.au
+Thomas Kwasnitza
 Jae-hyeon Park          jhpark@kias.re.kr
-Dominik Stöckinger      Dominik.Stoeckinger@tu-dresden.de
-Alexander Voigt         Alexander.Voigt@desy.de
-
-Contributors
-============
-
-Dylan Harries
 Tom Steudtner
+Dominik Stöckinger      Dominik.Stoeckinger@tu-dresden.de
+Alexander Voigt         Alexander.Voigt@physik.rwth-aachen.de
 Jobst Ziebell
diff --git a/ChangeLog b/ChangeLog
@@ -142,6 +142,11 @@ FlexibleSUSY 2.0.0 [not released yet]
 
    Thanks to Peter Athron and Pat Scott.
 
+ * Bugfix (commit 0cb4042): Matrix products of the form A*B in the
+   boundary conditions are now interpreted as element-wise products
+   (as in Mathematica).  For genuine matrix products use MatMul[A,B]
+   or A.B.
+
 FlexibleSUSY-1.7.4 [April, 12 2017]
 
  * Bugfix (commit f434e30): Rename internal IndexSum symbol which

diff --git a/addons/GM2Calc/MSSMNoFV_onshell.cpp b/addons/GM2Calc/MSSMNoFV_onshell.cpp
@@ -638,7 +638,7 @@ double MSSMNoFV_onshell::convert_me2_root(
                 << precision << " GeV\n";
    }
 
-   return precision;
+   return std::isfinite(precision) ? precision : 0.;
 }
 
 /**
@@ -722,7 +722,7 @@ double MSSMNoFV_onshell::convert_me2_fpi(
                 << precision << " GeV\n";
    }
 
-   return precision;
+   return std::isfinite(precision) ? precision : 0.;
 }
 
 std::ostream& operator<<(std::ostream& os, const MSSMNoFV_onshell& model)

diff --git a/doc/mainpage.dox.in b/doc/mainpage.dox.in
@@ -53,12 +53,23 @@ the shell
 
 \section References
 
- - FlexibleSUSY — A spectrum generator generator for supersymmetric
-   models
-   [<a href="http://arxiv.org/abs/1406.2319">CPC 190 (2015) 139-172</a>]
+ - Peter Athron, Jae-hyeon Park, Dominik Stöckinger, Alexander Voigt:
+   "FlexibleSUSY — A spectrum generator generator for supersymmetric
+   models" [<a href="https://arxiv.org/abs/1406.2319">CPC 190 (2015)
+   139-172</a>]
 
-If you use FlexibleSUSY in your work please cite the <a
-  href="http://inspirehep.net/record/1299998">FlexibleSUSY manual</a>.
+ - Peter Athron, Jae-hyeon Park, Tom Steudtner, Dominik Stöckinger,
+   Alexander Voigt: "Precise Higgs mass calculations in (non-)minimal
+   supersymmetry at both high and low scales" [<a
+   href="https://arxiv.org/abs/1609.00371">JHEP 1701 (2017) 079</a>]
+
+If you use FlexibleSUSY in your work please cite [<a
+href="http://inspirehep.net/record/1299998">CPC 190 (2015)
+139-172</a>].
+
+If you use the FlexibleEFTHiggs approach in your work please cite [<a
+href="https://inspirehep.net/record/1484857">JHEP 1701 (2017)
+079</a>].
 
 FlexibleSUSY depends on SARAH and contains components from
 Softsusy. Therefore, please also cite the following publications along
@@ -79,17 +90,13 @@ with FlexibleSUSY:
 FlexibleSUSY has been written by
 
 - Peter Athron
+- Markus Bach
+- Dylan Harries
+- Thomas Kwasnitza
 - Jae-hyeon Park
 - Dominik Stöckinger
 - Alexander Voigt
-
-The following people have made major contributions to FlexibleSUSY:
-
-- Dylan Harries (loop-induced effective vertices
-  \f$h\rightarrow\gamma\gamma\f$ and \f$h\rightarrow gg\f$)
-- Tom Steudtner (FlexibleEFTHiggs)
-- Jobst Ziebell (1-loop BSM contributions to the anomalous magnetic
-  moment of the muon, \f$a_\mu\f$)
+- Jobst Ziebell
 
 \section fs_resources Resources
 

diff --git a/meta/BetaFunction.m b/meta/BetaFunction.m
@@ -43,10 +43,16 @@
 
 GetAllBetaFunctions[BetaFunction[name_, type_, beta_List]] := beta;
 
-GuessType[sym_[Susyno`LieGroups`i1, SARAH`i2]] :=
+betaIndices = {
+    Susyno`LieGroups`i1 , SARAH`i2 , SARAH`i3 , SARAH`i4
+};
+
+IsBetaIdx[i_] := MemberQ[betaIndices, i];
+
+GuessType[sym_[_?IsBetaIdx, _?IsBetaIdx]] :=
     Parameters`GetType[sym];
 
-GuessType[sym_[Susyno`LieGroups`i1]] :=
+GuessType[sym_[_?IsBetaIdx]] :=
     Parameters`GetType[sym];
 
 GuessType[sym_] :=

diff --git a/meta/CConversion.m b/meta/CConversion.m
@@ -42,6 +42,8 @@
 AbsSqrt::usage="";
 FSKroneckerDelta::usage="";
 TensorProd::usage="";
+ElementwiseProd::usage = "represents an elemen-wise product of
+ vectors, matrices or tensors";
 
 HaveSameDimension::usage = "Checks if given types have same
 dimension";
@@ -768,6 +770,16 @@
     Format["(" <> ToString[CForm[HoldForm[x]]] <> ")*(" <>
            ToString[CForm[HoldForm[y]]] <> ").transpose()", OutputForm];
 
+Format[CConversion`ElementwiseProd[HoldPattern[fac_]],CForm] :=
+    Format[ToString[CForm[HoldForm[fac]]], OutputForm];
+
+Format[CConversion`ElementwiseProd[HoldPattern[fac__]],CForm] :=
+    Format[
+        "(" <> StringJoin[
+            Riffle[ToString[CForm[HoldForm[#]]]& /@ {fac}, ").cwiseProduct("]
+         ] <> ")"
+        , OutputForm];
+
 (* Finds all Greek symbols in an expression.
    Note: All arguments of Which and If are evaluated.
  *)
@@ -777,6 +789,21 @@
               Cases[expr, x_Symbol | x_Symbol[__] :> x, {0,Infinity}, Heads->True], GreekQ]
          ];
 
+CConversion`ElementwiseProd[] = 1;
+
+(* finds all vectors, matrices and tensors in a given list *)
+FindMatrices[syms_List] :=
+    Select[syms,
+           Parameters`IsParameter[#] &&
+           Length[Parameters`GetParameterDimensions[#]] > 1 &];
+
+FactorElementwiseProd[fac_] :=
+    Module[{facs, mat},
+           facs = If[Head[fac] === Times, List @@ fac, {fac}];
+           mat = FindMatrices[facs];
+           Times[Sequence @@ Complement[facs, mat]] CConversion`ElementwiseProd[Sequence @@ mat]
+          ];
+
 (* Converts an expression to CForm and expands SARAH symbols
  *
  *   MatMul[A]      ->   A
@@ -813,6 +840,7 @@
                     FlexibleSUSY`BETA[l_,p_]     :> FlexibleSUSY`BETA1[l,p] /.
                     Susyno`LieGroups`conj    -> SARAH`Conj //.
                     conjSimplification /.
+                    Times[fac__] :> FactorElementwiseProd[fac] /.
                     SARAH`Delta[a_,a_]       -> 1 /.
                     Power[a_?NumericQ,n_?NumericQ] :> N[Power[a,n]] /.
                     Sqrt[a_?NumericQ]        :> N[Sqrt[a]] /.

diff --git a/meta/FlexibleSUSY.m b/meta/FlexibleSUSY.m
@@ -34,17 +34,21 @@
 
 FS`Version = StringTrim[FSImportString[FileNameJoin[{$flexiblesusyConfigDir,"version"}]]];
 FS`GitCommit = StringTrim[FSImportString[FileNameJoin[{$flexiblesusyConfigDir,"git_commit"}]]];
-FS`Authors = {"P. Athron", "J.-h. Park", "D. Stöckinger", "A. Voigt"};
-FS`Contributors = {"D. Harries", "T. Steudtner", "J. Ziebell"};
+FS`Authors = {"P. Athron", "T. Kwasnitza", "D. Harries",
+              "J.-h. Park", "T. Steudtner", "D. Stöckinger",
+              "A. Voigt", "J. Ziebell"};
+FS`Contributors = {};
 FS`Years   = "2013-2017";
 FS`References = Get[FileNameJoin[{$flexiblesusyConfigDir,"references"}]];
 
 Print[""];
 Utils`FSFancyLine["="];
 Utils`FSFancyPrint["FlexibleSUSY " <> FS`Version, 0];
-Print["  by " <> Utils`StringJoinWithSeparator[FS`Authors, ", "] <> ", " <>
-      FS`Years];
-Print["  contributions by " <> Utils`StringJoinWithSeparator[FS`Contributors, ", "]];
+Print["  by " <> StringJoin[Riffle[Riffle[FS`Authors, ", "], "\n  ", 11]] <>
+      "\n  " <> FS`Years];
+If[FS`Contributors =!= {},
+   Print["  contributions by " <> Utils`StringJoinWithSeparator[FS`Contributors, ", "]];
+  ];
 Print[""];
 Utils`FSFancyPrint["References:"];
 Print["  " <> #]& /@ FS`References;
@@ -141,8 +145,16 @@ FlexibleSUSY model file (FlexibleSUSY.m).
 UseYukawa3LoopQCD = Automatic;
 FSRGELoopOrder = 2; (* RGE loop order (0, 1 or 2) *)
 PotentialLSPParticles = {};
-ExtraSLHAOutputBlocks = {};
-FSExtraInputParameters;
+ExtraSLHAOutputBlocks = {
+    {FlexibleSUSYLowEnergy,
+        {{1, FlexibleSUSYObservable`aMuon} } },
+    {EFFHIGGSCOUPLINGS, NoScale,
+        {{1, FlexibleSUSYObservable`CpHiggsPhotonPhoton},
+         {2, FlexibleSUSYObservable`CpHiggsGluonGluon},
+         {3, FlexibleSUSYObservable`CpPseudoScalarPhotonPhoton},
+         {4, FlexibleSUSYObservable`CpPseudoScalarGluonGluon} } },
+};
+FSExtraInputParameters = {};
 FSAuxiliaryParameterInfo = {};
 IMMINPAR = {};
 IMEXTPAR = {};

diff --git a/meta/LoopMasses.m b/meta/LoopMasses.m
@@ -425,7 +425,7 @@ be set to the (positive) tree-level mass.  M_tree
 
 (* ********** medium diagonalization routines ********** *)
 
-DoMediumDiagonalization[particle_Symbol /; IsScalar[particle], inputMomentum_, tadpole_List] :=
+DoMediumDiagonalization[particle_Symbol /; IsScalar[particle], inputMomentum_, tadpole_List, calc2L_:True] :=
     Module[{result, dim, dimStr, massName, particleName, mixingMatrix, selfEnergyFunction,
             momentum = inputMomentum, U, V, Utemp, Vtemp, tadpoleMatrix, diagSnippet,
             massMatrixStr,
@@ -480,23 +480,9 @@ be set to the (positive) tree-level mass.  M_tree
                   FlexibleSUSY`UseHiggs2LoopNMSSM === True) &&
                  MemberQ[{SARAH`HiggsBoson, SARAH`PseudoScalar}, particle],
                  addHigherLoopHiggsContributions = "self_energy += self_energy_2l;\n";
-                 calcHigherLoopHiggsContributions = "
-// two-loop Higgs self-energy contributions
-" <> selfEnergyMatrixCType <> " self_energy_2l(" <> selfEnergyMatrixCType <> "::Zero());
-
-if (pole_mass_loop_order > 1) {
-" <> IndentText["\
-self_energy_2l = self_energy_" <> CConversion`ToValidCSymbolString[particle] <> "_2loop();
-for (int i = 0; i < " <> dimStr <> "; i++) {
-   for (int k = 0; k < " <> dimStr <> "; k++) {
-      if (!std::isfinite(self_energy_2l(i,k))) {
-         self_energy_2l(i,k) = 0.;
-         problems.flag_bad_mass(" <> FlexibleSUSY`FSModelName <> "_info::" <> CConversion`ToValidCSymbolString[particle] <> ");
-      }
-   }
-}
-"] <> "}
-";
+                 If[calc2L,
+                    calcHigherLoopHiggsContributions = CalcEffPot2L[particle];
+                   ];
                 ];
               result = tadpoleMatrix <>
                        "const " <> selfEnergyMatrixCType <> " M_tree(" <> massMatrixStr <> "());\n" <>
@@ -520,7 +506,7 @@ be set to the (positive) tree-level mass.  M_tree
            Return[result];
           ];
 
-DoMediumDiagonalization[particle_Symbol /; IsFermion[particle], inputMomentum_, _] :=
+DoMediumDiagonalization[particle_Symbol /; IsFermion[particle], inputMomentum_, __] :=
     Module[{result, dim, dimStr, massName, mixingMatrix, U, V,
             selfEnergyFunctionS, selfEnergyFunctionPL, selfEnergyFunctionPR,
             momentum = inputMomentum, massMatrixStr,
@@ -675,7 +661,7 @@ be set to the (positive) tree-level mass.  M_tree
            Return[result];
           ];
 
-DoMediumDiagonalization[particle_Symbol /; IsVector[particle], inputMomentum_, _] :=
+DoMediumDiagonalization[particle_Symbol /; IsVector[particle], inputMomentum_, __] :=
     Module[{result, dim, dimStr, massName, particleName, mixingMatrix, selfEnergyFunction,
             momentum = inputMomentum, selfEnergyMatrixType, selfEnergyMatrixCType},
            dim = GetDimension[particle];
@@ -702,22 +688,52 @@ be set to the (positive) tree-level mass.  M_tree
            Return[result];
           ];
 
-DoMediumDiagonalization[particle_Symbol, inputMomentum_:"", _] :=
-    "ERROR(\"medium diagonalization of " <> ToString[particle] <> " not implemented\");\n";
+DoMediumDiagonalization[particle_Symbol, __] := (
+    Print["Error: medium diagonalization of ", particle, " not implemented"];
+    "ERROR(\"medium diagonalization of " <> ToString[particle] <> " not implemented\");\n"
+);
 
 
 (* ********** high precision diagonalization routines ********** *)
 
+CalcEffPot2L[particle_] :=
+    Module[{dim = GetDimension[particle], dimStr, selfEnergyMatrixCType},
+           dimStr = ToString[dim];
+           selfEnergyMatrixCType = CreateCType[TreeMasses`GetMassMatrixType[particle]];
+           "\
+// two-loop Higgs self-energy contributions
+" <> selfEnergyMatrixCType <> " self_energy_2l(" <> selfEnergyMatrixCType <> "::Zero());
+
+if (pole_mass_loop_order > 1) {
+" <> IndentText["\
+self_energy_2l = self_energy_" <> CConversion`ToValidCSymbolString[particle] <> "_2loop();
+for (int i = 0; i < " <> dimStr <> "; i++) {
+   for (int k = 0; k < " <> dimStr <> "; k++) {
+      if (!std::isfinite(self_energy_2l(i,k))) {
+         self_energy_2l(i,k) = 0.;
+         problems.flag_bad_mass(" <> FlexibleSUSY`FSModelName <> "_info::" <> CConversion`ToValidCSymbolString[particle] <> ");
+      }
+   }
+}
+"] <> "}
+"
+          ];
+
 DoSlowDiagonalization[particle_Symbol, tadpole_] :=
     Module[{result, dim, dimStr, massName, inputMomenta, outputMomenta,
-            body, particleStr},
+            body, particleStr, effPot2L = ""},
            dim = GetDimension[particle];
            dimStr = ToString[dim];
            particleStr = CConversion`ToValidCSymbolString[particle];
            massName = ToValidCSymbolString[FlexibleSUSY`M[particle]];
            inputMomenta = "old_" <> massName;
            outputMomenta = "new_" <> massName;
-           body = DoMediumDiagonalization[particle, inputMomenta, tadpole] <> "\n" <>
+           If[dim > 1 && (SARAH`UseHiggs2LoopMSSM === True ||
+                          FlexibleSUSY`UseHiggs2LoopNMSSM === True) &&
+              MemberQ[{SARAH`HiggsBoson, SARAH`PseudoScalar}, particle],
+              effPot2L = CalcEffPot2L[particle];
+             ];
+           body = DoMediumDiagonalization[particle, inputMomenta, tadpole, effPot2L === ""] <> "\n" <>
                   outputMomenta <> " = PHYSICAL(" <> massName <> ");\n" <>
                   "diff = MaxRelDiff(" <> outputMomenta <> ", " <> inputMomenta <> ");\n" <>
                   inputMomenta <> " = " <> outputMomenta <> ";\n" <>
@@ -728,6 +744,7 @@ be set to the (positive) tree-level mass.  M_tree
                     "decltype(" <> massName <> ") " <>
                     inputMomenta  <> "(" <> massName <> "), " <>
                     outputMomenta <> "(" <> massName <> ");\n\n" <>
+                    effPot2L <> "\n" <>
                     "do {\n" <>
                     IndentText[body] <>
                     "\
@@ -747,7 +764,7 @@ be set to the (positive) tree-level mass.  M_tree
     "// diagonalization with low precision\n" <> DoFastDiagonalization[particle, tadpole];
 
 DoDiagonalization[particle_Symbol, FlexibleSUSY`MediumPrecision, tadpole_] :=
-    "// diagonalization with medium precision\n" <> DoMediumDiagonalization[particle, "", tadpole];
+    "// diagonalization with medium precision\n" <> DoMediumDiagonalization[particle, "", tadpole, True];
 
 DoDiagonalization[particle_Symbol, FlexibleSUSY`HighPrecision, tadpole_] :=
     "// diagonalization with high precision\n" <> DoSlowDiagonalization[particle, tadpole];

diff --git a/meta/Parameters.m b/meta/Parameters.m
@@ -73,6 +73,7 @@
 IsMatrix::usage="returns True if parameter is a matrix";
 IsSymmetricMatrixParameter::usage="returns True if parameter is a matrix";
 IsTensor::usage="returns True if parameter is a matrix";
+IsParameter::usage="returns True if symbol is a model/input/output/phase parameter";
 IsModelParameter::usage="returns True if parameter is a model parameter";
 IsInputParameter::usage="returns True if parameter is an input parameter";
 IsOutputParameter::usage="returns True if parameter is a defined output parameter";
@@ -671,6 +672,12 @@
 IsExtraParameter[parameter_[indices__] /; And @@ (IsIndex /@ {indices})] :=
     IsExtraParameter[parameter];
 
+IsParameter[sym_] :=
+    IsModelParameter[sym] ||
+    IsInputParameter[sym] ||
+    IsExtraParameter[sym] ||
+    IsPhase[sym];
+
 IsRealParameter[Re[sym_]] := True;
 IsRealParameter[Im[sym_]] := True;
 IsRealParameter[FlexibleSUSY`M[_]] := True;
@@ -853,6 +860,9 @@
 GetComplexTypeFromDimension[{dims__} /; Length[{dims}] > 2 && (And @@ (NumberQ /@ {dims}))] :=
     CConversion`TensorType[CConversion`complexScalarCType, dims];
 
+GetType[sym_[indices__] /; And @@ (IsIndex /@ {indices})] :=
+    CConversion`GetScalarElementType[GetType[sym]];
+
 GetType[FlexibleSUSY`SCALE] := GetRealTypeFromDimension[{}];
 
 GetType[FlexibleSUSY`M[sym_]] :=
@@ -867,9 +877,6 @@
 GetType[sym_] :=
     GetTypeFromDimension[sym, SARAH`getDimParameters[sym]];
 
-GetType[sym_[indices__] /; And @@ (IsIndex /@ {indices})] :=
-    GetType[sym]
-
 GetParameterDimensions[sym_ /; (IsInputParameter[sym] || IsExtraParameter[sym])] :=
     Module[{type},
            type = GetType[sym];