Skip to content

Commit

Permalink
Begin refactoring of GMM2 into new module AMuon
Browse files Browse the repository at this point in the history
  • Loading branch information
iolojz committed Sep 3, 2017
1 parent 0c2ce7e commit bd3b586
Show file tree
Hide file tree
Showing 6 changed files with 703 additions and 11 deletions.
152 changes: 152 additions & 0 deletions meta/AMuon.m
@@ -0,0 +1,152 @@
(* ::Package:: *)

BeginPackage["AMuon`", {"SARAH`", "CConversion`", "TextFormatting`", "TreeMasses`", "LoopMasses`", "Vertices`"}];

(* The graphs that contribute to the EDM are precisely those with three
external lines given by the field in question, its Lorentz conjugate
and a photon.
They are given as a List of undirected adjacency matrices where
1 is the field itself
2 is its Lorentz conjugate
3 is the photon
and all other indices unspecified. *)
vertexCorrectionGraph = {{0,0,0,1,0,0},
{0,0,0,0,1,0},
{0,0,0,0,0,1},
{1,0,0,0,1,1},
{0,1,0,1,0,1},
{0,0,1,1,1,0}};
contributingGraphs = {vertexCorrectionGraph};

ContributingGraphs[] := contributingGraphs

GetPhoton[] := SARAH`Photon
GetMuon[] := If[TreeMasses`GetDimension[TreeMasses`GetSMMuonLeptonMultiplet[]] =!= 1,
TreeMasses`GetSMMuonLeptonMultiplet[],
Cases[SARAH`ParticleDefinitions[FlexibleSUSY`FSEigenstates],
{p_, {Description -> "Muon", ___}} -> p, 1][[1]]
]
GetMuonIndex[] := If[TreeMasses`GetDimension[TreeMasses`GetSMMuonLeptonMultiplet[]] =!= 1,
2,
Null]

ContributingDiagramsForGraph[graph_] :=
Module[{diagrams},
diagrams = CXXDiagrams`FeynmanDiagramsOfType[graph,
{1 -> GetMuonFamily[], 2 -> SARAH`AntiField[GetMuonFamily[]], 3 -> GetPhoton[]}];

Select[diagrams,IsDiagramSupported[graph,#] &]
]

IsDiagramSupported[vertexCorrectionGraph,diagram_] :=
Module[{photonEmitter,exchangeParticle},
photonEmitter = diagram[[4,3]]; (* Edge between vertices 4 and 6 (3rd edge of vertex 4) *)
exchangeParticle = diagram[[4,2]]; (* Edge between vertices 4 and 5 (2nd edge of vertex 4) *)

If[diagram[[6]] =!= {GetPhoton[],CXXDiagrams`LorentzConjugate[photonEmitter],photonEmitter},
Return[False]];
If[TreeMasses`IsFermion[photonEmitter] && TreeMasses`IsScalar[exchangeParticle],
Return[True]];
If[TreeMasses`IsFermion[exchangeParticle] && TreeMasses`IsScalar[photonEmitter],
Return[True]];

Return[False];
]

CreateInterfaceFunction[gTaggedDiagrams_List] :=
Module[{muon = GetMuon[], muonIndex = GetMuonIndex[],
prototype,definition,numberOfIndices},
numberOfIndices = NumberOfFieldIndices[muon];
prototype = "namespace " <> FlexibleSUSY`FSModelName <> "_a_muon {\n" <>
"double calculate_a_muon( " <>
"const " <> FlexibleSUSY`FSModelName <> "_mass_eigenstates& model );" <>
"\n}";

definition = "double " <> FlexibleSUSY`FSModelName <> "_edm::calculate_a_muon( " <>
"const " <> FlexibleSUSY`FSModelName <> "_mass_eigenstates& model_ )\n" <>
"{\n" <>
IndentText[
FlexibleSUSY`FSModelName <> "_mass_eigenstates model = model_;\n" <>
"EvaluationContext context{ model };\n" <>
"std::array<unsigned, " <> ToString @ numberOfIndices <>
"> indices = {" <>
If[muonIndex =!= Null,
ToString @ muonIndex <>
If[numberOfIndices =!= 1,
StringJoin @ Table[", 1", {numberOfIndices-1}],
""] <> " ",
If[numberOfIndices =!= 0,
StringJoin @ Riffle[Table[" 1", {numberOfIndices}], ","] <> " ",
""]
] <> "};\n\n" <>

"double val = 0.0;\n\n" <>

StringJoin @ Riffle[("val += " <> ToString @ # <> "::value(indices, context);") & /@
Flatten[CXXEvaluatorsForDiagramsFromGraph[#[[2]],#[[1]]] & /@ gTaggedDiagrams],
"\n"] <> "\n\n" <>

"return val;"
] <> "\n}";

{prototype, definition}
];

CXXEvaluatorsForDiagramsFromGraph[diagrams_,graph_] :=
CXXEvaluatorForDiagramFromGraph[#,graph] & /@ diagrams
CXXEvaluatorForDiagramFromGraph[diagram_,vertexCorrectionGraph] :=
Module[{photonEmitter,exchangeParticle},
photonEmitter = diagram[[4,3]]; (* Edge between vertices 4 and 6 (3rd edge of vertex 4) *)
exchangeParticle = diagram[[4,2]]; (* Edge between vertices 4 and 5 (2nd edge of vertex 4) *)

If[diagram[[6]] =!= {GetPhoton[],CXXDiagrams`LorentzConjugate[photonEmitter],photonEmitter},
Return["(unknown diagram)"]];
If[TreeMasses`IsFermion[photonEmitter] && TreeMasses`IsScalar[exchangeParticle],
Return[CXXEvaluatorFS[photonEmitter,exchangeParticle]]];
If[TreeMasses`IsFermion[exchangeParticle] && TreeMasses`IsScalar[photonEmitter],
Return[CXXEvaluatorSF[photonEmitter,exchangeParticle]]];

Return["(unknown diagram)"];
]

CXXEvaluatorFS[photonEmitter_,exchangeParticle_] :=
"AMuonVertexCorrectionFS<" <>
CXXDiagrams`CXXNameOfField[photonEmitter] <> ", " <>
CXXDiagrams`CXXNameOfField[exchangeParticle] <> ">"

CXXEvaluatorSF[photonEmitter_,exchangeParticle_] :=
"AMuonVertexCorrectionSF<" <>
CXXDiagrams`CXXNameOfField[photonEmitter] <> ", " <>
CXXDiagrams`CXXNameOfField[exchangeParticle] <> ">"

GetMinMass[particle_] :=
Module[{dim = TreeMasses`GetDimension[particle],
mStr = CConversion`ToValidCSymbolString[FlexibleSUSY`M[particle]],
tail},
If[dim == 1,
"model.get_" <> mStr <> "()",
tail = ToString[GetDimension[particle] - GetDimensionStartSkippingGoldstones[particle] + 1];
"model.get_" <> mStr <> "().tail<" <> tail <> ">().minCoeff()"
]
];

GetMSUSY[] :=
Module[{susyParticles},
susyParticles = Select[TreeMasses`GetSusyParticles[], IsElectricallyCharged];
If[susyParticles === {},
"return 0.;",
"return Min(" <>
StringJoin[Riffle[GetMinMass /@ susyParticles, ", "]] <>
");"
]
];

GetQED2L[] :=
"const double MSUSY = Abs(get_MSUSY(context.model));\n" <>
"const double m_muon = muonPhysicalMass(context);\n" <>
"const double alpha_em = Sqr(muonCharge(context))/(4*Pi);\n" <>
"const double qed_2L = alpha_em/(4*Pi) * 16 * FiniteLog(m_muon/MSUSY);\n\n" <>
"return qed_2L;";

EndPackage[];

37 changes: 35 additions & 2 deletions meta/FlexibleSUSY.m
Expand Up @@ -26,6 +26,7 @@
"Observables`",
"CXXDiagrams`",
"EDM2`",
"AMuon`",
"GMuonMinus2`",
"EDM`",
"EffectiveCouplings`",
Expand Down Expand Up @@ -1865,6 +1866,31 @@ corresponding tadpole is real or imaginary (only in models with CP
vertices
]

(* Write the AMuon c++ files *)
WriteAMuonClass[vertexRules_List, files_List] :=
Module[{graphs,diagrams,vertices,
interfacePrototype,interfaceDefinition,
getMSUSY, getQED2L},
graphs = AMuon`ContributingGraphs[];
diagrams = AMuon`ContributingDiagramsForGraph /@ graphs;

vertices = Flatten[CXXDiagrams`VerticesForDiagram /@ Flatten[diagrams,1],1];

{interfacePrototypes,interfaceDefinitions} =
AMuon`CreateInterfaceFunction @@@ Transpose[{graphs,#}] & /@ diagrams;

getMSUSY = GMuonMinus2`GetMSUSY[];
getQED2L = GMuonMinus2`GetQED2L[];

WriteOut`ReplaceInFiles[files,
{"@EDM2_InterfacePrototypes@" -> interfacePrototypes,
"@EDM2_InterfaceDefinitions@" -> interfaceDefinitions,
"@GMuonMinus2_GetMSUSY@" -> IndentText[WrapLines[getMSUSY]],
"@GMuonMinus2_QED_2L@" -> IndentText[WrapLines[getQED2L]],
Sequence @@ GeneralReplacementRules[]
}];
];

(* Write the GMM2 c++ files *)
WriteGMuonMinus2Class[vertexRules_List, files_List] :=
Module[{particles, muonFunctionPrototypes, diagrams, vertexFunctionData,
Expand Down Expand Up @@ -2941,7 +2967,7 @@ corresponding tadpole is real or imaginary (only in models with CP

MakeFlexibleSUSY[OptionsPattern[]] :=
Module[{nPointFunctions, runInputFile, initialGuesserInputFile,
edm2Vertices,edm2NPointFunctions,cxxVertexRules,
edm2Vertices,aMuonVertices,edm2NPointFunctions,cxxVertexRules,
gmm2Vertices = {}, edmFields,
susyBetaFunctions, susyBreakingBetaFunctions,
numberOfSusyParameters, anomDim,
Expand Down Expand Up @@ -3812,7 +3838,14 @@ corresponding tadpole is real or imaginary (only in models with CP
{FileNameJoin[{$flexiblesusyTemplateDir, "edm2.cpp.in"}],
FileNameJoin[{FSOutputDir, FlexibleSUSY`FSModelName <> "_edm2.cpp"}]}}];

WriteCXXDiagramClass[edm2Vertices,Lat$massMatrices,
aMuonVertices =
WriteAMuonClass[edmFields,
{{FileNameJoin[{$flexiblesusyTemplateDir, "a_muon2.hpp.in"}],
FileNameJoin[{FSOutputDir, FlexibleSUSY`FSModelName <> "_a_muon2.hpp"}]},
{FileNameJoin[{$flexiblesusyTemplateDir, "a_muon2.cpp.in"}],
FileNameJoin[{FSOutputDir, FlexibleSUSY`FSModelName <> "_a_muon2.cpp"}]}}];

WriteCXXDiagramClass[Join[edm2Vertices,aMuonVertices],Lat$massMatrices,
{{FileNameJoin[{$flexiblesusyTemplateDir, "cxx_diagrams.hpp.in"}],
FileNameJoin[{FSOutputDir, FlexibleSUSY`FSModelName <> "_cxx_diagrams.hpp"}]}}];

Expand Down

0 comments on commit bd3b586

Please sign in to comment.