Skip to content

Commit

Permalink
Merge branch 'temp' into vertexRulesGMM2
Browse files Browse the repository at this point in the history
  • Loading branch information
Jobst Ziebell authored and Jobst Ziebell committed Oct 1, 2015
2 parents a717320 + 0775b10 commit 515e517
Show file tree
Hide file tree
Showing 4 changed files with 131 additions and 139 deletions.
2 changes: 1 addition & 1 deletion meta/FlexibleSUSY.m
View file
Expand Up @@ -989,7 +989,7 @@ corresponding tadpole is real or imaginary (only in models with CP

(* Write the GMM2 c++ files *)
WriteGMuonMinus2Class[vertexRules_List, files_List] :=
Module[{particles, muonFunctions, diagrams, vertexFunctions,
Module[{particles, muonFunctionPrototypes, diagrams, vertexFunctionData,
definitions, calculationCode, threadedCalculationCode},
particles = GMuonMinus2`CreateParticles[];
muonFunctionPrototypes = GMuonMinus2`CreateMuonFunctions[vertexRules][[1]];
Expand Down
30 changes: 8 additions & 22 deletions meta/GMuonMinus2.m
View file
Expand Up @@ -78,15 +78,15 @@ If you add new kinds of vertices (e.g for new diagram types):
muonIndex = GetMuonIndex[];
muonFamily = GetMuonFamily[];

prototypes = ("static const unsigned int muonIndex();\n" <>
"static const double muonPhysicalMass( EvaluationContext &context );\n" <>
"static const double muonCharge( EvaluationContext &context );");
prototypes = ("unsigned int muonIndex();\n" <>
"double muonPhysicalMass( EvaluationContext &context );\n" <>
"double muonCharge( EvaluationContext &context );");

definitions = ("static const unsigned int muonIndex()\n" <>
definitions = (" unsigned int muonIndex()\n" <>
"{ unsigned int muonIndex" <>
If[muonIndex =!= Null, " = " <> ToString[muonIndex-1], ""] <>
"; return muonIndex; }\n" <>
"static const double muonPhysicalMass( EvaluationContext &context )\n" <>
"double muonPhysicalMass( EvaluationContext &context )\n" <>
"{\n" <>
IndentText @
("static double m_muon_pole = 0.0;\n\n" <>
Expand All @@ -110,7 +110,7 @@ If you add new kinds of vertices (e.g for new diagram types):

"return m_muon_pole;\n") <>
"}\n" <>
"static const double muonCharge( EvaluationContext &context )\n" <>
"double muonCharge( EvaluationContext &context )\n" <>
"{ return 1.0; }");

muonFunctions = {prototypes, definitions};
Expand Down Expand Up @@ -200,8 +200,8 @@ If you add new kinds of vertices (e.g for new diagram types):
particles = Select[particles, (! TreeMasses`IsGhost[#] &)];

code = (StringJoin @
Riffle[("template<> double EvaluationContext::mass<" <> ToString[#] <> ">( " <>
If[TreeMasses`GetDimension[#] === 1, "", "unsigned int index"] <> " ) const\n" <>
Riffle[("template<> double EvaluationContext::mass<" <> ToString[#] <> ">(" <>
If[TreeMasses`GetDimension[#] === 1, "", " unsigned int index "] <> ") const\n" <>
"{ return model.get_M" <> ParticleToCXXName[#] <>
If[TreeMasses`GetDimension[#] === 1, "()", "( index )"] <> "; }"
&) /@ particles, "\n\n"]);
Expand Down Expand Up @@ -541,20 +541,6 @@ If you add new kinds of vertices (e.g for new diagram types):
If[indexParameters =!= "", indexParameters = " " <> indexParameters <> " "];

vertexClassName = SymbolName[VertexTypeForParticles[particles]];
(* vertexFunctionBody = Switch[vertexClassName,
"SingleComponentedVertex",
"return vertex_type( context.model." <>
NameOfCouplingFunction[particles] <>
"(" <> indexParameters <> ") );",
"LeftAndRightComponentedVertex",
"std::complex<double> left = context.model." <>
NameOfCouplingFunction[particles] <> "PL" <>
"(" <> indexParameters <> ");\n" <>
"std::complex<double> right = context.model." <>
NameOfCouplingFunction[particles] <> "PR" <>
"(" <> indexParameters <> ");\n\n" <>
"return vertex_type( left, right );"]; *)
vertexFunctionBody = Switch[vertexClassName,
"SingleComponentedVertex",
"return vertex_type( 0.0 );",
Expand Down
237 changes: 122 additions & 115 deletions templates/g_muon_minus_2.cpp.in
View file
Expand Up @@ -47,8 +47,7 @@

using namespace flexiblesusy;

namespace flexiblesusy {
namespace @ModelName@_GMuonMinus2 {
namespace {
/**
* @class IndexBounds<N>
* @brief A class representing multiple (N) index ranges.
Expand Down Expand Up @@ -215,6 +214,11 @@ template<class P> double EvaluationContext::mass(unsigned int index) const
*/
template<class ...Args> struct DiagramEvaluator;

double OneLoopFunctionF1C(double x);
double OneLoopFunctionF2C(double x);
double OneLoopFunctionF1N(double x);
double OneLoopFunctionF2N(double x);

@GMuonMinus2_MuonFunctionPrototypes@

@GMuonMinus2_Diagrams@
Expand Down Expand Up @@ -411,132 +415,136 @@ public:
};

@GMuonMinus2_VertexFunctionData@
}

double @ModelName@_GMuonMinus2::calculate_amuon(const @ModelName@_mass_eigenstates& model_)
{
// make copy since model is modified via call to calculate_MFe_pole()
@ModelName@_mass_eigenstates model(model_);

#ifdef ENABLE_THREADS
@GMuonMinus2_ThreadedCalculation@
#else
@GMuonMinus2_Calculation@
#endif
}

namespace {
/**
* @defgroup LoopFunctions
* @brief The loop functions necessary for GMM2 one-loop calculations.
*
* These are OneLoopFunctionF1C(), OneLoopFunctionF2C(),
* OneLoopFunctionF1N() and OneLoopFunctionF2N()
* as specified in arXiv:1311.1775
* @defgroup LoopFunctions
* @brief The loop functions necessary for GMM2 one-loop calculations.
*
* These are OneLoopFunctionF1C(), OneLoopFunctionF2C(),
* OneLoopFunctionF1N() and OneLoopFunctionF2N()
* as specified in arXiv:1311.1775
*/

static inline double OneLoopFunctionF1C(double x)
double OneLoopFunctionF1C(double x)
{
if (is_zero(x))
return 4.0;

/* error around x=1 is <= 10^-12 on an intel i7 */
double y = x - 1.0;
if (std::abs(y) < 0.21) {
return (1.0000000000000000000 -
0.60000000000000000000 * y +
0.40000000000000000000 * y * y -
0.28571428571428571429 * y * y * y +
0.21428571428571428571 * y * y * y * y -
0.16666666666666666667 * y * y * y * y * y +
0.13333333333333333333 * y * y * y * y * y * y -
0.10909090909090909091 * y * y * y * y * y * y * y +
0.090909090909090909091 * y * y * y * y * y * y * y * y -
0.076923076923076923077 * y * y * y * y * y * y * y * y * y +
0.065934065934065934066 * y * y * y * y * y * y * y * y * y * y -
0.057142857142857142857 * y * y * y * y * y * y * y * y * y * y * y +
0.050000000000000000000 * y * y * y * y * y * y * y * y * y * y * y * y -
0.044117647058823529412 * y * y * y * y * y * y * y * y * y * y * y * y * y +
0.039215686274509803922 * y * y * y * y * y * y * y * y * y * y * y * y * y * y -
0.035087719298245614035 * y * y * y * y * y * y * y * y * y * y * y * y * y * y * y);
}

return 2.0 / (y * y * y * y) * (2.0 + 3.0 * x - 6.0 * x * x + x * x * x + 6.0 * x * std::log(x));
if (is_zero(x))
return 4.0;

/* error around x=1 is <= 10^-12 on an intel i7 */
double y = x - 1.0;
if (std::abs(y) < 0.21) {
return (1.0000000000000000000 -
0.60000000000000000000 * y +
0.40000000000000000000 * y * y -
0.28571428571428571429 * y * y * y +
0.21428571428571428571 * y * y * y * y -
0.16666666666666666667 * y * y * y * y * y +
0.13333333333333333333 * y * y * y * y * y * y -
0.10909090909090909091 * y * y * y * y * y * y * y +
0.090909090909090909091 * y * y * y * y * y * y * y * y -
0.076923076923076923077 * y * y * y * y * y * y * y * y * y +
0.065934065934065934066 * y * y * y * y * y * y * y * y * y * y -
0.057142857142857142857 * y * y * y * y * y * y * y * y * y * y * y +
0.050000000000000000000 * y * y * y * y * y * y * y * y * y * y * y * y -
0.044117647058823529412 * y * y * y * y * y * y * y * y * y * y * y * y * y +
0.039215686274509803922 * y * y * y * y * y * y * y * y * y * y * y * y * y * y -
0.035087719298245614035 * y * y * y * y * y * y * y * y * y * y * y * y * y * y * y);
}
static inline double OneLoopFunctionF2C(double x)
{
/* error around x=1 is <= 10^-13 on an intel i7 */
double y = x - 1.0;
if (std::abs(y) < 0.155)
return (1.0 - 0.75 * y + 0.6 * y * y -
0.50000000000000000000 * y * y * y +
0.4285714285714285714 * y * y * y * y -
0.37500000000000000000 * y * y * y * y * y +
0.33333333333333333333 * y * y * y * y * y * y -
0.3000000000000000000 * y * y * y * y * y * y * y +
0.2727272727272727273 * y * y * y * y * y * y * y * y -
0.2500000000000000000 * y * y * y * y * y * y * y * y * y +
0.23076923076923076923 * y * y * y * y * y * y * y * y * y * y -
0.21428571428571428571 * y * y * y * y * y * y * y * y * y * y * y +
0.2000000000000000000 * y * y * y * y * y * y * y * y * y * y * y * y -
0.1875000000000000000 * y * y * y * y * y * y * y * y * y * y * y * y * y +
0.1764705882352941176 * y * y * y * y * y * y * y * y * y * y * y * y * y * y -
0.16666666666666666667 * y * y * y * y * y * y * y * y * y * y * y * y * y * y * y);

return -3.0 / (2.0 * y * y * y) * (-3.0 + 4.0 * x - x * x - 2.0 * std::log(x));

return 2.0 / (y * y * y * y) * (2.0 + 3.0 * x - 6.0 * x * x + x * x * x + 6.0 * x * std::log(x));
}

static inline double OneLoopFunctionF1N(double x)
double OneLoopFunctionF2C(double x)
{
if (is_zero(x))
return 2.0;

/* error around x=1 is <= 10^-12 on an intel i7 */
double y = x - 1.0;
if (std::abs(y) < 0.23)
return (1.0000000000000000000 -
0.4000000000000000000 * y +
0.2000000000000000000 * y * y -
0.11428571428571428571 * y * y * y +
0.07142857142857142857 * y * y * y * y -
0.04761904761904761905 * y * y * y * y * y +
0.03333333333333333333 * y * y * y * y * y * y -
0.02424242424242424242 * y * y * y * y * y * y * y +
0.0181818181818181818 * y * y * y * y * y * y * y * y -
0.01398601398601398601 * y * y * y * y * y * y * y * y * y +
0.01098901098901098901 * y * y * y * y * y * y * y * y * y * y -
0.0087912087912087912 * y * y * y * y * y * y * y * y * y * y * y +
0.00714285714285714286 * y * y * y * y * y * y * y * y * y * y * y * y -
0.0058823529411764706 * y * y * y * y * y * y * y * y * y * y * y * y * y +
0.0049019607843137255 * y * y * y * y * y * y * y * y * y * y * y * y * y * y -
0.0041279669762641899 * y * y * y * y * y * y * y * y * y * y * y * y * y * y * y);

return 2.0 / (y * y * y * y) * (1.0 - 6.0 * x + 3.0 * x * x + 2.0 * x * x * x - 6.0 * x * x * std::log(x));
/* error around x=1 is <= 10^-13 on an intel i7 */
double y = x - 1.0;
if (std::abs(y) < 0.155)
return (1.0 - 0.75 * y + 0.6 * y * y -
0.50000000000000000000 * y * y * y +
0.4285714285714285714 * y * y * y * y -
0.37500000000000000000 * y * y * y * y * y +
0.33333333333333333333 * y * y * y * y * y * y -
0.3000000000000000000 * y * y * y * y * y * y * y +
0.2727272727272727273 * y * y * y * y * y * y * y * y -
0.2500000000000000000 * y * y * y * y * y * y * y * y * y +
0.23076923076923076923 * y * y * y * y * y * y * y * y * y * y -
0.21428571428571428571 * y * y * y * y * y * y * y * y * y * y * y +
0.2000000000000000000 * y * y * y * y * y * y * y * y * y * y * y * y -
0.1875000000000000000 * y * y * y * y * y * y * y * y * y * y * y * y * y +
0.1764705882352941176 * y * y * y * y * y * y * y * y * y * y * y * y * y * y -
0.16666666666666666667 * y * y * y * y * y * y * y * y * y * y * y * y * y * y * y);

return -3.0 / (2.0 * y * y * y) * (-3.0 + 4.0 * x - x * x - 2.0 * std::log(x));
}
static inline double OneLoopFunctionF2N(double x)

double OneLoopFunctionF1N(double x)
{
if (is_zero(x))
return 3.0;

/* error around x=1 is <= 10^-13 on an intel i7 */
double y = x - 1.0;
if (std::abs(y) < 0.185)
return (1.0000000000000000000 -
0.50000000000000000000 * y +
0.30000000000000000000 * y * y -
0.2000000000000000000 * y * y * y +
0.14285714285714285714 * y * y * y * y -
0.10714285714285714286 * y * y * y * y * y +
0.08333333333333333333 * y * y * y * y * y * y -
0.06666666666666666667 * y * y * y * y * y * y * y +
0.05454545454545454545 * y * y * y * y * y * y * y * y -
0.0454545454545454545 * y * y * y * y * y * y * y * y * y +
0.0384615384615384615 * y * y * y * y * y * y * y * y * y * y -
0.03296703296703296703 * y * y * y * y * y * y * y * y * y * y * y +
0.0285714285714285714 * y * y * y * y * y * y * y * y * y * y * y * y -
0.02500000000000000000 * y * y * y * y * y * y * y * y * y * y * y * y * y +
0.0220588235294117647 * y * y * y * y * y * y * y * y * y * y * y * y * y * y -
0.0196078431372549020 * y * y * y * y * y * y * y * y * y * y * y * y * y * y * y);

return -3.0 / (y * y * y) * (1.0 - x * x + 2.0 * x * std::log(x));
if (is_zero(x))
return 2.0;

/* error around x=1 is <= 10^-12 on an intel i7 */
double y = x - 1.0;
if (std::abs(y) < 0.23)
return (1.0000000000000000000 -
0.4000000000000000000 * y +
0.2000000000000000000 * y * y -
0.11428571428571428571 * y * y * y +
0.07142857142857142857 * y * y * y * y -
0.04761904761904761905 * y * y * y * y * y +
0.03333333333333333333 * y * y * y * y * y * y -
0.02424242424242424242 * y * y * y * y * y * y * y +
0.0181818181818181818 * y * y * y * y * y * y * y * y -
0.01398601398601398601 * y * y * y * y * y * y * y * y * y +
0.01098901098901098901 * y * y * y * y * y * y * y * y * y * y -
0.0087912087912087912 * y * y * y * y * y * y * y * y * y * y * y +
0.00714285714285714286 * y * y * y * y * y * y * y * y * y * y * y * y -
0.0058823529411764706 * y * y * y * y * y * y * y * y * y * y * y * y * y +
0.0049019607843137255 * y * y * y * y * y * y * y * y * y * y * y * y * y * y -
0.0041279669762641899 * y * y * y * y * y * y * y * y * y * y * y * y * y * y * y);

return 2.0 / (y * y * y * y) * (1.0 - 6.0 * x + 3.0 * x * x + 2.0 * x * x * x - 6.0 * x * x * std::log(x));
}

double calculate_amuon(const @ModelName@_mass_eigenstates& model_)
double OneLoopFunctionF2N(double x)
{
// make copy since model is modified via call to calculate_MFe_pole()
@ModelName@_mass_eigenstates model(model_);

#ifdef ENABLE_THREADS
@GMuonMinus2_ThreadedCalculation@
#else
@GMuonMinus2_Calculation@
#endif
if (is_zero(x))
return 3.0;

/* error around x=1 is <= 10^-13 on an intel i7 */
double y = x - 1.0;
if (std::abs(y) < 0.185)
return (1.0000000000000000000 -
0.50000000000000000000 * y +
0.30000000000000000000 * y * y -
0.2000000000000000000 * y * y * y +
0.14285714285714285714 * y * y * y * y -
0.10714285714285714286 * y * y * y * y * y +
0.08333333333333333333 * y * y * y * y * y * y -
0.06666666666666666667 * y * y * y * y * y * y * y +
0.05454545454545454545 * y * y * y * y * y * y * y * y -
0.0454545454545454545 * y * y * y * y * y * y * y * y * y +
0.0384615384615384615 * y * y * y * y * y * y * y * y * y * y -
0.03296703296703296703 * y * y * y * y * y * y * y * y * y * y * y +
0.0285714285714285714 * y * y * y * y * y * y * y * y * y * y * y * y -
0.02500000000000000000 * y * y * y * y * y * y * y * y * y * y * y * y * y +
0.0220588235294117647 * y * y * y * y * y * y * y * y * y * y * y * y * y * y -
0.0196078431372549020 * y * y * y * y * y * y * y * y * y * y * y * y * y * y * y);

return -3.0 / (y * y * y) * (1.0 - x * x + 2.0 * x * std::log(x));
}

@GMuonMinus2_Definitions@
Expand Down Expand Up @@ -709,4 +717,3 @@ double DiagramEvaluator<OneLoopDiagram<4>, PhotonEmitter, ExchangeParticle>::val
return res;
}
}
}
1 change: 0 additions & 1 deletion templates/g_muon_minus_2.hpp.in
View file
Expand Up @@ -38,7 +38,6 @@ namespace @ModelName@_GMuonMinus2 {
*/
double calculate_amuon(const @ModelName@_mass_eigenstates& model);
} // namespace @ModelName@_GMuonMinus2

} // namespace flexiblesusy

#endif

0 comments on commit 515e517

Please sign in to comment.